企业官网建设流程全解析

简介：这套程序面向地球物理勘探与电磁法学习者，基于石应骏《大地电磁测深》教材中的解析法，实现各向同性均匀多层层状介质的大地电磁一维正演计算，适用于MT1D理论响应模拟、电性层参数验证及教学演示等场景。压缩包共4个文件，包含2个Matlab脚本（主程序MT1D_FWD.m与辅助函数yxcx.m）及对应的2个asv自动保存版本，整体大小仅2KB，属于轻量级可直接运行的源码工具。已有468人学习下载，社区认可度较好。通过程序可快速获得层状模型下的视电阻率与相位响应曲线，便于理解正演原理、对比不同层参数的影响，也为后续反演算法开发提供清晰的基础实现。

1. MT一维正演程序：层状介质模型是反演绕不开的起点

拿到一条实测视电阻率曲线，最快的验证方式不是直接上二维反演，而是先给一组层状模型参数，用一维正演程序算出理论响应，曲线形态贴近了再往下走。这个MATLAB项目实现的就是大地电磁各向同性均匀多层层状介质一维正演，输入各层电阻率与厚度，输出对应频点的视电阻率曲线和阻抗相位。我参与某项目时，需要用正演结果给二维反演准备初始模型，就把这个核心代码只有约二十行的小程序当作唯一的验证手段，结果比预想稳定，也让我真正理解了递推阻抗的来源。整个程序不依赖第三方工具箱，适合刚接触MT的同学入门，也适合老手把正演当黑匣子复用在处理流程里。

2. 算法基石：从Maxwell方程到阻抗递推，再到视电阻率输出

2.1 物理模型与基本假设

一维正演的物理图景并不复杂：把地下空间看作若干水平导电层，每层内部电阻率恒定、各向同性，只存在层间电性差异，不关心横向变化。场源视为从高空垂直入射的平面电磁波，在地表附近满足准静态条件，因此可以忽略位移电流。

在时谐因子取 e^{-iωt} 的前提下，从Maxwell方程组出发：

∇×E = iωμH
∇×H = σE

消去磁场后，电场水平分量满足方程 d²E_x/dz² - k²E_x = 0，其中波数 k 由角频率、磁导率和电阻率共同决定。对第 j 层，我习惯写成：

k_j = sqrt(-1i * omega * mu0 / rho_j)

这里用MATLAB表达式的形式写出，是因为程序里直接就是这么算的。复数的开方存在符号约定问题，但只要整个程序从波数、固有阻抗到递推公式统一使用同一套约定，最终视电阻率和相位数值就是稳定的。

2.2 均匀层内的解析解与固有阻抗

在一层均匀介质内部，电磁场由上行波和下行波叠加而成，这是递推公式的基础。对某一层 j，定义该层的固有阻抗：

Z0_j = i * omega * mu0 / k_j

固有阻抗代表均匀半空间里平面波垂直入射时的地表阻抗。当层是有限厚度介质时，上下行波在层顶和层底来回反射，等效阻抗不再等于固有阻抗，而是受到层厚 h_j 和传播因子 k_j 的调制。双曲正切函数 tanh(k_j * h_j) 恰好描述了这种调制关系，它同时吸收了波的衰减和相位变化。

2.3 层间递推公式：从底部半空间向上迭代

多层介质的核心递推公式如下，这也是整个正演程序的灵魂：

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

物理含义很直接：若已知下方第 j+1 层顶部的阻抗 Z_{j+1}，就能推出第 j 层顶部的阻抗 Z_j。递推从最底部开始，底部半空间没有上行波，其顶部阻抗就是固有阻抗，即 Z_N = Z0_N；然后自下而上逐层套用公式，最终得到地表阻抗。整个递推过程是一个标准的“初值加迭代”问题，完全适合用循环实现。

2.4 视电阻率与相位换算

得到地表阻抗 Z_surf 后，视电阻率和阻抗相位分别为：

rho_a = abs(Z_surf)^2 / (omega * mu0)
phi = atan2(imag(Z_surf), real(Z_surf))

视电阻率公式可以简单验证：当介质是均匀半空间时，Z_surf = Z0，代入上述公式得到 rho_a 恰好等于该层真实电阻率，这是一个非常重要的自检性质。相位则反映了电场相对磁场的滞后程度，在反演中与视电阻率联合使用，能提供额外约束。

频率与探测深度的关系也值得说明。趋肤深度 delta = sqrt(2 * rho / (omega * mu0))，高频电磁波衰减快、探测浅，低频衰减慢、探测深。所以一维正演必须遍历一整串频率，每个频点独立完成一次递推，最后拼接成完整的响应曲线。

3. 程序结构与数据流：五个文件解决“模型参数到响应曲线”

3.1 文件组织与职责划分

这个程序采用“入口脚本 + 功能函数”的分层方式，文件不多，分工明确：

把计算和绘图拆开是有意的设计。正演函数不负责画图，只返回数组，这样在批量参数扫描或反演迭代中，可以连续调用上万次而不需要创建图形窗口，性能上不会翻车。

3.2 入口脚本设计：一次完整的模型试算

% demo_run.m % 四层模型: 100Ω·m / 30Ω·m / 200Ω·m / 10Ω·m(半空间) rho = [100, 30, 200, 10]; h = [500, 200, 800]; freqs = make_frequency(320, 0.001, 8); [rho_a, pha] = mt1d_forward(rho, h, freqs); plot_response(freqs, rho_a, pha, '4-layer demo model');

代码逻辑说明：rho 有 4 个元素，h 只有 3 个，因为最后一层是半空间，没有厚度。freqs 从 320 Hz 到 0.001 Hz，每个数量级取 8 个频点，一共得到约 50 个频率值。这种“厚度比层数少一”的约定贯穿整个程序，使用时很容易记。

3.3 调用约定与数据结构

约定只有三条：电阻率单位是 Ω·m，厚度单位是 m，频率单位是 Hz；输出向量与输入频率向量等长；正演函数只做数值计算，不包含任何界面逻辑。这三条约定让函数天然适合嵌入后续的反演循环。

有一点值得提前说：mt1d_forward 内部使用循环逐频点计算，这在MATLAB中已经足够快。如果以后要处理上万频点，可以把外层循环改写成矩阵运算，把不同频率的递推同时算，但代价是代码可读性明显下降。我倾向于先保证逻辑清晰，性能优化放到确实需要时再做。

3.4 安装调试前的基准测试：半空间模型自检

修改代码后首先要跑一个基准测试：把模型设成单一半空间，例如 rho = [150]，调用 mt1d_forward，这时程序里层的递推循环不会执行，Z_surf 等于固有阻抗，输出视电阻率应该严格等于 150 Ω·m。如果不等于，说明波数、固有阻抗或换算公式之间存在符号或量纲不一致。

这个测试虽然简单，却能把最常见的公式约定问题一次性暴露出来。我从第一次写正演程序起就保留了这段基准测试脚本，后来每次改代码都会先跑一遍，再去算多层模型。

4. 核心代码精读：递推循环、相位换算与绘图出图

4.1 频率序列生成函数 make_frequency.m

function freqs = make_frequency(fmin, fmax, pts_per_decade) % 生成对数等间隔频率序列 % 输入: fmin - 最低频率(Hz), fmax - 最高频率(Hz) % pts_per_decade - 每十倍频程采样点数 % 输出: freqs - 频率向量，对数均匀分布 left = log10(fmin); right = log10(fmax); ndec = right - left; npts = round(ndec * pts_per_decade) + 1; freqs = 10 .^ linspace(left, right, npts); end

代码逻辑说明：linspace 在线性指数空间生成等间隔点，10 的幂次方转回频率值，所以频点在双对数坐标系中是均匀排列的。npts 加 1 是为了让 fmax 严格落在序列末尾，避免绘图时曲线少掉最右侧一个点。pts_per_decade 常用取值是 6 到 12，取太少曲线不够平滑，取太多则计算量增加。

4.2 核心正演函数 mt1d_forward.m

function [rho_a, pha] = mt1d_forward(rho, h, freqs) % 大地电磁各向同性均匀多层层状介质一维正演 % 输入: rho - 各层电阻率(Ω·m), 长度 N % h - 各层厚度(m), 长度 N-1, 最后一层为半空间 % freqs - 频率序列(Hz) % 输出: rho_a - 视电阻率(Ω·m) % pha - 阻抗相位(°) mu0 = 4*pi*1e-7; nLayers = numel(rho); M = numel(freqs); rho_a = zeros(1, M); pha = zeros(1, M); for iF = 1:M omega = 2 * pi * freqs(iF); % 底部半空间固有阻抗作为递推起点 kN = sqrt(-1i * omega * mu0 / rho(nLayers)); Z_next = 1i * omega * mu0 / kN; % 自下而上递推，直到地表 for iL = nLayers-1 : -1 : 1 k = sqrt(-1i * omega * mu0 / rho(iL)); Z0 = 1i * omega * mu0 / k; th = tanh(k * h(iL)); Z_next = Z0 * (Z_next + Z0 * th) / (Z0 + Z_next * th); end Z_surf = Z_next; rho_a(iF) = abs(Z_surf)^2 / (omega * mu0); pha(iF) = atan2(imag(Z_surf), real(Z_surf)) * 180 / pi; end end

参数说明：rho 和 h 的维度约定与第 3 章一致；mu0 是真空磁导率常数；omega 为角频率；Z_next 在循环中反复更新，含义是“当前已计算到的下一层顶部阻抗”。内层循环从 nLayers-1 递减到 1，保证递推方向始终从深部向地表。

逻辑说明：外层循环逐频点计算，内层循环逐层递推。k 是复波数，tanh 的输出是复数，指数因子 exp(-2kh) 在这里被自动吸收进双曲正切，不需要显式计算。rho_a 的换算公式在半空间模型下会收敛到真实电阻率，这是程序自洽性的体现。相位输出直接用 atan2，但要注意相位跳变问题，第 5 章会专门讨论。

4.3 绘图函数 plot_response.m

function plot_response(freqs, rho_a, pha, title_str) % 双面板绘图: 左图视电阻率, 右图阻抗相位 figure('Position', [200 100 960 420]); subplot(1, 2, 1); loglog(freqs, rho_a, 'o-', 'LineWidth', 1.4, 'MarkerSize', 6); xlabel('Frequency (Hz)'); ylabel('Apparent Resistivity (\Omega\cdotm)'); set(gca, 'XDir', 'reverse'); grid on; subplot(1, 2, 2); semilogx(freqs, pha, 's-', 'LineWidth', 1.4, 'MarkerSize', 6); xlabel('Frequency (Hz)'); ylabel('Impedance Phase (deg)'); set(gca, 'XDir', 'reverse'); grid on; if nargin > 3 sgtitle(title_str); end end

代码逻辑说明：XDir 设置成 reverse 是为了让高频在左、低频在右，这符合大地电磁测深曲线的行业习惯。左图使用双对数坐标，因为视电阻率跨越多个数量级；右图相位使用半对数坐标，横轴对数、纵轴线性。两组曲线共享同一频率轴，方便观察高频段和低频段的变化关系。

4.4 用示例模型跑一遍：认识曲线的形态特征

以三章中的四层模型为例，可以预期曲线大致呈现“三段式”特征。中等频段视电阻率接近 30 Ω·m，对应低阻第二层；低频段曲线最终逼近 10 Ω·m，对应半空间的基底电阻率；高频段曲线在首层电阻率 100 Ω·m 附近波动。相位曲线则在层间电性差异大的地方出现明显极值。

这种形态与真实MT测深曲线一致。正演的意义不在于算出“正确答案”，而在于帮助判断某个地电模型会产生怎样的地表响应——比如低阻层埋深增加时曲线极值向低频移动，这一规律在反演初始模型设计时非常有用。

5. 避坑指南：从相位跳变到高频数值失真，四个最容易翻车的地方

5.1 相位曲线在过零点附近跳变

现象：视电阻率曲线平滑，但相位曲线在某几个频点突然从正值跳到负值，或者在 ±90° 附近来回跳。

原因：atan2 返回的是主值，范围被截断在 -180° 到 180°，当阻抗实部或虚部过零时，角度会跨越分支线，产生看似不连续的跳变。物理上的相位曲线应当是连续的。

解决：在输出相位前使用 unwrap（注意unwrap按弧度处理）。具体做法是把 atan2 的结果求角度，再乘以 180/pi 前先unwrap，或者直接对弧度结果执行 unwrap。

pha(iF) = unwrap(angle(Z_surf)) * 180 / pi;

这个方法能消除大部分非物理跳变，保留相位曲线的整体趋势。如果仍有个别跳变点，检查该频点是否接近模型固有的相位翻转位置——那种跳变是真实的，不需要强行抹平。

5.2 高频段视电阻率锯齿状振荡

现象：频率高到几十千赫兹以上时，视电阻率曲线出现锯齿，相位曲线同时变得杂乱。

原因：高频时首层的电厚度 k*h 很大，tanh 数值趋近于 1，两个接近的复数在递推中相减相除，数值精度损失被放大。同时电磁波在首层内已经衰减殆尽，深层信息无法传递到地表，底层对响应几乎没有贡献。

解决：控制频率范围，让最小趋肤深度不小于首层厚度的一定比例，我一般按五分之一到三分之一来取。用趋肤深度公式反推频率上限：

f_max ≈ 2 * rho_1 / (pi * mu0 * h_1^2)

超过这个上限，程序算出来的数值没有物理意义，只会增加曲线毛刺。这个高频限制在真实MT资料处理中本来就存在，并不是程序的缺陷。

5.3 薄层响应被“吞掉”，曲线对某些参数不敏感

现象：模型里加了一层很薄的低阻层，正演出来曲线和没有这层时几乎一样；或者修改薄层厚度，曲线完全不动。

原因：一维MT对深度的分辨能力受趋肤深度限制，通常当地层厚度远小于趋肤深度而且电导差异不大时，该层对地表阻抗的影响低于数值分辨率。这是物理解释的非唯一性问题，不是程序算错了。

解决：正演模拟阶段把薄层合并到相邻层，使用等效电阻率参与计算，避免在一个必然不可观的参数上浪费精力。这类问题在阅读反演结果时尤其常见——反演算法反出一个薄层，不要立刻认为它真实存在，先看它在目标频段内是否真的影响了响应。

5.4 与参考程序对比时整体差了一个相位偏移

现象：视电阻率曲线和对齐，但相位曲线整体平移了约 90° 或方向相反。

原因：MT正演存在符号约定差异。不同的时谐因子选取、不同的极化方向定义、不同的阻抗正反比，都会让相位差出一个固定偏移。这不是程序bug，而是约定不一致。

解决：程序里始终保留一个半空间基准测试，用已知解析解校准相位输出的初始偏移。我在把正演结果接入其他处理系统时，会先构造一个均匀半空间模型对比两套结果的相位，统一后再接入后续流程。调试这类顺手记过一句话：相位对不齐，先查约定，不要改系数。记住这个口诀，能少走很多弯路。

6. 二次开发实战：把正演内核接入参数扫描与反演迭代

6.1 用半空间自检和两层模型做回归验证

每一次改动正演代码，我都固定用两个测试用例验证：半空间模型必须严格还原输入电阻率；两层模型的曲线特征必须符合教科书上的定性规律。这两个用例通过后，才敢把新版本的正演函数用于反演迭代。

6.2 批量参数扫描：观察某一层电阻率的影响

% 扫描中间层电阻率从10到100 Ω·m rho_scan = logspace(log10(10), log10(100), 15); figure; hold on; legend_txt = cell(1, 15); for k = 1:15 rho_now = [100, rho_scan(k), 200]; h_now = [500, 300]; [rho_a, ~] = mt1d_forward(rho_now, h_now, freqs); plot(1 ./ freqs, rho_a); legend_txt{k} = sprintf('\\rho_2 = %.0f', rho_scan(k)); end set(gca, 'XScale', 'log', 'YScale', 'log'); xlabel('Period (s)'); ylabel('Apparent Resistivity (\Omega\cdotm)'); legend(legend_txt, 'Location', 'best');

这段代码把正演函数套进一个循环，一次性生成一族曲线，直接反映中间层电阻率变化对高频段曲线的控制程度。正演函数返回数组、不绘图的好处在这里体现出来——调用 15 次也不会有图形窗口抖动的问题。

6.3 嵌入反演框架：用 fminsearch 做一个最小二乘反演

% 用合成观测数据反演中间层电阻率 rho_true = [100, 30, 200]; h_true = [500, 300]; freqs_obs = make_frequency(100, 0.01, 6); [rho_obs, ~] = mt1d_forward(rho_true, h_true, freqs_obs); % 给观测数据加2%高斯噪声 rng(1); rho_obs = rho_obs .* (1 + 0.02 * randn(size(rho_obs))); % 目标函数: 对数域均方误差 obj_fun = @(x) mean((log(rho_obs) - ... log(mt1d_forward([100, x, 200], h_true, freqs_obs))).^2); x_init = [80]; x_opt = fminsearch(obj_fun, x_init); fprintf('反演得到的第二层电阻率: %.1f Ω·m\n', x_opt);

目标函数选择对数域残差，是因为视电阻率跨越多个数量级，直接用线性残差会被高阻层主导，低阻层信息被稀释。fminsearch 是零导数优化算法，适合参数少的一维反演。如果反演参数变成五六个甚至更多，就需要换成带梯度或全局优化的方法，但核心的正演函数无需改动。

6.4 从一维到二维：正演内核的复用方式

二维反演程序里，每个测点的初始模型通常都需要一维正演来计算。把 mt1d_forward 做成独立函数后，只需要在二维反演的代价函数里循环调用它，每个测点调用一次，即可完成初始响应计算。我自己做过这类改造，最大的收益不是代码量减少，而是正演部分被单独验证过，出现问题时可以确定是二维网格部分的问题，不必怀疑正演。

从那以后，我每次改动正演内核，都会强制先跑半空间和两层模型两个回归测试，再去做参数扫描或接入反演。这个小习惯拦下了好几次因为符号约定改动引入的错误，省下的调试时间远超写测试代码的投入。希望帮到你。

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