关于布格重力异常计算及资料处理与反演和解释的报告.doc

关于布格重力异常计算及资料处理与反演和解释的报告 姓名：林俊 班级：061084-27 学号：2081003195 指导老师：陈超 日期：2011.4.14目录前言2 目的2 任务要求2 工作过程2 成果2工作内容及步骤3 1-布格重力异常计算3 2-布格重力异常处理3 1.绘制平面等值线图3 2.异常处理（分离区域异常和局部异常）6 3-布格重力异常反演特征点法反演11 4-布格重力异常的解释13评述与结论13 评述13 结论14关于布格重力异常计算及资料处理与反演和解释的报告前言目的：熟悉并掌握布格重力异常计算及资料处理与反演和解释任务要求： 根据在一个地区重力测量的结果，计算出布格重力异常，并根据异常进行资料处理和解释，并完成一份工作报告。工作过程： (1)利用实测的相对重力值、相对高程值和X,Y坐标值，计算各种校正（地形校正除外），纬度校正用 计算，自由空间（或高度）校正用 计算，中间层校正用 计算，已知地表物质密度为2.50g/cm3，起算点纬度为45； (2)获得各点处的布格重力异常值后，绘出平面等值线图，等值线距为0.5mGal； (3)根据异常（平面或剖面）特征，选用适当的方法进行处理（如压制干扰、消除区域场等）进行处理，并对处理效果进行描述； (4)将处理后的异常进行反演；(5)写出全部过程和所采用的处理与反演方法之应用理由。成果：根据布格重力异常数据计算及资料处理与反演初步结果判断，该异常应由地区下一球体引起，球体埋深98.8m，剩余质量，球体中心在地面的投影点坐标为（248.8,248.8）m。工作内容及步骤1-布格重力异常计算 用excel先算出各项校正（除地形校正外），需要注意的是在纬度校正中 为测点到总基点间纬向距离，由于测点都位于总基点以北，故取正值； 为总基点纬度即45; 单位要划为km。中间层校正中 取值为2.50 。2-布格重力异常处理 1.绘制平面等值线图获得各点处的布格重力异常值后，用surfer绘图软件中的自然邻点网格法绘出平面等值线图，等值线距为0.5mGal，旁边为经过九点平滑法处理后的等值线图如下： 5mGal5mGal 图例： 等值线 图例： 等值线 可以看到经过平滑后的等值线光滑了很多，相当于虑去了高频成分。平滑处理的目的是消除异常数据中由观测引起的偶然误差以及由地表附近密度分布不均匀引起的杂乱无章的重力效应，以得到有意义的地质体引起的异常。 为了更好的了解该区域重力异常特性，再用surfer做出其一阶导平面等值线图和二阶导平面等值线图： 一阶导在（y=250）线上的剖面图 二阶导在（y=250）线上的剖面图 由上两图跟以下球形正演得出的平面等值线图和剖面图相比较，可以基本判定地下引起重力异常的物体形态似球，且根据二阶导平面等值线图还可以推断下面可能是由一个大球和四个小球叠加而产生的重力异常。这为后面的反演过程选用什么模型提供了依据。下面是球体重力异常一阶及二阶导正演剖面图： 2.异常处理（分离区域异常和局部异常）根据异常特征，我们知道该地区异常比较稳定，高频成分较少，在前面已摆出经过九点二次平滑法压制干扰后的等值线图，前后基本上变化不大，只是等值线更加光滑滑而已；而我们看到异常的导数等值线图把局部异常突出的比较好，因此我们可以用高阶导数法来提取局部重力异常。又我们很容易看出该地区区域重力异常呈线性递增，因此我们还可选用趋势分析（最小二乘多项式拟合）的方法进行处理消除区域场等，最后对处理效果进行描述。鉴于对编程能力的需要，高阶导数法难度较大，短期内无法实现，借助软件的话，我们不知道其演算过程，反演引起的误差也无法估计，因此我暂时选用趋势分析法（最小二乘多项式拟合法）。趋势分析法是用多项式拟合区域性背景场。其原理如下：选用一个曲面函数n阶多项式，用于拟合一个区域内的异常；通过若干点上异常值代入多项式可以得到一个方程组，求解方程组可得到多项式的系数，即得到描述区域异常的多项式函数；由于多项式的阶次一般视情况取14阶，去拟合相对简单的异常背景，以实现区域异常与局部异常的分离。eg.下面我先尝试着用matlab针对NE剖面（y=x）编了一个函数来做最小二乘多项式拟合法进行实验，下面是源代码：clear,clf,clcx=0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20;y=0.125527 0.7676005 1.2901735 1.789785 2.529205 3.127242 4.11251 5.200932 6.7870845 8.761621 9.9646955 9.618232 8.793035 8.1096465 8.027642 8.1744065 8.4864405 8.804667 9.265239 9.661196 10.209; cc=polyfit(x,y,1) %求出多项式前的各项系数，改变数字分别xx=0:1:20; 代表做不同次数拟合，此处是1为线性拟合yy=polyval(cc,xx);plot(xx,yy,-) xy=y-yy; hold on plot(x,xy,-)plot(x,y,x) axis(0,20,-2,11)hold on xlabel(x) ylabel(y)实验结果如下： 一次线性拟合 二次曲线（抛物线）拟合 五次多项式拟合 十次多项式拟合 二十次多项式拟合 用多项式拟合进行异常区分，存在两个方面的问题，一是选用多项式阶次问题，二是选择多少点参加拟合，不同的阶次和不同的拟合点数，其结果差异很大。（实验结果证实了不同的阶次对结果的影响） 1）多项式阶次越高，对原始异常拟合越好，异常区分效果往往不好。阶次太高，会造成趋势值包含过多的局部异常成分，因而计算出局部异常会偏小。此处明显线性拟合最好，得到了较好的局部异常值。 2）计算点数量越多，趋势异常越平缓简单，会造成对复杂背景拟合不到位的情况。此处剖面只用21个数据点，效果还好。 3）虚假异常问题及其消除的措施 用趋势分析方法区分异常，会产生局部虚假异常。此处两端出现负值就是最好的证明，我打算重新做一个线性变化的区域异常，即取两端点的连线。或用迭代计算法消除虚假异常。下面是处理结果： y=x向剖面图可以看到分离出来的局部异常曲线并不光滑，是因为没做平滑处理和点数太少的关系，接下来再对其剖面数据进行三点平均法并对分离出的局部异常数据进行三次样条插值：三点平均公式：三次样条插值matlab源代码：csfit.m文件：function s=csfit(x,y,dx0,dxn)n=length(x)-1;h=diff(x);d=diff(y)./h;a=h(2:n-1);b=2*(h(1:n-1)+h(2:n);c=h(2:n);u=6*diff(d);b(1)=b(1)-h(1)/2;u(1)=u(1)-3*(d(1);b(n-1)=b(n-1)-h(n)/2;u(n-1)=u(n-1)-3*(-d(n);for k=2:n-1 temp=a(k-1)/b(k-1); b(k)=b(k)-temp*c(k-1); u(k)=u(k)-temp*u(k-1);endm(n)=u(n-1)/b(n-1);for k=n-2:-1:1 m(k+1)=(u(k)-c(k)*m(k+2)/b(k);endm(1)=3*(d(1)-dx0)/h(1)-m(2)/2;m(n+1)=3*(dxn-d(n)/h(n)-m(n)/2;for k=0:n-1 s(k+1,1)=(m(k+2)-m(k+1)/(6*h(k+1); s(k+1,2)=m(k+1)/2; s(k+1,3)=d(k+1)-h(k+1)*(2*m(k+1)+m(k+2)/6; s(k+1,4)=y(k+1);end主函数：clear,clf,clcx=0*25*2(1/2) 1*25*2(1/2) 2*25*2(1/2) 3*25*2(1/2) 4*25*2(1/2) 5*25*2(1/2) 6*25*2(1/2) 7*25*2(1/2) 8*25*2(1/2) 9*25*2(1/2) 10*25*2(1/2) 11*25*2(1/2) 12*25*2(1/2) 13*25*2(1/2) 14*25*2(1/2) 15*25*2(1/2) 16*25*2(1/2) 17*25*2(1/2) 18*25*2(1/2) 19*25*2(1/2) 20*25*2(1/2);y=0.125527 0.727767 1.2825196667 1.8697211667 2.4820773333 3.256319 4.1468946667 5.200932 6.7870845 8.761621 9.9646955 9.618231 8.793035 8.1096465 8.1038983333 8.2294963333 8.4885046667 8.8521155 9.2437006667 9.7118116667 10.209; plot(x,y,x)hold on cc=(10.209-0.125527)/(20*25*2(1/2) 0.125527; xx=0:25*2(1/2):20*25*2(1/2); yy=polyval(cc,xx); plot(xx,yy,-) hold on xy=y-yy; dx0=0;dxn=0; s=csfit(x,xy,dx0,dxn) x1=0:0.01:1*25*2(1/2);y1=polyval(s(1,:),x1-x(1); x2=1*25*2(1/2):0.01:2*25*2(1/2);y2=polyval(s(2,:),x2-x(2); x3=2*25*2(1/2):0.01:3*25*2(1/2);y3=polyval(s(3,:),x3-x(3); x4=3*25*2(1/2):0.01:4*25*2(1/2);y4=polyval(s(4,:),x4-x(4); x5=4*25*2(1/2):0.01:5*25*2(1/2);y5=polyval(s(5,:),x5-x(5); x6=5*25*2(1/2):0.01:6*25*2(1/2);y6=polyval(s(6,:),x6-x(6); x7=6*25*2(1/2):0.01:7*25*2(1/2);y7=polyval(s(7,:),x7-x(7); x8=7*25*2(1/2):0.01:8*25*2(1/2);y8=polyval(s(8,:),x8-x(8); x9=8*25*2(1/2):0.01:9*25*2(1/2);y9=polyval(s(9,:),x9-x(9); x10=9*25*2(1/2):0.01:10*25*2(1/2);y10=polyval(s(10,:),x10-x(10); x11=10*25*2(1/2):0.01:11*25*2(1/2);y11=polyval(s(11,:),x11-x(11); x12=11*25*2(1/2):0.01:12*25*2(1/2);y12=polyval(s(12,:),x12-x(12); x13=12*25*2(1/2):0.01:13*25*2(1/2);y13=polyval(s(13,:),x13-x(13); x14=13*25*2(1/2):0.01:14*25*2(1/2);y14=polyval(s(14,:),x14-x(14); x15=14*25*2(1/2):0.01:15*25*2(1/2);y15=polyval(s(15,:),x15-x(15); x16=15*25*2(1/2):0.01:16*25*2(1/2);y16=polyval(s(16,:),x16-x(16); x17=16*25*2(1/2):0.01:17*25*2(1/2);y17=polyval(s(17,:),x17-x(17); x18=17*25*2(1/2):0.01:18*25*2(1/2);y18=polyval(s(18,:),x18-x(18); x19=18*25*2(1/2):0.01:19*25*2(1/2);y19=polyval(s(19,:),x19-x(19); x20=19*25*2(1/2):0.01:20*25*2(1/2);y20=polyval(s(20,:),x20-x(20);plot(x1,y1,x2,y2,x3,y3,x4,y4,x5,y5,x6,y6,x7,y7,x8,y8,x9,y9,x10,y10,x11,y11,x12,y12,x13,y13,x14,y14,x15,y15,x16,y16,x17,y17,x18,y18,x19,y19,x20,y20,x,xy,.) grid on xlabel(x)ylabel(y)处理结果如下： y=x向剖面图3-布格重力异常反演 根据重力异常的平面和剖面特征，初步选用球体模型做反演：特征点法反演g(Vz): 当x=351.8m时， 有 当，即时， 当，即 时， 当，即时， ， 球体中心点地面投影坐标：即（248.8,248.8）m 由于求出的几个D值比较接近，说明场源体确实近于球形，因此埋深D取其平均值。若求得D值彼此差别甚大或沿深度方向有规律的递增，则说明该异常不能当作球体引起的异常来解释，或者场源体沿深度还有一定的延伸。4-布格重力异常的解释由该工区的平面等值线图，可以看出区域性异常的走向为NE（SW）向，变化特征为从SW到NE异常呈线性递增，幅度大约为10mGal；局部重力异常等值线为圆形圈闭，呈等轴状，异常值中心高，四周低，有极大值，表现为重力高，幅值接近10mGal，初步推断相对应的几何形体为剩余密度为正值的均匀球体、铅直圆柱体等。通过，R3.24g/cm3（假设围岩密度=2.50g/cm3）。可能反映的地质因素为：透镜状的致密金属矿体，如铬铁矿、铁矿、铜矿等；中基性岩浆（密度较高）的侵入体，形成岩株状，穿插在较低密度的岩体或地层中；高密度岩层形成的穹窿、等轴背斜等；松散沉积物下面的基岩（密度较高）局部隆起；低密度岩层形成的向斜或凹陷内充填了高密度的岩体，如砾石等。评述与结论 评述： 1、重力测量工作程度 ：良好，总共441个数据点，21条测线，每条测线上为21个测点。工区面积为500m500m，很好的涵盖了有用