重力异常正演总结

一、正方体重力异常正演推导过程（蒋）直立长方体外任意一点的重力异常公式推导如下：p(x,y,z)：长方体外任一点：长方体密度长方体内任一点坐标222111222321112222221()()()()(,,)[()()()]abhabhabhabhGMVGdddrxyzVzgxyzGdddZxyz设：222()()()rxyz；ln[()]uyr；则1.()xdudyrr要求原式：222111222ln[()()()()]hbaabhGyxyzd在此先求：222ln[()()()()]yxyzd22.()11.()()...()()().ln[()].[()].udxuxxdyrrxxyrdyrr上式第二项：2().[()].xdyrr222222222().[()].[()()].()()().()()...()()[()()].xrydxzrzyyzdxdddxzrxzr上式最后一项：222().().[()()].yzdxzr2222222222222222222().().()2().()222().().[()()]..[()()].[()()]().().[()]..[().().()]().()..()().().yzxrxyryzyzdrxzyzyzrxdrxyrzyzdrzxy212222().[.()]()..[().().()]().()rzxrzdxyrzyx212().[.()]()..().()().()rzxrzrddyxyx2222.()().()1().().().()1().().arctan().()rzxyzrzdyxzrzyx由上述推导过程，可以得出：222ln[()()()()]yxyzd().ln[()]().arctanxxyrzz().().ln[()]().arctan().()zryxrzyx222111222(,,)ln[()()()()]hbaabhgxyzGyxyzd222111().().ln[()]().arctan().ln[()]().arctan().()hbaabhxzrGxyrzyxrzzyx将上式完全展开后，().arctanxzz可以消除。222111().().ln[()]().ln[()]().arctan().()hbaabhzrGxyryxrzyx由前面的推导结果，先假设一点该点在面上11326.67210/(.)Gmkgs3332.67/2.6710/gcmkgm假设长方体范围为：100100公里，高公里从岩石圈动力学图上，大概间隔即假设：12100,100100aa12100,100100bb12100,30100hh立方体体积：143310m由上述数据算得：20.00816061/gms沿轴平行方向两边同时压缩米后：121212100,100100101,10009999.75,30100.25aabbhh，，20.00816005/gms820.5610/0.56ggmsGal沿轴平行方向两边同时压缩米后：121212282101,100099100,10010099.75,30100.250.00815996/0.6510/0.65aabbhhgmsggmsGal，，（上述单位均为：米）对于剩余密度均匀的正方体来说，它与将其全部剩余质量集中在其质心处的点质量所产生的异常完全一样。设正方体中心埋深为，边长为。将坐标原点选在正方体中心在地面的投影点上。我们取正方体在上的中心剖面进行分析，其他与此剖面类似。g=G(为体积元的剩余密度。不同边长的正方体只影响异常的绝对值，并不影响其分布形态。Matlab画图程序：%长方体模型参数%a=2000;%长b=200;%宽c=100;%高x0=0;%质心坐标（x0,y0,z0)y0=0;H=1000;%质心埋深z0ph=2*10^3;%剩余密度x=(-2000:50:2000);y=(-2000:50:2000);z=0;G=6.67e-11;[x1,y1]=meshgrid(x,y);%生成网线节点矩阵r1=sqrt((x0+a/2-x1).^2+(y0+b/2-y1).^2+(H+c/2-z).^2);r2=sqrt((x0+a/2-x1).^2+(y0+b/2-y1).^2+(H-c/2-z).^2);r3=sqrt((x0+a/2-x1).^2+(y0-b/2-y1).^2+(H+c/2-z).^2);r4=sqrt((x0+a/2-x1).^2+(y0-b/2-y1).^2+(H-c/2-z).^2);r5=sqrt((x0-a/2-x1).^2+(y0+b/2-y1).^2+(H+c/2-z).^2);r6=sqrt((x0-a/2-x1).^2+(y0+b/2-y1).^2+(H-c/2-z).^2);r7=sqrt((x0-a/2-x1).^2+(y0-b/2-y1).^2+(H+c/2-z).^2);r8=sqrt((x0-a/2-x1).^2+(y0-b/2-y1).^2+(H-c/2-z).^2);g1=G*ph*((x0+a/2-x1).*log((y0+b/2-y1)+r1)+(y0+b/2-y1).*log((x0+a/2-x1)+r1)+(H+c/2-z).*atan((x0+a/2-x1).*(y0+b/2-y1))./((H+c/2-z).*r1));g2=G*ph*((x0+a/2-x1).*log((y0+b/2-y1)+r2)+(y0+b/2-y1).*log((x0+a/2-x1)+r2)+(H-c/2-z).*atan((x0+a/2-x1).*(y0+b/2-y1))./((H-c/2-z).*r2));g3=G*ph*((x0+a/2-x1).*log((y0-b/2-y1)+r3)+(y0-b/2-y1).*log((x0+a/2-x1)+r3)+(H+c/2-z).*atan((x0+a/2-x1).*(y0-b/2-y1))./((H+c/2-z).*r3));g4=G*ph*((x0+a/2-x1).*log((y0-b/2-y1)+r4)+(y0-b/2-y1).*log((x0+a/2-x1)+r4)+(H-c/2-z).*atan((x0+a/2-x1).*(y0-b/2-y1))./((H-c/2-z).*r4));g5=G*ph*((x0-a/2-x1).*log((y0+b/2-y1)+r5)+(y0+b/2-y1).*log((x0-a/2-x1)+r5)+(H+c/2-z).*atan((x0-a/2-x1).*(y0+b/2-y1))./((H+c/2-z).*r5));g6=G*ph*((x0-a/2-x1).*log((y0+b/2-y1)+r6)+(y0+b/2-y1).*log((x0-a/2-x1)+r6)+(H-c/2-z).*atan((x0-a/2-x1).*(y0+b/2-y1))./((H-c/2-z).*r6));g7=G*ph*((x0-a/2-x1).*log((y0-b/2-y1)+r7)+(y0-b/2-y1).*log((x0-a/2-x1)+r7)+(H+c/2-z).*atan((x0-a/2-x1).*(y0-b/2-y1))./((H+c/2-z).*r7));g8=G*ph*((x0-a/2-x1).*log((y0-b/2-y1)+r8)+(y0-b/2-y1).*log((x0-a/2-x1)+r8)+(H-c/2-z).*atan((x0-a/2-x1).*(y0-b/2-y1))./((H-c/2-z).*r8));gg=-(g1-g2-g3+g4-g5+g6+g7-g8)*10^5;%单位mGalfigure(1)mesh(x1,y1,gg)%三维xlabel('')ylabel('')title('均匀长方体重力异常')figure(2)contourf(x1,y1,gg)%二维title('均匀长方体重力异常')Figure1Figure2二、球体重力异常正演推导过程（李）对于剩余密度均匀的球体来说，他与将其全部质量集中在球心处的点质量所产生的异常完全一样。设球心的埋深为D，半径为R，则它的剩余质量34=R3M。为使计算简化，将坐标原点O选在球心于地面的投影点上。由对称性可知，只需研究过原点O的任意剖面上异常的分布即可。设该剖面与X轴中心重合，则剖面上任一点P,0,0x处的重力异常表达式为：3/222gGMDxD球体异常质量重力异常分布精确推导：球体半径为剩余密度为球外一点距离球心()xxR点处的重力为：22222223/2000cossin(2cos)vRGmGgvllGgdvlxrgGrdrddxrrx异常质量体半径大小对重力异常值得影响：3=0.35g/cm，500Dm，11226.67210/GNmkg。（）1150Rm，（）2175Rm，3200RmMatlab画图代码：G=6.672*10^(-11);R1=150;R2=175;R3=200;den=0.35;dM1=4/3*pi*R1^3*den;dM2=4/3*pi*R2^3*den;dM3=4/3*pi*R3^3*den;D=500;x=-2000:2:2000;g1=G.*dM1.*D./((x.^2+D^2).^(3/2));g2=G.*dM2.*D./((x.^2+D^2).^(3/2));g3=G.*dM3.*D./((x.^2+D^2).^(3/2));plot(x,g1,'-.r');xlabel('海拔高程(m)'),ylabel('重力异常值(m/s^2)'),title('重力异常随高度的变化'),legend('R=150m');figure(2);plot(x,g1,'r',x,g2,'b',x,g3,'k');xlabel('海拔高程(m)');yl