企业官网建设流程全解析

简介：这套MATLAB工具箱面向地质勘探与地球物理研究人员，用于对磁测和重力观测数据进行Euler反褶积处理，以推断地下异常体的位置、形态与规模。压缩包共9个文件，内含4个M脚本和5个MAT数据文件；M脚本分别承担参数设置、自动反褶积计算及多种地质模型构建，MAT文件则提供实测磁场分量和重力异常的样例数据，便于直接运行验证。包体仅约619KB，轻量实用。已有195人学习下载。使用者可获得一套可运行的Euler反褶积流程，既能通过样例数据快速上手，也可替换自备测量资料开展实际资料解释，适合有一定地球物理或MATLAB基础、希望掌握磁法/重力数据处理方法的研究者。

1. 磁测数据反演总被质疑“不唯一”？这套 MATLAB 的 Euler 反演代码给出了可复现的答案

做磁法数据处理的人都有一个共同的痛点：磁测数据反演结果解释起来太“软”了。同样一条剖面上的磁异常，有人说是岩墙，有人说是接触带，争论到最后往往靠经验拍板。Euler deconvolution（欧拉反褶积）是解决这个问题的经典方法，它不依赖初始模型，直接通过场源几何参数和异常衰减关系估算场源位置和深度，结果客观、可复核。我拆过不少 MATLAB 的磁法处理代码，这套 v29-08-02 版本的磁性体 Euler 反演工具包，代码结构紧凑，同时覆盖 magnetic 和 gravity 两套数据输入，对做重磁联合解释的人来说是个能直接用、也值得拆开改的底子。

这套资源的核心价值在于：它把 Euler 反演从“原理会说、代码不会写”的尴尬中解放出来，提供了窗口扫描、结构指数设定、深度估计、误差统计的完整流程。适合地球物理专业的研究生、物探队的数据处理工程师，以及想在自己项目里快速接入反演算法的 MATLAB 用户。接下来的内容，我会先说清楚这套代码在做什么、每个参数到底在控制什么，然后带你把数据准备好、把参数调对、把结果看懂，最后把我在实践中踩过的坑原原本本摆出来。

2. 把 Euler 反演从公式变成代码：结构指数、窗口扫描和四类场源模型

2.1 为什么选 Euler 反褶积而不是直接做反演拟合

磁异常反演的最大难题是“非唯一性”：不同的场源组合可能产生几乎一样的异常。传统做法是给定一个初始模型，然后用最优化方法去拟合观测数据，但初始模型本身就有主观性，收敛结果容易陷入局部极小。Euler 反褶积的思路完全不同——它不反演磁化强度，不猜场源形状，而是利用位场满足的齐次方程，直接从异常值和它的水平/垂向导数求解场源位置。

Euler 齐次方程的表达式很简洁：

(x - x0) * ∂T/∂x + (y - y0) * ∂T/∂y + (z - z0) * ∂T/∂z = -N * (T - B)

其中 T 是观测到的磁异常或重力异常，B 是背景场，N 是结构指数（Structural Index, SI），(x0, y0, z0) 就是要解的场源位置。方程的物理含义是：不同几何形状的场源，其异常随距离的衰减率不同，N 就是这个衰减率的表现。这套代码的核心功能，就是对网格数据开一个固定大小的窗口，在每个窗口内用最小二乘法解这个方程，得到一组 (x0, y0, z0)，然后把所有窗口的解汇总成深度估计图。

这里有个关键点：方程是线性的，未知数是 x0、y0、z0，但 N 是必须预先给定的。N 给错了，反演结果就是系统性偏差——这一点后面避坑章会重点展开。这套代码把 N 作为输入参数暴露给用户，而不是写死，这在工程上非常重要。

2.2 代码包的目录结构和主函数入口

拿到代码包之后，第一件事不是急着跑，而是先看懂文件组织。这套资源的文件结构大致如下：

v29-08-02_matlabmagnetic/ ├── euler_deconv_main.m # 主脚本：读数据 → 算导数 → 扫窗口 → 输出结果 ├── euler_deconv_func.m # 核心函数：单个窗口内的 Euler 方程求解 ├── derivatives.m # 水平/垂向导数计算（频率域算子） ├── window_scan.m # 窗口扫描与深度解汇总 ├── mag_data_demo.dat # 磁异常演示数据 ├── grav_data_demo.dat # 重力异常演示数据 └── README.txt # 参数说明与运行指引

主程序 euler_deconv_main.m 的流程是线性的，适合逐步断点调试。我建议你先用演示数据把流程跑通，再换成自己的网格数据。第一次跑的时候把 MATLAB 当前目录切到代码包根目录，然后在命令行执行：

% 加载主脚本，注意当前路径必须在代码包根目录 run('euler_deconv_main.m');

如果你用的是 MATLAB R2019b 之后的版本，run 命令没有问题；如果是更老的版本，建议直接用编辑器打开主脚本，按 Ctrl+Enter 分节执行。执行完成后，工作区应该出现 depth_solution、si_map、coverage_map 三个关键变量。它们分别代表深度解矩阵、每个网格点的结构指数设定值、以及有效解覆盖次数。

2.3 频率域导数计算：Euler 方程的输入质量决定一切

Euler 反褶积对导数的质量极其敏感。导数算不准，后面解出来的位置和深度全是噪声。这套代码的 derivatives.m 用的是频率域算子——先对异常做二维傅里叶变换，乘上导数因子，再做逆变换。这样做的好处是计算快、无方向偏差，坏处是边界效应明显。代码里已经做了边缘削边（taper），但我实测发现削边宽度默认只有 5 个网格点，对于 100×100 以上的网格来说偏小。

% derivatives.m 核心逻辑（节选） % 输入: anomaly_grid - 原始异常网格 % dx, dy - 网格间距（米） % 输出: dz - 垂向导数（磁异常通常取垂向分量） [ny, nx] = size(anomaly_grid); kx = (0:nx-1) * 2 * pi / (nx * dx); % 波数矢量 ky = (0:ny-1) * 2 * pi / (ny * dy); [KX, KY] = meshgrid(ifftshift(kx), ifftshift(ky)); K = sqrt(KX.^2 + KY.^2); % 径向波数 % 垂向导子：|k|，注意高频放大效应 DZ = K; ft_grid = fft2(anomaly_grid); dz_grid = real(ifft2(ft_grid .* DZ));

参数说明：dx 和 dy 必须是实际物理间距，单位要和坐标一致，比如都是米。很多新手把网格点距当成 1 来处理，结果反演深度全部偏小——因为导数幅值和间距的平方成反比。代码里用的是 2π 周期的波数定义，如果你原来的数据是按度（经纬度）存储的，需要先做投影转换，否则波数单位对不上。

另外一个容易忽略的点：频率域导数会把高频噪声放大。代码在 DZ 上做了一个简单的 cosine 衰减窗，但我建议你在这个步骤之前对原始异常做一次低通滤波。我一般会用向上延拓代替滤波——向上延拓 2 个点距后的数据再求导，比直接滤波的相位畸变更小。

2.4 窗口扫描逻辑：窗口大小与步长的配合是深度分辨率的命门

window_scan.m 是整个代码里最值得逐行读的部分。Euler 反演的窗口扫描策略是：用一个固定大小的矩形窗口，在网格上按一定步长滑动，每到一个位置，就把窗口内的所有网格点代入 Euler 方程，用最小二乘求解一个场源位置。窗口越小，空间分辨率越高，但每个窗口内的方程数量越少，抗噪能力越差；窗口越大则相反。

% window_scan.m 核心逻辑（节选） % 参数块：窗口尺寸和步长 win_size = 15; % 窗口边长（网格点数），推荐 7~20 step = 2; % 滑动步长（网格点数），通常取 win_size/5 ~ win_size/2 % 扫描主循环 for iy = 1:step:(ny - win_size) for ix = 1:step:(nx - win_size) % 取窗口内的异常值和导数值 T_win = anomaly_grid(iy:iy+win_size-1, ix:ix+win_size-1); Tx_win = dx_grid(iy:iy+win_size-1, ix:ix+win_size-1); Ty_win = dy_grid(iy:iy+win_size-1, ix:ix+win_size-1); Tz_win = dz_grid(iy:iy+win_size-1, ix:ix+win_size-1); % 构建 Euler 方程系数矩阵 A 和观测向量 b % A = [Tx, Ty, Tz]，b = -N*T（背景场 B 此处置零或另行估计） A = [Tx_win(:), Ty_win(:), Tz_win(:)]; b = -si * T_win(:); % 最小二乘解：x = (A'*A)\A'*b x_ls = (A'*A) \ (A'*b); % 记录解的位置（窗口中心 + 偏移解） sol_x(iy, ix) = ix + step/2 + x_ls(1)/dx; sol_y(iy, ix) = iy + step/2 + x_ls(2)/dy; sol_z(iy, ix) = x_ls(3); end end

这段代码的逻辑很清晰，但有一点需要特别说明：A'*A 的条件数决定了求解稳定性。当窗口内异常变化平缓时，Tx、Ty、Tz 的幅值都很小，A'*A 接近奇异，解出来的 (x0, y0, z0) 会飘到窗口外面甚至发散到几千公里。所以实际工程中必须加阻尼：

% 加阻尼最小二乘，防止矩阵奇异造成深度解发散 lambda = 0.001 * trace(A'*A) / 3; % 阻尼系数取迹的比例 x_ls = (A'*A + lambda * eye(3)) \ (A'*b);

这个阻尼项是我强烈建议你加上去的。原始代码里如果没有它，当窗口落在异常平缓区时，解的质量会非常差。你可以先不加阻尼跑一遍演示数据，再对比加阻尼的结果，差异会直观到让你吃惊。

3. 数据准备与参数配置：把一份实测网格数据喂给 Euler 反演的全流程

3.1 输入数据的格式要求与预处理清单

这套代码的输入是规则网格化的异常数据，ASCII 格式纯文本即可。常见格式是三列或四列：X 坐标、Y 坐标、异常值，以及可选的测点高程。我在实际使用中一般会把高程剔除，因为 Euler 反演默认把数据当水平面处理，如果地表起伏明显，而你又没有做地形校正，高程信息反而会干扰导数计算。

数据预处理的顺序很重要，顺序错了结果就废了。我惯用的流程是：

原始数据 → 网格化（最小曲率或克里金） → 去趋势（移除背景场） → 低通滤波 → 向上延拓一个点距 → 输入 Euler 反演

去趋势这一步尤其关键。Euler 方程里的 B（背景场）如果没去掉，反演得到的深度会出现系统性的偏向——异常梯度被背景场污染，方程右边的 N*T 和左边导数之间的均衡被破坏。代码的默认做法是把窗口均值当作背景场直接减掉，这个方法对区域性背景场有效，但如果背景场本身有线性梯度，窗口均值就镇不住了。我一般会在数据加载之后、反演之前，先做一次二阶多项式去趋势：

% 二阶多项式去趋势（推荐在进 Euler 前做一次） % 数据文件格式: X Y anomaly data = load('your_grid_data.txt'); x = data(:,1); y = data(:,2); t = data(:,3); % 构建设计矩阵：[1, x, y, x.^2, x.*y, y.^2] G = [ones(size(x)), x, y, x.^2, x.*y, y.^2]; m = (G'*G) \ (G'*t); % 多项式系数最小二乘解 trend = G * m; t_detrended = t - trend; % 去趋势后的异常 % 写回文件，供主脚本读取 out = [x, y, t_detrended]; save('your_grid_detrended.txt', '-ascii', 'out');

不管代码的 README 里怎么写的，把原始数据直接扔进 Euler 反演都是最不推荐的。你可以先用演示数据跑通，然后拿自己的数据做两次对比——一次直接跑，一次做去趋势后再跑，深度解的散差对比会让你意识到预处理不是可有可无的。

3.2 结构指数 SI 的选取原则与四种常见场源

结构指数是 Euler 反演里最重要的先验参数，它决定了方程右端的缩放因子。这套代码把 SI 作为一个变量放在主脚本里，你需要自己确定它的值。选错 SI 的后果不是解偏移一点，而是深度系统性失真——SI 偏大时深度被低估，SI 偏小时深度被高估，偏差比例可以达到 20%~40%。

下表是磁异常和重力异常对应的常见结构指数取值，我在实际项目中按这个表初选，然后用剖面拟合做交叉验证：

使用提示：SI 取整数值是教科书推荐，但实际数据处理中，SI 可以取小数。比如一个有限延伸的岩墙，它的磁异常衰减率介于薄板（SI=1）和圆柱（SI=2）之间，取 1.5 往往比取 1 或 2 都能得到更稳定的解。这套代码的 SI 是浮点变量，直接用小数即可，不需要改代码。

另一个经验：如果你在同一工区内有多个不同类型的场源（比如浅部有岩墙、深部有岩体），全局用单一 SI 会顾此失彼。处理办法是先做一次小窗口、SI=1 的反演，把浅部结果提取出来，再做一次大窗口、SI=2 的反演，把深部结果单独看。这套代码支持你多次运行主脚本改参数，所以这种分层反演策略不需要改任何代码——同一套数据跑两遍，设置不同的 SI 和窗口，然后把两套结果叠加到一张图上。

3.3 窗口大小与网格点距的匹配关系：两个示例参数组合

窗口大小和网格点距的关系，决定了反演结果能“看到”的最小场源规模。窗口的物理尺寸必须覆盖场源异常的主要波长范围。经验法则是：窗口边长至少要达到目标场源深度预估值的 2~3 倍。比如你预计场源深度在 50~100 米，网格点距是 10 米，那么窗口边长取 30~50 个点距比较合理。如果窗口太小，窗口内只有不到半个异常波峰，Euler 方程求解时信息量不够，解的位置会偏向窗口边缘。

我实际调试时常用两组参数组合，分别应对不同尺度的目标：

% 浅部小目标场景：窗口小、步长小、SI 取 1 win_size = 15; % 窗口边长 15 个点距 step = 3; % 步长 3 个点距，保证重叠率 si = 1.0; % 岩墙型场源 % 深部大目标场景：窗口大、步长放宽、SI 取 2 % 注意：窗口变大后，步长太小时相邻窗口解高度相关，不利于独立统计 win_size = 31; % 窗口边长 31 个点距 step = 8; % 步长 8 个点距 si = 2.0; % 圆柱型场源

这两组参数对应的深度分辨率差异很大。第一组能分辨横向距离 50 米左右的场源变化，第二组则更侧重深部趋势。把两组结果做差，你还能得到“深浅分离”的初步效果。不要指望一组参数打天下，这也是 Euler 反演实践里最常见的误区。

3.4 运行主脚本并导出深度解分布图

数据准备好、参数设好之后，运行主脚本，然后直接把 sol_z 矩阵用 imagesc 画出来，加上地形等高线叠加显示。这一步能快速判断反演是否合理：

% 反演结果可视化：深度解平面分布 figure('Color', 'w'); imagesc(sol_z); % sol_z 是窗口扫描得到的深度解矩阵 colorbar; caxis([0, 300]); % 根据实际深度范围调整色标上限 xlabel('X 网格点'); ylabel('Y 网格点'); title('Euler 反演深度解分布 (SI=1.5, win=15)'); % 叠加已知地质体的边界（如果有的话） hold on; contour(geology_boundary, 'k-', 'LineWidth', 1.2);

色标上限 caxis 很重要——深度解矩阵里经常有少数发散值（比如几千），如果不限制色标范围，一张图会被几个飞点拉成全红或全蓝，真实分布完全看不出。我一般会把深度解做一次 3σ 滤波再画图：超过 3σ 的直接置为 NaN，图上显示空白，宁可少一块区域可以看到的信息，也别被飞点污染视觉。

4. 磁法与重力数据联合反演：同一套框架处理 gravity 数据的差异点

4.1 重力数据走 Euler 反演时哪些环节要改

这套代码包同时带了 mag_data_demo.dat 和 grav_data_demo.dat 两组演示数据。重力异常和磁异常在物理本质上都是位场——都满足 Laplace 方程，所以 Euler 反演的数学框架完全通用。但切换数据源时，有三个环节必须调整，否则结果不可信：

第一是导数分量。磁异常通常用总磁场强度异常（TMI）或垂直分量，垂向导数对应 ∂T/∂z；重力异常则用布格重力异常，垂向导数是 ∂g/∂z（即重力梯度）。代码的 derivatives.m 对这二者没有本质区别，只要异常值单位正确即可。但要注意：磁异常的单位是 nT（纳特斯拉），重力异常的单位是 mGal（毫伽），数值量级差 3~6 个数量级——这不是问题，因为 Euler 方程里的 T 和导数项是线性齐次的，单位只会影响解的幅值尺度，不会影响位置和深度。

第二是结构指数。上面参数表里列了重力数据的 SI 对应关系：球体取 2，水平圆柱取 1，台阶取 0，接触带取 -1。注意 SI=-1 时方程右端为 -(-1)*T = +T，此时背景场估计的精度要求更高，建议在重力数据上优先做区域场剥离。

第三是背景场处理。磁异常的背景场往往可以通过 IGRF（国际地磁参考场）扣除，比较容易做干净；重力异常则受深部构造和区域密度变化的叠加影响，背景场经常是倾斜面而不是水平面。同样的去趋势代码，重力的多项式阶数建议比磁力高一阶——磁力用二阶，重力用三阶：

% 重力数据去趋势：用三阶多项式，更好地吸收区域场影响 G = [ones(size(x)), x, y, x.^2, x.*y, y.^2, x.^3, x.^2.*y, x.*y.^2, y.^3]; m = (G'*G) \ (G'*t); trend = G * m; t_detrended = t - trend;

4.2 重磁联合解的交叉验证策略

同时拥有磁和重力数据时，可以做联合约束：磁法 Euler 解和重力 Euler 解的深度分布如果在空间上吻合，场源位置的可信度就高；如果出现系统性的深浅差异，则说明某一种数据没有被处理好。最常见的偏差模式是：磁法结果浅、重力结果深。这通常意味着磁测数据里混入了浅层磁性干扰（比如地表铁磁性物质），或者重力数据没有完全剥离区域性深部场。

我处理重磁联合数据时的标准流程是：

先在磁异常上用 SI=1~2 扫一遍，标记浅部磁源集中的区块； 再在重力异常上用 SI=0~1 扫一遍，提取深部密度体分布； 把两者的解画在同一张剖面图上，按横向位置对齐； 如果同一地质体在两种数据上都出现，取两个深度的加权平均作为最终估计； 如果只在一种数据上出现，在报告中标注为“单方法解”，慎重采信。

这套策略可以在不修改代码的前提下实现——代码每次只处理一种数据，你只需要跑两遍，导出两个 sol_z 矩阵，然后在自己的脚本里做空间匹配。

4.3 单位与坐标系统的注意细节

磁法和重力数据的坐标投影不一致，是联合解释中最容易翻车的点。MATLAB 代码本身不关心坐标系，它只处理数值网格，但如果你把磁法结果（比如是经纬度网格）和重力结果（比如是 UTM 投影网格）直接叠加，差几公里是常有的事。建议在进入 MATLAB 之前，把所有数据统一投影到同一个坐标系（推荐 UTM），并在文件名后缀注明投影带。

另一个细节是网格间距的单位。Euler 方程里的深度解是“以网格间距为单位”的——代码里 x_ls(3)/dz 这种换算是把深度解从波数域拉回到物理距离。如果网格间距是 100 米，而你的坐标文件里用的是米，dx=100 没问题；但如果你的网格是 0.01 度的经纬网格，dx 是随纬度变化的——代码里 dx 是单个标量，这在高纬度地区会产生明显的水平畸变。处理办法：先把经纬网格投影到公里网格，再进 Euler 反演。千万别偷懒。

5. 避坑指南：Euler 反演的六个常见翻车现场与排查方法

5.1 深度解大面积发散，画出来全是红蓝杂点

现象：反演完成后，sol_z 矩阵里有大量值超过正常范围（比如负深度或大于网格对角线长度的深度），平面图上看起来像彩色噪点。

原因：最典型的是窗口落在异常平缓区。当窗口内 T、Tx、Ty、Tz 的幅值都接近零时，Euler 方程的系数矩阵 A 接近零矩阵，最小二乘解无约束，深度就随便飘了。其次是 SI 设置过小，方程右端太小，解对噪声的放大倍数大增。

解决：给 A'*A 加阻尼项是第一步（前面提过的 lambda 方法）；更治本的办法是设置解的筛选条件——只保留深度解在窗口范围内（或者窗口对角线的 1.5 倍以内）的点。代码通常会在后处理里做一遍筛选，如果没有，你自己加一行条件判断即可。还有一招：对 sol_z 做 3σ 截断，把飞点设为 NaN。

5.2 深度结果系统性偏浅或偏深

现象：反演得到的场源深度和钻孔验证深度对比，整体偏差超过 30%，且偏同一个方向。

原因：几乎可以确定是 SI 选错。SI 偏大，方程右端 N*T 偏大，为了让等式成立，左边导数的权重相对变小，解出的 z0 偏浅；SI 偏小则相反。如果 SI 用对了还是偏差，下一个疑点是背景场没做干净——区域背景场叠加在异常上等效于改变了 T 的直流分量，让深度解产生固定偏移。

解决：先做 SI 扫描测试——把 SI 从 0 到 3 每隔 0.5 跑一遍，看深度解的统计中位数随 SI 的变化曲线。曲线斜率最大的区间就是最敏感的 SI 区间，取曲线拐点附近的值作为最优 SI。然后对比原始数据和去趋势数据的反演结果，如果系统性偏差显著缩小，说明背景场是大问题。

5.3 磁化方向没校正导致场源水平位置偏移

现象：反演出的场源位置沿磁化方向平移，比如在低磁纬度区，解的 x 坐标整体向东或向西偏移数倍于预期误差。

原因：Euler 反演假设异常已经做了化极处理（RTP），或者假设磁化方向是垂直的。在中低纬度地区，斜磁化会使异常峰值位置偏离场源真实位置，直接转入 Euler 反演，场的几何关系就不对。

解决：在做导数计算之前，先对磁异常做化极。MATLAB 里可以用频率域化极算子，但化极在低磁纬度区是不稳定的，通常配合维纳滤波或倾斜校正来压制放大噪声。如果你没有化极工具，至少要做到在报告里说明未化极，并用剖面拟合来补偿位置偏移。

5.4 窗口边缘解质量差，形成一条“边框效应”

现象：深度解图像的边缘两圈出现规则的高值或低值条纹，内部区域正常。

原因：窗口扫描到网格边界附近时，窗口内的有效网格点不足——代码里没有对越界做反射或裁剪，而是直接截断，导致窗口内数据量骤减，方程欠定。再加上频率域导数的边缘 taper 效应，边界数据失真。

解决：最简单的方法是在反演前对网格做外扩——每边加宽 10~20 个点距，用边缘均值或向上延拓值填充，反演后只取中间区域的结果。这个做法会轻微平滑边界附近的细节，但能保证有效区域的解不受边缘效应干扰。我个人做法是：外扩 20 个点距，不参与最终结果解释。

5.5 窗口太大把多个场源当成一个解

现象：两个相距很近的异常体，反演解只落在一个中间位置，深度偏大，像是“把两个场源平均了”。

原因：窗口物理尺寸超过了两个场源的间距。Euler 方程在窗口内的最小二乘解是统计意义上的中心位置——当窗口内有两个不同深度的场源时，解很难同时满足两组场的梯度关系，结果会偏向信号更强的一个。

解决：缩小窗口尺寸，或者做分层解释。先用小窗口快速扫描，找到浅部解的横向定位；再判断是否需要深部场源——这通常靠对比不同窗口尺寸的结果图来实现。还有一个技巧：计算每个窗口的拟合残差，残差大的位置意味着单一场源假设不成立，那里需要更细的窗口划分。

5.6 网格点距不均匀造成解的横向漂移

现象：反演解的 x 坐标和 y 坐标误差不均匀——一个方向误差小，另一个方向误差大一倍以上。

原因：代码假设 dx=dy。如果你的网格是等经纬度划分的，x 方向的物理间距随纬度变化而 y 方向不变；即使是在投影坐标系下，数据网格化时也可能产生各向异性的点距。

解决：在进入反演前，把网格重采样为等间距方形网格。MATLAB 里用 interp2 做重采样即可：

% 不规则网格重采样为等间距网格（以 dx=dy=50m 为例） xq = min_x:50:max_x; yq = min_y:50:max_y; [XQ, YQ] = meshgrid(xq, yq); TQ = interp2(X, Y, T, XQ, YQ, 'cubic');

重采样后检查 dx 和 dy 是否在 1% 误差范围内，再进反演。这个坑非常隐蔽——因为代码不会报错，只会产出轻微畸变的结果，很多人因此错过正确结论。

6. 让反演结果真正拿得出手：参数验证、置信度评估和剖面出图的一条龙操作

反演跑完只是第一步，你还需要把结果验证到可以写进报告的程度。我最常用的是“合成数据回验 + 实测数据剖面交叉验证”的组合拳。合成回验的做法是：设计一个已知深度的场源模型（比如一个深 150 米的球体），正演得到理论异常，加入 2% 的随机噪声，然后跑同一套 Euler 代码，统计反演深度与真实深度的误差分布。这一步能同时校验你的 SI 设置、窗口大小和阻尼系数。代码不用另外写，MATLAB 自带的 sphere 模型磁异常正演 20 行就能搞定：

% 合成球体磁异常正演（用于验证反演参数） % 球体参数：埋深 150m，磁化率 0.01 SI，地磁场强度 50000nT % 网格范围：1000m x 1000m，点距 10m [xg, yg] = meshgrid(-500:10:500, -500:10:500); depth = 150; r = sqrt(xg.^2 + yg.^2 + depth^2); % 简化公式：垂直磁化球的 TMI 异常（仅作验证用） T_syn = 50000 * 0.01 * (4/3*pi*20^3) * (depth^2 - xg.^2 - yg.^2) ./ r.^5; % 加 2% 高斯噪声 T_noisy = T_syn + 0.02 * std(T_syn(:)) * randn(size(T_syn));

跑完回验，把合成数据的深度解统计出来——中位数、标准差、偏差百分比都记录到文本里。这套代码能承受的噪声水平，决定了你的实测数据能达到的精度上限。

剖面交叉验证是另一道必做的工序。Euler 反演的结果是离散点云，地质解释时很难直接使用。我的习惯是：在目标异常体走向上切一条剖面，把 Euler 解的 x 坐标投影到剖面线上，深度保留，用散点图叠加在实测异常剖面下方。这样做的价值是，你可以直观看到解点的深度分布是否集中，是否存在多解簇，是否和已知钻孔位置吻合。如果解点在剖面上呈现从浅到深的连续条带，说明这可能是接触带或断裂带；如果解点呈团簇状，说明是有限延伸的独立场源。

最后一个技巧和出图有关。Euler 反演结果直接输出为彩色平面图在学术论文里说服力不够，因为读者无法快速评估解的置信度。我通常会在深度图上叠加两个附加图层：一是每个窗口解的拟合残差（residual norm），二是解的覆盖次数。拟合残差低、覆盖次数高的区域，用高饱和色；反之则淡化。这样一张图同时传达位置、深度和可信度，审稿人或技术评审可以一眼看出哪些结论是硬的、哪些是软的。

素材准备好了之后，我强烈建议你把每一步参数——包括去趋势阶数、窗口尺寸、步长、SI 设定值、是否化极、是否加阻尼——都记录在一个独立脚本里，确保未来任何时刻能完整复现同一份结果。这一步被很多人忽略，等三个月后要回溯数据时追悔莫及。从那以后我每次跑完反演都会强制走一遍“脚本记录 + 合成回验 + 剖面出图”这套流程，不让任何一组参数成为无源之水。希望这套思路能帮你在自己的磁测数据上少走几个来回。

本文还有配套的精品资源，点击获取