简介：压缩包提供面向三维瞬变电磁法（3D TEM）正演模拟的SemiAirMultiSourc程序集合，包含Fortran 90源码、模块文件与可执行程序，适用于地球物理勘探、矿产资源勘查及地下水探测领域的科研人员和工程师，用于构建多场源三维地质模型并计算电磁响应，为反演解释提供正演基础。包内共83个文件，以f90源程序、mod模块、txt说明/数据文件、exe可执行程序和obj/pdb调试文件为主，整体约3.74MB，结构包含主程序、电偶源/线源正演模块、滤波系数与汉克尔系数等，便于直接运行或二次开发；该程序集还提供示例数据与说明文档，可快速验证正演效果。目前已有445人学习浏览。通过这份资源，读者可掌握SemiAirMultiSourc的源码组织与多源正演计算逻辑，借助示例数据和说明文档快速上手，有助于在复杂地质构造中提高目标体定位与识别精度，为三维瞬变电磁法研究提供实用工具。
1. SemiAirMultiSourc 与三维瞬变电磁正演的痛点
SemiAirMultiSourc 这个名字拆开看就很直白：Semi-Air（半航空）、Multi-Source（多源），加上"三维正演/三维瞬变电磁"这几个限定词，指的是一套面向半航空瞬变电磁（semi-airborne TEM）勘探场景的三维正演程序。所谓半航空，是发射源在地面、接收装置在空中（无人机或飞机拖曳），它比传统地面 TEM 更适合复杂地形，比全航空 TEM 又有更强的发射功率和更大的探测深度。而多源，指的是在一次勘探中布设多个发射源，用一次三维正演把这些源的响应都算出来。
做这个方向的人通常是物探或算法工程师：手里有工区电导率模型（来自 MT 或测井约束），想快速预演半航空 TEM 在这个模型下的响应特征，评估不同源位置、不同观测高度下的分辨能力。SemiAirMultiSourc 这类工具的核心价值在于，三维正演不再只是一两个点的试算，而是把多个源、多条测线的计算统一到一个流程里，省去逐个源重复建模和网格剖分的时间。这篇文章不介绍某个特定版本的程序细节，而是顺着"半航空多源三维瞬变电磁正演"这条技术路径，把控制方程、多源叠加策略、FDTD 实现和参数设计完整过一遍。读完你会知道这类正演程序内部在算什么、网格和时间窗怎么设、多源结果怎么组织和验证。
2. 三维瞬变电磁正演的控制方程与多源叠加原理
2.1 扩散方程才是 TEM 正演的物理基础
瞬变电磁和频率域电磁最大的差别在于：它不是解波动方程，而是解扩散方程。发射源在 t=0 时刻关断后，电磁场在地下以扩散形式传播，场的衰减速度由电导率决定。控制方程从 Maxwell 方程组出发，在忽略位移电流的低频假设下得到：
∇ × E = -μ ∂H/∂t
∇ × H = σE + J_s
把两式消元，可以得到纯磁场或纯电场的扩散方程。对 FDTD 类算法来说，一般保留一阶耦合形式，直接在交错网格上做时间步进。这里的 σ 是电导率，μ 取真空磁导率即可，J_s 是发射源电流密度。
注意：半航空场景里，发射源在地面、接收点在空中，模型必须包含空气层，空气电导率通常取 1e-8 S/m 甚至更低。空气层的作用是让"关断后的感应场"能传播到空中接收点，这一点和全空间假设的地下模型完全不同，网格剖分时不能省。
三维正演的离散方式有 FDTD、有限元（矢量棱边元）、有限体积和积分方程几种。SemiAirMultiSourc 这类半航空工具常见选型是 FDTD 或有限元，原因在于：半航空模型里常有范围不小的空气层和地形起伏，FDTD 的矩形网格对地形的处理粗糙，而有限元能贴合地形；但如果工区相对平缓、关注的是层状背景上的异常体，FDTD 的效率和内存控制更有优势。具体选哪个，取决于你的网格精度要求，而不是哪个"更高级"。
2.2 多源激励的本质是线性叠加
多源（Multi-Source）在正演里有两层含义。第一层是最朴素的：多个发射源位置不同，每个源单独算一次三维正演，然后把所有源的响应放在一个文件里存起来。SemiAirMultiSourc 名字里的 MultiSourc 更多是这个意思——把多源计算作为程序的一等公民，而不是用户自己写循环反复调单源正演。
第二层是物理上的叠加原理：TEM 正演方程是线性的（σ 不随场强变化），所以多个源同时激励的响应等于各源单独响应的和。也就是说，如果源 A 和源 B 在时间上错开，你可以分别算 A、B 的响应再叠加，也可以把两个源放进一个模型同时算。后者在 FDTD 里非常划算，因为空间网格和时间步进完全不变，只是一次迭代里要同时更新多组场值。
关键技巧在这里：半航空多源正演最省时间的做法，不是逐源循环跑 N 次，而是把 N 个源的场分量合并成一个维度，在同一次时间步进里做批量更新。对 FDTD 来说，这相当于把原来的单个场数组 H[3, nx, ny, nz] 变成 H[nsrc, 3, nx, ny, nz]，内存涨 N 倍，但循环开销、网格生成、吸收边界的计算只做一次。对隐式有限元来说，多右端项（multiple RHS）几乎白送——矩阵分解只做一次，回代 N 次就行。
2.2.1 多源叠加要求源之间无耦合
叠加原理成立的前提是各源激励之间没有非线性耦合。实际操作里要区分两种情况：如果所有源采用时间同步关断（同一时刻关断、波形一致），那直接叠加即可；如果源之间有时分或频分机制，就要按各自的时窗单独截取响应，再组装。半航空 TEM 实际施工里，多源通常沿测线铺开，发射源间距几百米到几公里，接收机在空中飞过各源上方。正演时不需要把多个源同时激励，只需要把每个源的响应独立算出来，按测线和源的位置归档。
所以 SemiAirMultiSourc 内部的正演流程通常长这样：读入模型和源列表 → 对每个源生成激励 → 时间步进计算出全时间道的场值 → 在接收点位置做空间插值提取响应 → 按源-测线-时间道组织输出。这个流水线不复杂，但每一步的参数错误都会让最终响应错得莫名其妙，下面几章逐个拆开讲。
3. FDTD 实现三维瞬变电磁多源正演的最小流程
3.1 交错网格与场值更新顺序
FDTD 用 Yee 网格交错放置电场和磁场分量，每个电分量位于棱边中心，每个磁分量位于面中心。三维情况下，Hx、Hy、Hz 和 Ex、Ey、Ez 六个分量交替更新。瞬变电磁用的是扩散方程，时间步进不能直接套用波动方程的 CFL 条件，而是用无条件稳定的 Du Fort-Frankel 格式或者向后 Euler（隐式）。
不过对大多数三维 TEM 正演来说，隐式求解每步要解一个大型稀疏线性方程组，单步代价高；显式 Du Fort-Frankel 因稳定性约束放宽，步长可以取得较大，实际使用较多的是基于磁场扩散方程的显式格式。以磁场扩散方程为例，每个网格点的更新只涉及自身和相邻点的前两个时间层，纯局部操作，天然适合并行。
3.2 一个最小可运行的多源正演骨架
下面给出一段 python + numpy 的最小骨架，展示多源半航空 TEM 正演的核心循环结构。这段代码不是工业级实现，但能把"多源批量更新、时间步进、接收点插值"三个关键动作说清楚。
import numpy as np def run_semi_air_multi_source(sigma, nsrc, dt, nt, rx_positions, source_positions): """ sigma: (nx, ny, nz) 电导率模型，已包含空气层 nsrc: 发射源个数 dt: 时间步长, nt: 总步数 rx_positions: 接收点坐标列表 source_positions: 每个源的起始位置 返回: responses[nsrc, n_rx, nt]，多源多接收点的全时间道响应 """ nx, ny, nz = sigma.shape # 场值数组，最后加一维 nsrc，实现多源并行携带 Hx = np.zeros((nsrc, nx, ny, nz)) Hy = np.zeros((nsrc, nx, ny, nz)) Hz = np.zeros((nsrc, nx, ny, nz)) Hx_prev = Hx.copy(); Hy_prev = Hy.copy(); Hz_prev = Hz.copy() responses = np.zeros((nsrc, len(rx_positions), nt)) for it in range(nt): # 注入激励：第 s 个源在 source_positions[s] 处沿指定方向加电流脉冲 # 实际程序应先用静态场求初始条件，这里简化为直接赋值 # 更新磁场：Du Fort-Frankel 中心差分格式 Hx_new = update_dufortfrankel(Hx, Hx_prev, sigma, dt) Hy_new = update_dufortfrankel(Hy, Hy_prev, sigma, dt) Hz_new = update_dufortfrankel(Hz, Hz_prev, sigma, dt) # 滚动时间层 Hx_prev, Hy_prev, Hz_prev = Hx, Hy, Hz Hx, Hy, Hz = Hx_new, Hy_new, Hz_new # 接收点插值提取（只取空中接收高度上的场） responses[:, :, it] = interpolate_rx(Hx, Hy, Hz, rx_positions) return responses这段骨架有三个值得注意的设计点。第一个是场值数组多了 nsrc 维度，这是多源 FDTD 的关键——每个源独立携带一组场值，在同一个循环里推进，网格、边界、时间步长完全共享，省去重复建模。第二个是 dt 的取值，FDTD 显式格式的 dt 受电导率和网格步长约束，实际取 dt ≤ μ·σ_min·h²/4（h 为最小网格步长），这里的 σ_min 要取整个模型里非空气的最小电导率。空气层电导率太低会收紧约束，通常把空气电导率裁剪到 1e-6 S/m 以上来换取合理步长。第三个是激励注入：阶跃关断的初始场应该是发射源关断前一瞬间的稳定静态场，严格做法是先解一次直流或低频问题求初始场，再用感应场做后续步进，骨架里简化为直接赋值，正式实现时不能省。
3.3 多源结果的装配与输出
多源正演跑完以后，输出组织直接决定后续反演和解释的效率。一个实用做法是把响应存成按"源"分块的文件，块内再按"测线（接收点）"和"时间道"排列：
# 常见做法：用 npz 或 h5 组织多源响应 # run_semi_air_multisource --model model.h5 --sources sources.txt \ # --rx rx_line1.txt --nt 40 --dt 2e-6 --output resp_all.npz输出文件里建议同时写入模型网格坐标、源位置、接收点坐标和时间道向量，原因是后处理时几乎一定会做"按源提取子集""按时间道切片对比"这类操作，元信息缺失会逼着你重新跑一次正演去查坐标对应关系。时间道选择参考发射源关断后的观测窗口，半航空 TEM 通常从几十微秒观测到几十毫秒，对数等间隔取 20~40 道即可覆盖浅部到深部的响应。
4. 网格剖分、时间窗和多源布置的参数设计
4.1 网格尺寸的三个约束：最小步长、最大步长与空气层厚度
网格设计是三维瞬变电磁正演里最影响成败的一步。最小网格步长要能分辨异常体：一般取异常体最小尺寸的 1/3 到 1/2；最大网格步长往边界方向逐渐拉大，拉大倍数控制在 1.2~1.5 之间，超过 1.5 会导致网格畸变、数值色散增大。空气层的高度要覆盖接收机飞行高度以上足够空间：接收机通常在几十到几百米高度飞行，空气层顶部高度一般取模型最大探测深度的 1~2 倍，再往上用吸收边界截断。
| 参数 | 经验取值 | 说明 |
|---|---|---|
| 最小网格步长 h_min | 异常体尺寸的 1/3~1/2 | 保证异常体跨 2~3 个网格以上 |
| 网格放大倍率 | 1.2~1.5 | 超过 1.5 数值色散明显增加 |
| 空气层电导率 | 1e-6 ~ 1e-8 S/m | 过低会收紧 dt 约束，建议裁剪到 1e-6 |
| 空气层高度 | 最大探测深度的 1~2 倍 | 覆盖接收高度并留出扩散空间 |
| 吸收边界厚度 | 10~20 个网格 | 用 PML 或向上延拓近似 |
网格数量直接决定内存：单精度下每网格存 6 个场分量加电导率，约 28 字节。一个 200×200×150 的网格大约 1.68 亿格点，多源并行后乘上源数，内存是主要瓶颈。所以网格设计的第一原则是"异常区加密、背景区放纵"：把最细网格限制在发射源和异常体附近，往外按 1.3 倍等比放大，源和异常体之间用过渡网格衔接。
4.2 时间窗设计：从关断时刻到最深响应
半航空 TEM 的时间窗要覆盖感应场从浅到深的整个扩散过程。早期道（几微秒到几十微秒）反映近地表，晚期道（几毫秒到几十毫秒）反映深部。时间窗的上下界受两个因素限制：上限是接收机噪声和仪器动态范围，下限是发射源关断时间（ramp time）。如果关断时间 tc 和早期道时间可比，就必须对关断波形做去卷积处理，否则早期道响应全是关断效应。
正演里时间道建议按对数均匀分布。dt 由稳定性条件决定，输出道之间可以跨多个时间步采样。一个常见配置是：dt 取 1~5 微秒，总步数几千到几万步，输出 30~40 个对数均匀时间道。要注意晚期道的响应幅值可能衰减到早期道的 1e-6 以下，单精度浮点会损失精度，必要时晚期道改用双精度累加，或在数据处理中对响应取对数后再存储。
4.3 多源布置怎么影响正演效率
多源布置对正演效率的影响经常被忽略。源与源之间的距离如果小于扩散半径，扩散后期两个源的响应会互相重叠——这在物理上是真实的（叠加原理），但在正演组织上会导致冗余：如果源间距远大于最大扩散半径，后期响应互不影响，批量更新时每个源还是占满整个网格的内存，浪费严重。
实际工程里常见的处理是"分块并行"：把测线分成若干段，每段包含 2~4 个源，段内共享网格和空气层，段与段之间并行提交给不同计算节点。这样既保留批量更新的效率，又不至于让内存随源数线性爆炸。对于源间距特别大（比如超过最大探测深度两倍）的情况，逐源算反而更划算，把省下的内存投入到更大的网格上，效果更好。
5. 验证、提速与常见踩坑的三个技巧
5.1 用半空间解析解做第一道验证
拿到多源正演结果后，第一件事不是看异常响应，而是先做正确性验证。最常用的验证手段是均匀半空间解析解：水平层状模型的解析解（汉克尔变换或滤波系数法）有成熟实现，把三维正演模型的电导率设成均匀半空间，源和接收点按实际位置布置，把三维结果和解析解画在同一张图上。早期道应当几乎重合，晚期道的偏差通常来自边界反射和网格离散过粗，允许范围是 5% 以内。
如果晚期道偏大且随网格细化不收敛，优先怀疑吸收边界；如果早期道抖动，优先怀疑初始场注入和关断处理。这个顺序性检查能省掉大量排错时间。
5.2 多源共享一次网格，用局部加密控制规模
一个实用技巧是把"多源"和"网格分层加密"结合起来。正演模型里，源位置不同意味着最细网格区要覆盖所有源和接收测线，这会白白加密大片无异常区域。常见做法是两级网格：全局背景网格较粗（覆盖所有源和接收区域），局部加密网格只挂在异常体周围，两组网格之间用插值传递场值。半航空多源场景特别适合这种做法，因为接收点在空中飞行，不在近地表，网格约束没有频率域地面方法那么严格。
5.3 一次正演跑太久，先检查这三个参数
三维瞬变电磁正演跑得慢，大多数情况出在三个参数上：网格步长过细（内存和每步耗时同时涨）、dt 过小（总步数膨胀）、时间道输出过多（I/O 占用过高）。习惯做法是先放大一倍步长跑通流程看趋势，再把 dt 放大到稳定性边界附近，最后把输出道减半。如果三步之后还是慢，那就是模型本身太大，该考虑多源并行了。
并行时注意场值更新的空间局部性：按坐标方向分块，块边界需要一层 ghost cell 交换，每个时间步通信一次。通信量是 6 个分量 × 界面面积 × 源数，源数多时通信量线性增长，这也是多源批量更新在分布式环境下的权衡点——如果单节点内存够放全部源的场，优先用共享内存并行而不是分布式。实际观测数据进来后，还可以用多源线性叠加做快速敏感性排序，挑出对目标异常响应最大的几个源，用它们做反演的数据子集，能进一步压缩计算成本。
本文还有配套的精品资源，点击获取
