简介这是面向地球物理勘探与电磁法学习者的MATLAB正演程序用于求解大地电磁在各向同性均匀多层层状介质中的一维响应。代码严格按石应骏《大地电磁测深》教材中的解析法编写结构直接、变量清晰便于对照理论公式验证结果适合初学者快速上手也可作为课堂教学或毕业设计的参考实现。压缩包共收录4个文件包括两个.m源码文件及两个.asv自动保存副本前者是可直接运行的核心脚本后者为编辑备份整体仅2KB轻量无依赖能快速部署运行。目前已有468人学习说明该程序在相关专业学生与科研人员中具有一定参考价值。通过阅读源码可掌握层状大地模型下视电阻率、阻抗相位的正演计算流程理解各向同性分层介质中的电磁场传播规律在此基础上还可进一步扩展频带、增加层数或向各向异性、二维/三维正演方向延伸是开展大地电磁数值模拟的良好基石。1. 大地电磁一维正演一份 Matlab 程序能解决什么做大地电磁MT解释的人几乎都绕不开一维正演。二维、三维反演出结果后我习惯先拿一维正演核对深部电性变化趋势确认初始模型是否合理再决定要不要继续加网格细化。这套流程听起来基础但真自己写程序跑起来就会发现表面阻抗递推式到处都是坑从符号约定到频段设置一个细节不对整条曲线就走样。这份使用 Matlab 实现的各向同性均匀多层一维正演资源核心解决的就是“给定电阻率与厚度剖面快速得到表面视电阻率和相位响应”这件事。适合三类人做 MT 反演前想验证初始模型的解释人员、刚接触电磁法想理解多层递推逻辑的学生以及想找一个可靠参照程序来校核自己代码的开发者。2. 从半空间向上递推表面阻抗的解析式与边界条件2.1 为什么只算表面阻抗平面波分解与边界连续条件大地电磁一维正演的理论起点并不复杂。自然界中的 MT 信号可以看成近似垂直入射的平面波电场和磁场在地下介质里传播时满足亥姆霍兹方程。均匀介质中电场 E、磁场 H 沿深度方向的分布可以用传播常数 k 来描述而更重要的是多层介质中任意深度处的 E 与 H 的比值在边界上必须连续。这个“比值连续”就是整条递推链能成立的根本原因。把每一层单独拿出来看层内电磁场可以写成上行波加下行波的组合。层与层交界面上切向电场和切向磁场两边相等。于是下层界面的阻抗值就能通过一层内部的波传播关系平移到上层界面再传给再往上的层。整个过程压缩完最终只需要关注地表处的比值 Z_surface E_x / H_y。视电阻率与相位全部由这个复阻抗直接算出不需要保留电磁场在每一层的具体数值。这一点决定了程序的结构不是解一个大方程组而是用一个从下到上的标量递推。用“阻抗反射”的眼光看问题比对着波动方程硬解要直观得多。我一般会建议刚接触正演的人先把“比值连续”这四个字写在注释第一行后面所有公式如果推不下去就回到这个物理条件自查。2.2 递推式的符号约定时谐因子与传播常数一起定递推式看起来只有一行但最容易翻车的恰恰是符号。程序采用 e^{-iωt} 时谐约定这是 MT 领域教科书里最常见的写法。此时电磁波沿深度方向的传播常数写成k sqrt(-i·ω·μ₀ / ρ)电阻率 ρ、角频率 ω、真空磁导率 μ₀。这里必须用复数开方而不是实数开方因为传播常数的实部决定衰减虚部决定相位二者都由同一个表达式给出。选平方根主分支后实部大于零对应波向深部衰减这个行为正好和趋肤深度公式吻合。为了确认符号没选反先用均匀半空间模型跑一遍视电阻率应当等于输入的 ρ相位应当等于 45°这是最廉价的自检。每一层的本征阻抗写为Z₀ sqrt(i·ω·μ₀·ρ)它是“如果这一层无限厚时地表看进去的阻抗”。多层情况下从最底下的半空间开始把它的本征阻抗当作初始值然后逐层向上更新。第 j 层的表面阻抗由下式给出Z_j Z₀j · (Z_{j1} Z₀j·tanh(k_j·h_j)) / (Z₀j Z_{j1}·tanh(k_j·h_j))其中 h_j 是这一层的厚度k_j 是这一层的传播常数。需要注意两个边界极限当 h_j 趋近于 0Z_j 应该退化为下面一层的阻抗当 h_j 非常大Z_j 应该收敛到本层自己的本征阻抗。用这两条极限去检验代码如果结果不满足说明递推方向或 tanh 参数写法有误。这里强调 tanh 而不是 tan是因为 tanh 在复数自变量下不会引入周期性极点数值上安静得多。2.3 为什么必须从最下面往上算半空间没有上行反射递推方向是一个容易被忽略的“逻辑定势”。初学者容易从上往下算因为地表是我们观测的位置。可实际上地表不是天然边界它上面是空气空气与地表的阻抗差会产生反射这个界面本身就是模型的一部分。真正确定性的边界在最底部即半空间。半空间的含义是这一层无限厚不存在来自更下方的反射波。因此在半空间顶面下行波依然存在上行波被物理地排除了E/H 的比值就恰好等于本征阻抗 Z₀。如果从上往下推地表处没有确定的初值只会得到一堆未定系数。所以程序里必然是先取 Z Z₀(end)再 for 循环从倒数第二层一直算到第一层。我见过把初始值写成第一层本征阻抗的改法结果是高频渐近线不对低频看起来还凑合整条曲线在中间频段出现莫名其妙的弯曲。遇到这种症状先检查初值再检查循环顺序。3. 跑通代码参数怎么定义、频段怎么选、输出看什么3.1 模型输入与频率范围先定厚度阵列再定频点这份程序的第一步是组织输入参数。电阻率数组 rho 按从上到下的顺序排列厚度数组 h 只给到倒数第二层最后一层是半空间没有厚度。频率数组通常用对数等间隔分布因为 MT 响应在周期对数轴上变化缓慢线性的频率网格不仅浪费点还会在高频段采样过密、低频段采样不足。%% 模型与频率定义 rho [100, 1000, 10]; % 电阻率单位 Ohm·m最后一层为半空间 h [500, 2000]; % 厚度单位 m长度比 rho 少一个 fmin 0.001; % 最低频率单位 Hz fmax 1000; % 最高频率单位 Hz nf 60; % 频点个数 f logspace(log10(fmin), log10(fmax), nf); % 对数等间隔频点这里的 rho 对应三层模型表层 100 Ω·m、中层 1000 Ω·m、底部 10 Ω·m。厚度 h 第一层 500 m、第二层 2000 m。选频率范围时不要凭感觉可以用趋肤深度 δ ≈ 503·sqrt(ρ·T) 估算T 是周期。想让探测深度覆盖到几十公里低频端至少要保证衰减深度大于目标深度的两倍。从工程经验看最低频对应的趋肤深度达不到目标层底部时曲线只是缓慢地向半空间渐近线爬根本看不出电阻率真值这种情况下反演初始模型也会跟着偏。3.2 主循环实现先算本征阻抗再向上递推程序的核心循环我一般这样组织。对每个频点逐次计算传播常数、层本征阻抗然后从最底层向上做递推。之所以不整体向量化是因为第 j 层的结果依赖第 j1 层的输出携带着明显的顺序依赖硬用矩阵运算反而把逻辑弄脏。mu0 4 * pi * 1e-7; % 真空磁导率 rho_a zeros(size(f)); % 视电阻率输出 phase zeros(size(f)); % 相位输出 for n 1:nf w 2 * pi * f(n); % 角频率 k sqrt(-1i * w * mu0 ./ rho); % 每层传播常数量纲 1/m Z0 sqrt( 1i * w * mu0 .* rho); % 每层本征阻抗单位 Ohm Zsurf Z0(end); % 初始值为半空间本征阻抗 for j length(rho)-1:-1:1 t tanh(k(j) * h(j)); % 层内波传播引起的阻抗变换 Zsurf Z0(j) * (Zsurf Z0(j) * t) / (Z0(j) Zsurf * t); end rho_a(n) abs(Zsurf).^2 / (w * mu0); % 视电阻率公式 phase(n) angle(Zsurf) * 180 / pi; % 相位单位度 end重点解释两处。一是复数开方Matlab 的 sqrt 对复数输入返回的是主分支结果所以sqrt(-1i * ...)能同时给出衰减常数和相位常数。二是递推中间用了 tanh 而不是手写双曲函数展开目的就是避免大参数时指数爆炸。如果在某些频段出现Inf或NaN优先检查 k*h 是否过大以及是否误用了实数 tan。输出时视电阻率公式里的 w 是当前角频率不是频率单位问题很容易在这里造成 2π 的倍数误差。3.3 输出检查高频和低频渐近线是第一个红绿灯程序跑完不要先急着反演。先看两条渐近线最高频段的视电阻率必须接近第一层真实电阻率最低频段的视电阻率必须接近半空间电阻率。高频时趋肤深度浅信号全在第一层内部看不到下面低频时波长太长所有薄层都被“平均”掉只剩半空间。我拿这份代码跑三层模型时高频端从 100 Ω·m 出发到 1000 Hz 基本稳定在 100105 Ω·m低频端在 0.001 Hz 处回落到 1012 Ω·m这个误差水平在反演初始化里完全够用。相位输出也要看一眼。在视电阻率随周期上升的区段相位通常小于 45°在下降的区段相位大于 45°。这个规律可以当快速粗检工具。如果相位曲线走向与视电阻率斜率完全脱节多半是 Z 的复角被多值性问题污染了。下一章集中写这些坑。4. 避坑指南五个容易让结果翻脸的细节4.1 高频渐近线不水平问题出在初始阻抗和符号约定现象跑均匀半空间模型时视电阻率高频段偏离输入电阻率比如输入 100 Ω·m 输出却随着频率上升缓慢爬升到 110 Ω·m 甚至更高。原因最常见是把初始阻抗设成了第一层的本征阻抗或者把 e^{-iωt} 约定写成了 e^{iωt}导致传播常数虚部的符号反了。还有一种情况是半空间本征阻抗取模后忘记保留复角直接把 abs(Z0) 当初始值高频段必然出错。解决先把初始值改成Z0(end)全代码只用一套时谐约定。写程序前在注释里固定“e^{-iωt}”后面所有 sqrt 和 tanh 参数都按照这套约定来。我习惯在代码里顺手加一行自查注释均匀半空间时视电阻率应恒等于 rho相位恒等于 45°。出现偏差就用这个模型二分定位哪里不对改哪里。4.2 曲线中段出现毛刺或尖峰tan 与 tanh 混用的后果现象三层模型跑出来的视电阻率在某个频段突然往上跳一个尖峰然后又恢复正常相位在同一个位置反向突变。原因递推公式里把 tanh(k·h) 写成了 tan(1i·k·h) 或者直接用了 tan。复数 tan 在某些频率下会遇到周期性极点k·h 的值一但跨越极点结果就剧烈跳变。虽然数学上两者可以互相转化但数值行为完全不同。解决统一使用tanh不要做人工展开。如果因为某些教材写的是正切形式也请先确认程序里最终化简为 tanh。遇到尖峰后回头查公式来源常见符号表里阻抗递推往往给出的是双曲正切版本正切版本只是等价改写不适合直接落到代码。4.3 低频段始终探不到底频段选窄了趋肤深度不够现象视电阻率曲线在低频端还在持续下降或上升看起来没有任何趋于渐近线的意思。比如我预设基底 10 Ω·m跑到 0.001 Hz 时视电阻率还有 20 Ω·m再往下算也能继续降却发现程序里低频端已经到边界了。原因最低频率对应的趋肤深度不够深。用 δ ≈ 503·sqrt(ρ·T) 估算电阻率越低穿透越浅。想看到 50 km 深处的高导层若路径上平均电阻率只有几 Ω·m需要的周期可能长达几千秒。解决先根据目标深度反推所需最低周期再把这个周期对应的频率打进代码。一般是最高频到最低频扩展 4 到 5 个数量级比如从 1000 Hz 到 0.0001 Hz。如果资源里的默认频段范围不够直接改 fmin程序结构不用动。低频段频点也可以单独加密避免尾部曲线拖着长长的“假尾巴”影响反演约束。4.4 相位输出像锯齿angle 的包装范围与符号分支现象视电阻率曲线光滑相位却从 -170° 突然跳到 20°或者正常应该在 45° 附近输出却是 -135°。原因angle 函数返回的是 -π 到 π 之间的主值但递推过程中阻抗的复角会因为符号约定或第 2 层厚度较大而越过 ±π 包装线。对 MT 来说相位本应落在 0° 到 90° 之间出现 -135° 说明阻抗整体被乘了负数因子或者分支取向与约定不一致。解决输出前统一做一次相位修正。若相位小于 -90°先加上 180°若大于 90°则减去 180°。实际代码里可以用wrapTo360再转成想要的范围。更彻底的办法是在递推结束后检查 Zsurf 的实部是否为负如果是就整体乘 -1让视电阻率不变的同时相位回到第一象限。这个操作经常是最后一公里很多程序卡在视电阻率对、相位不对的尴尬点上。4.5 厚度阵列比例失衡深层与薄层混排引发数值放大现象在某个中频段出现不应有的小凸起例如 5 km 深度附近多出一个假的高阻尖峰而模型里根本没有这个层。原因厚度数组如果按等差排列深层厚度动辄上万米浅层又只有几米层间厚度比超过两个数量级时tanh(k·h) 在部分频点会接近 ±1递推时浮点舍入误差被放大最终形成有规律的伪响应。解决厚度阵列按近似对数排列相邻层厚度比控制在 1.52 倍左右需要重点刻画的分界面可以局部加密但不要让相邻层厚度差超过 100 倍。我在构造模型时通常先写下目标深度再反向推算每层厚度保证每一层的趋肤深度与该层厚度之比落在合理区间。这也是为什么很多教科书推荐先用对数规则网格正演再考虑反演细化。5. 三步验证半空间自检、夹层对照与频散互检5.1 第一步均匀半空间把程序逼到最极端拿到这份资源后我第一件事不是直接跑多层模型而是把 rho 设成一个常数比如 100 Ω·mh 留空。此时整个递推只剩初始阻抗每层循环都跳过程序退化成最简单的解析解。如果这一步视电阻率不等于 100、相位不等于 45°后面什么都别谈。自检通过标准可以定得严格一点视电阻率相对偏差小于 0.1%相位偏差小于 0.1°。均匀半空间没有层间反射这是检验符号和常数最纯粹的环境。如果这块出错优先查 mu0 数值和复数开方写法。5.2 第二步低–高–低夹层模型看曲线形态用一个经典三层模型验证递推的层间响应。参数选择如下表模型ρ1 (Ω·m)h1 (m)ρ2 (Ω·m)h2 (m)ρ3 (Ω·m)低阻覆盖105001000200010按程序输出后我预期看到视电阻率从高频端 10 Ω·m 左右开始随周期增加逐渐抬升在中频段靠近 1000 Ω·m随后在低频端回落到约 10 Ω·m。相位应该在电阻率上升段低于 45°在下降段高于 45°。这个“低–高–低”结构是 MT 解释里最常见的高阻屏蔽场景它同时考验了三个层之间的远场耦合比两层模型更能暴露程序问题。5.3 第三步频散互检不依赖外部工具的内部体检当没有标准解析解可以对照时可以用视电阻率与相位之间的频散关系做粗检。MT 响应满足近似的 K-K 关系工程上简化为一条经验规则在双对数坐标下视电阻率随周期上升越快相位越低下降越快相位越高。如果某个频段视电阻率快速下降而相位也降到 45° 以下或电阻率上升而相位升到 45° 以上基本可以判定程序里有分支写反了或厚度逻辑错位。从那以后我每次拿到陌生的一维正演程序强制自己走完这三步半空间、三层夹层、频散互检。三步全过才敢拿结果去做反演初始模型。这套工作流与具体代码无关你今天下载这份 Matlab 资源后也可以直接套用。希望帮到你。本文还有配套的精品资源点击获取