[转载]二度体重力正演

原文地址：二度体重力正演作者：用户crust

教科书中公式有误，借助Mathematica5.0软件推导了新公式，教科书的公式多了个2倍系数。

二度体重力正演公式推导：

为了测试公式的正确性，利用matlab进行了编程验证。

分别用二度体，三度体的公式计算了圆柱体的重力异常，验证了公式和编程的正确性。

 

matlab编程验证：

clear all;
G=6.67e-11*1e5;
%单位mgal
%二度体重力正演

%模型尺度
lx=100;
dx=2;
lz=100;
dz=2;
xx=(0:lx-1)*dx;
zz=(1:lz)*dz;
[xx,zz]=meshgrid(xx,zz);
%密度模型
drho_model=zeros(lx,lz);
%均匀圆柱体
    %半径 22.5（m）
    radius=22.5;
    %埋深 60（m）
    depth=60;
    %距离x轴零点距离45m
    x0=45;
    %密度 1000（kg/m3）
    drho=1000;
    drho_model(((xx-x0).^2+(zz-depth).^2)
    subplot(2,1,2);
    hold off;
    contour(drho_model);
    set(gca,'ydir','reverse');
%观测点位置：
xo=(0:lx-1)*dx;
zo=0;

%1简化公式
g(1:lx)=0;
gg=g;
    for (i=1:lx)
        for(j=1:lz)           
            if (drho_model(i,j)~=0)              
                gg=drho_model(i,j)*(zz(i,j))*dx*dz./((xo-xx(i,j)).^2+(zz(i,j))^2);              
                g=g+gg;
            end
        end
    end
g=g*2*G;
subplot(2,1,1);
hold off;
set(gca,'ydir','normal');

plot(g);

%2二度体公式
g(1:lx)=0;
gg=g;
for (k=1:lx)
    for (i=1:lx)
        for(j=1:lz)
            x1=xx(i,j)-dx/2-xo(k);
            x2=xx(i,j)+dx/2-xo(k);
            z1=zz(i,j)-dz/2-zo;
            z2=zz(i,j)+dz/2-zo;
            if (drho_model(i,j)~=0)
                gg=x1*log(x1^2+z1^2)+2*z1*atan(x1/z1)...
                    -x1*log(x1^2+z2^2)-2*z2*atan(x1/z2)...
                    -x2*log(x2^2+z1^2)-2*z1*atan(x2/z1)...
                    +x2*log(x2^2+z2^2)+2*z2*atan(x2/z2);
                g(k)=g(k)+gg*drho_model(i,j);
            end
        end
    end
end
g=g*G;
hold on;

plot(g,'red');
%3圆柱体理论公式
for(i=1:lx)
    gt(i)=2*G*3.14*radius^2*1000*depth/((xo(i)-x0)^2+depth^2);
end
plot(gt,'black');
%4三度体近似解
g(1:lx)=0;
gg=0;
for (k=1:lx)
    for (i=1:lx)
        for(j=1:lz)
            x1=xx(i,j)-dx/2;
            x2=xx(i,j)+dx/2;
            xxx=[x1,x2];
            yyy=[-1e5,1e5];
            z1=zz(i,j)-dz/2;
            z2=zz(i,j)+dz/2;
            zzz=[z1,z2];
            if (drho_model(i,j)~=0)
                gg=GravityByPrism(xo(k),0,0,xxx,yyy,zzz,1000);
                g(k)=g(k)+gg;
            end
        end
    end
end
plot(g,'green');

 

计算结果：

4种公式计算结果基本一致。
 

上图中上部是地表的重力异常曲线，下部是重力异常体的位置。