1. 从二维到三维一个地球物理人的执念干了十几年地球物理勘探尤其是电磁法这块我最大的感受就是二维模型越来越不够用了。客户拿着复杂地形、多层矿体或者城市地下空间探测的需求过来你再用一个简单的二维剖面去解释自己心里都发虚。瞬变电磁法TEM响应快、探测深度大是找矿和工程勘查的利器但它的三维正演长期以来都是业内公认的“硬骨头”。为什么非得啃这块骨头因为现实世界是三维的。一个倾斜的板状体、一个不规则的采空区、几个相互靠近的异常体在二维近似下它们的响应会被严重扭曲甚至掩盖导致反演结果失真轻则漏掉矿体重则工程误判。所以研发一套自主可控、高效稳定的三维瞬变电磁正演程序成了我们团队必须攻下的山头。这条路我们走了好几年从最初的迷茫试错到中间的瓶颈攻坚再到最后的方案定型每一步都踩在理论和代码的刀刃上。今天我就把这几年做“正演”这部分的核心历程、技术选型的思考、以及那些教科书上不会写的“坑”和“爽点”掰开揉碎了和大家聊聊。无论你是刚入行的研究生还是正在寻找合适工具的一线工程师希望这些实打实的经验能给你一些启发。2. 技术路线抉择时域还是频域这是个问题启动三维正演研发第一个灵魂拷问就是走时域直接求解路线还是走频域转换路线这直接决定了后续整个算法框架的复杂度和计算效率的天花板。2.1 频域法的诱惑与门槛频域法的思路很直观先在频域下求解麦克斯韦方程组得到各个频率点的谐波场再通过正弦变换或余弦变换最常用的是余弦变换将频域响应转换到时域。它的巨大优势在于对于线性介质频域方程是线性的求解相对成熟。很多成熟的频域电磁法如CSAMT代码库可以借鉴开源的数值计算库如SimPEG、EM1DFM也提供了很好的频域求解器基础。我们最初也尝试过这条路线。利用有限差分或有限元在频域求解技术上确实有迹可循。但很快我们就遇到了频域法用于瞬变电磁的三个致命伤计算量指数增长瞬变电磁需要很宽的频率范围通常从零点几赫兹到几十万赫兹来捕捉早期到晚期的全过程响应。每个频率点都需要独立求解一次大型复数线性方程组。对于三维模型网格数动辄几十万甚至上百万每个频率点的求解成本都非常高。要获得一条光滑的时域衰减曲线可能需要计算上百个频率点总计算时间难以承受。变换的数值稳定性从频域到时域的数值变换是个精细活。特别是对于晚期信号对应低频成分变换核函数接近奇异对频域数据的精度要求极高。频域求解中微小的误差或截断在变换后会被急剧放大导致晚期道数据震荡甚至出现非物理的负值。为了稳定常常需要引入滤波、窗函数等技巧这又引入了人为干预和不确定性。处理关断时间效应复杂实际瞬变电磁发射电流是梯形波有关断时间。在频域处理关断效应需要将梯形波进行傅里叶分解对每个频率成分乘以相应的频谱因子这进一步增加了计算和处理的复杂度。注意虽然频域法在理论上是完备的但在追求高效、稳定、适用于大规模三维模型的正演生产中这些门槛让我们最终决定放弃。除非你的模型非常简单或者只关心早期响应否则不建议在三维瞬变电磁正演中首选纯频域路线。2.2 时域直接求解一条更艰难但更直接的路既然频域转换麻烦那能不能直接在时间维度上推进这就是时域有限差分FDTD或时域有限元FETD的思路。我们选择了基于有限差分的时域推进方案原因在于其概念清晰、易于并行且内存访问模式规则对现代计算架构友好。时域直接求解的核心是离散化麦克斯韦旋度方程在时间上一步步迭代。对于瞬变电磁控制方程通常是扩散方程形式的∇ × (∇ × e) μσ ∂e/∂t -μ ∂j_s/∂t其中e是电场σ是电导率μ是磁导率通常假定为真空磁导率j_s是源电流密度。直接显式时间推进如经典的Yee网格FDTD对电磁波传播很好但对扩散问题稳定性条件极其苛刻时间步长Δt需要小于μσ(Δx)^2的量级对于导电介质和细网格这意味着需要推进海量的时间步完全不现实。因此隐式时间积分方案成为唯一可行的选择。我们采用了后向欧拉Backward Euler方法将时间导数离散在每个时间步需要求解一个大型线性方程组[K] * {e}^{n1} {b}^n其中[K]是包含电导率和网格信息的系统矩阵{e}^{n1}是下一时间步待求的电场向量{b}^n由上一时间步的场和源项构成。这条路“更直接”是因为它一步到位得到时域解无需频域变换。但它也更“艰难”因为每个时间步都要解一个可能病态的大型稀疏线性系统对求解器的效率和稳定性要求极高。3. 核心引擎搭建从网格剖分到方程求解选定了时域隐式求解的路线接下来就是打造核心计算引擎。这就像造一辆车底盘网格、发动机矩阵组装、变速箱求解器必须协同工作。3.1 非结构化网格贴合复杂几何的必然选择对于三维复杂模型如起伏地形、任意形态地质体规则的结构化网格如长方体网格会带来两个问题一是用大量的小网格去填充空气区域浪费计算资源二是在模型边界处产生阶梯状近似影响计算精度。因此我们采用了非结构化四面体网格。我们使用了TetGen和Gmsh作为前处理工具。Gmsh强大的几何建模和尺寸场定义功能可以让我们根据模型复杂度和预期分辨率灵活控制不同区域的网格密度。例如在发射源附近、接收点位置、以及异常体边界我们会加密网格在远离源且物性均匀的区域则采用较粗的网格。一个典型的百万级单元网格生成流程如下几何模型导入/创建通过CAD文件或程序定义地形曲面和地质体。定义尺寸场根据到源的距离、到地质界面的距离等规则设定网格尺寸函数。生成网格调用Gmsh生成四面体网格并导出节点坐标和单元连接关系。质量检查与优化检查网格单元的纵横比、体积等质量指标对质量过差的单元进行优化或重新生成。实操心得网格质量直接决定求解的成败。特别要注意避免出现“银条”sliver单元即非常扁平的四面体它们会导致系统矩阵条件数恶化求解器难以收敛。在Gmsh中可以通过设置Mesh.Algorithm5(Delaunay) 或Mesh.Algorithm6(Frontal) 并配合合适的优化选项来改善网格质量。生成网格后务必用checkMesh之类的工具或自己写脚本做一遍质量筛查。3.2 有限元离散与矩阵组装精度与效率的平衡在非结构化四面体网格上我们采用矢量有限元Edge-based FEM来离散电场。与节点有限元相比矢量有限元能自然地满足电场在介质界面切向连续、法向可跃变的物理条件并且能消除伪解非常适合电磁场问题。具体来说我们在每个四面体边上定义基函数通常采用一阶Whitney边元将待求电场表示为这些基函数的线性组合。通过Galerkin加权余量法将扩散方程转化为线性系统。矩阵组装是计算密集环节。系统矩阵[K]是稀疏的其非零元素由单元刚度矩阵和与电导率相关的质量矩阵贡献。我们采用单元级组装策略遍历所有四面体单元计算每个单元对全局矩阵的贡献然后通过“散射”scatter操作累加到全局矩阵的相应位置。这个过程非常适合并行化我们使用OpenMP对单元循环进行了多线程加速。这里有一个关键细节由于采用了隐式时间积分系统矩阵[K]是随时间步变化的吗对于后向欧拉[K] [S] [M]/Δt其中[S]是刚度矩阵与电导率分布有关[M]是质量矩阵与电导率有关。如果电导率模型不随时间变化绝大多数正演场景那么[S]和[M]是常数。只有Δt可能变化。因此如果采用固定的时间步长[K]在整个时间推进过程中是不变的这带来了一个巨大的优化机会我们只需要在第一步组装并分解一次矩阵后续每一步都复用这个分解结果进行快速回代即可。3.3 线性求解器稳定计算的定海神针时域隐式求解的成败系于线性求解器一身。我们需要求解A*x b其中A就是我们的系统矩阵[K]。这个矩阵具有以下特点稀疏、对称正定SPD、但可能病态尤其是当网格质量差或电导率对比度极大时。基于上述特点我们放弃了直接求解器如LU分解因为对于百万量级未知数其内存消耗和计算时间都是灾难性的。迭代求解器是唯一选择而预处理共轭梯度法PCG因其对SPD矩阵的高效性成为我们的首选。关键在于预条件子Preconditioner。一个糟糕的预条件子会让PCG迭代数百上千次也不收敛。我们测试了多种方案预条件子类型构建成本应用成本效果评估适用场景雅可比对角极低极低效果差迭代次数多仅用于调试不用于生产SSOR低中有一定改善但对强异性问题效果有限中小规模、条件较好的问题不完全乔列斯基分解IC中中效果显著是常用选择通用性较好但需要调整填充级别代数多重网格AMG高低效果极佳迭代次数少且与问题规模弱相关大规模、强异性问题的首选我们最终选择了基于Hypre库的BoomerAMG作为预条件子。AMG通过构建一系列粗细网格在粗网格上快速消除低频误差其收敛速度几乎与网格规模无关。虽然其构建阶段setup耗时较长但正如前面提到的我们的系统矩阵[K]在整个时间推进中不变因此AMG只需要在第一步构建一次后续每个时间步的求解成本极低只需几次PCG迭代即可收敛。配置BoomerAMG参数是个经验活。我们经过大量测试找到了一组相对稳健的参数PCG求解器容差1e-8 AMG Cycle类型V-cycle 平滑器Gauss-Seidel或Chebyshev 粗网格最大层数25 强连接阈值0.25这些参数不是金科玉律需要根据具体模型调整。例如对于导电性极好的模型如良导体可能需要更强的平滑或调整阈值。4. 源与接收的细节魔鬼在这里正演的准确性不仅取决于求解器更取决于如何正确地“注入”源和“读取”接收点数据。这里面的细节差之毫厘谬以千里。4.1 发射源的处理从理想电流到有限尺寸导线很多教科书或简单代码里把发射源当作一个点电流元或电偶极子。但在实际三维正演中尤其是针对大定源回线或长导线源必须考虑源的有限尺寸和几何形态。我们的做法是将发射线圈或导线离散为一系列短线段。对于每个短线段根据其方向和电流大小计算它在所经过的网格单元中形成的源电流密度矢量j_s。在有限元框架下这个源项会出现在方程右端项{b}中。对于回线源需要确保所有线段首尾相连形成一个闭合回路。更关键的是关断时间的模拟。瞬变电磁观测的是发射电流关断后产生的二次场。我们采用梯形波模拟发射电流上升时间、稳定时间、关断时间。在时域推进中我们需要计算∂j_s/∂t。在电流上升和关断阶段这个导数不为零是激发二次场的源在稳定发射期间导数为零此时方程右端项仅由之前的场变化贡献。我们采用线性关断模型。假设关断时间为t_off在关断阶段[t0, t0t_off]内电流从I0线性下降到0。那么∂j_s/∂t就是一个负的常数。将这个常数项正确地集成到受影响的网格单元中是获得准确早期响应的关键。4.2 接收响应的计算与存储接收点通常测量的是垂直磁场分量Bz或其时间导数dBz/dt。我们的有限元解直接得到的是电场e。因此需要从电场推导出磁场。根据法拉第电磁感应定律∇ × e -∂b/∂t。在时域我们对求得的电场e进行旋度运算得到-∂b/∂t然后通过时间积分在时域推进中自然累积或解一个泊松方程来得到磁场b本身。我们选择在时域推进中同步计算磁场在每个时间步求解得到e^{n1}。在接收点所在的网格单元内利用单元形函数和e^{n1}插值得到该点的电场。对该点电场进行数值旋度运算得到(∂b/∂t)^{n1}。采用梯形积分公式更新该点的磁场b^{n1} b^n 0.5*Δt*[(∂b/∂t)^n (∂b/∂t)^{n1}]。为了节省存储我们并不保存每个时间步所有网格点的场值而是只保存接收点处的磁场时间序列。输出时根据用户需要输出b(t)或db/dt(t)。踩坑实录接收点定位。如果接收点恰好位于网格节点或边上插值很简单。但更多时候接收点位于网格单元内部。我们必须快速定位接收点位于哪个四面体内点定位问题。如果每次计算都全局搜索开销巨大。我们的优化方法是在预处理阶段利用网格的邻接关系为每个接收点建立一个“可能单元”列表在时间推进时从一个初始猜测单元开始利用重心坐标判断点是否在该单元内如果不是则根据重心坐标的符号跳转到相邻单元。这个过程通常几步内就能完成定位效率极高。5. 时间步进策略与整体流程时域推进需要选择时间步长Δt。瞬变电磁响应跨越多个数量级从微秒到秒早期变化快需要小步长晚期变化慢可以用大步长。采用固定的Δt要么早期计算冗余要么晚期精度不足。因此自适应时间步长是必须的。我们的策略基于局部截断误差估计。基本思想是先用一个步长Δt推进一步得到解u1再用两个Δt/2步推进一步得到解u2。比较u1和u2的差异可以估计出当前步长的误差。如果误差小于设定容差则接受该步并尝试增大下一步的步长如果误差过大则拒绝该步减小步长重新计算。然而对于瞬变电磁扩散问题场随时间衰减早期误差可能主要来自空间离散而非时间离散。因此我们结合了对数等间隔采样的启发式方法。预先设定一个时间道序列如从1e-7秒到1秒按对数等间隔取30个道。时间推进的目标就是“命中”这些时间道。算法会动态调整步长力求在到达每个目标时间道时恰好完成一个时间步并将该时刻的场值记录下来。这样既能保证关键时间点有输出又能让步长在早期小、晚期大。整个三维正演程序的流程可以总结如下前处理读入模型文件地形、物性分区。调用网格生成器Gmsh生成非结构化四面体网格。设置发射源参数类型、位置、波形、关断时间。设置接收点位置。设置时间道序列。初始化在网格上分配电导率σ。组装系统矩阵[K]刚度矩阵质量矩阵/初始Δt。构建AMG预条件子仅一次。初始化电场e和磁场b为零。时间步进循环While当前时间 最大时间根据当前时间和目标时间道确定尝试步长Δt_try。计算右端项{b}包含源项-μ ∂j_s/∂t。使用PCGAMG求解器求解[K] * {e}^{new} {b}。计算新时刻的∂b/∂t和b^{new}。进行误差估计。若误差可接受更新场值e e^{new},b b^{new}。如果当前时间接近某个目标时间道则记录接收点的磁场值。根据误差调整下一步的步长。时间前进。若误差不可接受减小Δt_try重新计算本步。后处理输出所有接收点在所有时间道上的磁场或磁场时间导数。可选项计算视电阻率、绘制衰减曲线图等。6. 验证、性能与那些“坑”程序写完了不代表就能用了。验证和优化是更漫长的过程。6.1 如何验证你的正演程序是对的我们采用“由简入繁交叉验证”的策略一维模型验证这是黄金标准。构建一个水平层状大地模型用我们成熟的一维正演程序解析解或半解析解计算出响应。然后用三维程序在一个很大的网格中模拟同样的层状模型比较两者的结果。在排除边界影响后两者的衰减曲线应该几乎完全重合。这是我们建立信心的第一步。与商业/开源软件对比寻找一些公开的模型算例或者用简单的三维模型如均匀半空间中的立方体与成熟的商业软件如EMIGMA、Maxwell或公认的开源代码如SimPEG的TDEM模块进行对比。注意对比时要确保所有参数网格、源、时间道设置一致。对称性检验对于具有对称性的模型如中心回线下的球体其响应也应具有对称性。计算出的多测点响应如果偏离对称性可能预示着网格、源或接收处理有bug。收敛性测试逐步加密网格观察计算结果是否趋于稳定。如果结果随网格加密剧烈变化说明离散化误差还很大网格不够密。6.2 性能优化实战当程序正确性得到保证后性能就成了关键。三维正演是计算密集型任务优化手段包括并行化我们实现了两级并行。线程级并行OpenMP用于矩阵组装、场量插值、旋度计算等循环密集操作。进程级并行MPI用于超大规模模型。采用区域分解将网格划分到不同进程每个进程负责一部分区域的矩阵组装和场计算进程间通过通信边界交换数据。这对求解器提出了更高要求我们使用了Hypre的并行AMG和PCG求解器。内存优化系统矩阵[K]采用CSR压缩稀疏行格式存储只存非零元。对于固定时间步长的情况[K]只存一份。场向量e,b是主要的存储开销。I/O优化时间推进过程中避免频繁写盘。只在到达目标时间道或每隔若干步时将接收点数据写入内存缓冲区最后循环结束时一次性写入文件。6.3 踩过的“坑”与经验边界条件设置我们采用狄利克雷边界条件边界处电场切向分量为零。这要求计算区域足够大使得边界处的场已衰减到可忽略不计。如何确定“足够大”一个经验法则是计算区域边界到源和异常体的距离应大于当前计算的最大时间所对应的扩散深度的3-5倍。扩散深度δ ≈ sqrt(2t/μσ)。你需要根据最晚时间和最低背景电导率来估算。初始场的处理时域扩散问题需要初始条件。我们通常设t0时场为零。但要注意在电流关断完成的瞬间t t_off一次场瞬间消失此时二次场从零开始衰减。我们的时间推进是从t0开始的包含了关断过程因此初始零场是合理的。晚期噪声与精度随着时间推移场值指数衰减可能接近机器精度。晚期道的计算容易出现数值噪声。除了保证求解器的高精度如PCG容差设为1e-10在计算db/dt时直接对衰减的b做数值微分会放大噪声。一个技巧是在记录b(t)的同时也记录∂b/∂t在时间步进中已算出后者通常更稳定。或者可以对b(t)的衰减曲线进行平滑后再求导。空气层处理空气的电导率为零这会导致系统矩阵[K]中对应空气单元的部分奇异因为σ0质量矩阵项为零。为了避免奇异性我们通常给空气赋予一个非常小的电导率值如1e-8 S/m。这个值小到对电磁响应的影响可忽略不计但能保证矩阵非奇异。研发三维瞬变电磁正演就像在理论和代码的迷宫中摸索前行。每一个技术选型的背后都是对计算量、精度、通用性和开发成本的一次权衡。从频域到时域的抉择从结构化网格到非结构化的跨越从直接求解器到迭代求解器加强大预条件子的演进每一步都伴随着无数次的测试、调试和优化。当程序第一次正确复现出一维层状模型的响应曲线时当第一次成功模拟出三维异常体的典型“双峰”响应时那种成就感是无可替代的。这套核心引擎也为我们后续更艰巨的“反演”研发打下了不可或缺的坚实基础。反演是正演的逆过程但一个快速、稳定、准确的正演是反演能否成功收敛的前提。关于反演的故事那又是另一段充满挑战的征程了。