
简介压缩包提供面向三维瞬变电磁法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 正演的物理基础瞬变电磁和频率域电磁最大的差别在于它不是解波动方程而是解扩散方程。发射源在 t0 时刻关断后电磁场在地下以扩散形式传播场的衰减速度由电导率决定。控制方程从 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²/4h 为最小网格步长这里的 σ_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 个分量 × 界面面积 × 源数源数多时通信量线性增长这也是多源批量更新在分布式环境下的权衡点——如果单节点内存够放全部源的场优先用共享内存并行而不是分布式。实际观测数据进来后还可以用多源线性叠加做快速敏感性排序挑出对目标异常响应最大的几个源用它们做反演的数据子集能进一步压缩计算成本。本文还有配套的精品资源点击获取