简介这份资源是面向地球物理电磁法学习者的MATLAB一维正演程序聚焦大地电磁各向同性均匀多层层状介质模型依据石应骏《大地电磁测深》教材中的解析法编写适合初学者理解正演原理、验证算法实现也可作为课程作业与科研入门的参考代码。压缩包共4个文件以m脚本与asv自动备份文件为主前者承载正演计算主逻辑后者为MATLAB编辑过程留存的备份整体仅约2KB体量轻巧、便于直接阅读与二次修改。目前已有468人学习下载说明其在同类教学代码中具有一定参考价值。代码结构简单明了读者可从中获取层状介质视电阻率与阻抗的计算流程、解析公式的编程落地方式以及正演结果的组织思路便于对照教材公式逐行理解并快速上手调试。1. 从一张视电阻率曲线说起这个一维正演程序到底能算什么做大地电磁MT数据处理的人大概都遇到过这种场景野外采回来一批测点视电阻率和相位曲线画出来形态看着像层状介质响应但真要判断地下到底是三层还是五层、每层电阻率大概多少光靠眼睛盯曲线是不够的。这时候需要的是一个能快速算出理论响应的正演工具——给定层数、每层电阻率和厚度算出对应频率下的视电阻率与相位再拿去和实测曲线比对。这份基于 MATLAB 的大地电磁各向同性均匀多层层状介质一维正演程序解决的就是这件事。它面向的是地球物理专业的学生、做 MT 数据反演前需要构造初始模型的科研人员以及想把正演环节嵌进自己处理流程的开发者。程序本身不复杂但它是整个 MT 反演链条的起点模型响应算得对不对后面反演出来的东西就全跟着跑偏。2. 程序结构与核心算法从麦克斯韦方程到递推阻抗2.1 一维层状介质的物理模型与递推关系大地电磁一维正演的基本假设是地下由若干水平层组成每层内部电阻率均匀各向同性层与层之间电阻率突变。在这种假设下电磁场在每一层内满足亥姆霍兹方程解的形式是指数函数组合。关键的处理技巧是把每一层顶面的波阻抗和底面波阻抗用一个递推关系连起来从最底层往上逐层推算最终得到地表波阻抗再由地表波阻抗算出视电阻率和相位。这个递推关系的核心是所谓的阻抗递推公式。对于第 j 层已知其厚度 h_j、电阻率 ρ_j以及从下方传上来的底阻抗 Z_{j1}该层顶阻抗 Z_j 可以写成Z_j Z_{0j} * (Z_{j1} Z_{0j} * tanh(k_j * h_j)) / (Z_{0j} Z_{j1} * tanh(k_j * h_j))其中 Z_{0j} 是该层的特征阻抗k_j 是复波数。最底层第 N 层因为向下无限延伸没有从下方来的反射所以底阻抗直接取该层的特征阻抗。整个计算从底层开始逐层向上递推到地表。在 MATLAB 里实现这个递推最自然的做法是用一个从 N 到 1 的循环。下面这段代码展示了核心递推过程function [Z_surface] mt1d_impedance(rho, h, freq, mu) % rho: 各层电阻率向量 (Ohm*m)长度 N % h: 各层厚度向量 (m)长度 N最后一层厚度任意半空间 % freq: 频率向量 (Hz) % mu: 磁导率 (H/m)通常取 4*pi*1e-7 N length(rho); omega 2 * pi * freq; Z_surface zeros(size(freq)); for ifreq 1:length(freq) w omega(ifreq); % 从最底层开始底阻抗取该层特征阻抗 k_N sqrt(-1i * w * mu / rho(N)); Z sqrt(-1i * w * mu * rho(N)); % 底层特征阻抗 % 从 N-1 层向上递推到第 1 层 for j N-1:-1:1 k_j sqrt(-1i * w * mu / rho(j)); Z0_j sqrt(-1i * w * mu * rho(j)); tanh_kh tanh(k_j * h(j)); Z Z0_j * (Z Z0_j * tanh_kh) / (Z0_j Z * tanh_kh); end Z_surface(ifreq) Z; end end这段代码的逻辑很直接外层循环遍历频率内层循环从最底层往上递推。k_j是复波数Z0_j是特征阻抗tanh(k_j * h_j)体现了层的厚度效应。注意sqrt(-1i * w * mu / rho(j))里用了-1i这是因为在 MT 的符号约定下波数取负虚部保证随深度衰减。如果符号搞反算出来的相位会整体偏移 180 度这是新手最容易翻车的地方之一。参数方面rho和h的长度必须一致最后一层的h虽然不参与递推因为半空间无限厚但为了向量对齐还是给一个值。mu一般取真空磁导率除非专门研究磁性层。频率向量freq决定了输出曲线的采样密度通常按对数等间隔取比如 0.001 到 1000 Hz 之间取 40 个点。2.2 视电阻率与相位的换算拿到地表波阻抗 Z_surface 之后视电阻率和相位的计算就是标准公式function [rho_a, phase] impedance_to_response(Z, freq, mu) % 由地表波阻抗计算视电阻率和相位 omega 2 * pi * freq; rho_a abs(Z).^2 ./ (omega * mu); phase angle(Z) * 180 / pi; endrho_a的单位是欧姆米phase的单位是度。这里abs(Z).^2取的是阻抗模的平方angle(Z)返回的是弧度乘180/pi转成度。有一点需要注意不同文献对相位的符号约定可能不同有的定义相位为atan2(imag(Z), real(Z))再取负有的直接取angle(Z)。这份程序采用后者和大多数 MT 教材一致。如果你拿它和别的软件对比发现相位差个负号先检查约定别急着改代码。2.3 主脚本的组织方式一个完整的正演程序除了核心递推函数还需要一个主脚本来组织输入、调用函数、画图。典型的组织方式是这样的% main_mt1d.m clear; clc; % 定义三层模型第一层 100 Ohm*m 厚 500 m第二层 10 Ohm*m 厚 2000 m底层 1000 Ohm*m rho [100, 10, 1000]; h [500, 2000, 0]; % 最后一层厚度给 0 占位 % 频率范围0.001 到 1000 Hz对数等间隔 40 点 freq logspace(-3, 3, 40); mu 4 * pi * 1e-7; % 计算地表阻抗再换算视电阻率和相位 Z mt1d_impedance(rho, h, freq, mu); [rho_a, phase] impedance_to_response(Z, freq, mu); % 双对数坐标画视电阻率线性坐标画相位 figure; subplot(2,1,1); loglog(freq, rho_a, b-, LineWidth, 1.5); xlabel(Frequency (Hz)); ylabel(Apparent Resistivity (\Omega\cdot m)); grid on; subplot(2,1,2); semilogx(freq, phase, r-, LineWidth, 1.5); xlabel(Frequency (Hz)); ylabel(Phase (deg)); grid on;主脚本把模型参数、频率设置、函数调用和绘图串起来。logspace(-3, 3, 40)生成 0.001 到 1000 Hz 的 40 个对数等间隔频率点这是 MT 正演里最常见的频率范围。画图时视电阻率用双对数坐标相位用半对数坐标这也是行业惯例。如果你要批量算多个模型把主脚本改成循环遍历模型参数即可核心函数不用动。3. 跑起来环境准备、参数设置与结果验证3.1 MATLAB 环境与文件组织这份程序对 MATLAB 版本要求不高R2015a 之后的版本基本都能跑因为它只用了基础数学函数和绘图函数没有依赖任何工具箱。如果你用的是 Octave大部分代码也能直接运行但logspace和semilogx的行为在 Octave 里略有差异画图部分可能需要微调。文件组织上建议把核心函数和主脚本分开放在同一个目录下mt1d_forward/ ├── main_mt1d.m % 主脚本定义模型并调用计算 ├── mt1d_impedance.m % 阻抗递推核心函数 ├── impedance_to_response.m % 视电阻率与相位换算 └── plot_mt_response.m % 绘图辅助函数可选把函数单独存成.m文件的好处是可以在其他脚本里复用。比如你后面要做反演反演循环里反复调用mt1d_impedance直接函数调用比复制粘贴代码干净得多。MATLAB 要求函数文件名和函数名一致所以mt1d_impedance.m里定义的函数必须叫mt1d_impedance否则调用会报错。3.2 模型参数怎么设从两层到五层参数设置是这份程序最需要经验的地方。层数 N 决定了rho和h向量的长度两层模型就是rho [rho1, rho2]、h [h1, 0]三层就是三个元素依此类推。厚度向量的最后一个元素不参与计算但必须存在否则向量长度不匹配会报错。一个常见的误区是把厚度单位搞混。程序里厚度用米电阻率用欧姆米频率用赫兹磁导率用 H/m。如果你从别的资料里拿到的厚度是千米记得乘 1000。另一个误区是频率范围设得太窄导致曲线只显示了一部分特征。比如一个四层模型中间有低阻层如果频率下限只到 0.1 Hz可能看不到深部高阻层的响应。我一般会先把频率范围放宽到 0.0001 到 10000 Hz 跑一遍看看整体形态再缩到实际需要的范围。下面是一个五层模型的参数示例% 五层模型表层低阻中间两个低阻层底部高阻 rho [50, 200, 5, 500, 1000]; h [200, 800, 1500, 3000, 0]; freq logspace(-4, 4, 60);这个模型里第二层是高阻第三层是低阻第四层又变高阻底层最高。跑出来的视电阻率曲线会呈现多个拐点相位曲线也会有对应的极值。层数越多曲线形态越复杂参数之间的相互影响也越强这也是为什么反演时通常从简单模型开始逐步加层。3.3 结果验证和解析解、极限情况对比程序跑通之后第一件事是验证结果对不对。最直接的办法是拿两层模型和解析解对比。两层模型的视电阻率有闭式解虽然表达式有点长但可以查手册或者用已知的数值表对照。另一个办法是看极限情况当频率趋于零时视电阻率应该趋近于底层电阻率当频率趋于无穷时视电阻率应该趋近于第一层电阻率。相位在低频和高频都应该趋近于零对均匀半空间而言。还有一个实用的验证手段把模型设成均匀半空间即rho [100, 100]、h [任意, 0]这时候算出来的视电阻率应该在所有频率上都等于 100相位应该恒为零。如果相位不是零说明符号约定或者公式实现有问题。这个测试我每次改完代码都会跑一遍算是血泪经验——有一次我把tanh的参数写成了k_j * h(j) * 2结果均匀半空间算出来相位不为零查了半天才发现是手滑多打了个 2。提示验证时优先用均匀半空间和两层模型这两种情况有明确的预期结果比直接上复杂模型更容易定位问题。4. 避坑与排查符号约定、数值溢出与频率采样4.1 相位符号反了曲线整体偏移 180 度现象算出来的相位曲线在低频段接近 180 度而不是 0 度或者视电阻率曲线形态对但相位完全对不上。原因MT 里波数的符号约定有两种一种取k sqrt(-i * omega * mu / rho)另一种取k sqrt(i * omega * mu / rho)。两种约定下阻抗的相位差一个负号。如果递推公式里的tanh参数和特征阻抗的符号不匹配相位就会整体偏移。解决统一用-1i的约定即k sqrt(-1i * omega * mu / rho)、Z0 sqrt(-1i * omega * mu * rho)。改完之后用均匀半空间测试相位应该恒为零。如果还是不对检查angle(Z)之后有没有额外取负。4.2 高频段数值溢出tanh 算出 NaN现象频率较高时tanh(k_j * h_j)返回 NaN 或者 Inf导致整个阻抗计算失败。原因高频时k_j的实部和虚部都很大k_j * h_j的实部可能达到几十甚至上百tanh在实部很大时趋近于 1本身不会溢出但如果h_j设得特别大比如几万米而频率又高k_j * h_j的实部可能超过 MATLAB 双精度的表示范围导致 Inf。另一种情况是电阻率设得极小比如 0.001 欧姆米k_j变得很大同样会溢出。解决检查厚度和电阻率的量级是否合理。实际 MT 里电阻率一般在 1 到 10000 欧姆米之间厚度从几十米到几千米。如果确实需要算极端参数可以在tanh之前判断real(k_j * h_j)是否大于某个阈值比如 20大于就直接取 1避免数值问题。4.3 频率采样太稀曲线拐点被漏掉现象画出来的视电阻率曲线看起来平滑但和理论形态对比发现某个拐点没出来或者相位极值偏小。原因logspace的点数太少频率间隔太大导致曲线在变化剧烈的频段采样不足。层状模型的响应在层参数对应的特征频率附近变化最快如果采样点正好跨过这个区域就会漏掉细节。解决把频率点数从 40 增加到 60 或 80或者在对数坐标下检查相邻频率点的比值是否超过 1.5。我一般先用 100 个点跑一遍看整体形态确认没问题后再降到 40 个点用于批量计算。批量反演时点数太多会拖慢速度但正演阶段多花点时间把曲线算准是值得的。4.4 厚度向量最后一个元素给错报错或结果异常现象运行时报「矩阵维度不匹配」或者「索引超出范围」或者结果明显不对但没报错。原因h向量的长度和rho不一致或者最后一层的厚度给了非零值但代码里又把它当半空间处理。有些实现会在递推时跳过最后一层但如果h的长度不对循环索引就会出错。解决强制检查length(rho) length(h)在函数开头加一句assert。最后一层的h给 0 或者任意值都行但必须存在。如果想让代码更健壮可以在函数里判断j N时直接取特征阻抗不参与tanh计算。5. 二次开发把正演嵌进反演流程与批量计算5.1 从单模型到批量计算向量化与并行单模型正演跑通之后下一步往往是批量计算——比如为反演生成雅可比矩阵或者做模型参数扫描。最直接的做法是外面套一层循环% 批量计算多个模型的响应 models {[100,10,1000], [50,200,5,500,1000], [200,20,2000]}; h_list {[500,2000,0], [200,800,1500,3000,0], [300,1000,0]}; freq logspace(-3, 3, 40); mu 4*pi*1e-7; responses cell(length(models), 1); for k 1:length(models) Z mt1d_impedance(models{k}, h_list{k}, freq, mu); [rho_a, phase] impedance_to_response(Z, freq, mu); responses{k} struct(rho_a, rho_a, phase, phase); end如果模型数量多可以用parfor替换for做并行。但要注意parfor里不能有依赖循环顺序的变量每个迭代必须独立。上面的代码里每个模型的计算互不依赖可以直接改成parfor。并行之前记得先parpool开池子否则parfor会退化成串行。5.2 嵌进反演正演函数的接口设计反演的核心是不断调整模型参数让正演响应逼近实测数据。所以正演函数最好设计成「参数向量进、响应向量出」的形式方便优化算法调用。一个常见的接口是这样的function [d_pred] forward_model(m, freq, mu) % m: 模型参数向量前 N 个是各层电阻率后 N-1 个是各层厚度 % 返回 d_pred: 拼接的响应向量 [log10(rho_a); phase] N (length(m) 1) / 2; rho m(1:N); h [m(N1:end), 0]; Z mt1d_impedance(rho, h, freq, mu); [rho_a, phase] impedance_to_response(Z, freq, mu); d_pred [log10(rho_a(:)); phase(:)]; end这个接口把模型参数压成一个向量电阻率在前、厚度在后返回的响应也压成一个向量视电阻率取对数、相位保持原样。这样优化算法比如lsqnonlin或自己写的 Gauss-Newton可以直接调用不用关心内部结构。注意log10(rho_a)是为了让视电阻率和相位在同一个量级上避免优化时被大数主导。5.3 改造技巧加一层缓存避免重复计算反演迭代时如果某次迭代的模型参数和之前算过的完全一样重复计算就浪费了。可以在正演函数外面包一层缓存function d_pred forward_cached(m, freq, mu, cache) key num2str(m(:), %.6f ); if isKey(cache, key) d_pred cache(key); return; end d_pred forward_model(m, freq, mu); cache(key) d_pred; end用containers.Map做缓存键是模型参数的字符串表示。这个技巧在参数扫描或者网格搜索时特别有用能省掉大量重复计算。不过要注意浮点数的字符串化精度%.6f可能不够参数变化很小时会误判为同一个模型我一般用%.10f。5.4 验证改造后的代码回归测试不能省每次改完正演函数我都会拿几个固定模型跑一遍和改动前的输出对比。最省事的办法是把均匀半空间、两层、三层模型的响应存成.mat文件改完代码后重新算一遍用max(abs(d_new - d_old))检查差异。如果差异超过 1e-10说明改动影响了数值结果得回去查。这个习惯帮我拦下过好几次「以为只是重构、其实改了公式」的事故。从那以后我每次动核心递推代码都强制走一遍回归测试哪怕只是改了个变量名。希望这套正演程序能帮你把 MT 反演的第一步走稳。本文还有配套的精品资源点击获取