MATLAB实现水平圆柱体重力异常正演：从原理到地质解释

1. 项目概述从“异常”中寻找规律搞地球物理勘探或者地质工程的朋友对“重力异常”这个词肯定不陌生。简单来说我们脚下的大地其密度分布并不是均匀的。一个埋藏在地下的矿体、一个地质构造甚至一个古代遗存都会因为其密度与周围岩土存在差异导致其所在位置的地球重力场产生微小的、局部的变化。这个变化就是重力异常。我们的工作就是在地面上布设测点用高精度的重力仪测量这些微小的变化然后像解谜一样反推地下到底是什么东西引起了这个异常。“正演”就是这个解谜过程的“前半部分”或者说是构建谜题的过程。假设我们已经知道了地下物体的形状、大小、埋深、密度差那么理论上它在地表会产生一个什么样的重力异常响应计算出这个理论响应值就是正演。它是反演即从实测异常数据推断地下物体属性的基础和前提。只有正演模型足够准确、计算足够高效我们后续的反演解释才靠谱。这次要聊的就是用 MATLAB 来模拟一个非常经典的地球物理模型水平圆柱体的重力异常正演。为什么是水平圆柱体因为它形状规则其重力异常有解析解也就是可以用一个数学公式精确计算出来是验证算法、理解物理概念、乃至模拟一些近似柱状的地质体如盐丘、矿脉、隧道、管线的绝佳起点。很多复杂的模型初期都可以用一系列不同参数的圆柱体来近似。所以这个项目的核心价值在于掌握用数值计算工具MATLAB实现地球物理基本正演模型的能力为后续更复杂的反演解释和实际数据处理打下坚实的基础。无论你是地质、测绘、地球物理专业的学生还是相关领域的工程师自己动手实现一遍远比只看教科书上的公式印象深得多。2. 核心原理与模型解析公式背后的物理图景在动手写代码之前我们必须把模型和公式吃透。一知半解地套公式出了错都不知道在哪。2.1 水平圆柱体模型假设我们首先明确模型的几何和物理假设这是所有计算的出发点形状与产状地下物体是一个无限长的水平圆柱体。注意“无限长”这个假设它意味着我们只考虑垂直于圆柱轴线方向的横截面重力异常在轴线方向上是没有变化的。这极大地简化了问题将三维问题降维成了二维问题。物性参数圆柱体与围岩之间存在一个恒定的密度差 Δρ。设圆柱体密度为 ρ_c围岩密度为 ρ_s则 Δρ ρ_c - ρ_s。Δρ 可正可负正值对应高密度体如铁矿引起正异常负值对应低密度体如盐丘、空洞引起负异常。空间位置圆柱体轴线平行于地面埋藏深度为 d指轴线到地面的垂直距离圆柱体的半径为 R。观测方式我们在地表沿一条垂直于圆柱体轴线的测线进行观测。测点坐标为 (x, 0)其中 x 是距测线中心点通常投影在圆柱体中心正上方的水平距离。把这些假设在脑子里或者纸上画出来形成一个清晰的二维剖面图地面是一条水平线y0地下深度 d 处有一个半径为 R 的圆圆的中心即圆柱体的轴线。2.2 重力异常公式推导理解性回顾对于横截面为任意形状的二度体无限长柱体其重力异常 Δg 可以通过计算截面面积 S 与密度差 Δρ 的乘积再结合距离积分来求得。对于圆形的横截面水平圆柱体这个积分有优美的解析解。设测点坐标为 (x, 0)圆柱体中心坐标为 (0, d)半径为 R密度差为 Δρ万有引力常数为 G≈ 6.674×10⁻¹¹ m³/kg/s²。则在该测点产生的重力异常垂直分量通常我们测量和计算的就是这个垂直分量为Δg(x) 2πG Δρ R² * d / (x² d²)这个公式是核心。我们来拆解一下它的物理意义2πG Δρ R²这一部分可以看作是一个“强度因子”。πR²是圆柱体的横截面积Δρ是密度差G是引力常数。它们的乘积再乘以 2π整体反映了异常源本身的“引力强度”。d / (x² d²)这一部分是“几何衰减因子”。它描述了异常强度随着观测点与源体相对位置水平距离 x 和埋深 d的变化而衰减的规律。当测点正好在圆柱体正上方 (x0) 时异常取得最大值Δg_max 2πG Δρ (R²/d)。看最大异常值与半径的平方成正比与埋深成反比。一个埋藏更浅或半径更大的圆柱体其异常峰值会更明显。随着 |x| 增大异常值对称地减小。当 |x| d 时异常衰减近似与 x² 成反比。注意公式中使用的单位是国际单位制SI。在实际地球物理勘探中重力异常通常非常小常用单位是毫伽mGal1 mGal 10⁻⁵ m/s²。计算时G 取 6.674×10⁻¹¹密度单位用 kg/m³长度单位用米m计算出的 Δg 单位是 m/s²乘以 10⁵ 即得到 mGal。为了直观我们编程时可以先按 SI 单位计算最后统一转换。2.3 模型参数的影响分析理解每个参数如何影响异常曲线形态对于后续的反演和解释至关重要。我们可以做一下“思想实验”密度差 Δρ线性缩放因子。Δρ 增大一倍整个异常曲线幅度就增大一倍。它影响异常的“振幅”。半径 R与异常幅值成平方关系影响显著。R 也略微影响异常的“宽度”半径越大异常曲线越宽缓。埋深 d是最关键的参数之一。它同时影响异常的“振幅”和“宽度”。埋深增加峰值异常减小反比关系同时异常曲线变得更加宽缓因为衰减因子分母中的 d 增大了。一个深部的大物体和一个浅部的小物体可能产生幅值相近但形态不同的异常。水平位置公式中的 x 是以圆柱体中心在地表投影为原点的。模型默认是对称的异常曲线关于 x0 对称。3. MATLAB 实现从公式到图形理论清晰了接下来就是用 MATLAB 把它实现出来。我们的目标是输入一组模型参数Δρ, R, d计算并绘制出一条测线上的重力异常曲线。3.1 环境准备与参数定义首先我们定义模型的基本参数。这里我建议用一个结构体model来管理所有参数这样代码更清晰也便于后续进行参数化研究。% 定义水平圆柱体模型参数 model struct(); model.delta_rho 500; % 密度差单位kg/m^3 (例如砂岩与灰岩的密度差) model.R 50; % 圆柱体半径单位米 model.d 100; % 圆柱体中心埋深单位米 model.G 6.67430e-11; % 万有引力常数单位m^3 kg^-1 s^-2 % 定义观测测线 x_min -300; % 测线最小x坐标单位米 x_max 300; % 测线最大x坐标单位米 num_points 601; % 测点数量建议为奇数便于中心对称 x_profile linspace(x_min, x_max, num_points); % 生成等间距测点坐标参数选择心得delta_rho常见岩石密度差在几百到一千多 kg/m³ 之间。500 是一个中等偏下的值计算出的异常大小比较适中便于观察。R和d通常埋深d会大于半径R。这里设d100m,R50m即埋深是半径的2倍是一个比较典型的场景。如果R接近或大于d异常会非常尖锐接近“未完全埋藏”的状态。x_profile测线范围要足够覆盖异常区域。一个经验法则是测线半宽至少取3*d到5*d以确保能捕捉到异常衰减到接近背景值的部分。这里从 -300m 到 300m是埋深的3倍是合理的。测点数量要足够多曲线才会光滑。3.2 核心正演计算函数我们将正演计算封装成一个函数这是代码的核心模块。function delta_g forward_gravity_horizontal_cylinder(x, delta_rho, R, d, G) % 计算水平圆柱体在各测点引起的重力异常 % 输入 % x : 测点水平坐标向量 (米) % delta_rho : 密度差 (kg/m^3) % R : 圆柱体半径 (米) % d : 圆柱体中心埋深 (米) % G : 万有引力常数 % 输出 % delta_g : 重力异常向量 (m/s^2) % 使用解析公式直接计算 % 公式: Δg(x) 2 * π * G * Δρ * R^2 * d / (x.^2 d^2) delta_g 2 * pi * G * delta_rho * R^2 * d ./ (x.^2 d.^2); % 注意这里使用了点除 (./) 和点幂 (.^)以便对向量x进行逐元素计算。 end这个函数极其简洁就是公式的直接翻译。在命令行或脚本中调用它% 计算重力异常 g_anomaly_si forward_gravity_horizontal_cylinder(x_profile, ... model.delta_rho, ... model.R, ... model.d, ... model.G); % 将单位从 m/s^2 转换为更常用的毫伽 (mGal) % 1 m/s^2 100,000 mGal (即 10^5 mGal) g_anomaly_mgal g_anomaly_si * 1e5;代码细节与陷阱向量化运算x是一个向量公式中的x.^2和除法./必须使用点运算符这样才能对x的每个元素独立计算最终输出一个同长度的异常向量delta_g。如果误写成/(x^2 d^2)MATLAB 会尝试矩阵运算导致错误或结果不对。单位换算地球表面的重力加速度约为 9.8 m/s²而一个地质体引起的异常可能只有其百万分之一即微伽量级。所以用 SI 单位计算出来的值会非常小例如 1e-6 量级。转换为毫伽mGal后数值更直观例如几个到几百 mGal。记住换算关系1 mGal 10⁻⁵ m/s²所以乘以 10⁵ 即可。3.3 结果可视化与初步分析计算出数据后可视化是关键。一张好的图能传达大量信息。% 创建图形窗口 figure(Position, [100, 100, 900, 600]); % 设置图形位置和大小 % 子图1重力异常剖面曲线 subplot(2, 2, [1, 3]); % 占据左半部分 plot(x_profile, g_anomaly_mgal, b-, LineWidth, 2); grid on; xlabel(测点水平位置 x (m)); ylabel(重力异常 \Deltag (mGal)); title(水平圆柱体重力异常剖面曲线); % 标记最大值点 [max_val, max_idx] max(g_anomaly_mgal); hold on; plot(x_profile(max_idx), max_val, ro, MarkerSize, 10, MarkerFaceColor, r); text(x_profile(max_idx), max_val*1.05, sprintf(最大值: %.2f mGal, max_val), ... HorizontalAlignment, center); % 添加模型参数标注 param_text sprintf(\\Delta\\rho %d kg/m^3\nR %d m\nd %d m, ... model.delta_rho, model.R, model.d); text(0.05, 0.95, param_text, Units, normalized, ... VerticalAlignment, top, BackgroundColor, w, EdgeColor, k); % 子图2地下模型示意图 subplot(2, 2, 2); % 绘制地面线 plot([x_min, x_max], [0, 0], k-, LineWidth, 2); hold on; % 绘制圆柱体横截面圆 theta linspace(0, 2*pi, 100); circle_x model.R * cos(theta); circle_y model.d model.R * sin(theta); % 注意MATLAB图形y轴向下为正这里d是正值 fill(circle_x, circle_y, [0.8, 0.8, 1], EdgeColor, b, LineWidth, 1.5); % 浅蓝色填充 % 标记圆心轴线 plot(0, model.d, k, MarkerSize, 12, LineWidth, 2); % 标注 xlabel(水平距离 (m)); ylabel(深度 (m)); title(地下模型示意图 (横截面)); axis equal; grid on; % 设置y轴方向使深度向下为正 set(gca, YDir, reverse); ylim([0, model.d model.R 20]); % 添加标注线 annotation(arrow, [0.5, 0.5], [0.6, 0.75], String, 埋深 d); annotation(arrow, [0.5, 0.55], [0.5, 0.5], String, 半径 R); % 子图3异常等值线图二维平面图 subplot(2, 2, 4); % 假设圆柱体沿y方向无限延伸我们计算x-y平面上的异常y是沿走向方向 [y_grid, x_grid] meshgrid(linspace(-150, 150, 60), linspace(x_min, x_max, 80)); % 计算网格上每点的异常此时公式中距离应为 sqrt(x^2 d^2)因为y方向无变化 % 注意对于二度体在垂直于走向的剖面上异常不随y变化。但为了画平面图我们假设在y方向有限范围内观测。 % 更严谨的二维平面图需要计算全空间重力位这里为简化展示剖面曲线在y方向的“拉伸”。 g_2d forward_gravity_horizontal_cylinder(x_grid, model.delta_rho, model.R, model.d, model.G) * 1e5; contourf(x_grid, y_grid, g_2d, 20, LineStyle, none); colorbar; xlabel(x (m)); ylabel(y (沿走向m)); title(重力异常平面等值线图 (示意)); axis equal tight;绘图技巧与解读多子图布局使用subplot将剖面曲线、模型示意图和平面图放在一起信息呈现非常完整。模型图y轴反转在地球物理和地质剖面中深度向下为正。使用set(gca, YDir, reverse)实现这一点更符合专业习惯。异常曲线特征生成的剖面曲线应该是一条关于 x0 对称的、光滑的“钟形”曲线。峰值位于圆柱体中心正上方。曲线宽度与埋深d密切相关。等值线图这里的等值线图是一个示意。对于真正的二度体其重力异常在走向y方向是无限延伸且不变的所以等值线图应该是一组平行直线。我们这里计算了一个小范围的y网格只是为了视觉上展示一个“平面分布”的概念。在实际中对于有限长度的三维物体等值线图会是封闭的椭圆形。4. 参数化研究与模型影响探究仅仅计算一个模型是不够的。我们需要通过改变参数系统地观察异常曲线如何响应这能培养我们的“地质直觉”。4.1 设计对比实验我们将分别改变密度差、半径和埋深观察异常曲线的变化。% 基础参数 base_delta_rho 500; % kg/m^3 base_R 50; % m base_d 100; % m x linspace(-300, 300, 601); % 创建新图形 figure(Position, [100, 100, 1200, 800]); % 实验1改变密度差 Δρ subplot(2, 3, 1); delta_rho_values [200, 500, 800]; % kg/m^3 colors lines(length(delta_rho_values)); % 获取区分度好的颜色 for i 1:length(delta_rho_values) g forward_gravity_horizontal_cylinder(x, delta_rho_values(i), base_R, base_d, model.G) * 1e5; plot(x, g, -, Color, colors(i, :), LineWidth, 2, ... DisplayName, sprintf(\\Delta\\rho %d, delta_rho_values(i))); hold on; end grid on; xlabel(x (m)); ylabel(\Deltag (mGal)); title((a) 不同密度差的影响); legend(show, Location, best); % 实验2改变半径 R subplot(2, 3, 2); R_values [30, 50, 70]; % m for i 1:length(R_values) g forward_gravity_horizontal_cylinder(x, base_delta_rho, R_values(i), base_d, model.G) * 1e5; plot(x, g, -, Color, colors(i, :), LineWidth, 2, ... DisplayName, sprintf(R %d m, R_values(i))); hold on; end grid on; xlabel(x (m)); ylabel(\Deltag (mGal)); title((b) 不同半径的影响); legend(show, Location, best); % 实验3改变埋深 d subplot(2, 3, 3); d_values [80, 100, 120]; % m for i 1:length(d_values) g forward_gravity_horizontal_cylinder(x, base_delta_rho, base_R, d_values(i), model.G) * 1e5; plot(x, g, -, Color, colors(i, :), LineWidth, 2, ... DisplayName, sprintf(d %d m, d_values(i))); hold on; end grid on; xlabel(x (m)); ylabel(\Deltag (mGal)); title((c) 不同埋深的影响); legend(show, Location, best); % 实验4综合对比 - 峰值异常与参数关系理论值 subplot(2, 3, 4); % 理论峰值公式Δg_max 2πG Δρ R^2 / d delta_rho_range 200:100:800; peak_vs_drho 2*pi*model.G * delta_rho_range * base_R^2 / base_d * 1e5; plot(delta_rho_range, peak_vs_drho, o-, LineWidth, 2); grid on; xlabel(密度差 \Delta\rho (kg/m^3)); ylabel(峰值异常 \Deltag_{max} (mGal)); title((d) 峰值异常 vs. 密度差 (线性)); subplot(2, 3, 5); R_range 20:10:80; peak_vs_R 2*pi*model.G * base_delta_rho * R_range.^2 / base_d * 1e5; plot(R_range, peak_vs_R, s-, LineWidth, 2); grid on; xlabel(半径 R (m)); ylabel(峰值异常 \Deltag_{max} (mGal)); title((e) 峰值异常 vs. 半径 (平方关系)); subplot(2, 3, 6); d_range 60:10:140; peak_vs_d 2*pi*model.G * base_delta_rho * base_R^2 ./ d_range * 1e5; plot(d_range, peak_vs_d, ^-, LineWidth, 2); grid on; xlabel(埋深 d (m)); ylabel(峰值异常 \Deltag_{max} (mGal)); title((f) 峰值异常 vs. 埋深 (反比关系));4.2 实验结果分析与地质解释运行上述代码后我们可以得到六张图前三张是异常曲线形态对比后三张是峰值异常与各参数的定量关系。(a) 不同密度差的影响三条曲线形态完全一致只是振幅按比例缩放。密度差从200增加到800 kg/m³异常峰值也几乎按相同比例约4倍增加。这验证了 Δρ 是一个线性缩放因子。地质意义在野外如果我们看到两个形态相似但幅值不同的异常可能意味着相似的地质体具有不同的密度差。(b) 不同半径的影响半径增大异常峰值显著增加注意是平方关系同时异常曲线也略微变宽。R70m 的曲线比 R30m 的曲线不仅峰值高很多而且“山脚”也更宽。地质意义异常幅值和宽度同时增大通常指示着异常源体积更大。(c) 不同埋深的影响这是最有意思的。埋深增加异常峰值急剧减小d120m 的峰值约为 d80m 的 (80/120)≈0.67倍同时异常曲线变得更加宽缓。d80m 的曲线又高又瘦d120m 的曲线又矮又胖。这是重力勘探中一个非常重要的现象称为“等效原理”一个埋深大、体积大的地质体其产生的异常可能与一个埋深浅、体积小的地质体异常形态相似。这给反演解释带来了多解性。(d, e, f) 定量关系图清晰地展示了理论公式揭示的关系峰值异常与 Δρ 成正比与 R² 成正比与 d 成反比。这些图是连接模型参数与观测数据的桥梁。实操心得做参数化研究时一次只改变一个参数其他参数保持不变这是控制变量法的基本思想。画图时使用不同的线型和颜色并添加清晰的图例能让对比结果一目了然。把这些图保存下来就是一份非常好的学习笔记或报告素材。5. 高级应用与扩展思考掌握了基础正演后我们可以尝试一些更贴近实际应用的扩展。5.1 叠加异常与复杂模型近似真实地下往往不止一个地质体。多个水平圆柱体的异常可以通过线性叠加来计算。这就是“复杂模型可以由简单模型组合”的思想。% 定义两个水平圆柱体模型 model1.delta_rho 600; model1.R 40; model1.d 80; model1.x_center -50; model2.delta_rho -300; model2.R 30; model2.d 120; model2.x_center 70; % 计算测线覆盖两个物体 x linspace(-200, 200, 401); % 计算单个异常 g1 forward_gravity_horizontal_cylinder(x - model1.x_center, ... % 注意坐标平移 model1.delta_rho, model1.R, model1.d, model.G); g2 forward_gravity_horizontal_cylinder(x - model2.x_center, ... model2.delta_rho, model2.R, model2.d, model.G); % 叠加总异常 g_total g1 g2; g_total_mgal g_total * 1e5; % 绘图 figure; plot(x, g1*1e5, b--, LineWidth, 1.5, DisplayName, 高密度体 (正异常)); hold on; plot(x, g2*1e5, r--, LineWidth, 1.5, DisplayName, 低密度体 (负异常)); plot(x, g_total_mgal, k-, LineWidth, 2.5, DisplayName, 叠加总异常); grid on; xlabel(测点位置 x (m)); ylabel(重力异常 \Deltag (mGal)); title(多个水平圆柱体重力异常叠加); legend(show, Location, best);结果分析你会看到总异常曲线不再是简单的钟形。它可能有两个峰值或者一个正异常旁边伴随一个负的“尾巴”形态变得复杂。这模拟了真实地下多个地质体共存的情况。反演解释时就需要设法将这样的复合异常分解成多个简单异常源。5.2 加入观测噪声与反演概念引入野外实测数据永远包含噪声。为了模拟更真实的数据我们可以给理论异常添加随机噪声。% 生成理论异常 g_theory forward_gravity_horizontal_cylinder(x_profile, model.delta_rho, model.R, model.d, model.G) * 1e5; % 添加高斯白噪声假设噪声水平为峰值异常的2% noise_level 0.02 * max(abs(g_theory)); g_noisy g_theory noise_level * randn(size(g_theory)); % randn生成标准正态分布噪声 % 绘图对比 figure; plot(x_profile, g_theory, b-, LineWidth, 2, DisplayName, 理论异常); hold on; plot(x_profile, g_noisy, r., MarkerSize, 8, DisplayName, 含噪声“观测”数据); grid on; xlabel(测点位置 x (m)); ylabel(重力异常 \Deltag (mGal)); title(理论异常与含噪声数据对比); legend(show);意义添加噪声后光滑的钟形曲线变成了上下波动的散点。这引出了地球物理反演的核心挑战如何从带有噪声的、有限的观测数据中稳定地估计出地下的模型参数Δρ, R, d, x_center这就需要进行反演。最简单的反演思路可能是最小二乘法寻找一组模型参数使得其正演结果与观测数据之间的误差平方和最小。你可以尝试用fminsearch或lsqnonlin这样的优化函数来实现一个简单的反演这将是这个项目极好的延伸。5.3 从二度体到三度体球体模型水平圆柱体是二度体而很多地质体如矿囊、溶洞更接近三度体。最简单的三度体模型是球体。球体重力异常的公式是Δg(x) (4/3)πG Δρ R³ * d / (x² d²)^(3/2)你可以仿照本文的流程用 MATLAB 实现球体的正演。对比球体和圆柱体的异常曲线你会发现球体的异常衰减得更快分母是3/2次方曲线更“瘦高”。这是区分物体延展度二度还是三度的重要标志。6. 常见问题、调试技巧与避坑指南在实现和调试过程中你肯定会遇到一些问题。这里总结一些常见坑点。6.1 公式输入错误这是最常犯的错误。请逐字检查公式检查π是pi。检查G的值是否正确6.67430e-11。检查指数和除法运算符R^2和./ (x.^2 d^2)。务必使用点运算符。检查括号匹配。调试技巧先计算一个点的值比如 x0。此时公式简化为Δg(0) 2πG Δρ R² / d。手动用计算器算一下再与程序输出对比。这是快速验证公式编码是否正确的好方法。6.2 单位混乱导致数量级离谱症状计算出的异常值要么极大如几万 mGal要么极小如10^-10 mGal。检查密度单位岩石密度通常是2-3 g/cm³即 2000-3000 kg/m³。密度差通常是几百 kg/m³。如果你误用了 g/cm³即 500 g/cm³结果会大1000倍。检查长度单位公式默认是米m。如果你的埋深 d 以为是米实际数据是公里km忘了换算结果会差1000倍。牢记最终单位公式直接算出的是 m/s²。乘以 10^5 得到 mGal。如果你期望看到的是几十个 mGal 的量级而程序输出是 0.000几那很可能忘了转换单位。6.3 图形显示异常曲线是一条水平线检查计算函数的输入参数是否传对了。特别是x向量是否正确生成delta_rho是否为正数。曲线形状奇怪不对称检查公式中分母是不是(x.^2 d^2)确保是x的平方。如果写成(x d^2)曲线就不对称了。模型示意图中物体位置不对记住 MATLAB 绘图坐标原点在左下角y轴向上为正。而地质剖面深度向下为正。使用set(gca, YDir, reverse)来翻转y轴。计算圆的y坐标时是d R*sin(theta)因为d是正值深度。6.4 性能与向量化我们的计算很简单向量化后效率很高。但如果未来做更复杂的计算如三度体积分或者对非常大的网格进行计算可能会遇到性能问题。预分配数组在循环前用zeros()或ones()函数预先分配好存储结果数组的空间这能显著提升循环速度。利用矩阵运算尽量避免多层嵌套循环。MATLAB 擅长矩阵运算思考能否将问题转化为矩阵乘法或数组运算。匿名函数对于简单的正演公式可以定义为匿名函数使代码更简洁calc_g (x, drho, R, d) 2*pi*G*drho*R^2*d./(x.^2d^2);6.5 从正演到反演的思维转变正演是“给定模型求响应”。反演是“给定响应求模型”。做完正演一定要思考反问题多解性如前所述不同参数组合可能产生相似的异常曲线。噪声影响噪声会掩盖异常细节使反演结果不稳定。约束的重要性在实际反演中必须加入先验地质信息作为约束如密度差的范围、埋深不可能为负等才能得到地质上合理的解。自己动手实现这个水平圆柱体的正演就像是拿到了地球物理勘探的一把“钥匙”。它虽然简单但蕴含了重力方法最核心的思想通过地表观测的微小引力变化去推测地下不可见世界的奥秘。当你看到自己写出的几行代码成功绘制出那条优美的、符合物理规律的异常曲线时那种将理论付诸实践的感觉是单纯看书无法比拟的。接下来你可以尝试球体模型可以尝试叠加多个物体甚至可以挑战一下最基础的网格化反演。每一步的扩展都会让你对地球物理数据的处理和解释有更深的理解。