重力反演理论例子.docxVIP

地球物理反演理论考核 姓名：陈高祥 学号：095211085 成绩: 一个有限延伸的二度板状体，我们假设板状体上顶面与下底面平行且水平， 板厚度一致并且内部均匀，定义上顶和下底宽度为 2b，延伸深度21，中心坐标 为(Xo,Zo)，密度差为c，倾角为a，如下图所示： 假设 1=10m，b=20m, c =500kg/m3, x0=0，z0=50m， 45。 请编制程序计算Z=0, x=-200m~200m范围内，点距为5m的重 力异常，并绘制重力异常曲线。 对于(1)中计算的重力异常，编制程序，加上5%的随机误差， 绘制含有误差的重力异常曲线。 把(2)中含有误差的重力异常曲线作为观测数据， 请用最小二 乘法进行反演，要求给出目标函数，收敛准则，程序框图，编 制程序，绘制观测数据和反演结果的对比曲线，并对反演结果 进行评价。 z=0）的重解：(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.a-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)).A2+xi.A2*(sin(a)).A4; p32=(h+xi*sin(a)*cos(a)).A2+xi.A2*(sin(a)).A4; p3=xi*(sin(a)).A2*log(p31/p32); p41=xi*(H-h)*(sin(a)).A2; p42=xi.A2*(sin(a)).A2+(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.A6;% 理论公式 %计算另一个台阶的重力异常 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)).A2+xi.A2*(sin(a)).A4; p32=(h+xi*sin(a)*cos(a)).A2+xi.A2*(sin(a)).A4; p3=xi*(sin(a)).A2*log(p31/p32); p41=xi*(H-h)*(sin(a)).A2; p42=xi.A2*(sin(a)).A2+(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.A6;% 理论公式 %求总异常 g(i)=g1(i)-g2(i); end 计算结果如图(1): 以n ■iw doo 3 o so割点横餌弼1i(n w XN?買由仲 ——国樹渥 r —?朝嗾 以n ■iw doo 3 o so 割点横餌弼1 i(n w XN ?買由仲 ——国樹渥 r —?朝嗾 _1?圖 a 0 50 側点橫处掾怖 roa