本文旨在介绍一篇发表于1996年1-2月期《Radio Science》期刊（第31卷第1期，第1-23页）上的研究论文。该论文题为“Three-dimensional wideband electromagnetic modeling on massively parallel computers”，作者为David L. Alumbaugh, Gregory A. Newman, Lydie Prevost, 以及 John N. Shadid，均来自美国桑迪亚国家实验室（Sandia National Laboratories, Albuquerque, New Mexico）。

一、 研究背景与目标

本研究属于地球物理电磁学（Geophysical Electromagnetics）与计算电磁学领域。其核心目标是开发一种能够在超大规模并行计算机上运行的三维（3-D）宽频带电磁（EM）建模方法。在20世纪90年代中期，随着计算机性能的提升，三维电磁建模已取得显著进展，但对于更复杂、更真实的模型，传统串行计算机在内存和计算速度上存在局限。积分方程法（Integral Equation, IE）在处理复杂大模型时计算量巨大，而各种近似方法则在频率范围和电导率对比度方面存在限制。此外，当时的方法尚未能有效应用于雷达频率范围（>10 MHz），在该频段波传播效应变得重要。

因此，本研究旨在提出一种基于微分方程的解决方案，能够精确模拟从电磁扩散主导（<100 kHz）到波传播主导（>10 MHz）的宽频带（100 Hz 至 100 MHz）电磁响应。为了实现对大尺度、高分辨率现实模型的模拟，研究团队将算法修改为可在超大规模并行（Massively Parallel, MP）计算机架构上运行，以期突破计算瓶颈。

二、 研究方法与工作流程

本研究提出的方法本质上是一种频域有限差分（Finite Difference, FD）方案。其详细工作流程可概括为以下几个关键步骤：

1. 控制方程与公式化： 研究从描述散射电场的矢量亥姆霍兹方程（Vector Helmholtz Equation）的修正形式出发。该修正形式通过引入复坐标拉伸（complex grid stretching）变量（e_j, h_j, j=x,y,z），为应用完美匹配层（Perfectly Matched Layer, PML）吸收边界条件（Absorbing Boundary Conditions, ABCs）奠定了基础。PML边界条件对于模拟高频（>10 MHz）波传播至关重要，它能有效吸收外向波，避免网格边界反射造成的误差。方程中还包含了可变电导率（σ）、介电常数（ε）和磁导率（μ）的参数，并采用了散射场（scattered field）公式，将总场分解为背景场（primary field）和散射场，这有助于在源点附近获得更平滑的源项行为。

2. 数值离散化： 采用Yee氏交错网格（staggered grid）对计算区域（包括地下和空气）进行离散化。电场分量定义在网格棱边中心，磁场分量定义在网格面中心。这种网格结构能自然满足麦克斯韦方程的微分形式。通过有限差分近似，将连续的修正亥姆霍兹方程离散化为一个大型、稀疏、复对称的线性方程组：A f = s。其中，A 是包含空间导数离散化和模型电磁参数的复对称矩阵，f 是待求的散射电场向量，s 是由背景场和异常体参数构成的等效源向量。

在离散化过程中，关键步骤包括对电导率和磁导率进行平均化处理。导纳率（σ + iωε）在网格棱边中点处计算，采用基于安培定律的加权平均法（如图2所示）。磁导率则在网格面中心处计算，采用几何平均法以确保磁感应强度法向分量的连续性。

3. 线性系统求解： 离散化后得到的线性方程组采用预处理的Krylov子空间迭代法求解。研究对比了双共轭梯度法（Biconjugate Gradient, BiCG）和拟最小残差法（Quasi-Minimum Residual, QMR）两种迭代求解器。由于QMR方法在稳定性方面表现更优（如图3所示，BiCG在迭代后期出现残差振荡，而QMR收敛更稳定），研究最终选择QMR作为主要求解器。 为了加速收敛，采用了雅可比缩放（Jacobi scaling）预处理器。该方法通过对角矩阵对原方程组进行预处理，降低了矩阵的条件数，且只需在求解器调用前应用一次，计算效率高。

4. 并行化实现： 为了处理大规模模型（高达数千万个未知数），研究将串行代码移植到MIMD（多指令多数据）型超大规模并行计算机（如1840处理器的Intel Paragon）上。并行化策略包括： * 数据分解： 将整个三维计算模型网格在物理上分解到多个处理器上，每个处理器负责一个子区域的计算，实现负载均衡。 * 输入/输出管理： 将输入数据分为全局数据（如源、接收器位置、频率等）和局部数据（各网格单元的电磁参数）。全局数据由“领导”处理器读取后广播，局部数据则从并行磁盘系统同时读入各处理器。 * 通信模式： 迭代求解过程中的矩阵-向量乘法和向量点积运算需要处理器间通信。根据有限差分模板（图1），每个处理器需要与相邻处理器交换其“子区域”边界上的未知数信息。研究定义了清晰的处理器通信模板（图4），包括通过立方体“面”的通信（交换2-3个未知数）和通过特定“边角”的通信（交换1个未知数）。 * 代码可移植性： 使用消息传递接口（Message-Passing Interface, MPI）库进行进程间通信，确保了代码在不同并行平台上的可移植性。

5. 验证与演示： 为了验证算法的正确性和通用性，研究模拟了三种不同的地球物理模型，并将结果与已知的解析解或数值解进行对比： * 航空电磁模拟： 模拟了位于空气中、覆盖在具有磁性的均匀半空间上方的垂直磁偶极子（VMD）和水平磁偶极子（HMD）响应（频率0.9 kHz和56 kHz）。有限差分结果与一维解析解（Lee的代码）高度吻合，证明了算法能正确处理空气-大地界面以及跨越不连续面的电磁场突变。 * 井间电磁模拟： 模拟了位于低阻层中或上方的垂直磁偶极子（VMD）和垂直电偶极子（VED）源，在另一口井中接收电磁场（频率0.1 kHz和10 kHz）。有限差分结果与三维积分方程解（Newman等人的代码）进行了对比。结果显示，对于源位于异常区域之外的情况，两者吻合极好。当源位于异常区域（低阻层）内部时，VMD结果依然良好，但VED结果对背景电导率的选择更敏感，需谨慎处理。研究还指出，积分方程解的精度受其自身网格离散化程度影响。 * 高频模拟与吸收边界验证： 模拟了一个多层地质模型在高达28.5 MHz频率下的响应，用于验证PML吸收边界条件的必要性。结果表明，在高于10 MHz的频率下，不使用吸收边界或仅使用实坐标拉伸（real grid stretching）会导致严重的边界反射误差。而采用复坐标拉伸（complex grid stretching）的PML边界条件，即使在网格尺寸较小的情况下，也能获得与一维解析解几乎完全一致的结果（图9d, e）。研究还发现，采用复坐标拉伸时，Krylov求解器的收敛速度对频率的依赖性较小，系统谱特性更稳定（图10）。

三、 主要研究结果

1. 算法有效性： 研究成功开发并验证了一种基于修正矢量亥姆霍兹方程和交错网格有限差分法的三维宽频带电磁建模算法。该算法能够精确模拟从100 Hz到30 MHz频率范围内，不同源类型（电/磁偶极子）和极化方式下的电磁响应，并能正确处理复杂的电磁参数变化（电导率、介电常数、磁导率）以及空气-大地界面。
2. PML边界条件的成功应用： 研究首次将基于复坐标拉伸的PML吸收边界条件应用于地球物理电磁建模的频域有限差分法中。结果表明，PML在高于10 MHz的频段至关重要，它能有效吸收外向波，允许使用更小的计算网格，从而显著减少计算量。研究还探讨了实部和虚部拉伸参数（a, b）的选择经验：低于100 kHz主要使用实拉伸（b=0），高于10 MHz需使用复拉伸（b≠0，a=0），中间频段可结合使用。
3. 并行计算性能： 算法在Intel Paragon并行计算机上成功实现。性能测试表明（图11），当每个处理器分配约10,000至24,000个未知数时，并行效率最高。使用1820个处理器时，理论上可求解的最大模型规模为280×260×200个网格单元（双精度存储下为210×195×150），对应数千万个未知数，最高浮点运算速率可达14.7 GFLOPS。相比于高端工作站（如IBM RS6000-590），计算速度提升了两个数量级。
4. 求解器与预处理器选择： 对比研究表明，QMR迭代求解器比BiCG更稳定，尽管每次迭代计算量稍大，但总体收敛时间更短。雅可比缩放预处理器简单有效，能显著加速收敛。
5. 源位置与背景场选择： 研究发现，当发射源位于电磁参数异常区域内部或附近时，散射场公式中的等效源项可能变化剧烈，导致数值困难。对于VED源位于异常区内的井间模型，需要仔细选择背景电导率以获得准确结果。

四、 研究结论与意义

本研究提出并验证了一种适用于超大规模并行计算机的三维宽频带电磁正演建模算法。该算法结合了修正的矢量亥姆霍兹方程、交错网格有限差分、PML吸收边界条件、QMR迭代求解器以及高效的并行化策略，能够精确模拟从低频扩散到高频传播的宽频带电磁响应。

其科学价值在于：1）将计算地球物理电磁建模的频率范围成功扩展至雷达频段；2）证明了PML吸收边界条件在频域地球物理电磁计算中的有效性和必要性；3）展示了超大规模并行计算技术在地球物理大规模数值模拟中的巨大潜力，使得模拟更复杂、更真实的地球模型成为可能。

其应用价值直接体现在多个领域：1）辅助地球物理仪器设计（如文中提到的VETEM系统）；2）模拟航空电磁勘探；3）井间电磁监测（如提高石油采收率过程监测）；4）为后续的三维电磁反演研究提供了高效、精确的正演引擎。

五、 研究亮点

1. 宽频带能力： 首次实现了在同一算法框架下，对从低频（100 Hz，扩散主导）到高频（>10 MHz，传播主导）的宽频带三维电磁响应的精确模拟。
2. 高效吸收边界： 成功引入并验证了基于复坐标拉伸的PML吸收边界条件，解决了高频模拟中的边界反射难题，并允许缩小计算网格。
3. 大规模并行实现： 系统性地描述了将三维有限差分电磁建模算法移植到MIMD超大规模并行计算机的完整流程，包括数据分解、通信模板设计和MPI实现，为后续研究提供了宝贵范例。
4. 算法稳定性优化： 通过对比实验，确定了QMR迭代求解器结合雅可比缩放预处理器是该类问题的高效稳定组合。
5. 综合验证： 通过多个具有代表性的地球物理模型（航空、井间、高频层状模型），并与独立的一维、三维数值解进行详细对比，全面验证了算法的正确性和鲁棒性。

六、 其他有价值的内容

论文在最后讨论了该频域方法与有限差分时域（FDTD）方法相比的优缺点。虽然频域方法因需处理复数运算、存储大型矩阵而内存需求更大，但对于频率域非线性反演问题（本研究算法的最终应用目标）以及需要对数频率采样或研究频散效应的情况，频域方法更具优势。此外，论文还指出了未来可能的改进方向，包括：开发更好的预处理器（如多重网格）、系统研究网格拉伸参数的选择、改进不连续面附近场的插值方法、处理低频模拟中的通道电流问题、以及更好地处理空气-大地界面等。这些讨论为后续研究指明了方向。