本文由Ying Zhang（张莹）、Jianhui Li（李建辉，通讯作者）、Shuang Liu（刘双）以及Kun Li（李昆）等人合作完成。其中，第一作者与通讯作者来自中国地质大学（武汉）地球物理与空间信息学院及地下多尺度成像湖北省重点实验室，同时该单位还隶属于地质过程与矿产资源国家重点实验室；李昆来自西南石油大学地球科学与技术学院。该研究发表于由Oxford University Press代表Royal Astronomical Society出版的期刊 *Geophysical Journal International (GJI)*，发表时间为2026年。
在高磁化率矿体中，当磁化率（magnetic susceptibility）超过0.1 SI时，自退磁场（self-demagnetization field）的强度可与背景地磁场（geomagnetic field）相当，从而显著削弱并重新分布内部磁化强度。与此同时，许多含磁铁矿和钛磁铁矿的岩石由于构造与变形作用，矿物往往具有优势取向，其磁性更适合用各向异性磁化率（anisotropic magnetic susceptibility, AMS）张量而非标量来描述。上述两种效应还会与剩磁（remanent magnetization）相互耦合，进一步增加三维磁法正演与反演的难度。
已有研究大致可归为两类：一类聚焦于强退磁与强剩磁条件下的磁场数值模拟，但大多假设介质为各向同性；另一类关注各向异性对磁异常的影响，却往往简化甚至忽略强退磁效应与复杂三维几何条件。因此，目前尚缺乏能够同时考虑各向异性磁化率张量、强自退磁和剩余磁化强度的统一高效数值方法。基于此，本文提出一种迭代积分方程算法，旨在同时耦合上述三种因素，为复杂强磁性各向异性地质体的高精度磁法正演提供可靠工具。
本研究的方法流程可归纳为以下六个主要步骤。
第一步：建立磁化强度方程。 研究从基本磁化关系出发，将总磁化强度向量m表达为磁化率张量χ、背景磁场h₀、自退磁场（异常场）hₐ以及剩余磁化强度mᵣ的函数。背景磁场与剩磁的方向由磁倾角（inclination）和磁偏角（declination）定义。
第二步：构建磁化率张量。 各向异性介质的磁化率由三个主磁化率（χ₁、χ₂、χ₃）和三个欧拉角（Euler angles，α、β、γ）唯一确定。通过旋转矩阵R及对角主磁化率张量χ₀，利用变换χ = Rᵀχ₀R得到全局坐标系下的对称二阶张量。该旋转矩阵参考了Takahashi & Oliveira Jr. (2017)中的定义。
第三步：推导磁位积分方程与空间-波数域变换。 异常磁势（anomalous magnetic potential）uₐ与磁化强度m满足空间域积分方程。通过对水平方向实施二维傅里叶变换（基于标准FFT的Gauss-FFT方案），将三维积分方程降至一维垂向形式，有效降低了计算维度。
第四步：垂向形函数离散。 在一维积分方程中，垂向区间被划分为若干积分单元。研究分别构建了线性插值（linear shape function）与二次插值（quadratic shape function）。线性插值在两相邻节点间构造一阶多项式；二次插值利用三个垂向节点构造二阶多项式。文中分别推导了两种插值格式下的积分系数表达式，并对观测点位于积分单元上方、下方和内部三种情况进行分类求解。最终，通过求和获得空间-波数域异常磁势。
第五步：计算异常磁场并构建迭代方案。 利用hₐ = −∇uₐ及傅里叶变换的微分性质，在空间-波数域求得异常磁场分量，再通过Gauss逆FFT转换回空间域。为解决强磁条件下传统Born迭代不收敛的问题，本文借鉴Zhdanov和Fang (1997)以及Ouyang和Chen (2020)的改进格林算子思想，将原有各向同性迭代格式推广为张量形式。改进后的迭代格式为h⁽ⁱ⁺¹⁾ = h₀ + b(h⁽ⁱ⁾) + hₐ⁽ⁱ⁾ + d h⁽ⁱ⁾，其中b与d为与磁化率张量相关的矩阵。理论分析表明，矩阵d的谱半径始终小于1，结合改进格林算子的非扩张性（non-expansive property），证明该迭代方案在各向异性强退磁条件下可稳定收敛。此外，剩磁项在相邻两次迭代差值中被抵消，因此不影响迭代收敛速度，只影响最终异常场的幅值与空间分布。收敛判据设定为两次迭代异常磁场的相对均方根差小于10⁻⁴。
第六步：验证与测试对象设计。 为验证算法，研究选取具有解析解的椭球模型（Takahashi & Oliveira Jr. 2017）作为基准。模型区域为1000 m × 1000 m × 1000 m，网格规模为201 × 201 × 201。背景磁场强度为50000 nT，倾角45.9°，偏角5.9°。设计了三个模型：模型1为主轴各向异性椭球，三个半轴分别为300 m、100 m、50 m，主磁化率为[5, 3, 1] SI，欧拉角均为0°；模型2为任意各向异性椭球，欧拉角设为60°、45°、30°；模型3在模型2基础上进一步将椭球空间姿态旋转，走向、倾伏角和侧伏角分别设定为45°、75°和30°。所有模型剩余磁化强度为5 A/m，剩磁倾角10°，偏角45°。通过与COMSOL Multiphysics有限元模拟的基准模型对比，进一步验证了算法的正确性。
算法精度验证。 对于模型1、模型2和模型3，算法分别在17次、11次和14次迭代后收敛。三个磁场分量相对解析解的RRMS误差分别为0.67%/0.31%/0.41%、0.68%/0.54%/0.48%和0.52%/0.82%/0.67%，均低于1%。对比结果表明，欧拉角与空间姿态的旋转会显著改变磁异常形态、梯度和幅值。与有限元比较时，当计算域扩展因子kb从0增加到3，有限元结果的RRMS从94.55%降至0.87%，而本文算法无需扩边即可在目标区域内获得高精度结果。
插值阶数影响。 通过对比线性与二次插值，发现两种格式在各模型下的RRMS均低于1%。但二次插值在每个迭代步的耗时和内存占用均低于线性插值。例如模型1中，二次插值每次迭代耗时37.67秒、内存6.61 GB，而线性插值为39.68秒和7.93 GB，说明二次插值在保持精度同时具有更高计算效率和更低内存开销。因此，研究最终选用二次插值作为垂向离散方案。
主磁化率敏感性。 在各向同性条件下，随磁化率从1 SI增加到20 SI，RRMS维持在1%左右，但迭代次数从7次增至35次。在线性各向异性条件（λ=0）下，误差保持稳定，而迭代次数几乎仅由最大主磁化率χ₁决定。例如当χ₁=20时，无论χ₂和χ₃如何变化，迭代次数均保持为36。在固定χ₁的条件下改变形状参数λ（即调整χ₂在χ₁和χ₃之间的位置），误差随λ增加而下降，但迭代次数保持不变。这表明收敛速度主要受最大主磁化率控制，数值误差则更多受各向异性结构影响。
欧拉角影响。 在[3, 1, 1] SI主磁化率条件下，单独改变α或γ会引起误差与迭代次数的周期性波动，而改变β则几乎无影响。将主磁化率改为[3, 2, 1] SI后，β的影响变得显著。原因在于[3,1,1] SI具有轴对称性，而[3,2,1] SI打破了这种对称性。
剩磁与几何参数。 剩磁强度从1 A/m增加到50 A/m时，迭代次数保持17不变，误差也无显著变化，验证了理论分析中剩磁不影响收敛特性的结论。埋深和几何尺寸对迭代次数影响微弱，但会改变不同分量的误差分布：随埋深增加，bax误差上升而bay、baz误差下降；异常体尺寸越大，边界离散误差越小，各分量误差均下降。
收敛稳定性。 在不同水平、垂向及三维同步网格细化条件下，算法在高达100 SI的磁化率时仍能稳定收敛。迭代次数在低磁化率下基本不变，在高磁化率下仅在粗网格上有轻微波动后迅速稳定，证明算法具备极强的稳定性。
本文提出并验证了一种基于二次形函数插值的迭代积分方程算法，可在强自退磁与剩磁耦合条件下对各向异性磁化率介质进行高精度三维磁法正演。该方法通过二维傅里叶变换将积分方程降至一维垂向形式，并将改进的格林算子迭代格式推广至张量形式，实现了稳定收敛。数值实验表明，最大主磁化率是控制收敛速度的主导因素，而剩磁仅作为固定源项改变异常场的幅值和分布。与有限元方法相比，该算法在相同离散条件下所需运行时间和内存更少。应用于复杂模型时，磁化率各向异性和欧拉角旋转可显著改变磁异常形态。该研究为复杂三维强磁性各向异性地质体的高精度磁法正演提供了高效可靠的工具。
本研究的突出亮点在于：其一，将改进的收敛迭代算子从各向同性推广至各向异性张量形式，并通过谱半径分析严格证明了收敛稳定性；其二，系统对比了线性与二次插值在精度、时间和内存方面的表现，确立了二次插值的综合优势；其三，全面揭示了最大主磁化率、欧拉角、剩磁强度、埋深和几何尺寸等因素对误差和迭代次数的差异化影响机制，为该类算法的实际应用提供了重要指导。