企业官网建设流程全解析

简介：面向大地电磁测深学习与研究人群，这份Matlab程序实现了各向同性均匀多层层状介质的一维正演计算，依据石应骏《大地电磁测深》教材中的解析法编写，代码简单明了，适合地球物理专业学生与入门研究者快速理解MT一维正演原理。资源包含4个文件，其中两个.m脚本分别承担主流程与核心计算，另有两个.asv自动保存备份文件，整套压缩包仅2KB，精炼无冗余。目前已有464人学习下载，得到不少学习者的实际验证。通过学习这份代码，读者既能直接运行并得到层状模型的视电阻率响应，也能对照教材公式逐行理解解析法推导；借助清晰的层参数输入与结果输出，还可修改各层电阻率、厚度等参数快速构建自己的地电模型，用于理论曲线分析、反演初始模型设计以及教学演示，是兼顾学习价值与实用性的正演小工具。 我做了十几年地球物理数据处理，各类正演程序算是吃饭的家伙。很多人刚接触大地电磁法（MT）时会碰上一个看起来特别“教科书”的任务：写一个各向同性、均匀多层层状介质的一维正演程序。标题听着又长又绕，实际干的事不难——把地下简化成一层一层的水平介质，每层电阻率恒定、方向同性，然后模拟天然电磁场在其中的传播，最终算出地表会观测到的视电阻率和阻抗相位。这个小程序能解决什么问题呢？反演之前的初始模型试算、野外测深方案设计、快速判断某个地电断面大概长什么样，全都离不开它。这篇博文我会从原理到代码、从验证到踩坑，把程序完整拆一遍，适合刚入门的同学照着复现，也适合做实测处理的老手拿去当标定工具。

1. 一维正演到底在算什么：先把这个模型在脑子里立起来

1.1 从麦克斯韦方程组到“地表阻抗”

大地电磁法的物理基础不复杂：天然电磁波从高空照射到地表，在地下介质中感应出电磁场，而我们在地表观测的是电场和磁场的正交分量。这些分量满足麦克斯韦方程组，对于一维水平层状介质来说，波可以看成垂直入射，问题退化为沿深度方向传播的一维Helmholtz方程。

真正关键的是引入地表波阻抗这个概念，通常写成Z = E_x / H_y（对应TE模式）。为什么说它关键？因为MT方法最终关心的视电阻率和阻抗相位，都可以直接从阻抗换算出来，不需要把所有深度上的电场磁场值全算出来。一维层状介质的处理思路，就是利用每层介质中电磁波的传播特性，把层与层之间的边界条件一层层衔接起来，最终把整条剖面的响应写成层层递推的形式。

我不建议一上来就盯着公式看。先建立物理图像：电磁波从地表向下传播，在每一层界面都会发生反射和透射，地表观测到的阻抗是所有这些反射波叠加的结果。递推的过程，本质上就是把这无数次的反射“压缩”成一个等效阻抗。

1.2 为什么只有各向同性、均匀层状介质才有解析解

标题里的限定词不是随便写的。“各向同性”意味着电阻率是标量而不是张量。碰上各向异性介质，比如页岩地层常见的水平电阻率和垂直电阻率不一致的情况，一维正演的公式马上就不够用了，得引入额外的张量参数。而“均匀多层层状”则保证了每一层内的电阻率ρ和厚度h都是常数，电磁波在每一层内部可以写成简单的解析表达式，在界面处通过电场切向连续、磁场切向连续来衔接。

这套逻辑成立的根本原因，是层状介质本身允许写出闭式解。对比一下二三维正演：那需要在空间上做网格剖分，用有限差分或有限元去求解整个空间的场分布。而一维层状模型，数学上就是一组复指数函数边界条件的连接问题，只需要做阻抗递推就能得到答案，速度快得可以忽略不计。一维正演虽然简单，但它是一切更复杂正反演手段的地基。很多二三维反演程序在初始模型、约束条件、标定检验等环节都要反复调用一维正演。所以把一维正演吃透，后面路会平顺很多。

2. 把公式写成人话：阻抗递推与视电阻率计算

2.1 传播常数、固有阻抗和递推初值

先约定符号。设第j层的电阻率为ρ_j，厚度为h_j，电磁波角频率为ω，真空磁导率μ0 = 4π × 10⁻⁷ H/m。电磁波在第j层内的传播常数为：

k_j = sqrt(i ω μ0 / ρ_j)

这个k_j是复数，实部描述衰减，虚部描述相位变化。我们常说的趋肤深度δ = sqrt(2ρ / (ω μ0))，就和k_j直接相关：电阻率越高、周期越长，电磁波能穿透的深度越大，这就是为什么大地电磁法能勘探深部结构。

每一层还有一个固有阻抗，物理上可以理解为电磁波在该层介质中传播时的“特征阻抗”：

Z_0j = i ω μ0 / k_j

对于最底层，它被当成向下无限延伸的半空间。在底部没有界面反射，所以该层的上表面阻抗就直接等于固有阻抗，这成为递推的初值：

Z_N = Z_0N

这一条是整个递推的起点，也是我建议所有人在写代码时最先确认的地方。

2.2 从最深一层往上推：递推公式与响应换算

确定了底层阻抗后，逐层向上递推。第j层上表面的阻抗，由下层传递上来的阻抗Z_{j+1}和本层参数共同决定：

Z_j = Z_0j × (Z_{j+1} + Z_0j × tanh(i k_j h_j)) / (Z_0j + Z_{j+1} × tanh(i k_j h_j))

物理上可以这样理解：括号里的每一项都代表电磁波在层内往返传播后的阻抗叠加，tanh函数就是描述反射波与透射波“叠加”的数学工具。当h_j远大于趋肤深度时，tanh趋近于1，下层的影响被屏蔽掉，上表面阻抗逼近本层固有阻抗；当h_j很薄时，tanh的幅值很小，递推公式平稳地将下层阻抗传递上来。

从底层一直递推到地表，得到Z_1就是地表波阻抗。之后换算视电阻率和相位：

ρ_a = |Z_1|² / (ω μ0)

φ = arctan(Im(Z_1) / Re(Z_1))

这里的相位需要注意符号约定。MT行业习惯上按阻抗相位的辐角来显示，有的程序输出正值，表示相位超前，有的输出负值。我个人在程序里统一用arctan2取辐角，再换算成角度制，具体正负看设定的右手定则，关键是整套流程保持一致。实际资料处理中，如果相位曲线整体反号，通常不是介质模型问题，而是坐标轴方向约定错了。

还有一个换算容易踩坑：MT习惯用周期T而不是频率f来描述测深点，角频率ω = 2π / T。比如周期1秒对应角频率约6.28 rad/s，没见过直接把1/T代入公式算的——那样算出来的结果会整体偏移一个2π因子，曲线形态对但数值完全不对。

3. 写一套能直接跑的程序：核心实现与验证

3.1 模型参数与频段怎么选

一维正演的输入参数非常精简：每一层的电阻率（单位Ω·m）、每一层的厚度（单位m）、要做正演的周期序列（单位s）。在野外采集和理论模拟中，周期范围一般取0.001秒到1000秒，覆盖高频浅部和低频深部。实际设计频点时，建议采用对数等间隔采样，每个数量级取4到8个点，太密会拖慢计算，太稀则曲线形态不够光滑。

层数和厚度的选择取决于你要模拟的对象。比如想模拟一个典型的低阻夹层（H型地电断面）：表层高阻覆盖、中间低阻层、底部高阻基底，可以用三层模型：

• 第1层：电阻率100 Ω·m，厚度100 m
• 第2层：电阻率10 Ω·m，厚度200 m
• 第3层：电阻率1000 Ω·m，半无限

这类模型在实际生产中很常见，对应的是浅部风化层、中部含水层或破碎带、深部完整基底的结构。

3.2 核心代码：Python实现一维正演

直接把正演函数封装好，方便反复调用。我用numpy实现，输入输出都走数组，批量计算周期序列很顺手：

import numpy as np def mt1d(resistivities, thicknesses, periods, mu0=4 * np.pi * 1e-7): """ 大地电磁各向同性均匀多层层状介质一维正演。 参数 ---- resistivities : array_like 各层电阻率，单位 ohm.m，按地表向下的顺序排列。 thicknesses : array_like 各层厚度，单位 m，长度比电阻率少一层即可，最后底层半空间可不传。 periods : array_like 周期序列，单位 s。 返回 ---- rho_app : ndarray 视电阻率，单位 ohm.m。 phase : ndarray 阻抗相位，单位度。 """ rho = np.asarray(resistivities, dtype=float) h = np.asarray(thicknesses, dtype=float) T = np.asarray(periods, dtype=float) nlayers = len(rho) omega = 2.0 * np.pi / T rho_app = np.zeros_like(T) phase = np.zeros_like(T) for i, w in enumerate(omega): k = np.sqrt(1j * w * mu0 / rho) # 每层传播常数 z0 = 1j * w * mu0 / k # 每层固有阻抗 z_below = z0[-1] # 底层半空间阻抗 for j in range(nlayers - 2, -1, -1): t = np.tanh(1j * k[j] * h[j]) z_below = z0[j] * (z_below + z0[j] * t) / (z0[j] + z_below * t) rho_app[i] = np.abs(z_below) ** 2 / (w * mu0) phase[i] = np.degrees(np.angle(z_below)) return rho_app, phase

使用示例：

periods = np.logspace(-3, 3, 49) rho_app, phase = mt1d([100, 10, 1000], [100, 200], periods)

注意代码里没有对层数做过多限制，理论上可以传任意多层。实测下来，几十层的模型算几千个频点也就几毫秒，性能完全不用操心。循环写成频率在最外层，是因为阻抗递推本身有严格的顺序依赖，不适合直接对层方向做向量化；强行向量化反而容易把逻辑绕混。

3.3 自检清单：怎么确认程序没写错

写归写，不验证等于白写。我的习惯是每步只加一点复杂度，逐层验证：

1. 半空间模型：只有一层时，程序输出的视电阻率应该是一条水平直线，数值等于输入的电阻率，相位恒定在45度（按默认符号约定）。
2. 两层模型：高频端视电阻率趋近第一层电阻率，低频端趋近第二层电阻率，中间是平滑过渡的S形曲线。
3. 三层H型模型：高频端接近第一层，随周期增大视电阻率先下降出现极小值，再回升向基底电阻率靠近。相位在极小值低频一侧会明显抬升。

用上面那个三层模型实测，高频端1000Hz附近视电阻率约在105 Ω·m上下浮动，周期增大到几秒时视电阻率达到谷底，只有几欧姆米，随后回升，在1000秒附近逼近1000 Ω·m。相位从高频的接近45度，经过低阻层区域时抬到60到70度，再随周期增大回到接近45度。如果你把自检模型的输出和这个趋势对不上，优先怀疑递推方向或符号约定。

4. 这些年踩过的坑：常见问题与排查经验

4.1 异常结果速查表

我把实际使用中遇到过的反常情况整理成一张表，按“症状-原因-处理方式”排列，排查时对照着找比从头看公式快得多。

表格：

4.2 几个值得反复强调的细节

复数tanh的数值稳定性问题，值得单独说。当某层很厚或电阻率很低时，i k_j h_j的实部很大，tanh趋近于1；当层很薄时，tanh趋近于小量。直接使用numpy的np.tanh是复数值稳定的，但如果你为了“看清公式”，手写e指数展开的版本，就可能出现指数项溢出。我第一次写的时候吃过这个亏，高频段曲线像锯齿一样，排查半天发现是e的指数项超过了浮点上限。

另一个高频踩坑点是厚度列表的长度。底层是半空间，不需要也不可能给定厚度。代码里循环是从nlayers-2倒推，也就是第1层到倒数第2层，这个边界条件对应着厚度列表只传前nlayers-1层。如果你习惯把最后一层厚度传成0或者一个很大的数，结果会出现虚假的低频异常。我的建议是：接口定义里强制要求厚度列表比电阻率列表少一个元素，传入即报错，从源头杜绝这类问题。

相位换算还有个细节：np.angle返回的辐角范围在(-π, π]之间，如果你需要给结果做插值或绘图，碰到相位跨越±180度边界时会出现跳变。虽然一维正演本身很少碰到这种情况，但后续做反演目标函数时要小心，最好先做相位解缠或用连续化的表示方法。

最后再分享两个实用经验

根据我的实际经验，这个程序最大的价值不只是“算出一条曲线”，而是帮你建立地电模型和远区响应之间的直觉。我建议你写一个批量扫参脚本，把中间层电阻率从1扫描到10000，厚度从10扫描到1000，所有曲线叠在一张图上看。多扫几组，以后看到实测曲线时，脑子里能立刻浮现出大概的地下结构。

另外一个小技巧：把正演函数封装好之后，顺手给它加一个辅助函数，能输出各频点的趋肤深度以及等效探测深度估算值。这样设计的野外测深方案时，能快速判断当前频段到底能探测到多深的目标层，特别实用。这个一维正演程序很小，但做MT的人几乎天天都要用到它。

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