Scalable parallel 3-D TEM inversion via rational approximation of the matrix exponential

技术摘要：基于矩阵指数有理逼近的可扩展并行三维瞬变电磁反演

问题陈述
瞬变电磁法（TEM）是应用地球物理中用于矿产勘探、水文地质和环境研究的关键技术。然而，在非平凡环境中反演三维电导率模型在计算上极具挑战性。传统方法依赖隐式时间步进方案（例如向后欧拉法）在时域求解麦克斯韦方程组。这些方法围绕时间推进组织计算流程，其中正演建模和反演的成本与解析瞬变响应所需的时间步数紧密耦合。由于摩尔定律的限制导致单核性能提升停滞，这些固有的串行算法在现代并行硬件上难以有效扩展。

方法论
作者提出了一种新颖的反演框架，通过利用矩阵指数作用重构正演问题，将计算成本与观测时间点的数量解耦。

1. 矩阵指数的有理逼近：
   不再采用时间步进，瞬态解 $u(t)$ 被表示为矩阵指数 $\exp(-t M^{-1}K)$ 对初始向量的作用。作者采用基于 RKFIT 算法构建的该指数的有理近优逼近（RBA）。该逼近将解表示为移位线性系统的线性组合：
   $$u(t_j) \approx 2 \text{Re} \left( \sum_{i=1}^m \alpha_i^{(j)} (K - \xi_i M)^{-1} f \right)$$
   其中 $\{ \xi_i \}$ 是所有时间点共享的极点，$\{ \alpha_i^{(j)} \}$ 是预计算的系数。关键在于，移位系统矩阵 $(K - \xi_i M)$ 与时间无关。

2. 并行正演与灵敏度计算：
   该重构使得 $m$ 个移位线性系统能够独立并行求解。时间依赖性完全包含在残差（系数 $\alpha$）中，意味着正演响应和灵敏度的计算实际上与所需的观测时间数量 $K$ 无关。

   • 正演求解器： 使用并行直接求解器（MUMPS）对移位系统进行因子分解和求解。
   • 灵敏度（雅可比矩阵）： 雅可比矩阵利用同一组移位系统隐式计算。质量矩阵的导数被应用于预计算的解向量 $g_i = (K - \xi_i M)^{-1}f$。这避免了时间步进方法中典型的串行前向/伴随递归。

3. 反演框架：

   • 优化： 反演问题使用高斯 - 牛顿框架求解。线性化最小二乘问题通过 LSQR 算法求解，该算法将雅可比矩阵作为隐式算子（无矩阵）处理。
   • 正则化： 光滑性约束使用拉维雅特 - 托马斯（RT0）单元构建，以处理非结构化网格上的分段常数模型参数，确保梯度的离散化一致性。
   • 实现： 代码使用 Julia 编写，利用 Gridap.jl 框架进行有限元组装，并采用自定义的 MPI 接口连接 MUMPS，以管理共享内存架构上的并行因子分解。

主要贡献

• 算法设计范式的转变： 本文展示了从将串行时间推进方法适配到并行硬件，转向设计能够固有地暴露并发性的数学公式（RBA）的转变。
• 时间通道的解耦： 计算成本由极点数量 $m$（通常较小，例如 16–21）决定，而不是由时间通道数量 $K$（通常为数百计）决定。这使得在不增加线性代数工作量的情况下实现高分辨率时间采样成为可能。
• 可扩展的反演工作流： 通过在数学模型层面嵌入并行性，作者实现了可扩展的反演，其中主导成本（求解移位系统）可以分布在 MPI 进程上。
• 开源发布： 该实现作为开源软件发布，以支持地球物理界的可重复研究。

结果
该方法在一个具有约 70 万自由度的合成三维模型上进行了验证，该模型在半空间中包含导电和电阻异常体。

• 重建精度： 反演成功恢复了异常体的位置和大致形状。导电块得到了良好分辨，而电阻块显示出对比度降低，这与 TEM 灵敏度和光滑性正则化的已知局限性一致。
• 性能：
  • 扩展性： 强扩展实验表明，对于因子分解和回代，在 8 个 MPI 进程下均表现出接近理想的性能。由于共享内存 NUMA（非一致性内存访问）效应和内存带宽争用，在 16 个进程时效率略有下降。
  • 内存： 峰值内存消耗约为 620 GB，用于同时存储 21 个移位系统的因子分解。
  • 比较： 与隐式欧拉时间步进法（对于相同的时间窗口需要 31 次实值因子分解）相比，RBA 方法需要 21 次复值因子分解。虽然复数运算成本更高，但 RBA 方法消除了串行时间推进的关键路径，在分布式架构上提供了降低挂钟时间的优越潜力。

意义与主张
作者声称，这项工作提供了一种可扩展的替代方案，用于标准的三维 TEM 反演，专门解决了时间递归的瓶颈。主要优势不一定在于串行算术工作的减少，而在于并行关键路径的缩短以及成本与时间通道数量的独立性。

论文承认当前的局限性，指出同时存储多个稀疏因子分解会产生巨大的内存占用，且共享内存架构面临带宽瓶颈。作者建议未来的工作应侧重于转向分布式高性能计算（HPC）集群，并可能用预条件迭代求解器替代直接求解器，以进一步提高可扩展性并减少内存使用。研究结论认为，所提出的框架成功恢复了三维结构，并为更高效的大规模电磁反演奠定了基础。