简介根据石应骏《大地电磁测深》教材解析法编写的大地电磁各向同性均匀多层层状介质一维正演程序面向地球物理、勘探工程等专业本科生及初入MT正演领域的研究者可用于快速计算任意层数水平层状模型的理论视电阻率与阻抗相位响应辅助理解一维大地电磁正演原理与编程实现。资源压缩包共4个文件以2个MATLAB脚本.m为主另有2个.asv备份文件整体大小仅2KB代码精简、无冗余依赖适合直接阅读和二次修改核心计算与调用逻辑分离便于逐段理解递推过程。目前已有464人学习下载尤其适合希望从教材公式过渡到实际数值计算的学习者。通过该程序可掌握多层介质中电磁场递推解析法的编程思路结合简单的输入参数即可得到不同地电断面的正演结果并可通过修改层参数快速对比不同模型响应为进一步开展反演或更复杂正演研究打下基础。 搞大地电磁MT的人不管以后做正演还是反演都绕不开一维正演这个基本功。我最早写“大地电磁各向同性均匀多层层状介质一维正演程序”是因为手里的实测曲线特别乱需要在现场快速判断测点下方大致几层、每层大概什么电阻率拿几组模型响应和实测曲线叠着看。这个本来只是“临时工具”的程序后来被我反复打磨变成了做反演、设计观测方案、给新人解释MT响应都离不开的东西。这篇文章不谈教材里的公式推导怎么严谨只聊一个能直接复现的版本递推公式从哪来、代码怎么写、哪些细节会让你的结果和参考值对不上以及我踩过的几个坑。适合刚开始写MT正演程序的学生、需要自研反演内核的工程师以及任何想验证别人正演代码是否靠谱的人。1. 一维正演在MT解释流程中的位置和能力边界1.1 输入输出与核心用途这个程序解决的是这样一类问题假设地下介质由n层水平均匀层组成每层电阻率恒定且各向同性已知各层电阻率和厚度计算在地表观测到的视电阻率ρ_a和阻抗相位φ随频率或周期的变化关系。输出就是常说的“正演曲线”它是一切后续解释工作的参照系。我把它当做一个“最小可用”工具在三个场景里反复使用。第一个场景是现场快速判断实测曲线出来后先用简单的层状模型去拟合曲线拐点能迅速估计基底埋深和高低阻层的数量级比直接上二维反演高效得多。第二个场景是作为反演内核任何一维反演算法都需要成千上万次正演调用正演必须可靠而且快。第三个场景是验证别的软件二维三维程序能不能退化成层状介质解拿一维结果去对立刻见分晓。还有个很容易被忽略的作用给新人讲“为什么会出现极值”“为什么相位会低于45度”这类问题时一维正演是最直观的教具。它能让你在一个小时之内理解MT响应的基本物理图像比直接看二维剖面上的花样要清晰得多。1.2 为什么一维没有被二维三维替代可能有人会说现在二维三维电磁正演已经很成熟了一维正演还有必要自己写吗我的观点是不仅必要而且是一切的起点。二维三维反演通常需要一维模型来构建初始模型解释人员也需要一维正演来理解曲线的基本形态。三维正演的代码里最简单的测试就是“把三维网格整体退化为层状介质”而检验退化是否正确靠的还是一维正演。当然也要讲清楚边界。一维正演只适用于地下横向变化不大、测点间距远大于构造尺度的情况。当地形起伏、断层、块状异常体占主导时一维结果只能作为参考不能直接下结论。另外在层状各向同性假设下一维正演结果对TE和TM模式没有极化差异因为垂直入射平面波在两个主极化方向上对应的是同一个表面阻抗。实际资料里如果发现两种模式曲线差异明显那本身就是横向非均匀或各向异性介质存在的信号一维解释框架就解释不动了。这个能力边界在解释时一定要清楚。2. 从Maxwell方程到递归公式的核心推导2.1 模型假设与MT平面波前提一维正演建立在“地表入射波近似为垂直入射平面波”的基础上。天然电磁场源在电离层和远方的雷暴传播到地表时曲率已经很小在测点局部可以当作平面波处理。将地下按电性差异分成若干水平层每一层内部电阻率恒定、各向同性界面无限延伸这就是标题里说的“各向同性均匀多层层状介质”。在导电介质中控制方程不是普通波动方程而是扩散型方程。取垂直向下为正z轴在某个角频率ω下电场满足二阶常微分方程解的形式包含沿z方向衰减的波和反射波。这里不需要去解完整的全域解只需要关心地表处电场和磁场的比值也就是表面阻抗。2.2 均匀半空间阻抗递推的锚点先看最简单的情况整个地下是均匀半空间电阻率为ρ。此时电场和磁场通过表面阻抗Z0关联。定义Z0 E_x / H_y并取时间因子为e^{iωt}的约定本征阻抗为Z0 sqrt(i ω μ ρ)对应的视电阻率就是ρ_a |Z0|² / (ω μ) ρ均匀半空间相位为45度。这个结果简单到近乎废话但它有两个大用处一是标定代码是否正确二是作为递推的起始值。只有半空间情况下波场通解中只有一个衰减方向才能直接写出表面阻抗不需要处理反射波。2.3 递推公式与时间因子约定对第j层设电阻率ρ_j、厚度h_j最后一层为半空间不设厚度。本层的传播常数为k_j sqrt(i ω μ / ρ_j)本征阻抗为Z0j sqrt(i ω μ ρ_j)给定下层顶面阻抗Z_{j1}本层顶面阻抗由阻抗传输公式求出Z_j Z0j * (Z_{j1} Z0j * tanh(k_j h_j)) / (Z0j Z_{j1} * tanh(k_j h_j))从最底层开始向上逐层递推最终得到地表阻抗Z_1进而计算所有频率点的视电阻率和相位。这里有一个几乎所有参考资料都会出现的“坑”时间因子用e^{iωt}还是e^{-iωt}会让传播常数前面差一个i公式可能写成tanh(k_j h_j)也可能写成tanh(i k_j h_j)。两种写法在各自约定下都是对的但混着用会让高频极限完全错误。我的建议是全程序固定统一约定并把这个约定写进注释别靠脑子记。3. 程序实现中的关键细节与完整代码3.1 一份能直接跑的Python实现下面是我日常使用的精简版函数。只做计算不做绘图。import numpy as np def mt1d(freqs, rho, h): 大地电磁各向同性均匀多层层状介质一维正演 参数 ---- freqs : array_like 频率数组单位 Hz rho : array_like 各层电阻率单位 Ω·m长度 n h : array_like 各层厚度单位 m长度 n-1最后一层不设厚度 返回 ---- rho_a : ndarray 视电阻率单位 Ω·m phase : ndarray 阻抗相位单位 deg mu0 4 * np.pi * 1e-7 omega 2 * np.pi * np.asarray(freqs, dtypenp.float64) # 兼容标量频率输入 if omega.ndim 0: omega omega.reshape(1) # 半空间的地表阻抗采用 e^{iωt} 约定 Z np.sqrt(1j * omega * mu0 * rho[-1]) # 从倒数第二层向上递推 for j in range(len(rho) - 2, -1, -1): Z0j np.sqrt(1j * omega * mu0 * rho[j]) kh np.sqrt(1j * omega * mu0 / rho[j]) * h[j] t np.tanh(kh) Z Z0j * (Z Z0j * t) / (Z0j Z * t) rho_a np.abs(Z) ** 2 / (omega * mu0) phase np.rad2deg(np.angle(Z)) return rho_a, phase这段代码很简洁但包含了所有核心逻辑频率数组向量化、最底层阻抗初始化、从倒数第二层向上递推、最后返回视电阻率和相位。rho和h分别按“欧姆米”和“米”传入这是最容易出单位错误的地方。3.2 为什么从最底层向上递推初始阻抗来自半空间本质是因为半空间里没有来自下方的反射波电场的通解只有衰减波一支阻抗直接等于介质本征阻抗。从下方开始递推每向上走一层都把下一层的阻抗“折算”到当前层顶面。这个过程和传输线理论里的输入阻抗计算完全同构理解了这一点公式就不会记错方向。能不能反过来从地表往下推地表阻抗正是我们要求的未知数如果从地表往下等于手里只有一个“待求量”却要以它为起点物理上不自然。从底层向上是因为底层边界条件是已知的这是递推能成立的根本原因。3.3 高频数值稳定性、复数精度与频率范围高频段当k_j h_j的实部很大时tanh会趋近于1递推公式自然退化到Z0j意思是趋肤深度远小于第一层厚度高频只能看到第一层。数值上直接用numpy.tanh不会溢出但有些人为了省事自己写(e^{2x}-1)/(e^{2x}1)这种展开形式x一大就会溢出或得到错误结果。这个坑我见过不止一次用库函数就行。低频段k_j h_j趋于0tanh近似等于k_j h_j曲线会逐渐逼近最后一层电阻率。当模型包含高阻厚层且频率低到0.0001 Hz以下时ω μ 这个因子变得非常小视电阻率计算对复数精度比较敏感。普通双精度复数complex128在绝大多数场景下够用但如果你做极低频的深探测研究建议先拿均匀半空间模型校验过精度再信任结果。频率范围选择要和目标深度匹配。用趋肤深度公式快速估算δ ≈ 503 * sqrt(ρ / f)单位是米。假设第一层电阻率100 Ω·mf1 Hz时δ约5030 m如果想探测到20 km深频率得到0.001 Hz量级。选频段之前算一笔这种估算账能避免算完发现曲线什么都没反映的尴尬。3.4 接口设计的小习惯函数只做计算、不做绘图。返回的相位用角度制单位明确写在docstring里。绘图、写文件、反演调用都在函数外面做这样反演循环里调用几千次也不会有I/O开销。另一个习惯是把时间因子约定和阻抗定义注释写在开头三个月后回来看代码会发现这个注释非常救命。4. 验证程序的三种方法缺一不可4.1 均匀半空间测试把rho设为[100]h设为空列表或长度0数组跑一遍全频段。理论结果是视电阻率恒等于100、相位恒等于45度。如果程序返回的不是这个结果说明阻抗定义、视电阻率公式或时间因子约定有问题。这一步简单得近乎无聊但我在帮别人调代码时发现超过一半的“正演结果不对”都能在这一步暴露。别因为简单就跳过这是最快排除整体性约定错误的办法。4.2 高频和低频极限测试用三层模型做极限测试ρ [100, 1000, 10] Ω·mh [500, 2000] m。高频极限比如1000 Hz下第一层趋肤深度约160 m和第一层厚度500 m相比明显偏小所以视电阻率应当接近第一层的100 Ω·m偏差不会太大。低频极限比如0.0001 Hz下穿透深度远超整套层系曲线应当逼近最后一层的10 Ω·m。如果递推方向写反大概率会在这种极限测试中立刻暴露高频段表现得像最后一层低频段却像第一层。这种“反直觉”的特征一眼就能看出来比对着整条曲线找差异要快得多。4.3 典型三层模型的响应特征继续用同一个三层模型频率从1000 Hz到0.001 Hz对数间隔。曲线在三个频段有明显特征频率范围主导层视电阻率趋势相位特征高频段约1000 Hz以上第一层接近第一层约100 Ω·m接近45度中频段约1~100 Hz第二层高阻出现明显峰值低于45度的谷值低频段约0.001~0.01 Hz第三层低阻逐步下降至接近10 Ω·m界面附近出现极值后向45度恢复这里的高阻层会让视电阻率出现峰值是因为高阻层阻碍了电流向下传递电流重新分配造成地表场变化的“虚假抬升”相位在视电阻率上升段低于45度下降段高于45度这个对应关系是MT解释里最常用的定性规则之一。你的程序算出来的趋势如果和这个规则冲突说明递推公式的符号约定出了问题而不是模型本身有什么特殊。5. 调试经验与向各向异性、二维三维的扩展准备5.1 常见Bug排查顺序按我的经验排错顺序可以固定下来先跑均匀半空间测试再跑高频低频极限最后对照典型三层模型。具体排查点时优先级从高到低厚度数组长度是不是len(rho)-1漏传长度会导致曲线末端突然异常。电阻率和厚度单位是否统一km当m用结果相差几千倍。频率和周期是否混用MT曲线常按周期画但正演接口大多按频率调用时一定统一。时间因子约定是否一致结果整体相位翻转或者高频段像最后一层多半是这里出了问题。是不是错误地用了实数类型存复数Z必须是复数数组用float存会被直接截断结果出来全是乱码。5.2 从正演到反演的性能优化一维正演单次计算量极小哪怕层数只有几层反演中反复调用时主要优化点就是把频率向量化。上面代码已经是向量化实现层数循环内部处理的是整个频率数组。配套反演时我一般还会限制频率点数量在几十个以内因为MT曲线的有效信息密度本身有限点数再多对拟合结果提升也很微弱。如果做马尔可夫链蒙特卡洛这类需要几千万次采样的统计反演可以考虑用numba加jit或者C扩展。但坦白说大部分场景下性能瓶颈不在正演本身而在于反演里的模型扰动、似然计算和结果存储。过早优化正演循环意义不大。5.3 各向异性层状介质的扩展方向“各向同性”意味着每层只有一个电阻率标量方向无关。实际岩石由于层理、裂隙定向发育宏观上往往呈现各向异性水平电阻率和垂直电阻率不同。扩展成各向异性之后事情会变复杂TE模式和TM模式不再拥有同一个表面阻抗两个模式需要分开递推视电阻率随极化方向出现差异。单层内部的本征阻抗和传播常数要和各向异性电导率张量耦合递推公式的形式仍然像阻抗传输但里面每个量都要换成张量推导出来的等效值。一维各向异性正演是理解“双模式曲线不一致”现象的好工具也是向更复杂介质过渡的中间台阶。从实用角度先把各向同性一维正演吃透再扩展各向异性比直接上二维三维要稳得多。二维三维的正演代码里一维结果是用来做退化验证的基准所以这一步基础打得越扎实后面碰到的疑难问题就越少。最后再分享一条我的个人习惯。我写完这个一维正演程序以后把均匀半空间测试写成了一个最短的独立脚本每次改动代码都先跑一遍。正演程序是后面所有反演和解释工作的地基地基上任何一个小误差都会被放大成曲线形态的“灵异事件”。还有在代码注释里写清楚自己采用了哪种时间因子约定和阻抗定义这行注释看起来不起眼但三个月后回来看代码时真的会救你一命。本文还有配套的精品资源点击获取 大地电磁一维正演程序实现与调试：从递推公式到Python代码
简介根据石应骏《大地电磁测深》教材解析法编写的大地电磁各向同性均匀多层层状介质一维正演程序面向地球物理、勘探工程等专业本科生及初入MT正演领域的研究者可用于快速计算任意层数水平层状模型的理论视电阻率与阻抗相位响应辅助理解一维大地电磁正演原理与编程实现。资源压缩包共4个文件以2个MATLAB脚本.m为主另有2个.asv备份文件整体大小仅2KB代码精简、无冗余依赖适合直接阅读和二次修改核心计算与调用逻辑分离便于逐段理解递推过程。目前已有464人学习下载尤其适合希望从教材公式过渡到实际数值计算的学习者。通过该程序可掌握多层介质中电磁场递推解析法的编程思路结合简单的输入参数即可得到不同地电断面的正演结果并可通过修改层参数快速对比不同模型响应为进一步开展反演或更复杂正演研究打下基础。 搞大地电磁MT的人不管以后做正演还是反演都绕不开一维正演这个基本功。我最早写“大地电磁各向同性均匀多层层状介质一维正演程序”是因为手里的实测曲线特别乱需要在现场快速判断测点下方大致几层、每层大概什么电阻率拿几组模型响应和实测曲线叠着看。这个本来只是“临时工具”的程序后来被我反复打磨变成了做反演、设计观测方案、给新人解释MT响应都离不开的东西。这篇文章不谈教材里的公式推导怎么严谨只聊一个能直接复现的版本递推公式从哪来、代码怎么写、哪些细节会让你的结果和参考值对不上以及我踩过的几个坑。适合刚开始写MT正演程序的学生、需要自研反演内核的工程师以及任何想验证别人正演代码是否靠谱的人。1. 一维正演在MT解释流程中的位置和能力边界1.1 输入输出与核心用途这个程序解决的是这样一类问题假设地下介质由n层水平均匀层组成每层电阻率恒定且各向同性已知各层电阻率和厚度计算在地表观测到的视电阻率ρ_a和阻抗相位φ随频率或周期的变化关系。输出就是常说的“正演曲线”它是一切后续解释工作的参照系。我把它当做一个“最小可用”工具在三个场景里反复使用。第一个场景是现场快速判断实测曲线出来后先用简单的层状模型去拟合曲线拐点能迅速估计基底埋深和高低阻层的数量级比直接上二维反演高效得多。第二个场景是作为反演内核任何一维反演算法都需要成千上万次正演调用正演必须可靠而且快。第三个场景是验证别的软件二维三维程序能不能退化成层状介质解拿一维结果去对立刻见分晓。还有个很容易被忽略的作用给新人讲“为什么会出现极值”“为什么相位会低于45度”这类问题时一维正演是最直观的教具。它能让你在一个小时之内理解MT响应的基本物理图像比直接看二维剖面上的花样要清晰得多。1.2 为什么一维没有被二维三维替代可能有人会说现在二维三维电磁正演已经很成熟了一维正演还有必要自己写吗我的观点是不仅必要而且是一切的起点。二维三维反演通常需要一维模型来构建初始模型解释人员也需要一维正演来理解曲线的基本形态。三维正演的代码里最简单的测试就是“把三维网格整体退化为层状介质”而检验退化是否正确靠的还是一维正演。当然也要讲清楚边界。一维正演只适用于地下横向变化不大、测点间距远大于构造尺度的情况。当地形起伏、断层、块状异常体占主导时一维结果只能作为参考不能直接下结论。另外在层状各向同性假设下一维正演结果对TE和TM模式没有极化差异因为垂直入射平面波在两个主极化方向上对应的是同一个表面阻抗。实际资料里如果发现两种模式曲线差异明显那本身就是横向非均匀或各向异性介质存在的信号一维解释框架就解释不动了。这个能力边界在解释时一定要清楚。2. 从Maxwell方程到递归公式的核心推导2.1 模型假设与MT平面波前提一维正演建立在“地表入射波近似为垂直入射平面波”的基础上。天然电磁场源在电离层和远方的雷暴传播到地表时曲率已经很小在测点局部可以当作平面波处理。将地下按电性差异分成若干水平层每一层内部电阻率恒定、各向同性界面无限延伸这就是标题里说的“各向同性均匀多层层状介质”。在导电介质中控制方程不是普通波动方程而是扩散型方程。取垂直向下为正z轴在某个角频率ω下电场满足二阶常微分方程解的形式包含沿z方向衰减的波和反射波。这里不需要去解完整的全域解只需要关心地表处电场和磁场的比值也就是表面阻抗。2.2 均匀半空间阻抗递推的锚点先看最简单的情况整个地下是均匀半空间电阻率为ρ。此时电场和磁场通过表面阻抗Z0关联。定义Z0 E_x / H_y并取时间因子为e^{iωt}的约定本征阻抗为Z0 sqrt(i ω μ ρ)对应的视电阻率就是ρ_a |Z0|² / (ω μ) ρ均匀半空间相位为45度。这个结果简单到近乎废话但它有两个大用处一是标定代码是否正确二是作为递推的起始值。只有半空间情况下波场通解中只有一个衰减方向才能直接写出表面阻抗不需要处理反射波。2.3 递推公式与时间因子约定对第j层设电阻率ρ_j、厚度h_j最后一层为半空间不设厚度。本层的传播常数为k_j sqrt(i ω μ / ρ_j)本征阻抗为Z0j sqrt(i ω μ ρ_j)给定下层顶面阻抗Z_{j1}本层顶面阻抗由阻抗传输公式求出Z_j Z0j * (Z_{j1} Z0j * tanh(k_j h_j)) / (Z0j Z_{j1} * tanh(k_j h_j))从最底层开始向上逐层递推最终得到地表阻抗Z_1进而计算所有频率点的视电阻率和相位。这里有一个几乎所有参考资料都会出现的“坑”时间因子用e^{iωt}还是e^{-iωt}会让传播常数前面差一个i公式可能写成tanh(k_j h_j)也可能写成tanh(i k_j h_j)。两种写法在各自约定下都是对的但混着用会让高频极限完全错误。我的建议是全程序固定统一约定并把这个约定写进注释别靠脑子记。3. 程序实现中的关键细节与完整代码3.1 一份能直接跑的Python实现下面是我日常使用的精简版函数。只做计算不做绘图。import numpy as np def mt1d(freqs, rho, h): 大地电磁各向同性均匀多层层状介质一维正演 参数 ---- freqs : array_like 频率数组单位 Hz rho : array_like 各层电阻率单位 Ω·m长度 n h : array_like 各层厚度单位 m长度 n-1最后一层不设厚度 返回 ---- rho_a : ndarray 视电阻率单位 Ω·m phase : ndarray 阻抗相位单位 deg mu0 4 * np.pi * 1e-7 omega 2 * np.pi * np.asarray(freqs, dtypenp.float64) # 兼容标量频率输入 if omega.ndim 0: omega omega.reshape(1) # 半空间的地表阻抗采用 e^{iωt} 约定 Z np.sqrt(1j * omega * mu0 * rho[-1]) # 从倒数第二层向上递推 for j in range(len(rho) - 2, -1, -1): Z0j np.sqrt(1j * omega * mu0 * rho[j]) kh np.sqrt(1j * omega * mu0 / rho[j]) * h[j] t np.tanh(kh) Z Z0j * (Z Z0j * t) / (Z0j Z * t) rho_a np.abs(Z) ** 2 / (omega * mu0) phase np.rad2deg(np.angle(Z)) return rho_a, phase这段代码很简洁但包含了所有核心逻辑频率数组向量化、最底层阻抗初始化、从倒数第二层向上递推、最后返回视电阻率和相位。rho和h分别按“欧姆米”和“米”传入这是最容易出单位错误的地方。3.2 为什么从最底层向上递推初始阻抗来自半空间本质是因为半空间里没有来自下方的反射波电场的通解只有衰减波一支阻抗直接等于介质本征阻抗。从下方开始递推每向上走一层都把下一层的阻抗“折算”到当前层顶面。这个过程和传输线理论里的输入阻抗计算完全同构理解了这一点公式就不会记错方向。能不能反过来从地表往下推地表阻抗正是我们要求的未知数如果从地表往下等于手里只有一个“待求量”却要以它为起点物理上不自然。从底层向上是因为底层边界条件是已知的这是递推能成立的根本原因。3.3 高频数值稳定性、复数精度与频率范围高频段当k_j h_j的实部很大时tanh会趋近于1递推公式自然退化到Z0j意思是趋肤深度远小于第一层厚度高频只能看到第一层。数值上直接用numpy.tanh不会溢出但有些人为了省事自己写(e^{2x}-1)/(e^{2x}1)这种展开形式x一大就会溢出或得到错误结果。这个坑我见过不止一次用库函数就行。低频段k_j h_j趋于0tanh近似等于k_j h_j曲线会逐渐逼近最后一层电阻率。当模型包含高阻厚层且频率低到0.0001 Hz以下时ω μ 这个因子变得非常小视电阻率计算对复数精度比较敏感。普通双精度复数complex128在绝大多数场景下够用但如果你做极低频的深探测研究建议先拿均匀半空间模型校验过精度再信任结果。频率范围选择要和目标深度匹配。用趋肤深度公式快速估算δ ≈ 503 * sqrt(ρ / f)单位是米。假设第一层电阻率100 Ω·mf1 Hz时δ约5030 m如果想探测到20 km深频率得到0.001 Hz量级。选频段之前算一笔这种估算账能避免算完发现曲线什么都没反映的尴尬。3.4 接口设计的小习惯函数只做计算、不做绘图。返回的相位用角度制单位明确写在docstring里。绘图、写文件、反演调用都在函数外面做这样反演循环里调用几千次也不会有I/O开销。另一个习惯是把时间因子约定和阻抗定义注释写在开头三个月后回来看代码会发现这个注释非常救命。4. 验证程序的三种方法缺一不可4.1 均匀半空间测试把rho设为[100]h设为空列表或长度0数组跑一遍全频段。理论结果是视电阻率恒等于100、相位恒等于45度。如果程序返回的不是这个结果说明阻抗定义、视电阻率公式或时间因子约定有问题。这一步简单得近乎无聊但我在帮别人调代码时发现超过一半的“正演结果不对”都能在这一步暴露。别因为简单就跳过这是最快排除整体性约定错误的办法。4.2 高频和低频极限测试用三层模型做极限测试ρ [100, 1000, 10] Ω·mh [500, 2000] m。高频极限比如1000 Hz下第一层趋肤深度约160 m和第一层厚度500 m相比明显偏小所以视电阻率应当接近第一层的100 Ω·m偏差不会太大。低频极限比如0.0001 Hz下穿透深度远超整套层系曲线应当逼近最后一层的10 Ω·m。如果递推方向写反大概率会在这种极限测试中立刻暴露高频段表现得像最后一层低频段却像第一层。这种“反直觉”的特征一眼就能看出来比对着整条曲线找差异要快得多。4.3 典型三层模型的响应特征继续用同一个三层模型频率从1000 Hz到0.001 Hz对数间隔。曲线在三个频段有明显特征频率范围主导层视电阻率趋势相位特征高频段约1000 Hz以上第一层接近第一层约100 Ω·m接近45度中频段约1~100 Hz第二层高阻出现明显峰值低于45度的谷值低频段约0.001~0.01 Hz第三层低阻逐步下降至接近10 Ω·m界面附近出现极值后向45度恢复这里的高阻层会让视电阻率出现峰值是因为高阻层阻碍了电流向下传递电流重新分配造成地表场变化的“虚假抬升”相位在视电阻率上升段低于45度下降段高于45度这个对应关系是MT解释里最常用的定性规则之一。你的程序算出来的趋势如果和这个规则冲突说明递推公式的符号约定出了问题而不是模型本身有什么特殊。5. 调试经验与向各向异性、二维三维的扩展准备5.1 常见Bug排查顺序按我的经验排错顺序可以固定下来先跑均匀半空间测试再跑高频低频极限最后对照典型三层模型。具体排查点时优先级从高到低厚度数组长度是不是len(rho)-1漏传长度会导致曲线末端突然异常。电阻率和厚度单位是否统一km当m用结果相差几千倍。频率和周期是否混用MT曲线常按周期画但正演接口大多按频率调用时一定统一。时间因子约定是否一致结果整体相位翻转或者高频段像最后一层多半是这里出了问题。是不是错误地用了实数类型存复数Z必须是复数数组用float存会被直接截断结果出来全是乱码。5.2 从正演到反演的性能优化一维正演单次计算量极小哪怕层数只有几层反演中反复调用时主要优化点就是把频率向量化。上面代码已经是向量化实现层数循环内部处理的是整个频率数组。配套反演时我一般还会限制频率点数量在几十个以内因为MT曲线的有效信息密度本身有限点数再多对拟合结果提升也很微弱。如果做马尔可夫链蒙特卡洛这类需要几千万次采样的统计反演可以考虑用numba加jit或者C扩展。但坦白说大部分场景下性能瓶颈不在正演本身而在于反演里的模型扰动、似然计算和结果存储。过早优化正演循环意义不大。5.3 各向异性层状介质的扩展方向“各向同性”意味着每层只有一个电阻率标量方向无关。实际岩石由于层理、裂隙定向发育宏观上往往呈现各向异性水平电阻率和垂直电阻率不同。扩展成各向异性之后事情会变复杂TE模式和TM模式不再拥有同一个表面阻抗两个模式需要分开递推视电阻率随极化方向出现差异。单层内部的本征阻抗和传播常数要和各向异性电导率张量耦合递推公式的形式仍然像阻抗传输但里面每个量都要换成张量推导出来的等效值。一维各向异性正演是理解“双模式曲线不一致”现象的好工具也是向更复杂介质过渡的中间台阶。从实用角度先把各向同性一维正演吃透再扩展各向异性比直接上二维三维要稳得多。二维三维的正演代码里一维结果是用来做退化验证的基准所以这一步基础打得越扎实后面碰到的疑难问题就越少。最后再分享一条我的个人习惯。我写完这个一维正演程序以后把均匀半空间测试写成了一个最短的独立脚本每次改动代码都先跑一遍。正演程序是后面所有反演和解释工作的地基地基上任何一个小误差都会被放大成曲线形态的“灵异事件”。还有在代码注释里写清楚自己采用了哪种时间因子约定和阻抗定义这行注释看起来不起眼但三个月后回来看代码时真的会救你一命。本文还有配套的精品资源点击获取