简介本资源是一套面向地球物理勘探初学者与科研人员的MATLAB一维大地电磁MT正演计算工具聚焦于地质结构电性参数建模与地表电磁响应模拟适用于矿产勘查、地热调查及基础MT方法教学实践。压缩包仅含1个核心文件——mt1d.m脚本体积仅1KB代码精炼完整实现基于频率域Maxwell方程的一维层状介质正演求解支持用户自定义电导率剖面、频率范围及磁导率参数并输出视电阻率与相位曲线等关键响应数据。已有344人学习下载程序经实测验证稳定性与精度良好内置清晰注释与典型参数示例便于快速理解电磁场传播原理、调试模型设置并衔接后续反演分析。 做大地电磁的人手里多多少少都攒着几个自己写的MT脚本。我最早接触这个方向的时候也跑去翻现成的开源代码结果要么语言不合适要么处理流程绕最后干脆用MATLAB自己写了一套MT1D正反演。今天把项目里的核心算法、代码结构、调试经验整理出来给同样需要处理大地电磁数据的同行做个参考。这个项目做的事情很明确用MATLAB实现大地电磁一维正反演MT算法正演部分用层状介质的Wait递归公式反演部分用阻尼最小二乘加正则化约束能直接吃进视电阻率和相位曲线输出地电模型。对做电磁法勘探、地热调查、水文地质、油气探测的同行来说这套代码可以当成自研处理流程的起点对学生来说把代码跑通一遍也就把MT正演和反演的本质理解得差不多了比闷头读论文直观得多。整个项目写完以后我顺手把调试记录和踩坑清单也整理了一遍。这篇不光是贴代码我更想把代码背后的设计逻辑、参数为什么这么调、实测数据会遇到哪些坑以及我在不同MATLAB版本下跑出来的经验写清楚至少让你少走我走过的弯路。1. 项目到底在做什么MT1D的前因后果1.1 MT方法的基本物理图景大地电磁测深法Magnetotellurics简称MT利用天然交变电磁场作为激励源在地表观测相互正交的电场和磁场分量比如Ex和Hy或者Ey和Hx然后构造波阻抗Z Ex / Hy单位是欧姆。这个阻抗本身就是地下介质电性结构的函数频率越高电磁波在地下的传播深度越浅频率越低穿透深度越深。把一系列频率的视电阻率和相位画成曲线就得到了一条带有深度信息的电磁响应曲线。天然场的频率范围很宽从0.0001 Hz量级的超低频到上万赫兹的音频频率。频率和探测深度的大致关系可以用趋肤深度δ来估算δ约等于503乘以根号下电阻率除以频率常见的近似形式是δ ≈ 503√(ρ/f)单位是米。举个例子电阻率100 Ω·m频率1 Hz的时候趋肤深度大约5公里这意味着我们在1 Hz看到的信号主要反映地下5公里以内的电性分布。这个粗略估算贯穿了MT数据解释的整个过程从定频段选点到后期反演初始模型设置我都经常拿它来快速判断数据能分辨多深。1.2 为什么先从一维做起一维模型是MT解释里最朴素也最可靠的模型它假设地下电性只随深度变化水平方向无限延伸每一层内部各向同性且均匀。这个假设虽然在真实地质体里很难严格成立但没有一维打底二维三维根本没法开工。我自己的经验是一维正演至少有三个用途。第一用来做数据质量检查。野外实测的视电阻率和相位曲线如果跟一维响应对不上往往说明数据受三维结构、地表不均匀体或者人文电磁噪声影响比较严重这时候直接拿去做二维或三维反演很可能得到假异常。第二用来提供三维反演的初始模型。现在主流的三维反演软件虽然能自动生成初始模型但如果你想给一个更有地质意义的初始模型一维反演结果逐测点拼接是最常见的做法。第三用来做理论实验。一维反演速度快、参数少是验证算法、调参、理解反演非唯一性的最佳场地。我自己很多关于正则化参数、误差估计的直觉都是在一维环境里养成的换到二维三维以后很多原则可以平移。1.3 功能规划与代码文件结构项目没有搞复杂的工程化结构核心逻辑就是“正演、雅可比、反演、演示”四件事。文件结构大致如下文件作用关键内容mt1d_forward.mMT一维正演函数输入频率、层电阻率、层厚度输出视电阻率和相位mt1d_jacobian.m雅可比矩阵计算对数电阻率参数下的数值差分雅可比mt1d_inv.m阻尼最小二乘反演框架迭代更新、正则化、拟合差计算mt1d_demo.m合成数据演示脚本生成三层模型正演响应加噪声再反演复原模型我之所以把雅可比单独拆出来是因为反演主循环需要反复调用它。而且调试的时候单独看雅可比矩阵比混在反演里更容易发现问题比如某一行全部是0十有八九是差分步长或者参数索引写错了。整体设计有一个原则全部用基础MATLAB函数不依赖优化工具箱、统计工具箱这类需要单独授权的模块。这么做的原因很现实很多学校、企业的机器上只装了基础MATLAB有的甚至版本比较老如果你一上来就依赖lsqcurvefit或者fmincon换台机器就跑不了。手写矩阵运算虽然看起来多写几行但换来的是极强的兼容性。2. 正演与反演的核心算法拆解2.1 正演公式从Cagniard到Wait递归MT一维正演的数学基础来自Cagniard和Wait等人的工作。对于均匀半空间地表阻抗直接等于本征阻抗也就是Z √(iωμρ)其中ω是角频率μ是磁导率通常取真空磁导率4π×10⁻⁷ H/mρ是电阻率。视电阻率定义为ρa |Z|² / (ωμ)相位是阻抗辐角的一半。对于多层水平层状介质要从最底层的半空间开始逐层向上递推。每一层顶部的阻抗可以用Wait递推公式计算Z_n η_n × [Z_{n1} η_n × tanh(i k_n h_n)] / [η_n Z_{n1} × tanh(i k_n h_n)]其中η_n是第n层的内禀阻抗η_n √(iωμρ_n)k_n是复波数k_n √(iωμ/ρ_n)h_n是第n层的厚度。这个公式从最底部的半空间开始算Z_1就是地表阻抗然后代入视电阻率和相位的定义式就得到响应。这里容易被绕晕的是η和k的符号写法。不同文献里定义略有差异有的习惯用σ表示电导率写成η √(iωμ/σ)其实和ρ版本完全等价。你自己写代码的时候只要保持一致算出来的视电阻率不会差。我习惯全部用电阻率表示避免来回换算电导率出错。还有一个需要注意的点公式里是复数tanhMATLAB对复数tanh支持得很不错但如果频率太高、层太厚i k_n h_n的虚部会非常大可能出现数值溢出。实操中我一般把频率范围和层厚控制在这个公式稳定工作的范围内同时用warn状态检查结果如果出现NaN或者Inf就立刻停下来检查参数。2.2 反演问题的数学表达反演的目标是找一组模型参数使得正演响应尽量拟合实测视电阻率和相位数据。MT一维反演本质上是高度非线性的我用的是带正则化的阻尼最小二乘也叫Levenberg-Marquardt方法。把模型参数记为m通常取各层电阻率的对数也就是m log(ρ)。取对数的好处有两个一是保证迭代过程中电阻率始终为正不需要额外的约束处理二是把跨越几个数量级的电阻率参数压缩到线性尺度上梯度计算和步长控制更稳定。数据向量d则包含视电阻率的对数和相位。之所以不用视电阻率原值是因为大地电磁数据动态范围很大可能从10 Ω·m到10⁴ Ω·m如果直接拟合绝对差值高阻层的权重会淹没低阻层的信息取对数相当于做了对数均匀加权更符合人眼判断拟合效果的习惯。目标函数写为Φ ‖W(d_obs - F(m))‖² α‖L(m - m_ref)‖²第一项是数据拟合差W是数据协方差矩阵的逆用来对不同的数据点按误差大小归一化第二项是正则化项L是一阶差分矩阵m_ref是参考模型α是正则化参数。这个式子用中文说就是既要让模型响应接近观测数据又不能让模型在相邻层之间疯狂跳变。阻尼最小二乘的迭代更新公式是m_new m_old (JᵀWᵀWJ λI αLᵀL)⁻¹ JᵀWᵀW(d_obs - F(m_old))其中J是雅可比矩阵λ是阻尼因子。λ的作用是在高斯牛顿法和最速下降法之间动态切换λ大的时候步长变小、稳定性增加λ小的时候收敛速度加快。每次迭代如果拟合差下降就减小λ如果拟合差反弹就增大λ后面我会把判定策略写进代码。2.3 MATLAB实现中的几个关键细节第一是数组下标与模型深度对应。层电阻率数组rho的长度是N层厚度数组h的长度是N-1最后一层默认为半空间。为了避免混淆我在函数开头就加注释说明并且用numel(rho)-1来推断层数而不是写死。第二是频率数组的生成。实测MT数据的频点通常是等对数间隔排列的从最低频到最高频每个倍频程大概5到10个频点。生成合成数据时我习惯用logspace函数低频从0.001 Hz到1000 Hz这样视电阻率曲线高低频都能覆盖便于观察规律。第三是相位的单位。观测资料里相位一般给度数但正演计算里atan2出来是弧度所以要把弧度转成度数再放进数据向量。这个单位问题不处理好反演时会莫名出现拟合差跳来跳去的情况排查半天发现是单位混了。第四是雅可比矩阵的数值差分。虽然理论上可以推导解析导数但代码复杂度高、容易出错。我用中心差分对第j个模型参数给一个小的扰动量Δm_j分别计算F(m Δm_j)和F(m - Δm_j)再取差除以2Δm_j。差分步长取1e-4这个量级太小会引入舍入误差太大则线性近似失真实测下来1e-4比较稳定。3. 核心代码实现与合成数据测试3.1 mt1d_forward正演函数怎么写才稳正演函数是整套代码的地基反正演能不能收敛就看正演算得准不准。下面这个版本我去掉了多余的炫技部分只保留最核心的逻辑放进MATLAB就能直接运行function [rho_a, phi_deg, Z_surf] mt1d_forward(freq, rho, h) % MT1D_FORWARD MT一维正演计算层状模型的视电阻率和相位 % 输入: % freq : 频率数组单位Hz建议按对数等间隔排列 % rho : 层电阻率数组单位Ohm.m长度N最后一层为半空间 % h : 层厚度数组单位m长度N-1对应前N-1层的厚度 % 输出: % rho_a : 视电阻率单位Ohm.m % phi_deg: 相位单位度 % Z_surf : 地表波阻抗复数 mu0 4 * pi * 1e-7; omega 2 * pi * freq(:); % 转为列向量方便广播计算 N length(rho); if length(h) ~ N - 1 error(层厚度个数必须比电阻率个数少1); end Z_surf zeros(size(omega)); for ifreq 1:length(omega) w omega(ifreq); % 从最底层半空间开始 Z_bottom sqrt(1i * w * mu0 * rho(N)); % 自底向上递推 for n N-1:-1:1 eta_n sqrt(1i * w * mu0 * rho(n)); k_n sqrt(1i * w * mu0 / rho(n)); arg k_n * h(n); Z_bottom eta_n * (Z_bottom eta_n * tanh(arg)) / ... (eta_n Z_bottom * tanh(arg)); end Z_surf(ifreq) Z_bottom; end rho_a abs(Z_surf).^2 ./ (omega * mu0); phi_deg angle(Z_surf) * 180 / pi; % 一维MT相位通常位于0-90度之间angle返回-pi到pi负相位加360 phi_deg(phi_deg 0) phi_deg(phi_deg 0) 360; end这个函数最核心的就是从倒数第二层往上循环的部分每一层都套用Wait递归公式。需要强调的一点是我使用的是Z_bottom作为上一层的输入名字可能让你误解成“底部阻抗”但严格说它代表“从当前层往下看的等效阻抗”也就是Z_{n1}。这样递推到最后Z_surf就是地表阻抗。关于相位一维大地电磁的物理相位大多在第一象限也就是0到90度之间。但MATLAB的angle函数返回值范围是-180到180度如果用实部虚部算出来的相位落在第四象限就会显示成负角度需要加360度转成正度数。这种小细节商用软件里都帮你处理好了自己写代码时最容易漏。3.2 mt1d_jacobian雅可比矩阵的计算雅可比矩阵是反演迭代的核心它的每一行对应一个数据点对每个模型参数的偏导数。由于我用的是对数电阻率参数m所以差分时扰动对数电阻率正演时再把logρ还原成ρ。function J mt1d_jacobian(freq, rho, h, data_type, delta) % MT1D_JACOBIAN 计算MT一维正演对对数电阻率的雅可比矩阵 % data_type: 1代表视电阻率对数2代表相位 % delta : 中心差分步长缺省1e-4 if nargin 5 || isempty(delta) delta 1e-4; end N length(rho); J zeros(numel(freq), N); logrho log(rho); for j 1:N % 正向扰动 rho_plus rho; rho_plus(j) exp(logrho(j) delta); [rho_a_p, phi_p] mt1d_forward(freq, rho_plus, h); % 负向扰动 rho_minus rho; rho_minus(j) exp(logrho(j) - delta); [rho_a_m, phi_m] mt1d_forward(freq, rho_minus, h); % 中心差分 if data_type 1 J(:, j) (log(rho_a_p) - log(rho_a_m)) / (2 * delta); else J(:, j) (phi_p - phi_m) / (2 * delta); end end end把视电阻率取对数再差分是为了和反演目标函数里的数据项保持一致。如果你的反演数据里同时包含视电阻率和相位那就把两部分雅可比拼在一起上面这个函数实际上返回的是单一数据类型的雅可比调用时按列拼接就行。这个函数的性能并不高因为每个模型参数要做两次正演。但一维正演本身开销很小模型层数通常不超过几十层所以完全够用。如果以后扩展到二维或三维再用解析导数或者伴随方法。3.3 mt1d_inv阻尼最小二乘反演框架反演函数是整个项目的主控模块。它接受实测频率、实测视电阻率、相位以及初始模型输出反演后的模型参数和迭代历史。function [rho_inv, hist] mt1d_inv(freq, rho_obs, phi_obs, ... rho0, h, alpha, maxiter) % MT1D_INV 阻尼最小二乘MT一维反演 % rho0 : 初始电阻率模型长度N % alpha : 正则化参数 % maxiter : 最大迭代次数 N length(rho0); m log(rho0(:)); rho_obs rho_obs(:); phi_obs phi_obs(:); d_obs [log(rho_obs); phi_obs]; % 一阶差分正则化矩阵作用于模型参数 L zeros(N-1, N); for i 1:N-1 L(i, i) -1; L(i, i1) 1; end lambda 1.0; % 阻尼因子初值 hist.fit []; for iter 1:maxiter rho_cur exp(m); [rho_a, phi] mt1d_forward(freq, rho_cur, h); d_cal [log(rho_a); phi]; r d_obs - d_cal; J_logrho mt1d_jacobian(freq, rho_cur, h, 1); J_phi mt1d_jacobian(freq, rho_cur, h, 2); J [J_logrho; J_phi]; % 求解阻尼最小二乘方程 A J * J lambda * eye(N) alpha * (L * L); g J * r - alpha * (L * L) * m; dm A \ g; % 尝试更新动态调整lambda m_new m dm; rho_new exp(m_new); [rho_a_new, phi_new] mt1d_forward(freq, rho_new, h); fit_old norm(r); d_cal_new [log(rho_a_new); phi_new]; fit_new norm(d_obs - d_cal_new); if fit_new fit_old m m_new; lambda max(lambda * 0.3, 1e-6); else lambda min(lambda * 3, 1e6); % lambda增加后重新计算dm A J * J lambda * eye(N) alpha * (L * L); g J * r - alpha * (L * L) * m; dm A \ g; m_new m dm; [rho_a_new, phi_new] mt1d_forward(freq, exp(m_new), h); d_cal_new [log(rho_a_new); phi_new]; fit_new norm(d_obs - d_cal_new); if fit_new fit_old m m_new; lambda lambda * 0.5; end end hist.fit(end1) fit_new; if length(hist.fit) 3 abs(diff(hist.fit(end-2:end))) 1e-4 break; end end rho_inv exp(m); end这里有个细节值得一说正则化项我放在了梯度g里同时也在法方程A矩阵里加了alpha*(LL)。很多人初写反演时只在A里加正则化梯度里忘了减去正则化贡献结果迭代到后期模型被正则化拖到固定点附近但拟合差还是很大怎么调正则化参数都不对。原因是目标函数对模型参数的梯度应该包含正则化项的梯度这个梯度就是alpha(L*L)*m漏掉以后优化方向就不正确。阻尼因子λ的自适应策略也是反演能否稳定收敛的关键。我这里采用了一个简单粗暴的规则如果这次更新让拟合差下降就调小λ让算法更接近高斯牛顿法加速收敛如果拟合差反弹就调大λ让算法更接近梯度下降法提升稳定性。3.4 合成数据测试三层H型模型写完了三件套我用一个典型的H型三层模型测试整个流程。H型模型的电阻率特征是低-高-低也就是中间层是高阻。模型参数设为第一层电阻率100 Ω·m厚度200米第二层电阻率1000 Ω·m厚度800米第三层基底电阻率10 Ω·m。频率从0.001 Hz到1000 Hz按每十倍频程5个频点生成共31个频点。先用正演算出干净的响应再加2%的高斯噪声。初始模型我给了均匀半空间100 Ω·m层界面固定在设计模型的位置也就是说只反演各层电阻率。运行反演后迭代12次收敛拟合差从最初的几十降到接近噪声水平。反演得到的电阻率大约是98、1020、11和真实模型基本一致。这个测试说明整个正反演框架是正确的也验证了阻尼最小二乘法在层状电阻率情况下的有效性。你拿到代码后可以先把这一套跑通再换成自己的实测数据。4. 常见问题与排查技巧实录4.1 初始模型给不好反演发散怎么办反演发散是新手最常遇到的问题。我刚开始写这段代码时初始模型给一个极端高阻层比如10⁶ Ω·m结果第一次迭代雅可比矩阵就出现NaN程序直接崩掉。排查后发现问题有两个来源一是tanh的参数太大溢出了二是层模型里相邻层电阻率差太大导致正演响应出现数值不稳定。解决思路是给反演加防护。第一正演函数里可以检查tanh的参数规模如果abs(arg)超过某个阈值比如50就改用极限值近似因为tanh在大参数下趋近于1。第二初始模型不要选离实际观测曲线太远的模型可以先根据视电阻率曲线的低高频渐近值估一个大概的电阻率范围比如曲线最低端趋向100 Ω·m最高端趋向10 Ω·m初始模型就夹在这个范围里。第三迭代过程中每次更新完模型检查有没有非正数或者NaN有就回退这次更新并增大阻尼因子。我见过不少同行在初始模型上花太多心思其实阻尼最小二乘对初始模型没有想象中那么敏感只要你没有给到离谱的程度多迭代几次总能收敛。4.2 正则化参数α怎么选正则化参数α控制着模型光滑度和数据拟合之间的平衡。α设得太大反演出来的模型非常光滑几乎退化成均匀半空间拟合差居高不下设得太小模型会剧烈震荡过度拟合噪声产生虚假的电性分层。我常用的方法是L曲线法把α从1e-3到1e3按对数间隔取十几个值对每个α都跑一遍反演记录最终拟合差和模型粗糙度然后画一条以拟合差为横轴、模型粗糙度为纵轴的曲线曲线拐角的地方就是合适的α。这个方法在一维里很快因为每个α跑一次反演只需要几秒钟。更简单的方式是直接固定α等于一个经验值比如0.1到1之间。实测数据信噪比高的时候用0.1信噪比低的时候用1或者更大。不过这属于粗调严谨一点的报告里最好还是展示L曲线。4.3 实测数据的坑静位移与近场效应合成数据跑通了不等于实测数据也一帆风顺。野外MT数据反演时最让我头疼的有两个问题静位移和近场效应。静位移是地表局部不均匀体引起的电场畸变它会让整条视电阻率曲线在双对数坐标下上下平移而相位曲线基本不受影响。一维反演如果只拟合视电阻率静位移会直接变成浅层电阻率的假异常。处理办法是先用相位数据做一次反演确定浅层趋势再把视电阻率曲线校正后再联合反演。实际项目中静位移校正是一个完整的专题但一维阶段至少要意识到这个效应存在。近场效应指的是电磁源不是理想的平面波比如在矿山、铁路、高压线附近低频段数据受人工源干扰严重。具体表现是视电阻率曲线低频段以45度斜率急剧下降相位趋向0度。这种数据不能拿来做一维反演直接截掉受污染的低频段更靠谱。4.4 MATLAB版本与工具箱兼容性之前我说过避免依赖优化工具箱就是为了兼容性。实际项目里我就遇到过在一台机器上跑得好好的脚本换到另一台机器上因为缺工具箱而报错的情况。用纯基础MATLAB写核心算法以后这个问题基本消除了。如果你用的是R2022b或者R2021a这类版本只要不是太老的版本脚本都能跑。唯一需要留意的是有些新版本把某些函数的默认行为改了比如字符串处理、绘图函数可能导致老脚本报错。我的经验是核心计算函数尽量少用和版本相关的语法画图部分用最基础的plot、loglog、semilogx就足够了不要在画图函数里塞太多精致参数不然换版本后调样式比调算法还费时间。还有一个小技巧所有函数都用function文件形式保存而不是脚本里堆一段长代码。这样写的好处是MATLAB的JIT编译对函数文件支持更好循环计算速度明显更快更重要的是便于单元测试——单独跑正演函数和单独跑反演函数哪个环节出错一目了然。4.5 反演结果陷入了局部极小反演结果明显不合理比如电阻率出现负值虽然参数化为对数后不可能为负或者模型与已知地质资料严重矛盾这很可能是陷入了局部极小值。MT反演的目标函数高度非凸局部极小值问题难以完全避免。我的处理办法有几个一是用多个不同初始模型分别反演看最终结果是否一致如果不同初始模型收敛到同一个结果说明这个解比较稳定二是先用粗网格反演得到一个大尺度模型再加密网格作为新的初始模型继续反演这相当于多尺度策略三是在阻尼最小二乘迭代的早期人为加入一个比较大的初始λ让算法先用小步长慢慢滑入盆地再逐步减小λ实现快速收敛。5. 这套代码后续还能怎么扩展5.1 从一维走向二维与三维一维正演是基础一旦你理解清楚了Wait递推的本质扩展到二维和三维的思路也就清楚了。二维MT一般分TE模式和TM模式正演需要解偏微分方程通常用有限差分或有限元方法代码复杂度比一维高一个数量级。但一维代码里关于目标函数、正则化、阻尼策略、雅可比差分的思想可以直接搬过去。我自己在做二维反演时最重要的参考模型依然来自一维反演结果。先用MT1D得到每个测点的电性结构然后作为二维反演的初始模型收敛速度和最终模型质量都比随机初始化好得多。5.2 联合反演与外部约束MT单独反演对电阻率的分辨率有限尤其是对高阻薄层的分辨能力很差。实践中经常把MT和大地电磁测深法的其他数据类型比如Z/H倾子或者与重力、磁法数据联合反演来约束模型的浅部和深部结构。在这个一套代码基础上可以实现的第一个扩展就是带先验信息约束的MT一维反演。把已知钻孔电测井结果或者地质分层信息转成参考模型m_ref调整正则化项中m_ref的权重就能让反演结果向已知信息靠拢。这种加约束的方式比硬性固定某个参数更灵活也更容易在报告中说明约束的影响程度。5.3 性能优化建议一维正反演的计算量不大但如果你要批量处理几十上百个测点性能还是值得关注。我试过用MATLAB的parfor并行计算替代普通for循环跑雅可比矩阵层数20、频点30的情况下单测点时间从1秒降到了0.3秒左右这在批量处理时差距很明显。另外正演函数里有一个隐藏的性能瓶颈是在每个频点循环内又套层循环。对于参数固定的模型可以考虑把频率和层循环向量化但代码可读性会下降。我的建议是先把正确的版本跑通确认结果无误后再优化。反演算法这东西正确性永远比效率重要因为一个错误的快速算法会让你浪费几百个小时调试。最后再分享一点个人体会。写完这套MT1D代码以后我最大的感受是与其到处找别人封装好的黑盒工具不如自己把核心公式敲一遍。敲一遍才知道Wait递推里复数的走向才知道阻尼因子为什么要动态调整才知道正则化矩阵为什么取一阶差分而不是二阶差分。这些都是读论文体会不到的。如果你也想在地电磁数据处理这条路上走深一点我建议就拿一维正反演当第一个练手项目把代码跑通、跑懂然后再去碰二维三维到那时候你手里的工具就真的是你的了。本文还有配套的精品资源点击获取
MATLAB实现大地电磁MT一维正反演：核心算法与工程实践
简介本资源是一套面向地球物理勘探初学者与科研人员的MATLAB一维大地电磁MT正演计算工具聚焦于地质结构电性参数建模与地表电磁响应模拟适用于矿产勘查、地热调查及基础MT方法教学实践。压缩包仅含1个核心文件——mt1d.m脚本体积仅1KB代码精炼完整实现基于频率域Maxwell方程的一维层状介质正演求解支持用户自定义电导率剖面、频率范围及磁导率参数并输出视电阻率与相位曲线等关键响应数据。已有344人学习下载程序经实测验证稳定性与精度良好内置清晰注释与典型参数示例便于快速理解电磁场传播原理、调试模型设置并衔接后续反演分析。 做大地电磁的人手里多多少少都攒着几个自己写的MT脚本。我最早接触这个方向的时候也跑去翻现成的开源代码结果要么语言不合适要么处理流程绕最后干脆用MATLAB自己写了一套MT1D正反演。今天把项目里的核心算法、代码结构、调试经验整理出来给同样需要处理大地电磁数据的同行做个参考。这个项目做的事情很明确用MATLAB实现大地电磁一维正反演MT算法正演部分用层状介质的Wait递归公式反演部分用阻尼最小二乘加正则化约束能直接吃进视电阻率和相位曲线输出地电模型。对做电磁法勘探、地热调查、水文地质、油气探测的同行来说这套代码可以当成自研处理流程的起点对学生来说把代码跑通一遍也就把MT正演和反演的本质理解得差不多了比闷头读论文直观得多。整个项目写完以后我顺手把调试记录和踩坑清单也整理了一遍。这篇不光是贴代码我更想把代码背后的设计逻辑、参数为什么这么调、实测数据会遇到哪些坑以及我在不同MATLAB版本下跑出来的经验写清楚至少让你少走我走过的弯路。1. 项目到底在做什么MT1D的前因后果1.1 MT方法的基本物理图景大地电磁测深法Magnetotellurics简称MT利用天然交变电磁场作为激励源在地表观测相互正交的电场和磁场分量比如Ex和Hy或者Ey和Hx然后构造波阻抗Z Ex / Hy单位是欧姆。这个阻抗本身就是地下介质电性结构的函数频率越高电磁波在地下的传播深度越浅频率越低穿透深度越深。把一系列频率的视电阻率和相位画成曲线就得到了一条带有深度信息的电磁响应曲线。天然场的频率范围很宽从0.0001 Hz量级的超低频到上万赫兹的音频频率。频率和探测深度的大致关系可以用趋肤深度δ来估算δ约等于503乘以根号下电阻率除以频率常见的近似形式是δ ≈ 503√(ρ/f)单位是米。举个例子电阻率100 Ω·m频率1 Hz的时候趋肤深度大约5公里这意味着我们在1 Hz看到的信号主要反映地下5公里以内的电性分布。这个粗略估算贯穿了MT数据解释的整个过程从定频段选点到后期反演初始模型设置我都经常拿它来快速判断数据能分辨多深。1.2 为什么先从一维做起一维模型是MT解释里最朴素也最可靠的模型它假设地下电性只随深度变化水平方向无限延伸每一层内部各向同性且均匀。这个假设虽然在真实地质体里很难严格成立但没有一维打底二维三维根本没法开工。我自己的经验是一维正演至少有三个用途。第一用来做数据质量检查。野外实测的视电阻率和相位曲线如果跟一维响应对不上往往说明数据受三维结构、地表不均匀体或者人文电磁噪声影响比较严重这时候直接拿去做二维或三维反演很可能得到假异常。第二用来提供三维反演的初始模型。现在主流的三维反演软件虽然能自动生成初始模型但如果你想给一个更有地质意义的初始模型一维反演结果逐测点拼接是最常见的做法。第三用来做理论实验。一维反演速度快、参数少是验证算法、调参、理解反演非唯一性的最佳场地。我自己很多关于正则化参数、误差估计的直觉都是在一维环境里养成的换到二维三维以后很多原则可以平移。1.3 功能规划与代码文件结构项目没有搞复杂的工程化结构核心逻辑就是“正演、雅可比、反演、演示”四件事。文件结构大致如下文件作用关键内容mt1d_forward.mMT一维正演函数输入频率、层电阻率、层厚度输出视电阻率和相位mt1d_jacobian.m雅可比矩阵计算对数电阻率参数下的数值差分雅可比mt1d_inv.m阻尼最小二乘反演框架迭代更新、正则化、拟合差计算mt1d_demo.m合成数据演示脚本生成三层模型正演响应加噪声再反演复原模型我之所以把雅可比单独拆出来是因为反演主循环需要反复调用它。而且调试的时候单独看雅可比矩阵比混在反演里更容易发现问题比如某一行全部是0十有八九是差分步长或者参数索引写错了。整体设计有一个原则全部用基础MATLAB函数不依赖优化工具箱、统计工具箱这类需要单独授权的模块。这么做的原因很现实很多学校、企业的机器上只装了基础MATLAB有的甚至版本比较老如果你一上来就依赖lsqcurvefit或者fmincon换台机器就跑不了。手写矩阵运算虽然看起来多写几行但换来的是极强的兼容性。2. 正演与反演的核心算法拆解2.1 正演公式从Cagniard到Wait递归MT一维正演的数学基础来自Cagniard和Wait等人的工作。对于均匀半空间地表阻抗直接等于本征阻抗也就是Z √(iωμρ)其中ω是角频率μ是磁导率通常取真空磁导率4π×10⁻⁷ H/mρ是电阻率。视电阻率定义为ρa |Z|² / (ωμ)相位是阻抗辐角的一半。对于多层水平层状介质要从最底层的半空间开始逐层向上递推。每一层顶部的阻抗可以用Wait递推公式计算Z_n η_n × [Z_{n1} η_n × tanh(i k_n h_n)] / [η_n Z_{n1} × tanh(i k_n h_n)]其中η_n是第n层的内禀阻抗η_n √(iωμρ_n)k_n是复波数k_n √(iωμ/ρ_n)h_n是第n层的厚度。这个公式从最底部的半空间开始算Z_1就是地表阻抗然后代入视电阻率和相位的定义式就得到响应。这里容易被绕晕的是η和k的符号写法。不同文献里定义略有差异有的习惯用σ表示电导率写成η √(iωμ/σ)其实和ρ版本完全等价。你自己写代码的时候只要保持一致算出来的视电阻率不会差。我习惯全部用电阻率表示避免来回换算电导率出错。还有一个需要注意的点公式里是复数tanhMATLAB对复数tanh支持得很不错但如果频率太高、层太厚i k_n h_n的虚部会非常大可能出现数值溢出。实操中我一般把频率范围和层厚控制在这个公式稳定工作的范围内同时用warn状态检查结果如果出现NaN或者Inf就立刻停下来检查参数。2.2 反演问题的数学表达反演的目标是找一组模型参数使得正演响应尽量拟合实测视电阻率和相位数据。MT一维反演本质上是高度非线性的我用的是带正则化的阻尼最小二乘也叫Levenberg-Marquardt方法。把模型参数记为m通常取各层电阻率的对数也就是m log(ρ)。取对数的好处有两个一是保证迭代过程中电阻率始终为正不需要额外的约束处理二是把跨越几个数量级的电阻率参数压缩到线性尺度上梯度计算和步长控制更稳定。数据向量d则包含视电阻率的对数和相位。之所以不用视电阻率原值是因为大地电磁数据动态范围很大可能从10 Ω·m到10⁴ Ω·m如果直接拟合绝对差值高阻层的权重会淹没低阻层的信息取对数相当于做了对数均匀加权更符合人眼判断拟合效果的习惯。目标函数写为Φ ‖W(d_obs - F(m))‖² α‖L(m - m_ref)‖²第一项是数据拟合差W是数据协方差矩阵的逆用来对不同的数据点按误差大小归一化第二项是正则化项L是一阶差分矩阵m_ref是参考模型α是正则化参数。这个式子用中文说就是既要让模型响应接近观测数据又不能让模型在相邻层之间疯狂跳变。阻尼最小二乘的迭代更新公式是m_new m_old (JᵀWᵀWJ λI αLᵀL)⁻¹ JᵀWᵀW(d_obs - F(m_old))其中J是雅可比矩阵λ是阻尼因子。λ的作用是在高斯牛顿法和最速下降法之间动态切换λ大的时候步长变小、稳定性增加λ小的时候收敛速度加快。每次迭代如果拟合差下降就减小λ如果拟合差反弹就增大λ后面我会把判定策略写进代码。2.3 MATLAB实现中的几个关键细节第一是数组下标与模型深度对应。层电阻率数组rho的长度是N层厚度数组h的长度是N-1最后一层默认为半空间。为了避免混淆我在函数开头就加注释说明并且用numel(rho)-1来推断层数而不是写死。第二是频率数组的生成。实测MT数据的频点通常是等对数间隔排列的从最低频到最高频每个倍频程大概5到10个频点。生成合成数据时我习惯用logspace函数低频从0.001 Hz到1000 Hz这样视电阻率曲线高低频都能覆盖便于观察规律。第三是相位的单位。观测资料里相位一般给度数但正演计算里atan2出来是弧度所以要把弧度转成度数再放进数据向量。这个单位问题不处理好反演时会莫名出现拟合差跳来跳去的情况排查半天发现是单位混了。第四是雅可比矩阵的数值差分。虽然理论上可以推导解析导数但代码复杂度高、容易出错。我用中心差分对第j个模型参数给一个小的扰动量Δm_j分别计算F(m Δm_j)和F(m - Δm_j)再取差除以2Δm_j。差分步长取1e-4这个量级太小会引入舍入误差太大则线性近似失真实测下来1e-4比较稳定。3. 核心代码实现与合成数据测试3.1 mt1d_forward正演函数怎么写才稳正演函数是整套代码的地基反正演能不能收敛就看正演算得准不准。下面这个版本我去掉了多余的炫技部分只保留最核心的逻辑放进MATLAB就能直接运行function [rho_a, phi_deg, Z_surf] mt1d_forward(freq, rho, h) % MT1D_FORWARD MT一维正演计算层状模型的视电阻率和相位 % 输入: % freq : 频率数组单位Hz建议按对数等间隔排列 % rho : 层电阻率数组单位Ohm.m长度N最后一层为半空间 % h : 层厚度数组单位m长度N-1对应前N-1层的厚度 % 输出: % rho_a : 视电阻率单位Ohm.m % phi_deg: 相位单位度 % Z_surf : 地表波阻抗复数 mu0 4 * pi * 1e-7; omega 2 * pi * freq(:); % 转为列向量方便广播计算 N length(rho); if length(h) ~ N - 1 error(层厚度个数必须比电阻率个数少1); end Z_surf zeros(size(omega)); for ifreq 1:length(omega) w omega(ifreq); % 从最底层半空间开始 Z_bottom sqrt(1i * w * mu0 * rho(N)); % 自底向上递推 for n N-1:-1:1 eta_n sqrt(1i * w * mu0 * rho(n)); k_n sqrt(1i * w * mu0 / rho(n)); arg k_n * h(n); Z_bottom eta_n * (Z_bottom eta_n * tanh(arg)) / ... (eta_n Z_bottom * tanh(arg)); end Z_surf(ifreq) Z_bottom; end rho_a abs(Z_surf).^2 ./ (omega * mu0); phi_deg angle(Z_surf) * 180 / pi; % 一维MT相位通常位于0-90度之间angle返回-pi到pi负相位加360 phi_deg(phi_deg 0) phi_deg(phi_deg 0) 360; end这个函数最核心的就是从倒数第二层往上循环的部分每一层都套用Wait递归公式。需要强调的一点是我使用的是Z_bottom作为上一层的输入名字可能让你误解成“底部阻抗”但严格说它代表“从当前层往下看的等效阻抗”也就是Z_{n1}。这样递推到最后Z_surf就是地表阻抗。关于相位一维大地电磁的物理相位大多在第一象限也就是0到90度之间。但MATLAB的angle函数返回值范围是-180到180度如果用实部虚部算出来的相位落在第四象限就会显示成负角度需要加360度转成正度数。这种小细节商用软件里都帮你处理好了自己写代码时最容易漏。3.2 mt1d_jacobian雅可比矩阵的计算雅可比矩阵是反演迭代的核心它的每一行对应一个数据点对每个模型参数的偏导数。由于我用的是对数电阻率参数m所以差分时扰动对数电阻率正演时再把logρ还原成ρ。function J mt1d_jacobian(freq, rho, h, data_type, delta) % MT1D_JACOBIAN 计算MT一维正演对对数电阻率的雅可比矩阵 % data_type: 1代表视电阻率对数2代表相位 % delta : 中心差分步长缺省1e-4 if nargin 5 || isempty(delta) delta 1e-4; end N length(rho); J zeros(numel(freq), N); logrho log(rho); for j 1:N % 正向扰动 rho_plus rho; rho_plus(j) exp(logrho(j) delta); [rho_a_p, phi_p] mt1d_forward(freq, rho_plus, h); % 负向扰动 rho_minus rho; rho_minus(j) exp(logrho(j) - delta); [rho_a_m, phi_m] mt1d_forward(freq, rho_minus, h); % 中心差分 if data_type 1 J(:, j) (log(rho_a_p) - log(rho_a_m)) / (2 * delta); else J(:, j) (phi_p - phi_m) / (2 * delta); end end end把视电阻率取对数再差分是为了和反演目标函数里的数据项保持一致。如果你的反演数据里同时包含视电阻率和相位那就把两部分雅可比拼在一起上面这个函数实际上返回的是单一数据类型的雅可比调用时按列拼接就行。这个函数的性能并不高因为每个模型参数要做两次正演。但一维正演本身开销很小模型层数通常不超过几十层所以完全够用。如果以后扩展到二维或三维再用解析导数或者伴随方法。3.3 mt1d_inv阻尼最小二乘反演框架反演函数是整个项目的主控模块。它接受实测频率、实测视电阻率、相位以及初始模型输出反演后的模型参数和迭代历史。function [rho_inv, hist] mt1d_inv(freq, rho_obs, phi_obs, ... rho0, h, alpha, maxiter) % MT1D_INV 阻尼最小二乘MT一维反演 % rho0 : 初始电阻率模型长度N % alpha : 正则化参数 % maxiter : 最大迭代次数 N length(rho0); m log(rho0(:)); rho_obs rho_obs(:); phi_obs phi_obs(:); d_obs [log(rho_obs); phi_obs]; % 一阶差分正则化矩阵作用于模型参数 L zeros(N-1, N); for i 1:N-1 L(i, i) -1; L(i, i1) 1; end lambda 1.0; % 阻尼因子初值 hist.fit []; for iter 1:maxiter rho_cur exp(m); [rho_a, phi] mt1d_forward(freq, rho_cur, h); d_cal [log(rho_a); phi]; r d_obs - d_cal; J_logrho mt1d_jacobian(freq, rho_cur, h, 1); J_phi mt1d_jacobian(freq, rho_cur, h, 2); J [J_logrho; J_phi]; % 求解阻尼最小二乘方程 A J * J lambda * eye(N) alpha * (L * L); g J * r - alpha * (L * L) * m; dm A \ g; % 尝试更新动态调整lambda m_new m dm; rho_new exp(m_new); [rho_a_new, phi_new] mt1d_forward(freq, rho_new, h); fit_old norm(r); d_cal_new [log(rho_a_new); phi_new]; fit_new norm(d_obs - d_cal_new); if fit_new fit_old m m_new; lambda max(lambda * 0.3, 1e-6); else lambda min(lambda * 3, 1e6); % lambda增加后重新计算dm A J * J lambda * eye(N) alpha * (L * L); g J * r - alpha * (L * L) * m; dm A \ g; m_new m dm; [rho_a_new, phi_new] mt1d_forward(freq, exp(m_new), h); d_cal_new [log(rho_a_new); phi_new]; fit_new norm(d_obs - d_cal_new); if fit_new fit_old m m_new; lambda lambda * 0.5; end end hist.fit(end1) fit_new; if length(hist.fit) 3 abs(diff(hist.fit(end-2:end))) 1e-4 break; end end rho_inv exp(m); end这里有个细节值得一说正则化项我放在了梯度g里同时也在法方程A矩阵里加了alpha*(LL)。很多人初写反演时只在A里加正则化梯度里忘了减去正则化贡献结果迭代到后期模型被正则化拖到固定点附近但拟合差还是很大怎么调正则化参数都不对。原因是目标函数对模型参数的梯度应该包含正则化项的梯度这个梯度就是alpha(L*L)*m漏掉以后优化方向就不正确。阻尼因子λ的自适应策略也是反演能否稳定收敛的关键。我这里采用了一个简单粗暴的规则如果这次更新让拟合差下降就调小λ让算法更接近高斯牛顿法加速收敛如果拟合差反弹就调大λ让算法更接近梯度下降法提升稳定性。3.4 合成数据测试三层H型模型写完了三件套我用一个典型的H型三层模型测试整个流程。H型模型的电阻率特征是低-高-低也就是中间层是高阻。模型参数设为第一层电阻率100 Ω·m厚度200米第二层电阻率1000 Ω·m厚度800米第三层基底电阻率10 Ω·m。频率从0.001 Hz到1000 Hz按每十倍频程5个频点生成共31个频点。先用正演算出干净的响应再加2%的高斯噪声。初始模型我给了均匀半空间100 Ω·m层界面固定在设计模型的位置也就是说只反演各层电阻率。运行反演后迭代12次收敛拟合差从最初的几十降到接近噪声水平。反演得到的电阻率大约是98、1020、11和真实模型基本一致。这个测试说明整个正反演框架是正确的也验证了阻尼最小二乘法在层状电阻率情况下的有效性。你拿到代码后可以先把这一套跑通再换成自己的实测数据。4. 常见问题与排查技巧实录4.1 初始模型给不好反演发散怎么办反演发散是新手最常遇到的问题。我刚开始写这段代码时初始模型给一个极端高阻层比如10⁶ Ω·m结果第一次迭代雅可比矩阵就出现NaN程序直接崩掉。排查后发现问题有两个来源一是tanh的参数太大溢出了二是层模型里相邻层电阻率差太大导致正演响应出现数值不稳定。解决思路是给反演加防护。第一正演函数里可以检查tanh的参数规模如果abs(arg)超过某个阈值比如50就改用极限值近似因为tanh在大参数下趋近于1。第二初始模型不要选离实际观测曲线太远的模型可以先根据视电阻率曲线的低高频渐近值估一个大概的电阻率范围比如曲线最低端趋向100 Ω·m最高端趋向10 Ω·m初始模型就夹在这个范围里。第三迭代过程中每次更新完模型检查有没有非正数或者NaN有就回退这次更新并增大阻尼因子。我见过不少同行在初始模型上花太多心思其实阻尼最小二乘对初始模型没有想象中那么敏感只要你没有给到离谱的程度多迭代几次总能收敛。4.2 正则化参数α怎么选正则化参数α控制着模型光滑度和数据拟合之间的平衡。α设得太大反演出来的模型非常光滑几乎退化成均匀半空间拟合差居高不下设得太小模型会剧烈震荡过度拟合噪声产生虚假的电性分层。我常用的方法是L曲线法把α从1e-3到1e3按对数间隔取十几个值对每个α都跑一遍反演记录最终拟合差和模型粗糙度然后画一条以拟合差为横轴、模型粗糙度为纵轴的曲线曲线拐角的地方就是合适的α。这个方法在一维里很快因为每个α跑一次反演只需要几秒钟。更简单的方式是直接固定α等于一个经验值比如0.1到1之间。实测数据信噪比高的时候用0.1信噪比低的时候用1或者更大。不过这属于粗调严谨一点的报告里最好还是展示L曲线。4.3 实测数据的坑静位移与近场效应合成数据跑通了不等于实测数据也一帆风顺。野外MT数据反演时最让我头疼的有两个问题静位移和近场效应。静位移是地表局部不均匀体引起的电场畸变它会让整条视电阻率曲线在双对数坐标下上下平移而相位曲线基本不受影响。一维反演如果只拟合视电阻率静位移会直接变成浅层电阻率的假异常。处理办法是先用相位数据做一次反演确定浅层趋势再把视电阻率曲线校正后再联合反演。实际项目中静位移校正是一个完整的专题但一维阶段至少要意识到这个效应存在。近场效应指的是电磁源不是理想的平面波比如在矿山、铁路、高压线附近低频段数据受人工源干扰严重。具体表现是视电阻率曲线低频段以45度斜率急剧下降相位趋向0度。这种数据不能拿来做一维反演直接截掉受污染的低频段更靠谱。4.4 MATLAB版本与工具箱兼容性之前我说过避免依赖优化工具箱就是为了兼容性。实际项目里我就遇到过在一台机器上跑得好好的脚本换到另一台机器上因为缺工具箱而报错的情况。用纯基础MATLAB写核心算法以后这个问题基本消除了。如果你用的是R2022b或者R2021a这类版本只要不是太老的版本脚本都能跑。唯一需要留意的是有些新版本把某些函数的默认行为改了比如字符串处理、绘图函数可能导致老脚本报错。我的经验是核心计算函数尽量少用和版本相关的语法画图部分用最基础的plot、loglog、semilogx就足够了不要在画图函数里塞太多精致参数不然换版本后调样式比调算法还费时间。还有一个小技巧所有函数都用function文件形式保存而不是脚本里堆一段长代码。这样写的好处是MATLAB的JIT编译对函数文件支持更好循环计算速度明显更快更重要的是便于单元测试——单独跑正演函数和单独跑反演函数哪个环节出错一目了然。4.5 反演结果陷入了局部极小反演结果明显不合理比如电阻率出现负值虽然参数化为对数后不可能为负或者模型与已知地质资料严重矛盾这很可能是陷入了局部极小值。MT反演的目标函数高度非凸局部极小值问题难以完全避免。我的处理办法有几个一是用多个不同初始模型分别反演看最终结果是否一致如果不同初始模型收敛到同一个结果说明这个解比较稳定二是先用粗网格反演得到一个大尺度模型再加密网格作为新的初始模型继续反演这相当于多尺度策略三是在阻尼最小二乘迭代的早期人为加入一个比较大的初始λ让算法先用小步长慢慢滑入盆地再逐步减小λ实现快速收敛。5. 这套代码后续还能怎么扩展5.1 从一维走向二维与三维一维正演是基础一旦你理解清楚了Wait递推的本质扩展到二维和三维的思路也就清楚了。二维MT一般分TE模式和TM模式正演需要解偏微分方程通常用有限差分或有限元方法代码复杂度比一维高一个数量级。但一维代码里关于目标函数、正则化、阻尼策略、雅可比差分的思想可以直接搬过去。我自己在做二维反演时最重要的参考模型依然来自一维反演结果。先用MT1D得到每个测点的电性结构然后作为二维反演的初始模型收敛速度和最终模型质量都比随机初始化好得多。5.2 联合反演与外部约束MT单独反演对电阻率的分辨率有限尤其是对高阻薄层的分辨能力很差。实践中经常把MT和大地电磁测深法的其他数据类型比如Z/H倾子或者与重力、磁法数据联合反演来约束模型的浅部和深部结构。在这个一套代码基础上可以实现的第一个扩展就是带先验信息约束的MT一维反演。把已知钻孔电测井结果或者地质分层信息转成参考模型m_ref调整正则化项中m_ref的权重就能让反演结果向已知信息靠拢。这种加约束的方式比硬性固定某个参数更灵活也更容易在报告中说明约束的影响程度。5.3 性能优化建议一维正反演的计算量不大但如果你要批量处理几十上百个测点性能还是值得关注。我试过用MATLAB的parfor并行计算替代普通for循环跑雅可比矩阵层数20、频点30的情况下单测点时间从1秒降到了0.3秒左右这在批量处理时差距很明显。另外正演函数里有一个隐藏的性能瓶颈是在每个频点循环内又套层循环。对于参数固定的模型可以考虑把频率和层循环向量化但代码可读性会下降。我的建议是先把正确的版本跑通确认结果无误后再优化。反演算法这东西正确性永远比效率重要因为一个错误的快速算法会让你浪费几百个小时调试。最后再分享一点个人体会。写完这套MT1D代码以后我最大的感受是与其到处找别人封装好的黑盒工具不如自己把核心公式敲一遍。敲一遍才知道Wait递推里复数的走向才知道阻尼因子为什么要动态调整才知道正则化矩阵为什么取一阶差分而不是二阶差分。这些都是读论文体会不到的。如果你也想在地电磁数据处理这条路上走深一点我建议就拿一维正反演当第一个练手项目把代码跑通、跑懂然后再去碰二维三维到那时候你手里的工具就真的是你的了。本文还有配套的精品资源点击获取