企业官网建设流程全解析

简介：这是面向地球物理勘探与电磁法学习者的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_{j+1} + Z₀j·tanh(k_j·h_j)) / (Z₀j + Z_{j+1}·tanh(k_j·h_j))

其中 h_j 是这一层的厚度，k_j 是这一层的传播常数。需要注意两个边界极限：当 h_j 趋近于 0，Z_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 层的结果依赖第 j+1 层的输出，携带着明显的顺序依赖，硬用矩阵运算反而把逻辑弄脏。

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 基本稳定在 100～105 Ω·m，低频端在 0.001 Hz 处回落到 10～12 Ω·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.5～2 倍左右，需要重点刻画的分界面可以局部加密，但不要让相邻层厚度差超过 100 倍。我在构造模型时通常先写下目标深度，再反向推算每层厚度，保证每一层的趋肤深度与该层厚度之比落在合理区间。这也是为什么很多教科书推荐先用对数规则网格正演，再考虑反演细化。

5. 三步验证：半空间自检、夹层对照与频散互检

5.1 第一步：均匀半空间，把程序逼到最极端

拿到这份资源后，我第一件事不是直接跑多层模型，而是把 rho 设成一个常数，比如 100 Ω·m，h 留空。此时整个递推只剩初始阻抗，每层循环都跳过，程序退化成最简单的解析解。如果这一步视电阻率不等于 100、相位不等于 45°，后面什么都别谈。

自检通过标准可以定得严格一点：视电阻率相对偏差小于 0.1%，相位偏差小于 0.1°。均匀半空间没有层间反射，这是检验符号和常数最纯粹的环境。如果这块出错，优先查 mu0 数值和复数开方写法。

5.2 第二步：低–高–低夹层模型，看曲线形态

用一个经典三层模型验证递推的层间响应。参数选择如下表：

按程序输出后，我预期看到视电阻率从高频端 10 Ω·m 左右开始，随周期增加逐渐抬升，在中频段靠近 1000 Ω·m，随后在低频端回落到约 10 Ω·m。相位应该在电阻率上升段低于 45°，在下降段高于 45°。这个“低–高–低”结构是 MT 解释里最常见的高阻屏蔽场景，它同时考验了三个层之间的远场耦合，比两层模型更能暴露程序问题。

5.3 第三步：频散互检，不依赖外部工具的内部体检

当没有标准解析解可以对照时，可以用视电阻率与相位之间的频散关系做粗检。MT 响应满足近似的 K-K 关系，工程上简化为一条经验规则：在双对数坐标下视电阻率随周期上升越快，相位越低；下降越快，相位越高。如果某个频段视电阻率快速下降而相位也降到 45° 以下，或电阻率上升而相位升到 45° 以上，基本可以判定程序里有分支写反了或厚度逻辑错位。

从那以后，我每次拿到陌生的一维正演程序，强制自己走完这三步：半空间、三层夹层、频散互检。三步全过，才敢拿结果去做反演初始模型。这套工作流与具体代码无关，你今天下载这份 Matlab 资源后也可以直接套用。希望帮到你。

本文还有配套的精品资源，点击获取