简介本资源是一份面向地球物理勘探方向研究生、科研人员及MATLAB初学者的一维大地电磁MT正演建模工具聚焦于无损地质探测中电导率结构的理论模拟与数据生成。它解决了MT方法教学与研究中正演计算门槛高、代码实现复杂的问题适用于矿产勘查、地壳电性结构分析等实际场景。压缩包为RAR格式仅含1个核心MATLAB脚本文件.m体积仅1KB轻量简洁可直接运行或嵌入更大反演流程该脚本实现了基于一维层状介质假设的频率域电磁场正演涵盖Maxwell方程离散、层间边界条件处理及阻抗张量计算等关键环节。已有344人学习下载代码经实测验证具备良好稳定性与计算精度附带清晰参数接口与物理量注释便于理解电磁场传播机制、调试模型参数并为后续反演提供可靠正演引擎。 搞地球物理的人尤其是做电磁法的十有八九绕不过大地电磁测深MT的一维正反演。不管你是本科生做毕业设计、研究生处理实测剖面还是工程师想快速验证一个盆地的深部电性结构MT1D都是最常用的一块敲门砖。配合Matlab这种随手就能跑的脚本环境从正演到反演的整套流程甚至可以在一个下午搭出雏形。这篇博文把我自己在 MT1D 上积累的“能直接抄”的经验写下来从物理原理、递推公式到 Matlab 代码从反演套路到实测数据踩坑。内容偏实用数学推导我会尽量讲清楚但不会堆满公式适合刚接触这个方向、又想快速上手的人。1. 从野外数据到一维地电模型MT1D在做什么1.1 大地电磁测深的基本物理图像大地电磁测深英文叫 Magnetotelluric简写 MT。它的物理图像可以用一句话说清楚在地表测量天然交变电磁场利用电磁波在导电介质中传播时的趋肤效应反推地下电阻率随深度的变化。天然场源主要是太阳风和雷电活动频带范围很宽从零点几毫赫兹到几千赫兹都能用。频率越低穿透深度越大频率越高探测深度越浅。所以把不同频率下的视电阻率和相位测出来就等于对地下从浅到深做了一次“电性CT”。MT1D 是其中最简单的一维模型假设地下只存在电阻率随深度变化、在水平方向无限延伸的层状介质。这个假设听着很理想实际上在盆地、平原、层状沉积岩地区相当可靠。很多情况下一维反演结果虽然不是最终答案但能给出非常合理的初始模型为二维、三维反演打底。那么一维正演到底算什么它计算的是给定一个层状地电模型每层电阻率厚度在某个频率的天然电磁场激发下地表能够观测到的视电阻率和相位。反过来反演就是给定一组实测视电阻率-相位曲线去寻找一个层状模型让它的理论响应尽可能接近实测数据。一正一反构成了 MT 资料解释最基础的两个操作。1.2 一维假设什么时候成立什么时候不成立这里经常有新手把一维结果直接拿去解释断层。要注意MT1D 只适合在横向电性变化比较弱的测点使用。如果测量剖面上有明显的大型断裂、岩体接触带、地下河道那么单测点的视电阻率曲线会受到三维效应影响一维反演得到的层厚度和电阻率会发生系统偏差。我自己的经验是两条判断标准区域构造背景清楚目的层基本呈层状同一构造单元内的相邻测点曲线形态一致没有突变。如果满足这两条一维反演往往比二维三维折腾半天还靠谱。如果不满足至少也要把一维结果当成约束或者初始模型来用别直接下地质结论。很多 MT 反演软件默认生成的一维反演剖面图其实只是在“画”每个测点单独反演的层状柱体它并不能代表真实的地下结构这一点在写报告时一定要小心。1.3 为什么用Matlab做MT1D正反演这个话题我经常被问现在 Python 那么多库为什么还用 Matlab我的回答是Matlab 在矩阵运算、快速可视化、交互调试上的体验确实好尤其是地球物理这种需要反复试探数据的场景。MT1D 正演本身不复杂用 Matlab 几十行就能写清楚反演里涉及的雅可比矩阵计算、线性方程组求解、误差统计Matlab 自带函数足够用。而且很多学校、地调单位的老代码都是 Matlab 写的兼容性是个现实问题。如果你所在团队已经迁移到了 Python那当然可以用 numpy/scipy 重写。但如果是自己学习 MT 算法Matlab 不是一个坏选择最关键的是它让你把精力放在算法本身而不是环境配置上。下面这部分我会从一个最简单的正演函数开始逐步扩展到反演。2. 一维正演的数学骨架阻抗递推与视电阻率2.1 从Maxwell方程组到Helmholtz方程先做一点必要的数学铺垫。在 MT 方法里我们假设介质是线性、各向同性、非磁性的磁导率取真空磁导率 (\mu_0 4\pi \times 10^{-7} , \mathrm{H/m})。在谐变时间因子 (e^{-i\omega t}) 约定下由 Maxwell 方程组可以得到电场 (E) 和磁场 (H) 满足的 Helmholtz 方程[ \nabla^2 E - i\omega \mu_0 \sigma E 0 ]其中 (\sigma 1/\rho) 是电导率。在一维层状介质中电磁场只沿深度方向变化水平方向的波数由介质本身决定。定义复波数[ k \sqrt{-i\omega \mu_0 / \rho} \sqrt{-i\omega \mu_0 \sigma} ]这个 (k) 的实部对应电磁波的衰减虚部对应相位变化。它的倒数就是通常说的趋肤深度。对于均匀半空间电场振幅从地表向下衰减到 (1/e) 的深度为[ \delta \approx 503 \sqrt{\rho / f} \quad (\mathrm{m}) ]这个公式非常有用。比如电阻率 100 Ω·m频率 1 Hz趋肤深度大约是 5 公里频率降到 0.001 Hz趋肤深度就接近 160 公里。你可以用这个公式快速估算每个频点大概能反映多深的电性结构比看反演结果还直观。2.2 层状介质中的阻抗递推公式一维 MT 正演的核心是算地面阻抗也就是地表电场水平分量与磁场水平分量的比值。这个比值满足从下往上的递推关系。假设模型一共有 (N) 层第 (j) 层的电阻率为 (\rho_j)厚度为 (h_j)其中第 (N) 层为无限厚半空间。第 (j) 层的固有阻抗也叫特征阻抗为[ Z_{0j} \sqrt{i\omega \mu_0 \rho_j} ]在最底层半空间由于没有下伏反射顶面阻抗就等于该层固有阻抗[ Z_N Z_{0N} ]向上递推在第 (j) 层顶面的阻抗为[ Z_j Z_{0j} \frac{Z_{j1} Z_{0j} \tanh(-i k_j h_j)}{Z_{0j} Z_{j1} \tanh(-i k_j h_j)} ]其中 (k_j \sqrt{-i\omega \mu_0 / \rho_j})。一直递推到第 1 层顶面就得到地表阻抗 (Z_1)。视电阻率和相位定义为[ \rho_a \frac{1}{\omega \mu_0} |Z_1|^2 ][ \phi \arg(Z_1) ]注意不同文献的相位符号可能不同取决于时间因子到底是 (e^{-i\omega t}) 还是 (e^{i\omega t})。Matlab 默认的复数计算没有这类约定问题但你在对比代码结果时一定要确认正负号规则。2.3 为什么用阻抗递推而不是直接解微分方程很多新手会问既然知道每一层都是均匀介质为什么不能直接写出解析表达式其实上面的递推公式就是由每一层内部的解析解和边界条件导出的。在实际编程中递推比直接求矩阵逆更稳定、更快。另外值得注意的是(\tanh) 的参数是复数而且随着频率降低、层厚度增加参数虚部的模会变得很大。Matlab 内置的复数双曲正切函数能够处理这种大数值不会像某些简单的近似公式那样产生数值溢出。这也是我建议直接用库函数、不要自己去展开 (\tanh) 的原因之一。3. Matlab实现MT1D正演代码有细节才稳定3.1 主函数设计一个完整的 MT1D 正演函数应该接收三个输入各层电阻率数组rho、各层厚度数组h、频率数组f。厚度数组长度比电阻率少一层因为最后一层是半空间不需要厚度。输出是视电阻率和相位。下面是一个可以直接运行的 Matlab 函数function [rhoa, phase] mt1d_fwd(rho, h, f) % MT1D_FWD 一维大地电磁正演 % rho : 各层电阻率单位 Ohm-m长度 N % h : 前 N-1 层层厚度单位 m长度 N-1 % f : 频率数组单位 Hz长度 M % rhoa : 视电阻率单位 Ohm-m长度 M % phase : 相位单位 degree长度 M mu0 4 * pi * 1e-7; w 2 * pi * f(:); % 转为列向量 N length(rho); Nf length(w); % 从最底层半空间开始Z 是每个频率下的地表阻抗 Z sqrt(1i * w * mu0 * rho(N)); % 最底层顶面阻抗 for j N-1:-1:1 k sqrt(-1i * w * mu0 / rho(j)); Z0 sqrt(1i * w * mu0 * rho(j)); % 递推公式注意 h(j) 是标量 Z Z0 .* (Z Z0 .* tanh(-1i * k * h(j))) ./ ... (Z0 Z .* tanh(-1i * k * h(j))); end rhoa abs(Z).^2 ./ (w * mu0); phase atan2(imag(Z), real(Z)) * 180 / pi; end这个函数的核心是for j N-1:-1:1的向上递推。每次循环里Z是所有频率在当前层顶面的阻抗向量既包含当前层的影响也一直携带下方所有层的信息。由于是逐层向上最后得到的就是地表阻抗。3.2 代码里有几个容易忽略的细节第一rhoa abs(Z).^2 ./ (w * mu0)这里的w是角频率不是频率。如果拿频率f去算视电阻率会差 (2\pi) 倍。第二f(:)把行向量转成列向量是为了避免w是行而Z是列时出现维度问题。第三phase atan2(imag(Z), real(Z)) * 180/pi用atan2而不是angle也完全可以但要注意输出范围是 -180° 到 180°。实际 MT 相位通常在 0° 到 90° 之间如果看到负相位先检查坐标轴方向和时间因子约定。另外这个函数默认最底层是半空间因此h的长度必须是length(rho)-1。如果少传一层厚度循环里的h(j)会越界。这类低级错误在调试中最浪费时间所以我会在函数开头加两个检查assert(length(h) N-1, h length must be N-1); assert(all(f 0), frequency must be positive);3.3 标准算例验证两层模型与三层模型写完函数第一件事不是上实测数据而是用解析程度很高的模型验证。最简单的两层模型上层电阻率 100 Ω·m、厚度 1000 m基底电阻率 1000 Ω·m。用正演函数算一条频带从 0.001 到 1000 Hz 的曲线rho [100, 1000]; h [1000]; f logspace(-3, 3, 60); [rhoa, phase] mt1d_fwd(rho, h, f); figure; subplot(2,1,1); loglog(f, rhoa, k-o); ylabel(视电阻率 / \Omega.m); set(gca, XDir, reverse); subplot(2,1,2); semilogx(f, phase, r-o); xlabel(频率 / Hz); ylabel(相位 / deg); set(gca, XDir, reverse);从物理上预期高频时趋肤深度远小于覆盖层厚度能量基本在上层传播视电阻率应该接近 100 Ω·m低频时电磁波穿透到基底视电阻率应该接近 1000 Ω·m。相位曲线则在两层电性差异对应的频段出现一个偏移高阻基底顶部往往对应相位谷或峰具体形态取决于界面两侧电阻率比值。这个算例跑通后再试三层模型比如“低阻盖层—高阻储层—半空间基底”。你会发现视电阻率曲线在某些频段出现“平台”对应不同的电性层。这一步能帮你建立曲线形态和地电模型之间的直觉比看十个理论图件都管用。4. 一维MT反演的常见套路在Matlab里怎么落地4.1 反演目标函数与阻尼最小二乘正演解决了反演才是真的难点。一维反演的目标是找一组模型参数使正演响应和实测数据在误差范围内一致。我建议把模型参数设置为对数电阻率和对数厚度而不是直接使用电阻率和厚度的线性值。原因有两条电阻率跨度常常有几个数量级线性最小二乘会把高阻层权重放得过大厚度必须为正如果用线性厚度迭代时可能出现负值导致正演崩溃。把数据也做对数处理尤其是视电阻率能压缩动态范围让反演更加稳定。相位本身已经是角度一般直接使用但也要换成弧度或者保持角度制都行。反演目标函数常见形式是[ \Phi(m) | W_d (d - F(m)) |^2 \lambda | C m |^2 ]其中 (d) 是实测数据向量可以是对数视电阻率相位(F(m)) 是正演响应(W_d) 是对角权重矩阵每个元素是数据误差的倒数(\lambda) 是正则化参数(C) 是粗糙度矩阵。第一项衡量拟合程度第二项压制模型的过度振荡。没有第二项这个反演多半是病态或不稳定的。最常用的求解方法是 Gauss-Newton 迭代[ m_{k1} m_k \left[J^T W_d^T W_d J \lambda C^T C\right]^{-1} J^T W_d^T W_d (d - F(m_k)) ]这里 (J) 是雅可比矩阵也叫灵敏度矩阵它的第 (i) 行第 (j) 列表示第 (i) 个数据对第 (j) 个模型参数的导数。在 Matlab 里最省事的方式是用有限差分近似[ J_{ij} \approx \frac{F_i(m \Delta m_j) - F_i(m - \Delta m_j)}{2 \Delta m_j} ]虽然比解析导数慢但对 MT1D 这种数据量不大的问题完全够用。4.2 Occam型反演与粗糙度约束如果你不想从零写最小二乘我强烈推荐实现一个 Occam 型反演。Constable 等人 1987 年提出的经典流程到目前为止仍是一维 MT 反演里最稳的算法。它的特点是在每一步迭代里显式地搜索正则化参数 (\lambda)使得模型足够平滑的同时满足拟合差要求。Occam 反演的迭代公式和阻尼最小二乘很像但有一个关键区别它用当前模型线性化后重新构造了一个“目标数据”[ \hat{d} d - F(m_k) J m_k ]然后解[ m_{k1} \left[J^T W_d^T W_d J \lambda C^T C\right]^{-1} J^T W_d^T W_d \hat{d} ]对于每一轮迭代给定一个 (\lambda)就得到一个候选模型 (m_{k1})。Occam 的做法是把 (\lambda) 在一个对数范围里扫一遍对每个 (\lambda) 都算一遍正演响应和模型粗糙度然后选择一条“拟合差下降最多、模型又足够平滑”的曲线。实际实现时经常画一张 (\lambda)-RMS 折中曲线找到拐点。粗糙度矩阵 (C) 的定义很灵活。最简单的是一阶差分C zeros(N-1, N); for j 1:N-1 C(j, j) -1; C(j, j1) 1; end这会把相邻两层电阻率之差作为惩罚项。如果想要更平滑的结果可以用二阶差分如果希望保留突变界面可以改用总变分TV正则化。后者实现稍麻烦但很多地质场景下效果更好。4.3 固定层厚度反演新手最该先掌握的方式在你自己写一维反演的第一版代码时我不建议把层电阻率和层厚度都当成未知数。厚度和电阻率存在强烈的“等效性”也就是增厚一层高阻层的同时减小其电阻率正演响应几乎不变。这个现象会让反演矩阵接近奇异产生非常夸张的锯齿模型。更稳的做法是预先设计一个固定厚度网格层厚度从浅到深逐层等比增加比如浅地表 50 m、100 m、200 m、500 m、1000 m……一直到底部半空间。反演只更新每层电阻率厚度全部固定。这样一来未知数从 (2N-1) 个降为 (N) 个稳定性大大提高最后得到的是一个“随深度变化的光滑电阻率剖面”。这个结果足够用于绝大多数解释工作。Matlab 实现时模型参数 (m_j \log_{10}\rho_j)层厚度h_grid固定不变。正演函数需要稍作修改因为最后一层是半空间h长度为 (N-1)。如果网格设计为 20 层那么h包含前 19 个厚度。固定厚度反演唯一需要担心的就是表层高阻薄层网格太粗会完全抹掉。不过这是 MT 本身分辨率的问题不是反演算法能解决的这部分后面再说细节。4.4 初始模型、正则化参数与收敛判断固定厚度反演的初始模型一般用均匀半空间比如所有层给 100 Ω·m。如果数据动态范围极大可以先求平均视电阻率作为初始值让第一轮正演响应不至于偏太远。关于正则化参数 (\lambda)我习惯从一个较大的初值开始比如 (\lambda10)然后按 0.5 倍或 0.2 倍递减。每一轮迭代后计算 RMS misfit[ RMS \sqrt{\frac{1}{M} \sum_{i1}^M \left( \frac{d_i^{obs} - d_i^{syn}}{\epsilon_i} \right)^2} ]其中 (M) 是视电阻率和相位所有数据点的总数(\epsilon_i) 是第 (i) 个数据的误差。当 RMS 接近 1 时说明数据基本被拟合在误差范围内如果 RMS 长期停在 5 以上不再下降先查数据质量再查正演代码不要盲目加大迭代次数。收敛判断除了 RMS 之外还要看模型更新的范数d_m norm(m_new - m_old) / norm(m_old); if d_m 1e-3 break; end如果模型参数更新已经很小但 RMS 还没降到合理范围说明当前模型已经卡在局部极小需要换一组初始模型或者调整粗糙度矩阵。5. 实测数据进反演前必须做的事5.1 数据筛选哪些点不能进入反演我见过太多人拿原始功率谱输出的视电阻率曲线直接反演结果模型一团糟。其实 MT 处理软件已经给出了质量控制参数关键是你会不会用。删除相干度coherency低的频点一般低于 0.8 的不要留删除“死频点”也就是由已知人文电磁干扰导致的 50 Hz 及其谐波附近的频点检查相位是否落在合理的物理范围内。一维正演的相位通常在 0° 到 90° 之间如果出现大量负相位或大于 90°要么是参考道方向选择错误要么是数据质量太差检查视电阻率曲线是否出现“飞点”也就是相邻几个频点连续变化但某个点突然跳出一个尖峰。飞点先删除不要直接平均。误差估计也不能随便给。如果处理软件给出了误差棒就用真实误差如果没给我会给视电阻率 10% 的相对误差、相位 5° 的绝对误差。这样反演里权重分配相对合理不会让个别低质量和高质量数据混在一锅端。5.2 TE、TM模式选择和联合一维正演本身不区分 TE 和 TM 模式因为在地下水平均匀时两者响应是一样的。但实测曲线的 TE 和 TM 可能差异很大这是地下电性不均匀造成的不能靠一维反演消除。我的处理习惯是如果只有单测点曲线优先用 TE 模式做一维反演因为 TE 模式受局部体充电效应也就是静位移影响通常小于 TM 模式如果有剖面数据把 TE 和 TM 分别做一维反演对比两条电阻率深度曲线不要在一维反演中把 TE 和 TM 硬塞进同一个数据向量里。两者差异不是一维能解释的强行联合只会让反演产生不真实的光滑模型甚至完全扭曲浅层结构。当然有一种情况例外如果两个模式曲线经过静位移校正后形态一致那么可以联合反演增强数据约束。但这种情况在实测中并不常见。5.3 静位移和近场畸变怎么识别静位移是近地表局部不均匀体引起的电流集中或稀疏效应它的特征是视电阻率曲线整体乘以一个频率无关的常数相位曲线几乎不变。在双对数坐标上静位移表现为整条视电阻率曲线相对于相邻测点上下平移。识别方法很简单把同一测点的 TE 和 TM 视电阻率曲线画在一起如果两条曲线形态平行但有一定垂直间距而相位曲线几乎重合那么大概率是静位移。一维反演遇到这种情况应该先做校正。校正方法很多用相位反演、用 AMT 高频数据约束浅层电阻率、或者做空间滤波。你至少应该知道不校正直接反演反演出来的浅层电阻率会被整体抬高或压低深度界面倒不一定失真太多。近场畸变则通常来自长周期段比如磁暴扰动、人工供电信号不满足平面波假设。表现为视电阻率曲线在最低频段斜率接近 -1相位趋近 0° 或 180°这是 MT 数据里很典型的“近场源”特征。这部分数据如果在反演前不做剔除反演会为了拟合它强行把深部电阻率改得极其反常结果完全不可信。6. 我踩过的坑和调试技巧6.1 高频段视电阻率震荡有次我写正演代码时在 1000 Hz 以上视电阻率总出小幅锯齿一开始以为是递推公式有问题。后来发现是因为我在定义频率数组时用了0.001:0.01:1000这种线性间隔高频段频点太密、低频段频点太少而且没做对数均匀分布。把频率改成logspace(-3, 3, 81)之后曲线立刻光滑了。如果你已经用了logspace还是震荡那大概率是表层网格太薄。在高频段电磁波趋肤深度可能只有几十米如果第一层厚度设成 500 m那么高频响应被第一层完全控制后续层的影响极小这时理论上不应该震荡。真正出现震荡通常是某个厚度或者电阻率输入了 NaN。我会在正演函数里打印每一层的k*h模值检查是否出现异常大或无穷。6.2 单位与坐标轴方向MT 的单位坑我踩了不止一次。一维正演里电阻率用 Ω·m厚度用 m频率用 Hz角频率 (ω2πf)。如果从老论文里抄公式可能用的是角频率而不是频率也可能用 CGS 单位。CGS 下 (\mu_01)电阻率单位是 cm·s换算极其痛苦。我的建议是程序中全部用国际单位制只在输入输出层做单位转换。还有坐标轴方向。MT 画图惯例是高频在左、低频在右也就是说横轴是“频率倒数”或“周期”低频对应深部。很多新手拿正演结果一画横轴默认从小到大低频在左侧视觉上“深部”反而在左边解释时特别容易混淆。我在验算时习惯把set(gca,XDir,reverse)让高频在左、低频在右之后看文献对应起来更舒服。6.3 薄层分辨率和等效模型MT 对低阻薄层有比较高的灵敏度但对高阻薄层灵敏度很差。一维反演想恢复一个几十米厚的火成岩高阻岩床基本不可能从反演曲线里得到准确厚度和电阻率。这并不一定是你算法有问题而是 MT 本身的分辨率限制。高阻薄层的等效性可以用一个小实验证明rho1 [100, 1000, 100]; h1 [500, 50, 5000]; rho2 [100, 500, 100]; h2 [500, 200, 5000]; [ra1, ph1] mt1d_fwd(rho1, h1, f); [ra2, ph2] mt1d_fwd(rho2, h2, f); max(abs(ra1 - ra2))你会发现两个模型的正演响应差异非常小。如果地质上已经知道薄高阻层存在反演时应该把该层厚度固定成已知厚度只反演电阻率否则反演会倾向于把它抹平或者加厚。6.4 反演不收敛先查雅可比矩阵我最早写反演时怎么迭代 RMS 都降不下去最后发现是雅可比矩阵的数值差分步长选错了。对对数电阻率参数步长一般取 0.01 到 0.1 之间。太大差分结果偏离真实导数太小两个正演响应之间的差异被浮点舍入误差淹没雅可比矩阵会变得全是噪声。调试办法随便选一个模型做一次小扰动扫描。比如把第一层电阻率从 (10^{-0.3}) 到 (10^{0.3}) 倍连续变化看视电阻率响应的变化是否平滑。如果曲线出现毛毛躁躁的跳动说明步长有问题或者rhoa计算里某个参数在边缘情况不稳定。另外反演每步迭代后不要只看总 RMS把视电阻率和相位的分项 RMS 分别打印。我经常遇到的情况是视电阻率拟合得很好但相位很差这说明模型参数可能陷入了一个等效模型——电阻率整体对但界面位置不对。此时应该适当调整相位误差权重或者加入相位对应的粗糙度约束。最后再分享一个小习惯任何正演和反演代码都要在“异常输入”场景下做测试。比如输入一个超高阻层 1e8 Ω·m或者一个超薄层 0.1 m看看代码是崩溃、返回 NaN 还是给你一个合理响应。把这些边界情况提前搞定野外数据来了你才能放心大胆地跑而不是一边看反演曲线一边心里打鼓。MT1D 虽然简单但真正让它“可复现、可信”的恰恰是这些不起眼的细节。本文还有配套的精品资源点击获取