企业官网建设流程全解析

简介：大地电磁各向同性均匀多层层状介质一维正演程序，是一份面向地球物理专业学生与研究人员的 MATLAB 工具，用于实现大地电磁法对多层介质的一维正演计算。程序依据石应骏《大地电磁测深》教材中的解析法编写，代码结构简洁、注释清晰，适合正在学习电磁测深原理、需要验证不同层状模型响应特征的初学者与科研人员使用，也可作为教学演示或后续反演研究的基础模块。压缩包共含四个文件，以两个 m 脚本为运行主体，另有两个 asv 自动备份文件，整体体积仅两 KB，轻量便携，解压后即可在 MATLAB 中直接运行，无需额外安装依赖。目前已有四百六十八人学习下载，在同类小工具中具备一定参考热度。借助该程序，读者能够快速获得不同层状介质的视电阻率与相位响应曲线，通过修改层参数即可扩展为任意多层模型，并进一步结合实测数据进行对比分析，是入门大地电磁一维正演的实用脚本。

1. 一个绕不开的起点：均匀层状介质里的电磁扩散方程

第一次拿到这份程序，我对着输出图发了一会儿呆：三层模型视电阻率曲线中间并没有出现教科书里那种平滑的“V”字形，而是带着明显的拐点和相位跳变，我以为代码写错了。后来才意识到，这就是大地电磁各向同性均匀多层层状介质一维正演该有的样子——地表观测到的MT响应，是所有地层电性、厚度和频率之间加权叠加的结果，不是简单把各层电阻率画出来。说得直白点，这份MATLAB程序就是给给定一维地电模型计算理论视电阻率和相位曲线的工具，是反演、灵敏度分析和电性结构解释绕不开的起点。适合刚开始做电磁法正反演的学生，也适合想快速验证层状介质响应形态的工程师。

2. 项目拆解：三个MATLAB文件组成的完整正演链路

拿到这份程序的第一件事，通常是先扫一遍目录。这个版本的核心部分用三个文件就能跑通：一个主脚本demo_MT1D.m，一个正演核心函数MT1D_forward.m，还有一个用来整理绘图参数的子函数plot_MT_resp.m。很多二次开发版本还会再带一个模型生成器，但最精简的骨架就是这三个文件，正演逻辑全部收敛在核心函数里，主脚本只负责给参数、调用、画图。

2.1 文件构成与各自职责

主脚本demo_MT1D.m是入口，它做的事情很直白：定义地电模型参数、定义频率采样向量、调用正演函数、把结果画成上下两联的视电阻率和相位曲线。我在阅读这类项目时习惯先看入口脚本，因为它能最快告诉我数据的流向和最终输出长什么样。

核心函数MT1D_forward.m是整个程序的心脏，它接收电阻率、厚度和频率三个输入，返回视电阻率和相位。但这个函数内部没有任何绘图代码，这是很好的习惯——计算和展示分离，后面你想接反演或者批量扫描，可以原封不动复用这个函数。

plot_MT_resp.m是绘图辅助函数，负责把频率轴设为对数坐标，把视电阻率轴设为对数，相位用线性坐标，并自动标出不同层的大致影响频段。严格来说它不是正演程序，但它决定了你第一眼看到的数据形状是否直观。

2.2 参数组织方式：电阻率、厚度、频率怎么传

这个程序的参数组织方式非常符合一维MT问题的数据结构。模型被定义为若干水平均匀层，每一层有电阻率和厚度，最后一层是无限厚的半空间，不需要厚度。

rho是一个 N 行 1 列的列向量，单位是欧姆米，表示各层电阻率，其中第 N 层是基底半空间。h是一个 N-1 行 1 列的列向量，单位是米，表示从地表到第 N-1 层的每层厚度。这两者的索引关系需要特别留意：rho(j)对应第 j 层，h(j)对应第 j 层的厚度，但第 N 层没有h(N)，因为它无限厚。如果你不小心把h的长度写成 N，程序会在循环里索引错误，这是最常见的入门报错。

频率f是一个任意长度的行向量或列向量，单位是赫兹。MT 测深的探测深度大约和根号下电阻率除以频率成正比，所以频率跨越几个数量级是很常见的。比如从 0.001 Hz 到 100 Hz，用对数均匀方式取 41 个点，既能覆盖浅部也能看到深部响应。

2.3 主流程串讲：从参数到输出到底发生了什么

主脚本里通常会有类似下面这一段：

% 三层模型：高阻浅层 - 低阻中间层 - 高阻基底 rho = [100; 10; 100]; % 电阻率：1层 100，2层 10，3层 100 h = [1000; 1000]; % 厚度：1层 1000m，2层 1000m，3层为半空间 f = logspace(-3, 2, 41); % 频点：0.001Hz ~ 100Hz，共41个对数间隔点 [rho_a, phase] = MT1D_forward(rho, h, f); subplot(2,1,1) loglog(f, rho_a) ylabel('视电阻率 / Ω·m') grid on subplot(2,1,2) semilogx(f, phase) ylabel('相位 / deg') xlabel('频率 / Hz') grid on

这里logspace(-3, 2, 41)生成的是对数等间隔频率向量，从 0.001 Hz 到 100 Hz。为什么用对数间隔？因为电磁波的衰减随指数变化，频率轴拉长后低频段的曲线细节会被压缩，对数采样能保证每个数量级都有足够多的计算点。rho_a和phase的长度与f一致，画图时直接一一对应。

在主流程里，正演函数会对每个频率单独跑一遍从底层到表层的阻抗递推，所以频率点越多计算耗时越长。不过一维问题的计算量本身很小，41 个频点在一台普通笔记本上的耗时基本在肉眼无法察觉的范围内，这为后面做反演迭代提供了很好的性能基础。

3. 核心算法拆解：递归阻抗法在一维层状介质中的落地

正演程序的算法价值不在于复杂，而在于用稳定形式实现理论公式。大地电磁一维正演的标准做法是从最底层半空间开始，逐层向上递推表层波阻抗，最后通过阻抗与视电阻率的换算关系得到结果。我要重点拆解这个递推过程，因为在反演程序里，这一段的导数计算和稳定性处理是决定迭代能否收敛的基石。

3.1 理论公式：从波数到本征阻抗

在忽略位移电流的前提下，电磁波在电阻率为rho的均匀介质中传播，波数可写成：

k = sqrt(-1i * omega * mu0 / rho)

其中omega = 2*pi*f是角频率，mu0 = 4*pi*1e-7是真空磁导率。对大地电磁而言，电阻率一般远大于 1，所以k是一个复数，其实部和虚部的值几乎相等（在简单均匀介质中相位固定为 0° 或 -90° 方向）。复数波数的物理意义是电磁波在介质中传播时同时发生衰减和相位偏移，这是大地电磁曲线拥有频率依赖性的根本原因。

每一层的本征阻抗用zeta表示：

zeta = 1i * omega * mu0 / k

经过恒等变换，等价于sqrt(1i * omega * mu0 * rho)。这个值代表无限厚均匀介质中电场和磁场之间的比例关系。

对包含基底在内的多层模型，从最底层向上递推。基底半空间的表面阻抗就是它的本征阻抗；对上方每一层，层顶阻抗由下一层顶阻抗和本层固有参数共同决定，表达式为：

Z_j = zeta_j * (Z_{j+1} + zeta_j * tanh(k_j * h_j)) / (zeta_j + Z_{j+1} * tanh(k_j * h_j))

这里h_j是第 j 层的厚度，Z_{j+1}是第 j+1 层顶部已推得的阻抗。这个递推公式的来源是电磁场在层界面处的切向电场和切向磁场连续条件，只要层与层之间没有横向变化，这个递推就能精确成立。

3.2 代码对照：核心函数逐段拆解

下面是最精简版本的正演函数。为了让你能直接读懂，我保留了变量名与公式的对应关系，并加了注释：

function [rho_a, phase] = MT1D_forward(rho, h, f) mu0 = 4 * pi * 1e-7; % 真空磁导率，固定常量 omega = 2 * pi * f; % 角频率向量 nf = numel(omega); % 频率点个数 nlayer = numel(rho); % 地层数（含半空间基底） rho_a = zeros(nf, 1); % 预分配输出 phase = zeros(nf, 1); for ii = 1:nf w = omega(ii); % 计算每一层的波数 k 和本征阻抗 zeta k = sqrt(-1i * w * mu0 ./ rho); % 注意 ./ 是逐元素除法 zeta = 1i * w * mu0 ./ k; % 从基底开始向上递推 Z = zeta(end); % 半空间表面阻抗 for jj = nlayer - 1 : -1 : 1 kh = k(jj) * h(jj); Z_next = zeta(jj) * (Z + zeta(jj) * tanh(kh)) ... / (zeta(jj) + Z * tanh(kh)); Z = Z_next; end % 由表层阻抗换算视电阻率和相位 rho_a(ii) = abs(Z).^2 / (w * mu0); % 卡尼亚视电阻率 phase(ii) = angle(Z) * 180 / pi; % 相位角，单位度 end end

这里有三个容易被忽略的参数细节。

第一，k和zeta是向量，长度等于层数nlayer，计算时用./ rho而不是/ rho，因为rho是列向量，如果不加点是矩阵除法，绝大多数情况会报维度不匹配。

第二，递推循环jj = nlayer-1:-1:1从倒数第二层开始向上走，每循环一次就更新Z。这个顺序不能反过来，否则先用浅层阻抗去算深层就会出错。

第三，abs(Z).^2 / (w*mu0)用的完全平方形式来自卡尼亚视电阻率定义。如果有人在某些教材里看到多加了一个系数 5 或者别的因子，多半是单位不同，这里统一用国际单位就能和实测仪器读数对应。

3.3 数值稳定性：为什么用 tanh 而不是 exp

很多初学版本会把递推公式改写为含exp(-2*k*h)的形式。数学上两者完全等价，但数值表现差异巨大。当层厚很大、电阻率很大或频率偏高时，k*h的实部可能超过几十甚至上百，这时exp(-2*k*h)会迅速衰减到接近零。如果在分母中出现这种极小数，除以它就会产生无穷大或极不稳定的数值，最终导致视电阻率曲线在某个频点突然出现 NaN 或者不合理的尖峰。

tanh(kh)在自变量实部很大时趋于 1，在实部接近 0 时趋于 0，它天然把这个比值限制在有限范围内，不会产生溢出。MATLAB 内置的tanh支持复数输入，内部会处理实部和虚部的周期性问题，因此直接用这个函数是最省心的。

不过tanh也不是万能的。如果频率极低，比如 1e-5 Hz，k*h的实部会小到接近机器精度，tanh的计算结果趋近于kh本身，此时分子分母同时出现小量的减法，不特别注意也可能丢失精度。常见的做法是对每个频点单独判断abs(real(kh))的阈值，小于 1e-8 时用kh直接替代tanh，大于 20 时直接用 1。这段逻辑会增加一点代码量，但能让程序在极宽频率范围内都稳定。我在实际工程代码里通常保留这个保护分支，因为反演时频率范围往往会跨越五个数量级以上。

4. 运行与验证：在MATLAB里跑通三层介质模型

这一章解决“怎么跑”的问题。一维MT正演程序不依赖任何外部工具箱，只需要 MATLAB 基础环境。为了让你能够立即验证程序是否工作正常，我给出一个可复现的三层模型测试流程，并解释输出曲线应该长什么样。

4.1 运行前置条件：一个干净的MATLAB环境

这个程序只需要 MATLAB 本身，不需要 Signal Processing Toolbox 或者 Optimization Toolbox。我自己在 R2021a 和 R2023b 上都跑过，脚本没有调用任何高级工具箱函数，所以理论上从 R2016b 之后的主流版本都能运行。

为了减少干扰，建议把三个文件放在同一个文件夹，并把当前工作目录切换到这个文件夹。在命令行窗口输入以下命令确认路径：

addpath(pwd)

如果你的文件名没有放在当前目录下，MATLAB 会找不到函数，但addpath(pwd)能确保当前目录被加入搜索路径，这是一个很实用的习惯。另外，脚本里用到subplot、loglog、semilogx，这些是 MATLAB 基础绘图命令，无需额外安装。

4.2 三层模型测试：输出曲线与物理解读

使用前面主脚本给出的三层模型：上半空间为 100 欧姆米的高阻层，中间为 10 欧姆米的低阻层，底层为 100 欧姆米的半空间。两层厚度均为 1 千米。这个模型相当于一个典型的“低阻异常夹层”场景。

运行后，视电阻率曲线在高频段会接近 100 欧姆米，因为高频趋肤效应使得电磁能量只集中在第一层；随着频率降低，电磁波穿透到低阻层，视电阻率下降，曲线出现凹陷；频率继续降到 0.01 Hz 以下时，电磁波穿透到基底高阻层，曲线又回升到 100 欧姆米附近。这个先降后升的形状就是这份程序的核心可复现输出之一。

相位曲线同样有意义。高频段相位接近 0 度，进入低阻层后相位逐渐向 -45 度甚至更低偏移，进入基底后回升。相位曲线的跳变位置和视电阻率曲线的拐点并不完全一致，两者一起使用才能更准确约束地层结构。只画视电阻率不看相位，等于浪费了一半信息。

4.3 批量计算：让频率扫描和结果导出自带

单次运行 41 个频点没问题，但在实际对比实验里我们常常需要计算多组模型。我习惯把模型循环包在外部，先定义一组模型矩阵，再逐列调用正演函数。下面是一个简单的批量示例：

model_set = [100 1000 100; % 模型1：三层 50 100 500]; % 模型2：另一种三层 h_set = [1000 1000]; f = logspace(-3, 2, 41); colors = {'r', 'b'}; figure for mi = 1:size(model_set,1) rho = model_set(mi,:)'; [rho_a, phase] = MT1D_forward(rho, h_set, f); loglog(f, rho_a, colors{mi}, 'LineWidth', 1.5); hold on end legend('模型1','模型2')

这段代码里的model_set是一个 2 行 3 列的矩阵，每一行是一组完整电阻率参数。循环里通过转置操作把行向量变成列向量再传给正演函数，避免维度不匹配。如果你需要把结果保存成文本文件，可以用writematrix在循环结束后一次性写出：

writematrix([f(:), rho_a(:), phase(:)], 'MT_result.txt')

这样生成的三列文件可以直接导入其他绘图软件或作为反演程序的初始数据。批量计算时要注意，MT1D_forward内部针对每个频率点都做了一次层循环，如果频率点很多且模型也很多，总耗时线性增长，但一维模型通常都在毫秒量级，除非你把频率点取到上千个，否则无需担心性能。

5. 避坑经验：一维MT正演最容易翻车的四个瞬间

任何正演程序在第一次投入使用时都会遇到几个典型问题。我踩过这些坑，所以把它们整理成可直接对照的排查清单。

5.1 高频段视电阻率出现锯齿

现象：在 1 Hz 到 100 Hz 的高频段，视电阻率曲线不是光滑下降，而是出现小幅度抖动或锯齿，偶尔还伴随负值。

原因：频率过高时，k*h的实部很大，tanh(kh)的结果会以指数衰减的速度逼近 1。此时，如果两个相邻频率点计算出的阻抗实部和虚部在浮点精度上差异极小，而递推公式中的分子分母又是两个接近相等的复数做减法，会触发灾难性抵消，导致高频端数值噪声被放大。

解决：最直接的办法是让频率采样点不要过密，保持每十倍频程约 10 到 20 个点即可。同时，对tanh的处理添加阈值判断：当real(kh) > 15时，直接令tanh(kh) = 1，避免不必要的浮点运算。在我自己的程序中加入这个判断后，高频锯齿彻底消失。

5.2 低频段出现NaN或Inf

现象：频率低于 1e-4 Hz 时，某些频点返回的rho_a是NaN或者Inf，把整条曲线打断。

原因：低频对应感知深度很大，但此时k*h的实部趋近于零，tanh(kh)约等于kh。如果计算过程中出现分母非常接近零的情况，浮点数就会溢出。另一种常见原因是初始化Z = zeta(end)时，如果基底电阻率极大或极小，会让zeta(end)的初值本身就导致后续分母异常。

解决：在递推前增加模型参数的合法性检查，确保所有电阻率大于零且厚度大于零。对超低频，手动设定一个截止频率，低于该频率的数据点直接跳过计算并赋值为半空间本征阻抗。更好的做法是在递推公式中改用coth形式，让它在高对比度层位时更稳定，但这通常会牺牲一部分计算速度，我这里仅作提示。

5.3 相位超出约定范围

现象：计算出的相位值落在-180°到180°之间，但实测资料或商用软件展示的相位通常在-90°到0°之间，导致曲线垂直翻转或者出现跳变。

原因：MATLAB 的angle函数返回复数阻抗的辐角主值，范围是(-π, π]。如果高频段阻抗实部为正、虚部为负，相位在第四象限，显示为负值；但如果数值噪声让实部变成负数，相位会跳到第三象限，显示为负 135 度或正 135 度，造成视觉上的大跳变。

解决：在程序里统一相位定义。我通常将相位按如下方式规范化：

phase = phase - 360 * floor((phase + 90) / 360);

这样可以把相位限制在-90°到270°之间，再根据电磁场方向约定取负值。具体实现取决于你的数据标准，但关键是必须在输入输出接口处写清楚相位定义，否则后续反演时很容易张冠李戴。

5.4 函数名与工具箱冲突

现象：正演函数在自己电脑上运行正常，换一台机器后报错 “Undefined function 'tanh'” 或者提示找不到MT1D_forward。

原因：很多用户的 MATLAB 安装路径中包含了其他同名函数，尤其是某些第三方工具箱里自带一个MT1D_forward.m，导致 MATLAB 路径搜索优先找到其他版本。还有人在函数文件里写了同名变量rho，与脚本变量产生冲突，造成随机性错误。

解决：在脚本开头使用which MT1D_forward查看当前被调用的文件路径，确认指向的是你的源码目录。如果路径不对，调整搜索路径优先级。另一个稳妥做法是把核心函数文件重命名为带前缀的名称，比如mt1d_forward_rp，降低命名冲突概率。我从那以后所有正演函数都会在文件头部加一段注释，写明版本日期和功能摘要，再用which验证，再开始跑数据，这个习惯帮我避开了很多莫名其妙的路径污染问题。

6. 二次开发：从正演程序到灵敏度矩阵的改造技巧

当你能稳定跑通正演，下一步自然是把它接进反演或者做参数影响分析。一维正演函数最大的优点是计算量小，适合做有限差分扰动求灵敏度矩阵。我提供两个具体改造方向。

第一个改造是把原来面向单模型的函数升级成面向模型向量的调用方式。反演迭代里，模型参数通常被压缩成一个长度为2*N-1的向量m = [log(rho1); ...; log(rhoN); log(h1); ...; log(h_{N-1})]。取对数是因为电阻率和厚度都跨越多个数量级，用对数参数能提高反演稳定性。在正演函数里把rho = exp(m(1:N)); h = exp(m(N+1:end));还原后调用原函数即可复用所有逻辑。

第二个改造是加入解析灵敏度计算。对阻抗递推公式中的Z_j关于rho_j和h_j求偏导数，可以用链式法则逐层回代。写代码时我会构造一个复变量递推，同时算出阻抗及其导数。由于推导过程长，这里给一个更通用的有限差分替代方案：

% 直接利用已有的正演函数计算当前模型的响应 [rho_a0, phase0] = MT1D_forward(rho, h, f); % 扰动第2层电阻率 rho2 = rho; rho2(2) = rho2(2) * 1.05; [rho_a1, ~] = MT1D_forward(rho2, h, f); % 灵敏度即一阶差分 d_rho_a_d_rho2 = (rho_a1 - rho_a0) / (rho2(2) - rho(2));

扰动比例取 5% 是经验值，太小会导致差分为零噪声主导，太大会让非线性误差增大。在批量实现时，用parfor把每个模型扰动计算分配到不同工作进程，能够显著压缩时间。我自己的做法是先跑通单点差分，再扩展到全模型扰动矩阵，最后用半空间解析解对比验证一列灵敏度值是否与理论吻合。

说到底，这个程序的终点不是画出一张漂亮的曲线，而是成为你理解观测数据与地下介质的桥梁。从那以后，我每拿到一组新模型参数，都会先跑一遍均匀半空间模型确认阻抗初值选取无误，再对目标模型做一次快速正演，确认曲线形态符合物理直觉，然后才放心进入反演流程。这个检查比任何公式推导都更早帮我发现参数输入错误，希望帮到你。

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