简介面向地球物理勘探与地质构造研究的重力异常正演工具包基于MATLAB实现专为需要定量模拟球体模型重力响应的研究人员与学习者设计。该程序可针对不同半径、密度和埋深参数快速计算理论重力异常生成可视化结果用于验证实测数据、辅助圈定矿藏或推断地下密度分布。压缩包共2个文件均为.m脚本总大小仅1KB结构轻量两个脚本分别负责构建球形地质体模型与完成重力异常数值计算和图形输出二者配合即可在MATLAB中执行完整正演流程。已有479人浏览学习适合课程实验、科研验证与重力正演算法的快速入门。通过这套脚本使用者既能深入理解球形体重力异常的计算原理也能直接修改参数开展多场景地质模拟为资源勘探与构造解释提供科学参考。1. 球体重力异常正演从解析公式到可复现的 gravity.zip 工作流如果一个刚入行重磁勘探的工程师手里只有一个gravity.zip程序包他最该先读懂的不是某个 GUI 界面而是一个半径为 100 m 的均质球体在埋深 200 m 时地表测得的重力异常曲线长什么样——因为球体是唯一一个能用初等函数完整描述重力异常的三维体模型。无论是矿体估算还是地热构造分析几乎所有三维异常体的正演都可以退化为「若干个球体叠加」的近似处理。本文从球体重力正演的解析基础出发把公式写成可直接运行的python代码给出剖面参数表和网格化参数设置再讨论埋深、半径、剩余密度三组参数的实际调节方向。整篇文章围绕gravity.zip里的核心任务——球体重力异常正演——展开读者拿到的将是一套能用于实际数据比对的工作底稿而不是教科书公式抄录。2. 球体重力正演的数学内核与单位体系2.1 从万有引力到垂直分量的重力异常表达式重力测量得到的不是引力场全向量而是其垂直于大地水准面方向的分量。一个密度均匀的球体在地表某观测点产生的重力异常其解析表达式为[ \Delta g \frac{G \cdot M \cdot z}{(x^2 z^2)^{3/2}} ]其中 (G) 是万有引力常量取值 (6.674 \times 10^{-11} , \mathrm{m^3 \cdot kg^{-1} \cdot s^{-2}})(M) 是球体总质量单位为 kg(z) 是球心到观测面的垂直距离即中心埋深单位为 m(x) 是观测点到球心在地表投影点的水平距离单位为 m。实际计算中(M) 用密度和体积写出(M \frac{4}{3}\pi a^3 \Delta\rho)(a) 是球体半径(\Delta\rho) 是球体与围岩的剩余密度差单位为 kg/m³。这里有一个必须强调的前提上述公式的成立条件是「观测点到球心的距离远大于球体半径」。当剖面测点非常靠近球体正上方且半径较大时公式会引入不可忽略的误差此时需要把球体剖分为多个小块分别计算再叠加。gravity.zip这类程序通常内嵌了这个近似判断逻辑但使用者在设置参数时需要自行确认模型尺度是否满足条件。一般来说如果 (z/a 3)点质量近似的误差可以控制在 5% 以内可以直接用解析公式。2.2 单位换算g.u.、mGal 与 μm/s² 的三重换算陷阱重力异常常用单位有三个m/s² 的国际单位制导出形式、g.u.gravity unit1 g.u. 1 μm/s²、以及历史遗留的 mGal1 mGal 10 μm/s²。在编写正演程序或解读结果时单位混用是最高频的「程序没错但数据全错」的根源。单位与 μm/s² 的关系使用场景m/s²1 m/s² 10⁶ μm/s²物理理论推导g.u.1 g.u. 1 μm/s²现代重力勘探成图mGal1 mGal 10 μm/s² 10 g.u.重磁资料解释报告gravity.zip中输入球体半径和埋深时若沿用米制输出异常若以 g.u. 为单位那么上述公式乘以 (10^6) 即可。以半径 100 m、埋深 200 m、剩余密度 1000 kg/m³ 的球体为例用 g.u. 表示的异常峰值约为地表直接叠加量级具体数值由下一章的代码给出。实际项目中布格重力异常的量级通常为数十到数百 g.u.一个地下矿体产生的局部异常往往只有几个到几十个 g.u.这要求正演程序的输出至少保留一位小数。2.3 剖面正演与平面网格正演的选择逻辑球体重力正演的两种任务是沿着一条测线计算剖面曲线或在平面网格上计算等值线图。剖面计算适合勘探初期的目标体埋深估算平面网格则适合与实测布格重力异常图直接对比。两者的数学本质相同只是观测点的坐标维度不同。剖面正演要指定测线方向与起止坐标、点距平面网格正演则要指定 X 方向和 Y 方向的起始点、终止点、间距。gravity.zip通常把这两种模式作为独立的功能模块参数文件在内容上有明显区别。3. 用 Python 重新实现球体重力正演的最小程序3.1 单球体剖面正演代码import numpy as np import matplotlib.pyplot as plt # 物理常量万有引力常数单位 m^3/(kg*s^2) G 6.674e-11 # 模型参数球体半径(m)球心埋深(m)剩余密度(kg/m^3) radius 100.0 depth 200.0 drho 1000.0 # 剖面参数起点、终点(m)点距(m) x_start -1000.0 x_end 1000.0 step 10.0 # 球体质量剩余质量 体积 * 剩余密度 volume 4.0 / 3.0 * np.pi * radius**3 mass volume * drho # 剖面坐标数组 x np.arange(x_start, x_end step, step) # 计算垂直方向重力异常结果为 g.u. (1 g.u. 1e-6 m/s^2) # 公式delta_g G * M * depth / (x^2 depth^2)^(3/2) delta_g G * mass * depth / (x**2 depth**2)**1.5 * 1e6 # 输出最大值及其位置 peak_index np.argmax(delta_g) print(f剖面长度: {len(x)} 个点) print(f最大重力异常: {delta_g[peak_index]:.2f} g.u.) print(f最大异常位于 x {x[peak_index]:.1f} m) # 绘制剖面曲线 plt.figure(figsize(10, 4)) plt.plot(x, delta_g, b-, linewidth1.5) plt.xlabel(水平距离 (m)) plt.ylabel(重力异常 (g.u.)) plt.title(球体正演剖面曲线) plt.grid(True, linestyle--, alpha0.6) plt.show()这段代码先算球体剩余质量再套用垂直分量公式计算每个观测点的异常值最后乘以 (10^6) 实现了从 m/s² 到 g.u. 的单位换算。需要注意的是np.argmax找出的峰值位置对应球心在地表的投影点这个位置理论上应该是剖面的中点如果异常峰值出现偏移说明输入的测线起止坐标没有把球心投影包含在内。剖面总长度设为埋深的 5 倍以上是因为球体重力异常在水平方向延伸较远截断边界会导致两端异常未归零影响后续的半幅宽度量取。3.2 参数配置文件gravity.ini 的设计思路gravity.zip这类程序通常以配置文件作为输入而不是把参数写死在代码里。一个实用的配置格式是 INI 风格分为[model]、[survey]、[output]三段。[model]段存放球体数量、每个球体的圆心坐标、半径、剩余密度[survey]段存放测线或网格的参数[output]段控制输出文件的格式和单位。多个球体的情况需要给每个球体编号程序按顺序遍历并累加所有球体在同一观测点的贡献。[model] sphere_count 1 sphere_1_center_x 0.0 sphere_1_center_y 0.0 sphere_1_center_z 200.0 sphere_1_radius 100.0 sphere_1_density 1000.0 [survey] mode profile x_start -1000.0 x_end 1000.0 point_spacing 10.0 [output] unit gu delimiter comma配置文件的优点是参数调整不需要重新编译或打开源代码适合批量试算。sphere_1_center_z是球心埋深不是球顶埋深这是常见的理解偏差——如果地质报告给出的是球顶埋深需要自行加上半径得到中心埋深。unit gu指定输出为 g.u.部分老程序支持mgal选项读取结果时要注意这一点。3.3 三维网格正演的数组化实现剖面正演只覆盖一条测线实际布格重力异常图是二维的因此网格正演是更常用的模式。将观测点坐标写成二维网格利用 numpy 的广播机制可以一次性算完全部点位的异常值。# 网格参数X方向起止与间距Y方向起止与间距 x_line np.arange(-1000.0, 1000.0 20.0, 20.0) y_line np.arange(-1000.0, 1000.0 20.0, 20.0) X, Y np.meshgrid(x_line, y_line) # 观测点到球心投影点的水平距离 r_horizontal np.sqrt(X**2 Y**2) # 代入球体重力异常公式 delta_g_grid G * mass * depth / (r_horizontal**2 depth**2)**1.5 * 1e6 # 输出峰值与网格形状 print(f网格规模: {X.shape[0]} x {X.shape[1]}) print(f网格峰值: {np.max(delta_g_grid):.2f} g.u.)网格正演的计算量随测点数量线性增长单球体在 100×100 的网格上计算只需要毫秒级的时间但多球体叠加时耗时按球体数量翻倍。如果球体数量超过几十个通常的做法是把球体坐标和半径向量化一次性用 numpy 完成所有球体对所有观测点的贡献计算而不是用 Python 的for循环逐球体累加。4. 参数敏感性分析与常见陷阱排查4.1 埋深、半径、剩余密度对异常曲线的独立影响球体重力异常的三个核心参数中埋深 (z) 的影响力最大。将公式简化为峰值表达式[ \Delta g_{\text{max}} \frac{G \cdot \frac{4}{3}\pi a^3 \Delta\rho}{z^2} ]峰值异常与埋深的平方成反比与半径的立方、剩余密度的一次方成正比。这意味着深部球体即使体积很大地表异常也可能被浅部小体积高密度体完全掩盖。半径与剩余密度的乘积等效于质量项单独调整半径或单独调整密度只要保持质量不变峰值异常不变——但曲线形态会变。半径越大而埋深不变意味着「观测点到球心的距离 / 球半径」之比变小异常曲线会在峰值附近变得更圆钝两侧衰减变缓。幅度特性参数 (A) 由球体质量与埋深决定形态参数由深度与点距的相对关系决定。实际操作中解释人员先根据异常峰值和半幅宽度估算深度范围再根据峰值反推质量项最后结合地质先验信息分配半径与密度——这个过程需要参数扫描来验证多解性范围。4.2 半幅宽度法估算埋深的经验规则球体重力异常剖面曲线的半幅宽度异常峰值一半处的水平宽度与埋深有近似的线性关系。对球体而言半幅宽度 (x_{1/2}) 与埋深 (z) 的关系满足[ z \approx 1.53 \times x_{1/2} ]这个系数是解析推导的结果适用于中心埋深与坐标原点的情形。实际操作时先在正演曲线上量取最大异常值的一半找到对应的两个横坐标计算两点之间的距离即为 (x_{1/2})再乘以 1.53 得到埋深估计。这个技巧在gravity.zip输出的剖面图场景中非常实用不需要额外写代码肉眼读图即可完成。需要特别注意半幅宽度法对噪声非常敏感。实测数据中叠加了区域场和随机噪声半幅宽度量取误差会被放大 1.53 倍。正演程序生成的理论曲线没有噪声问题但在与实测曲线对比时半幅宽度应该从「扣除区域场后的剩余异常」中量取而不是直接从原始布格异常上量取。4.3 剖面长度与网格范围的截断效应球体重力异常的理论衰减是 (1/x^2) 级别的与水准圆柱的 (1/x) 衰减不同衰减速度中等。如果剖面长度不够曲线两端不为零导致半幅宽度的基准线被抬高估算出的埋深偏大。设测线长度为 (L)球心埋深为 (z)当 (L/z 10) 时就必须考虑截断效应。推荐的剖面长度至少是埋深的 810 倍网格范围至少是埋深的 10 倍以上这保证边界处异常值衰减到峰值 (1%) 以下。如果程序输出的异常在边界处仍然有显著幅值解决方案是扩大测线范围而不是调整公式——因为公式本身没有错是观测范围没有覆盖球体的完整影响域。此外网格正演还要留意探测范围与球体投影位置的关系球心投影落在网格范围之外时整张网格都会表现为一个缓变的斜平面此时区间内看起来像线性背景的异常特征容易让人误判为存在区域场叠加球体的局部异常反而漏看。双球体模型验证了截断效应的影响力两球体间距 300 m、埋深分别为 150 m 和 250 m当测线覆盖长度达到两球体影响范围的边缘正演得到的曲线能明显看到两个叠加的峰值以及两个球体之间的鞍部。鞍部的深度与两球体的埋深差、水平距离比有关——正是这类合成模型的数值实验让解释人员建立起定性识别多源异常的经验。4.4 单位转换与球心/球顶埋深的两大致命错误球体重力正演常见的错误集中在两个环节。第一是防空区 —— 球心埋深和球顶埋深混用。若将球顶埋深 100 m、半径 100 m 的球体错把球心埋深设为 100 m峰值异常会被低估约 4 倍因为埋深从 100 m 变成 200 m 后峰值按平方反比衰减。第二是单位转换遗漏。若不乘以 (10^6)输出结果以 m/s² 为单位数值量级在 (10^{-6}) 左右绘制曲线时如果纵轴网格仍按整数刻度曲线就像一条水平直线容易被误判为「程序失效」。排查建议是先在球心投影正上方取一个点手算理论峰值并用计算器核对与代码输出比对无误后再进行批量计算。这个方法能在 5 分钟内定位绝大多数单位与深度的问题。5. 球体叠加正演与基于最小二乘的球心深度搜索5.1 多球体叠加正演应当使用向量化的计算机程序实际地质模型往往由多个球体组成例如两个相邻矿体、主矿体附带围岩蚀变体。多球体正演在代码层面只是单球体的循环叠加但球体之间的相互作用并不在重力异常计算中出现——重力满足叠加原理各球体的贡献独立相加即可。这一步有两个设计要点第一所有球体贡献的累加方向为正方向指向下的垂直分量不需要考虑球体间的引力影响第二当球体间距较小时叠加曲线可能只有一个峰峰值位置偏向质量较大的球体这就是通常说的「质量中心偏移效应」。# 多球体参数数组每行 [x中心, y中心, z中心(埋深), 半径, 剩余密度] spheres np.array([ [0.0, 0.0, 200.0, 100.0, 1000.0], [300.0, 0.0, 150.0, 80.0, 800.0], ]) x_line np.linspace(-500.0, 800.0, 131) total_anomaly np.zeros_like(x_line) # 遍历每个球体累加贡献 for sx, sy, sz, sr, sd in spheres: G_const 6.674e-11 s_mass 4/3 * np.pi * sr**3 * sd # 观测点到球心投影的水平距离剖面场景下y坐标取0 dist_h np.sqrt((x_line - sx)**2 (0.0 - sy)**2) # 累加重力异常 total_anomaly G_const * s_mass * sz / (dist_h**2 sz**2)**1.5 * 1e6在多球体场景下sz的取值直接决定各球体的相对贡献。深部大球体与浅部小球体的组合会产生明显的非对称曲线这种曲线与单一球体曲线存在显著差别前者可以通过左右不对称性和峰值的缓陡变化来识别。当两个球体的水平投影距离接近埋深之和时叠加曲线的主峰两侧会产生肩部或次级峰这些形态特征是判别多源的重要依据。5.2 基于网格搜索的最小二乘球心反演正演程序的进阶用途是反演即给定实测曲线搜索一组球体参数让正演结果最接近观测数据。最小二乘目标函数如下[ \chi^2 \sum_{i1}^N \left(\Delta g_{\text{obs}}^{(i)} - \Delta g_{\text{calc}}^{(i)}\right)^2 ]一个简单可靠的做法是网格搜索将埋深和半径的可能取值范围各自划分为若干候选值逐一正演计算并与观测曲线计算 (\chi^2)选择最小者作为估计结果。网格搜索虽然计算量大但不会陷入局部极小值适合球体数量较少的情形。通常设定的搜索范围要结合先验约束埋深搜索区间下限取半幅宽度法估算值的 0.5 倍上限取 2 倍剩余密度上下界根据岩石标本测量结果或地层资料确定。搜索步长的选取直接影响结果分辨率埋深方向的步长取 510 m半径方向的步长取 1020 m 比较合理。单球体在 50 个埋深候选值和 30 个半径候选值下总计算次数为 1500 次耗时极短完全可以在本地完成。5.3 用合成模型验证正演程序可靠性的三步流程拿到gravity.zip或自己写的正演程序第一步不是直接套用实测数据而是先完成三项验证。第一用解析峰值公式计算理论最大异常与程序输出对比相对误差应在 1% 以内。第二将输入的平均半径缩小到接近零的极限情况程序输出应接近点质量近似结果。第三两个完全相同球体叠加输出曲线应是单球体的两倍且峰值位置不变。这三步可以在十分钟内完成能排除绝大多数单位、公式和坐标系方向问题。之后再用球体正演去拟合实测剩余异常结合搜索得到的深度、质量与已知地质资料交叉验证推断的可靠性。正演的价值不在计算本身而是把模型假设变成可对照实测数据的定量预测从而约束解释人员对地下结构的想象空间。本文还有配套的精品资源点击获取