简介面向大地电磁测深学习与研究人群这份Matlab程序实现了各向同性均匀多层层状介质的一维正演计算依据石应骏《大地电磁测深》教材中的解析法编写代码简单明了适合地球物理专业学生与入门研究者快速理解MT一维正演原理。资源包含4个文件其中两个.m脚本分别承担主流程与核心计算另有两个.asv自动保存备份文件整套压缩包仅2KB精炼无冗余。目前已有464人学习下载得到不少学习者的实际验证。通过学习这份代码读者既能直接运行并得到层状模型的视电阻率响应也能对照教材公式逐行理解解析法推导借助清晰的层参数输入与结果输出还可修改各层电阻率、厚度等参数快速构建自己的地电模型用于理论曲线分析、反演初始模型设计以及教学演示是兼顾学习价值与实用性的正演小工具。 我做了十几年地球物理数据处理各类正演程序算是吃饭的家伙。很多人刚接触大地电磁法MT时会碰上一个看起来特别“教科书”的任务写一个各向同性、均匀多层层状介质的一维正演程序。标题听着又长又绕实际干的事不难——把地下简化成一层一层的水平介质每层电阻率恒定、方向同性然后模拟天然电磁场在其中的传播最终算出地表会观测到的视电阻率和阻抗相位。这个小程序能解决什么问题呢反演之前的初始模型试算、野外测深方案设计、快速判断某个地电断面大概长什么样全都离不开它。这篇博文我会从原理到代码、从验证到踩坑把程序完整拆一遍适合刚入门的同学照着复现也适合做实测处理的老手拿去当标定工具。1. 一维正演到底在算什么先把这个模型在脑子里立起来1.1 从麦克斯韦方程组到“地表阻抗”大地电磁法的物理基础不复杂天然电磁波从高空照射到地表在地下介质中感应出电磁场而我们在地表观测的是电场和磁场的正交分量。这些分量满足麦克斯韦方程组对于一维水平层状介质来说波可以看成垂直入射问题退化为沿深度方向传播的一维Helmholtz方程。真正关键的是引入地表波阻抗这个概念通常写成Z E_x / H_y对应TE模式。为什么说它关键因为MT方法最终关心的视电阻率和阻抗相位都可以直接从阻抗换算出来不需要把所有深度上的电场磁场值全算出来。一维层状介质的处理思路就是利用每层介质中电磁波的传播特性把层与层之间的边界条件一层层衔接起来最终把整条剖面的响应写成层层递推的形式。我不建议一上来就盯着公式看。先建立物理图像电磁波从地表向下传播在每一层界面都会发生反射和透射地表观测到的阻抗是所有这些反射波叠加的结果。递推的过程本质上就是把这无数次的反射“压缩”成一个等效阻抗。1.2 为什么只有各向同性、均匀层状介质才有解析解标题里的限定词不是随便写的。“各向同性”意味着电阻率是标量而不是张量。碰上各向异性介质比如页岩地层常见的水平电阻率和垂直电阻率不一致的情况一维正演的公式马上就不够用了得引入额外的张量参数。而“均匀多层层状”则保证了每一层内的电阻率ρ和厚度h都是常数电磁波在每一层内部可以写成简单的解析表达式在界面处通过电场切向连续、磁场切向连续来衔接。这套逻辑成立的根本原因是层状介质本身允许写出闭式解。对比一下二三维正演那需要在空间上做网格剖分用有限差分或有限元去求解整个空间的场分布。而一维层状模型数学上就是一组复指数函数边界条件的连接问题只需要做阻抗递推就能得到答案速度快得可以忽略不计。一维正演虽然简单但它是一切更复杂正反演手段的地基。很多二三维反演程序在初始模型、约束条件、标定检验等环节都要反复调用一维正演。所以把一维正演吃透后面路会平顺很多。2. 把公式写成人话阻抗递推与视电阻率计算2.1 传播常数、固有阻抗和递推初值先约定符号。设第j层的电阻率为ρ_j厚度为h_j电磁波角频率为ω真空磁导率μ0 4π × 10⁻⁷ H/m。电磁波在第j层内的传播常数为k_j sqrt(i ω μ0 / ρ_j)这个k_j是复数实部描述衰减虚部描述相位变化。我们常说的趋肤深度δ sqrt(2ρ / (ω μ0))就和k_j直接相关电阻率越高、周期越长电磁波能穿透的深度越大这就是为什么大地电磁法能勘探深部结构。每一层还有一个固有阻抗物理上可以理解为电磁波在该层介质中传播时的“特征阻抗”Z_0j i ω μ0 / k_j对于最底层它被当成向下无限延伸的半空间。在底部没有界面反射所以该层的上表面阻抗就直接等于固有阻抗这成为递推的初值Z_N Z_0N这一条是整个递推的起点也是我建议所有人在写代码时最先确认的地方。2.2 从最深一层往上推递推公式与响应换算确定了底层阻抗后逐层向上递推。第j层上表面的阻抗由下层传递上来的阻抗Z_{j1}和本层参数共同决定Z_j Z_0j × (Z_{j1} Z_0j × tanh(i k_j h_j)) / (Z_0j Z_{j1} × tanh(i k_j h_j))物理上可以这样理解括号里的每一项都代表电磁波在层内往返传播后的阻抗叠加tanh函数就是描述反射波与透射波“叠加”的数学工具。当h_j远大于趋肤深度时tanh趋近于1下层的影响被屏蔽掉上表面阻抗逼近本层固有阻抗当h_j很薄时tanh的幅值很小递推公式平稳地将下层阻抗传递上来。从底层一直递推到地表得到Z_1就是地表波阻抗。之后换算视电阻率和相位ρ_a |Z_1|² / (ω μ0)φ arctan(Im(Z_1) / Re(Z_1))这里的相位需要注意符号约定。MT行业习惯上按阻抗相位的辐角来显示有的程序输出正值表示相位超前有的输出负值。我个人在程序里统一用arctan2取辐角再换算成角度制具体正负看设定的右手定则关键是整套流程保持一致。实际资料处理中如果相位曲线整体反号通常不是介质模型问题而是坐标轴方向约定错了。还有一个换算容易踩坑MT习惯用周期T而不是频率f来描述测深点角频率ω 2π / T。比如周期1秒对应角频率约6.28 rad/s没见过直接把1/T代入公式算的——那样算出来的结果会整体偏移一个2π因子曲线形态对但数值完全不对。3. 写一套能直接跑的程序核心实现与验证3.1 模型参数与频段怎么选一维正演的输入参数非常精简每一层的电阻率单位Ω·m、每一层的厚度单位m、要做正演的周期序列单位s。在野外采集和理论模拟中周期范围一般取0.001秒到1000秒覆盖高频浅部和低频深部。实际设计频点时建议采用对数等间隔采样每个数量级取4到8个点太密会拖慢计算太稀则曲线形态不够光滑。层数和厚度的选择取决于你要模拟的对象。比如想模拟一个典型的低阻夹层H型地电断面表层高阻覆盖、中间低阻层、底部高阻基底可以用三层模型第1层电阻率100 Ω·m厚度100 m第2层电阻率10 Ω·m厚度200 m第3层电阻率1000 Ω·m半无限这类模型在实际生产中很常见对应的是浅部风化层、中部含水层或破碎带、深部完整基底的结构。3.2 核心代码Python实现一维正演直接把正演函数封装好方便反复调用。我用numpy实现输入输出都走数组批量计算周期序列很顺手import numpy as np def mt1d(resistivities, thicknesses, periods, mu04 * np.pi * 1e-7): 大地电磁各向同性均匀多层层状介质一维正演。 参数 ---- resistivities : array_like 各层电阻率单位 ohm.m按地表向下的顺序排列。 thicknesses : array_like 各层厚度单位 m长度比电阻率少一层即可最后底层半空间可不传。 periods : array_like 周期序列单位 s。 返回 ---- rho_app : ndarray 视电阻率单位 ohm.m。 phase : ndarray 阻抗相位单位度。 rho np.asarray(resistivities, dtypefloat) h np.asarray(thicknesses, dtypefloat) T np.asarray(periods, dtypefloat) nlayers len(rho) omega 2.0 * np.pi / T rho_app np.zeros_like(T) phase np.zeros_like(T) for i, w in enumerate(omega): k np.sqrt(1j * w * mu0 / rho) # 每层传播常数 z0 1j * w * mu0 / k # 每层固有阻抗 z_below z0[-1] # 底层半空间阻抗 for j in range(nlayers - 2, -1, -1): t np.tanh(1j * k[j] * h[j]) z_below z0[j] * (z_below z0[j] * t) / (z0[j] z_below * t) rho_app[i] np.abs(z_below) ** 2 / (w * mu0) phase[i] np.degrees(np.angle(z_below)) return rho_app, phase使用示例periods np.logspace(-3, 3, 49) rho_app, phase mt1d([100, 10, 1000], [100, 200], periods)注意代码里没有对层数做过多限制理论上可以传任意多层。实测下来几十层的模型算几千个频点也就几毫秒性能完全不用操心。循环写成频率在最外层是因为阻抗递推本身有严格的顺序依赖不适合直接对层方向做向量化强行向量化反而容易把逻辑绕混。3.3 自检清单怎么确认程序没写错写归写不验证等于白写。我的习惯是每步只加一点复杂度逐层验证半空间模型只有一层时程序输出的视电阻率应该是一条水平直线数值等于输入的电阻率相位恒定在45度按默认符号约定。两层模型高频端视电阻率趋近第一层电阻率低频端趋近第二层电阻率中间是平滑过渡的S形曲线。三层H型模型高频端接近第一层随周期增大视电阻率先下降出现极小值再回升向基底电阻率靠近。相位在极小值低频一侧会明显抬升。用上面那个三层模型实测高频端1000Hz附近视电阻率约在105 Ω·m上下浮动周期增大到几秒时视电阻率达到谷底只有几欧姆米随后回升在1000秒附近逼近1000 Ω·m。相位从高频的接近45度经过低阻层区域时抬到60到70度再随周期增大回到接近45度。如果你把自检模型的输出和这个趋势对不上优先怀疑递推方向或符号约定。4. 这些年踩过的坑常见问题与排查经验4.1 异常结果速查表我把实际使用中遇到过的反常情况整理成一张表按“症状-原因-处理方式”排列排查时对照着找比从头看公式快得多。表格症状常见原因处理思路输出NaN或inf电阻率传了0、负值厚度传了空值检查输入参数合法性电阻率必须为正数厚度必须大于0曲线形态对但数值整体偏移角频率算错把周期当频率直接用统一用ω 2π / T相位符号全部反号阻抗或坐标方向约定不一致检查E和H的方向定义保持符号约定统一高频段曲线剧烈震荡表层厚度太薄或者频点过密趋肤深度远小于层厚度tanh数值趋于饱和增加最浅层厚度或减少高频段采样密度结果不随周期变化递推初值写错底层阻抗用了0或无穷确认Z_N Z_0N从最底层开始再推一遍低频端没有趋近基底电阻率厚度列表长度比电阻率多传了或是层数循环范围写错检查thicknesses长度应为len(resistivities)-14.2 几个值得反复强调的细节复数tanh的数值稳定性问题值得单独说。当某层很厚或电阻率很低时i k_j h_j的实部很大tanh趋近于1当层很薄时tanh趋近于小量。直接使用numpy的np.tanh是复数值稳定的但如果你为了“看清公式”手写e指数展开的版本就可能出现指数项溢出。我第一次写的时候吃过这个亏高频段曲线像锯齿一样排查半天发现是e的指数项超过了浮点上限。另一个高频踩坑点是厚度列表的长度。底层是半空间不需要也不可能给定厚度。代码里循环是从nlayers-2倒推也就是第1层到倒数第2层这个边界条件对应着厚度列表只传前nlayers-1层。如果你习惯把最后一层厚度传成0或者一个很大的数结果会出现虚假的低频异常。我的建议是接口定义里强制要求厚度列表比电阻率列表少一个元素传入即报错从源头杜绝这类问题。相位换算还有个细节np.angle返回的辐角范围在(-π, π]之间如果你需要给结果做插值或绘图碰到相位跨越±180度边界时会出现跳变。虽然一维正演本身很少碰到这种情况但后续做反演目标函数时要小心最好先做相位解缠或用连续化的表示方法。最后再分享两个实用经验根据我的实际经验这个程序最大的价值不只是“算出一条曲线”而是帮你建立地电模型和远区响应之间的直觉。我建议你写一个批量扫参脚本把中间层电阻率从1扫描到10000厚度从10扫描到1000所有曲线叠在一张图上看。多扫几组以后看到实测曲线时脑子里能立刻浮现出大概的地下结构。另外一个小技巧把正演函数封装好之后顺手给它加一个辅助函数能输出各频点的趋肤深度以及等效探测深度估算值。这样设计的野外测深方案时能快速判断当前频段到底能探测到多深的目标层特别实用。这个一维正演程序很小但做MT的人几乎天天都要用到它。本文还有配套的精品资源点击获取