曲面上重力异常的计算之matlab程序

葛粲

在三维重力异常正演问题中，最常见的就是在直角坐标系下将异常体划分为一个个均匀密度的立方体，计算各个立方体的在测点产生的重力异常，然后累加起来。这种方法是最直接和显而易见的方法，但是在编程计算时往往效率很低，主要是涉及到多层循环（5层循环），既要对空间中的异常体循环（3层循环），又要对观测平面循环（2层循环）。利用matlab强大的矩阵计算功能，我们可以一次性获得一个密度异常立方体在整个观测平面上产生的观测异常。这样只需要对密度异常体进行循环并累加就可以完成三维重力异常正演，大大提高了计算效率。

一个密度异常立方体在观测点P产生的重力异常公式推导如下：

 http://s2/mw690/004f6djqgy6KJiWdaTvd1&690

http://s3/mw690/004f6djqgy6KJiWJ8jw12&690

http://s1/bmiddle/004f6djqgy6KJiWFTrid0&690

matlab子程序：

function [ gs ] = SurfaceGravityByPrism( xx,yy,zz,x,y,z,drho )

%LAYERCORRECTION Summary of this function goes here

% calculate the gravity on surf0 generating from a prism   

%[xx|yy|zz] define the depth of observation surface

%[x|y|z] define the prism position x[1]-x[2]; y[1]-y[2];z[1]-z[2]

%rho define the corrective density of the prism

%All input units must be referred to SI (m,kg/m,m/s).

%Output is in mgal (1 mgal=1e-5 m/s)

%example:

%% define prism parameter

%clear;

%x0=50;y0=50;z0=8.5;

%dx=40/2;dy=20/2;dz=16/2;

%rho=2.5e3;

%xi10=x0-dx;   xi20=x0+dx;

%eta10=y0-dy;  eta20=y0+dy;

%zeta10=z0-dz; zeta20=z0+dz;

%xi=[xi10,xi20];

%yi=[eta10,eta20];

%zi=[zeta10,zeta20];

%% define observation surface

%grid_dx=1;grid_dy=1;

%x=10:grid_dx:90;y=10:grid_dy:90;z=0;

%[xx,yy]=meshgrid(x,y);

%zz=rand(size(xx))-2*sin(xx/45);

%% calculate gravity at surface

% dg=SurfaceGravityByPrism(xx,yy,zz,xi,yi,zi,rho);

%%draw picture

%[c,h]=contour(xx,yy,dg*10,[0.5,1,2,2.5:1:10]);

% xlabel('x (m)')

% ylabel('y (m)')

% clabel(c,h)

G=6.67e-11*1e5;

v=0*zz;

for i=1:2

    for j=1:2

        for k=1:2

            sgn=(-1)^(i+j+k+1);

            r=sqrt((xx-x(i)).^2+(yy-y(j)).^2+(zz-z(k)).^2);                        

            gs=(  (x(i)-xx) .*  log((y(j)-yy)+r) ...

                 +(y(j)-yy) .*  log((x(i)-xx)+r) ...

                 -(z(k)-zz) .* atan((x(i)-xx).*(y(j)-yy)./(z(k)-zz)./r))*sgn;          

            v=v+gs;       

        end

    end

end

 gs=v*drho*G;

end

参考文献：

1，Richard J. Blakely,1995, POTENTIAL THEORY IN GRAVITY AND MAGNETIC APPLICATIONS CAMBRIDGE UNIVERSITY PRESS,P186-187

2，懿之的博客：http://hi.baidu.com/imheaventian/item/fdd743ea284b66e2fb42ba12