重力反演理论例子.doc

地球物理反演理论考核 姓名：陈高祥 学号：095211085 成绩： 一个有限延伸的二度板状体，我们假设板状体上顶面与下底面平行且水平，板厚度一致并且内部均匀，定义上顶和下底宽度为2b，延伸深度2l，中心坐标为（X0,Z0），密度差为σ，倾角为α，如下图所示： 假设l=10m，b=20m, σ=500kg/m3, x0=0，z0=50m，。请编制程序计算Z=0，x=-200m~200m范围内，点距为5m的重力异常，并绘制重力异常曲线。 对于（1）中计算的重力异常，编制程序，加上5%的随机误差，绘制含有误差的重力异常曲线。 把（2）中含有误差的重力异常曲线作为观测数据，请用最小二乘法进行反演，要求给出目标函数，收敛准则，程序框图，编制程序，绘制观测数据和反演结果的对比曲线，并对反演结果进行评价。 解：（1）对于台阶的重力异常，其存在解析解。此二度体在测线上（z=0）的重力异常计算可以看作两个倾斜台阶的重力异常差。 用MATLAB编写的程序及其计算结果如下 ⅰ）重力异常正演模拟程序 %计算平行四边形板状体的重力异常 clear tic;%计算运行时间 Xc=[-200:5:200;]; %测点横坐标 %给出正演模型参数T R=500; %剩余密度/ kg/m3 B=20; %上顶和下底的宽度的一半/m L=10; %延伸深度的一半/m X0=0;Z0=50; %中心坐标/m a=45; %倾角/度 T(1)=R; T(2)=B; T(3)=L; T(4)=X0; T(5)=Z0; T(6)=a; g=gravity(Xc,T); %调用子函数 plot(Xc,g) ylabel(\Deltag/g.u.,FontSize,15) xlabel(测点横坐标/m,FontSize,15) toc ⅱ）正演模拟调用子函数 %平行四边形板状体的重力异常调用子函数 function g=gravity(Xc,T) %g.u. %Xc为测点横坐标，R为剩余密度，B为上顶和下底的宽度的一半， %L为延伸深度的一半，X0、Y0为中心坐标，a为倾角（度） R=T(1); B=T(2); L=T(3); X0=T(4); Z0=T(5); a=(180-T(6))*pi/180; G=6.720*10.^-11; %引力常量 h=Z0-L;H=Z0+L; for i=1:length(Xc) %计算一个台阶的重力异常 xi=Xc(i)-(Z0*cot(a)+(X0-B));%坐标变换 p1=pi*(H-h)+2*H*atan((xi+H*cot(a))/H); p2=-2*h*atan((xi+h*cot(a))/h); p31=(H+xi*sin(a)*cos(a)).^2+xi.^2*(sin(a)).^4; p32=(h+xi*sin(a)*cos(a)).^2+xi.^2*(sin(a)).^4; p3=xi*(sin(a)).^2*log(p31/p32); p41=xi*(H-h)*(sin(a)).^2; p42=xi.^2*(sin(a)).^2+(H+h)*xi*sin(a)*cos(a)+H*h; p4=-2*xi*sin(a)*cos(a)*atan(p41/p42); g1(i)=G*R*(p1+p2+p3+p4)*10.^6;%理论公式 %计算另一个台阶的重力异常 xi=Xc(i)-(Z0*cot(a)+(X0+B));%坐标变换 p1=pi*(H-h)+2*H*atan((xi+H*cot(a))/H); p2=-2*h*atan((xi+h*cot(a))/h); p31=(H+xi*sin(a)*cos(a)).^2+xi.^2*(sin(a)).^4; p32=(h+xi*sin(a)*cos(a)).^2+xi.^2*(sin(a)).^4; p3=xi*(sin(a)).^2*log(p31/p32); p41=xi*(H-h)*(sin(a)).^2; p42=xi.^2*(sin(a)).^2+(H+h)*xi*sin(a)*cos(a)+H*h; p4=-2*xi*sin(a)*cos(a)*atan(p41/p42); g2(i)=G*R*(p1+p2+p3+p4)*10.^6;%理论公式 %求总异常 g(i)=g1(i)-g2(i); end 计算结果如图（1）： （a） 倾角不同 (b)倾角为45度 图（1）模拟重力异常曲线图 （