【MATLAB】球体重力异常正演

PourRevenir

已于 2023-12-01 12:48:07 修改

阅读量1.6k

点赞数 3

文章标签： matlab

于 2023-05-24 16:29:23 首次发布

原文链接：https://faculty.csu.edu.cn/xztong/zh_CN/zzcg/8518/content/1922.htm#zzcg

版权

在地球物理勘探工作中，一些近似等轴状的地质体可当作球体计算重力异常，这里为计算简便认为球体密度均匀，将地面当作水平面。

一、方程推导

设均匀球体埋深为D，半径为R，剩余密度为σ，剩余质量

$M= \frac{4}{3}{\pi}R^{3}\sigma$ 。

选球心在地面的投影位置为坐标原点，可以计算二度球体重力异常：

${\Delta}g=\frac{GMD}{{(x^{2}+D^{2})}^{\frac{3}{2}}}$

同理计算三度球体重力异常：

${\Delta}g=\frac{GMD}{{(x^{2}+y^{2}+D^{2})}^{\frac{3}{2}}}$

其中G为引力常数:

$G= 6.67\times 10^{-11} m^{3}/kg\cdot s^{2} = 6.67\times 10^{-2}(g.u.)m^{2}/t$

上式中D的单位为m（米），M的单位为t（吨），下同。

在原点处重力异常取得最大值为

$\Delta g_{max}=\frac{GM}{D^{2}}=6.67\times10^{-2}\frac{\left \{ M \right \}_{t}}{\left \{D \right \}^{2}_{m}}g.u.$

二、二度球体重力异常计算

前提埋深为100米，球体半径为50米，计算程序如下：

clear

% 定义参数
G = 6.67 * 1e-2;
R = 50;
D = 100;
sigma = 1;

% 定义绘图范围
x = -200:10:200;

% 计算重力异常
M = (4/3)*pi*(R^3)*sigma;
g = (G*M*D)./((x.^2+D^2).^1.5);

% 绘图
plot(x,g,'*-','MarkerEdgeColor','r');
box on;
xlabel('X/m');
ylabel('\Delta g/g.u.');

运行结果如下：

三、三度球体重力异常计算

计算程序如下：

clear

% 定义参数
G = 6.67 * 1e-2;
R = 50;
D = 100;
sigma = 1;

% 绘图范围
x = -200:10:200;
y = -200:10:200;

% 计算重力异常
M = (4/3)*pi*(R^3)*sigma;
for i = 1:size(x,2)
    for j = 1:size(x,2)
        g(i,j) = (G * M * D)/((x(i)^2 + y(j)^2 + D^2)^1.5);
    end
end

% 绘制等值线图
subplot(1,2,1)
contourf(x,y,g,12);
colorbar;
box on;
xlabel('X/m');
ylabel('Y/m');

% 绘制三维曲面图
subplot(1,2,2)
surf(x,y,g);
colorbar;
box on;
xlabel('X/m');
ylabel('Y/m');
zlabel('\Delta g/g.u.');

程序运行结果如下：