MATLAB实现水平圆柱体重力异常正演：从原理到地质解释-开发者社区

1. 项目概述：从“异常”中寻找规律

搞地球物理勘探或者地质工程的朋友，对“重力异常”这个词肯定不陌生。简单来说，我们脚下的大地，其密度分布并不是均匀的。一个埋藏在地下的矿体、一个地质构造，甚至一个古代遗存，都会因为其密度与周围岩土存在差异，导致其所在位置的地球重力场产生微小的、局部的变化。这个变化，就是重力异常。我们的工作，就是在地面上布设测点，用高精度的重力仪测量这些微小的变化，然后像解谜一样，反推地下到底是什么东西引起了这个异常。

“正演”，就是这个解谜过程的“前半部分”，或者说，是构建谜题的过程。假设我们已经知道了地下物体的形状、大小、埋深、密度差，那么理论上，它在地表会产生一个什么样的重力异常响应？计算出这个理论响应值，就是正演。它是反演（即从实测异常数据推断地下物体属性）的基础和前提。只有正演模型足够准确、计算足够高效，我们后续的反演解释才靠谱。

这次要聊的，就是用 MATLAB 来模拟一个非常经典的地球物理模型：水平圆柱体的重力异常正演。为什么是水平圆柱体？因为它形状规则，其重力异常有解析解（也就是可以用一个数学公式精确计算出来），是验证算法、理解物理概念、乃至模拟一些近似柱状的地质体（如盐丘、矿脉、隧道、管线）的绝佳起点。很多复杂的模型，初期都可以用一系列不同参数的圆柱体来近似。

所以，这个项目的核心价值在于：掌握用数值计算工具（MATLAB）实现地球物理基本正演模型的能力，为后续更复杂的反演解释和实际数据处理打下坚实的基础。无论你是地质、测绘、地球物理专业的学生，还是相关领域的工程师，自己动手实现一遍，远比只看教科书上的公式印象深得多。

2. 核心原理与模型解析：公式背后的物理图景

在动手写代码之前，我们必须把模型和公式吃透。一知半解地套公式，出了错都不知道在哪。

2.1 水平圆柱体模型假设

我们首先明确模型的几何和物理假设，这是所有计算的出发点：

形状与产状：地下物体是一个无限长的水平圆柱体。注意“无限长”这个假设，它意味着我们只考虑垂直于圆柱轴线方向的横截面，重力异常在轴线方向上是没有变化的。这极大地简化了问题，将三维问题降维成了二维问题。
物性参数：圆柱体与围岩之间存在一个恒定的密度差 Δρ。设圆柱体密度为 ρ_c，围岩密度为 ρ_s，则 Δρ = ρ_c - ρ_s。Δρ 可正可负，正值对应高密度体（如铁矿），引起正异常；负值对应低密度体（如盐丘、空洞），引起负异常。
空间位置：圆柱体轴线平行于地面，埋藏深度为 d（指轴线到地面的垂直距离），圆柱体的半径为 R。
观测方式：我们在地表沿一条垂直于圆柱体轴线的测线进行观测。测点坐标为 (x, 0)，其中 x 是距测线中心点（通常投影在圆柱体中心正上方）的水平距离。

把这些假设在脑子里或者纸上画出来，形成一个清晰的二维剖面图：地面是一条水平线（y=0），地下深度 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）的变化而衰减的规律。
- 当测点正好在圆柱体正上方 (x=0) 时，异常取得最大值：Δg_max = 2πG Δρ (R²/d)。看，最大异常值与半径的平方成正比，与埋深成反比。一个埋藏更浅或半径更大的圆柱体，其异常峰值会更明显。
- 随着 |x| 增大，异常值对称地减小。当 |x| >> d 时，异常衰减近似与 x² 成反比。

注意：公式中使用的单位是国际单位制（SI）。在实际地球物理勘探中，重力异常通常非常小，常用单位是毫伽（mGal），1 mGal = 10⁻⁵ m/s²。计算时，G 取 6.674×10⁻¹¹，密度单位用 kg/m³，长度单位用米（m），计算出的 Δg 单位是 m/s²，乘以 10⁵ 即得到 mGal。为了直观，我们编程时可以先按 SI 单位计算，最后统一转换。

2.3 模型参数的影响分析

理解每个参数如何影响异常曲线形态，对于后续的反演和解释至关重要。我们可以做一下“思想实验”：

密度差 Δρ：线性缩放因子。Δρ 增大一倍，整个异常曲线幅度就增大一倍。它影响异常的“振幅”。
半径 R：与异常幅值成平方关系，影响显著。R 也略微影响异常的“宽度”，半径越大，异常曲线越宽缓。
埋深 d：是最关键的参数之一。它同时影响异常的“振幅”和“宽度”。埋深增加，峰值异常减小（反比关系），同时异常曲线变得更加宽缓（因为衰减因子分母中的 d 增大了）。一个深部的大物体和一个浅部的小物体，可能产生幅值相近但形态不同的异常。
水平位置：公式中的 x 是以圆柱体中心在地表投影为原点的。模型默认是对称的，异常曲线关于 x=0 对称。

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。这里设d=100m,R=50m，即埋深是半径的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')实现这一点，更符合专业习惯。
异常曲线特征：生成的剖面曲线应该是一条关于 x=0 对称的、光滑的“钟形”曲线。峰值位于圆柱体中心正上方。曲线宽度与埋深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) 不同半径的影响：半径增大，异常峰值显著增加（注意是平方关系），同时异常曲线也略微变宽。R=70m 的曲线比 R=30m 的曲线不仅峰值高很多，而且“山脚”也更宽。地质意义：异常幅值和宽度同时增大，通常指示着异常源体积更大。

(c) 不同埋深的影响：这是最有意思的。埋深增加，异常峰值急剧减小（d=120m 的峰值约为 d=80m 的 (80/120)≈0.67倍），同时异常曲线变得更加宽缓。d=80m 的曲线又高又瘦，d=120m 的曲线又矮又胖。这是重力勘探中一个非常重要的现象，称为“等效原理”：一个埋深大、体积大的地质体，其产生的异常可能与一个埋深浅、体积小的地质体异常形态相似。这给反演解释带来了多解性。

(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)。务必使用点运算符。
检查括号匹配。

调试技巧：先计算一个点的值，比如 x=0。此时公式简化为Δ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.^2+d^2);

6.5 从正演到反演的思维转变

正演是“给定模型，求响应”。反演是“给定响应，求模型”。做完正演，一定要思考反问题：

多解性：如前所述，不同参数组合可能产生相似的异常曲线。
噪声影响：噪声会掩盖异常细节，使反演结果不稳定。
约束的重要性：在实际反演中，必须加入先验地质信息作为约束（如密度差的范围、埋深不可能为负等），才能得到地质上合理的解。

自己动手实现这个水平圆柱体的正演，就像是拿到了地球物理勘探的一把“钥匙”。它虽然简单，但蕴含了重力方法最核心的思想：通过地表观测的微小引力变化，去推测地下不可见世界的奥秘。当你看到自己写出的几行代码，成功绘制出那条优美的、符合物理规律的异常曲线时，那种将理论付诸实践的感觉，是单纯看书无法比拟的。接下来，你可以尝试球体模型，可以尝试叠加多个物体，甚至可以挑战一下最基础的网格化反演。每一步的扩展，都会让你对地球物理数据的处理和解释有更深的理解。