二维地质模型的表面重力值和重力异常计算

二维地质模型的表面重力值和重力异常计算,分为正演和反演两个核心步骤。

---

一、正演计算（已知模型求重力异常）

核心公式（多边形截面二度体）：

其中：

• \(G=6.674×10^{−11}m^3/kg/s^2\)（引力常数）
• ρ为密度差（单位：kg/m³）
• (\(x_i,z_i\))为多边形顶点坐标（逆时针顺序）

MATLAB实现：

function grav = poly_gravity(x_obs, z_obs, x, z, rho)
    G = 6.67408e-11;
    n_poly = length(x);
    x(end+1) = x(1); z(end+1) = z(1); % 闭合多边形
    grav = zeros(size(x_obs));
    for i = 1:length(x_obs)
        total = 0;
        for j = 1:n_poly
            ax1 = x(j)*(1-t) + x(j+1)*t - x_obs(i);
            ax2 = z(j)*(1-t) + z(j+1)*t - z_obs;
            total = total - 2 * rho * G * (atan(ax1./ax2)) .* (z(j+1)-z(j));
        end
        grav(i) = 1e5 * sum(total); % 转换为mGal
    end
end

调用示例：

x_obs = -500:10:500;       % 观测点位置（m）
z_obs = 0;                 % 地表高程（m）
x = [-200,-100,0,100,200]; % 多边形顶点X坐标
z = [0,200,200,0,0];       % 多边形顶点Z坐标
rho = 500;                 % 密度差（kg/m³）
gravity_anomaly = poly_gravity(x_obs, z_obs, x, z, rho);
plot(x_obs, gravity_anomaly);

---

二、反演计算（从异常反推模型）

粒子群优化（PSO）算法流程：

1. 参数编码：粒子位置表示断层顶点坐标 (\(x_l,y_l,x_h,y_h\))

2. 正演调用：计算当前模型的重力异常 model_gravity = poly_gravity(...)

3. 适应度函数：最小化残差平方和

   \(Cost=∑(observed\_gravity−model\_gravity)^2+λ⋅smoothness\_penalty\)

4. 迭代更新：通过PSO更新粒子位置直至收敛

MATLAB代码片段：

options = optimoptions('particleswarm', 'SwarmSize', 40, 'MaxIterations', 300);
cost_func = @(params) norm(observed_gravity - poly_gravity(x_obs, z_obs, params(1:2), params(3:4), rho));
[optimal_params, ~] = particleswarm(cost_func, 4, lb, ub, options);

---

三、关键技术要点

1. 模型离散化：复杂地质体需剖分为三角形/多边形单元
2. 数值稳定性：采用高斯-勒让德求积避免奇点
3. 多解性处理：添加光滑约束项 λ∥∇ρ∥2正则化

---

四、工具推荐

1. 代码： 反演计算二维地质模型的表面重力值和重力异常 www.youwenfan.com/contentcno/96047.html
2. 专业软件： IGMAS（Linux平台，多边形建模） MATLAB地质建模工具箱（Delaunay三角剖分） www.docin.com/p-1316125714.html