本文由中南大学的张莹、戴世坤、张钱江与桂林航天工业学院的凌家煊合作完成，于2024年11月7日以开放获取形式提前在线发表于国际地学期刊《Geophysical Journal International》（简称Geophys. J. Int.）2025年第240卷。论文标题为“一种考虑剩磁与退磁效应的磁场计算迭代算法”。

一、 学术背景

本研究属于地球物理勘探领域，具体涉及磁法勘探的正演模拟。磁法勘探是矿产资源勘查的重要手段。传统的磁异常解释方法通常忽略剩磁（remanent magnetization）和自退磁（self-demagnetization）效应。这在磁化率较低的情况下通常是可接受的。然而，在高磁化率条件下，磁性体被磁化后会产生一个与磁化场方向相反的自身退磁场，导致磁性体内部的有效磁化强度降低和方向改变，即自退磁效应。同时，高磁化率的岩石常伴随较强的剩磁，二者共同作用会显著改变总磁化强度的幅度和方向。忽略这些因素，会导致磁测数据处理与反演解释产生严重误差，例如低估矿体的磁化率和剩磁强度，甚至错误定位矿体，无法满足高精度勘探的需求。

在高磁化率条件下进行大规模、高精度的数值计算是磁法勘探应用的关键。现有的数值计算方法主要分为空间域和频率域两大类。空间域方法（如体积分、表面积分、有限体积法）在计算复杂形态或考虑两种效应时，往往计算量大、效率低。频率域方法（如快速傅里叶变换）能提高效率，但通常仅适用于低磁化率情况，且存在谱混叠和截断效应等问题，对于同时考虑强剩磁和强退磁的情况仍存在局限。

因此，本研究旨在提出一种能够同时、高效、高精度地计算强剩磁与强退磁条件下三维磁异常场及磁梯度张量的迭代算法，以弥补现有方法的不足，为复杂地质模型的高精度磁法正演与反演提供有力工具。

二、 研究详细工作流程

本研究是一个计算地球物理方法的开发与验证过程，不涉及传统意义上的实验对象和样本，其“研究流程”核心是算法的构建、实现与测试。具体流程如下：

1. 算法理论基础构建 该流程基于静磁场的麦克斯韦方程，从磁位满足的三维泊松方程出发。方程中，磁化强度包含了由背景场、异常场（即退磁场）诱导的感应磁化部分以及剩磁部分，完整地表达了退磁与剩磁效应。这是本算法的物理基础。

2. 维度缩减与方程变换 这是算法的核心创新步骤之一。为了高效求解三维偏微分方程，研究团队对水平方向（x, y）进行了二维傅里叶变换，将空间域的方程转换到空间-波数混合域。这一变换巧妙地将三维泊松方程转化为一系列相互独立的一维常微分方程，每个方程对应一个特定的水平波数对（kx, ky）。这一操作极大地降低了问题的计算维度，减少了计算量和存储需求，并具有高度的并行潜力。垂直方向（z）保持在空间域，便于根据实际模型灵活调整网格尺寸。

3. 边界条件设定与方程求解 在笛卡尔坐标系（z轴向下为正）下，为上述一维常微分方程设定了上、下边界条件。边界被设定为无源区域，分别应用了反映向上衰减波和向下衰减波的诺伊曼边界条件。 接下来，采用有限元法求解该边值问题。具体而言，在垂直方向对计算域进行离散，每个单元内采用二次插值形函数。通过单元分析和总体合成，将问题最终转化为一个五对角线性方程组。为了高效、稳定地求解此方程组，研究采用了“追赶法”。追赶法对严格对角占优的五对角方程组具有良好的数值稳定性，且计算和存储需求小，避免了中间结果阶数剧增和舍入误差累积的问题。求解后，即可得到空间-波数混合域中的异常磁位。

4. 磁场与磁梯度张量计算 基于磁位与磁场（磁场强度负梯度）以及磁梯度张量之间的数学关系，在波数域中可直接通过代数运算，由求得的异常磁位快速导出异常磁场三分量以及完整的磁梯度张量张量。最后，对水平方向进行二维逆傅里叶变换，即可得到空间域中的最终结果。

5. 迭代格式引入以确保收敛 由于在强磁条件下，未知的异常场（退磁场）对总磁化有不可忽略的贡献，方程无法直接求解。本研究采用了迭代方法：首先假设异常场为零进行初始计算，得到第一次逼近的异常场；然后利用该异常场更新总场，再代入方程进行下一次计算。为了确保迭代在高磁化率下仍能稳定收敛（避免传统逐次逼近法可能出现的发散），论文引入了一种基于电磁积分方程方法的收缩算子迭代格式。该格式源自修正的玻恩级数，已被证明在电磁场计算中能保证稳定收敛。研究将其等效地应用于此微分方程的迭代求解过程中。迭代的终止条件是相邻两次迭代总场绝对值之和的相对变化小于一个预设的极小值（本文设为10^-4）。

6. 算法实现与验证测试 研究团队使用Fortran 95语言实现了上述完整算法。验证与测试部分构成了研究的另一核心流程： * 正确性验证：设计了一个高磁化率（2 SI）且带有剩磁的球体模型，将算法的数值解与已知的解析解进行对比。计算区域为1000m立方体，网格规模为201x201x201。通过比较磁场三分量和磁梯度张量九个分量的数值解与解析解，并计算相对均方根误差来评估精度。 * 参数影响分析：使用球体模型，系统改变磁化率（从0.01 SI到100 SI）、剩磁强度、球体埋深、球体半径、网格节点数等参数，分析各因素对算法收敛所需迭代次数和最终计算精度的影响。 * 效率与内存分析：在相同计算机配置下，将本算法的单次迭代平均时间与一篇近期发表的基于积分方程的迭代算法（Ouyang et al., 2020）进行对比。同时统计了不同规模网格下的内存占用量。 * 边界适应性测试：设计模型，将高磁化率、强剩磁的异常体置于计算网格边界附近，验证算法在边界处的计算精度是否依然可靠。 * 复杂模型与地形适应性验证：构建包含多个不同磁化率和剩磁强度的异常体的复杂组合模型，分别计算不考虑退磁剩磁、仅考虑退磁、同时考虑退磁和剩磁三种情况下的磁场响应，以展示算法的适应性和同时考虑两种效应的必要性。最后，结合真实数字高程模型（DEM）数据，将起伏地形融入模型，通过设置多个不同深度的水平面进行计算和插值，验证算法在起伏地形条件下的适用性。

三、 主要研究结果

1. 算法正确性得到证实 对于高磁化率（2 SI）和强剩磁（20 A/m）的球体模型，算法经过7次迭代后收敛。数值解与解析解在磁场和磁梯度张量的所有分量上几乎完全一致，绝对误差比场值低两个数量级以上。各分量的相对均方根误差（RRMS）均低于0.5%，最高为0.35%（Baz分量）。这充分证明了该算法在强磁条件下的正确性和高精度。

2. 明确了自退磁效应的影响程度 通过改变球体磁化率（0.01至20 SI）并对比考虑与不考虑退磁的结果发现：当磁化率为0.01 SI时，自退磁效应可忽略不计；当磁化率增大到0.1 SI时，不考虑退磁的解与参考解（考虑退磁）的RRMS误差已达3.63%；当磁化率增至20 SI时，误差高达461.88%。这一结果量化地证明了，在磁化率大于0.1 SI的区域进行勘探时，必须考虑自退磁效应，否则将导致矿体储量估计严重偏低等解释错误。

3. 揭示了影响算法收敛的关键因素 参数敏感性分析显示：在测试的诸多因素（磁化率、剩磁强度、异常体埋深、半径、网格节点数）中，只有磁化率影响算法收敛所需的迭代次数。磁化率越高，达到收敛所需的迭代次数越多（从0.01 SI时的3次到100 SI时的51次）。而其他因素改变时，收敛迭代次数保持不变。这一发现明确了算法性能的核心制约参数。

4. 阐明了算法精度与网格、磁化率的关系 * 在相同网格配置下，磁化率越高，算法最终能达到的精度越低，且需要更多迭代次数收敛。例如，磁化率100 SI时，RRMS误差约为10%，而1 SI时则可低于1%。 * 在相同磁化率下，增加垂直方向或水平方向的网格节点数，都能有效提高算法的最终计算精度。 * 这意味着对于极高磁化率的模型，需要通过加密网格来获取更高精度的解。

5. 证明了算法的高效率与低内存占用 与基于积分方程的迭代算法相比，本算法在相同网格规模下显示出更高的计算效率。特别是在垂直节点数增加时，优势更加明显（当nx=ny=201, nz从51增至401时，本算法单次迭代时间优势从1.26倍提升至3.45倍）。在水平网格节点数很大时，优势虽有所减小，但仍优于对比算法。在201x201x401的大规模网格下，算法内存占用约为30 GB，远低于传统有限元法，使得在普通配置计算机上进行大规模数值模拟成为可能。

6. 验证了算法对边界异常和复杂模型的良好适应性 * 对于放置在计算区域边界附近的高磁化率强剩磁异常体（包括对称和非对称模型），算法依然能保持高精度（RRMS < 1%），边界效应可忽略。这对实际反演中常需在边界设置区域场源的情况具有重要意义。 * 复杂组合模型测试表明，仅考虑退磁或忽略两者都会导致与实际情况（同时考虑两者）显著不同的磁异常响应。同时考虑退磁和剩磁时，在强剩磁区域磁场值可增大数十倍，磁梯度张量能更清晰地反映异常体的分布细节。这证明了新算法处理复杂地质情况的能力。 * 起伏地形测试结果表明，地形效应会改变磁场的分布形态和强度，不可忽略。本算法通过分层计算与插值，能够有效地处理地形数据，获得起伏观测面上的磁场分布。

四、 研究结论与意义

本研究成功提出并验证了一种适用于强剩磁与强退磁条件的高效、高精度三维磁场及磁梯度张量迭代计算算法。主要结论如下： 1. 算法基于泊松方程，通过二维傅里叶变换降维、有限元离散、追赶法求解及引入稳定收敛的迭代格式，实现了复杂介质模型下磁场的高效稳定求解。 2. 算法正确性经解析解验证，精度高。 3. 磁化率是影响算法迭代收敛次数的唯一关键因素。 4. 算法在效率和内存占用上优于现有的积分方程类算法，适合在普通计算机上进行大规模计算。 5. 算法对边界附近的异常体、复杂地质结构以及起伏地形均有良好的适应性和高精度。

科学价值与应用价值：该研究解决了高磁化率、强剩磁条件下磁法正演模拟的数值计算难题，提供了一种兼具效率与精度的可靠工具。其科学价值在于深化了对退磁和剩磁效应耦合机制数值模拟的理解，推动了磁法勘探正演理论的发展。应用价值体现在：为高精度磁测数据的定量解释与反演提供了关键的技术支撑，有助于更准确地圈定矿体、估算资源储量，在金属矿产勘探，特别是在寻找强磁性铁矿、磁黄铁矿型矿床等方面具有广阔的应用前景。

五、 研究亮点

1. 创新性的算法框架：创造性结合了空间-波数域变换（降维）、一维有限元法（灵活离散）、追赶法（高效求解）和基于积分方程的收缩迭代格式（保证强磁收敛），形成了一套新颖、高效、稳定的混合域正演算法。
2. 全面考虑实际复杂因素：首次在一个算法中同时、有效地解决了强退磁、强剩磁、复杂模型形态、边界效应及起伏地形等多个在实际勘探中常遇到的难题。
3. 系统的性能评估：不仅验证了正确性，还深入、系统地分析了磁化率、网格等各种因素对算法收敛性和精度的影响规律，为算法使用者提供了清晰的参数选择指导。
4. 实际应用导向明确：所有测试模型的设计均紧密围绕实际勘探问题（如边界反演、复杂矿体、地形影响），证明了算法的实用性和工程应用潜力。

六、 其他有价值内容

论文的附录详细给出了算法涉及的变分问题推导、有限元单元矩阵的详细形式以及利用形函数求取磁位垂直方向导数的方法，为其他研究者理解和复现该算法提供了完整的数学与实现细节。作者同时表示，未来工作可将此正演算法进一步发展，通过计算灵敏度矩阵并结合正则化技术，将其应用于高磁条件下的磁化率与剩磁强度联合反演，从而进一步提升其在地球科学应用中的价值。