0 引言

随着我国经济的快速增长，矿产资源的消耗逐渐增多，地质条件有利地区的矿产资源已近枯竭，目标转向广大无人区。由于这些地区施工困难，常规的地面电磁法很难满足勘探需求，航空电磁法勘探因此成为重要的勘查技术手段^[1]。三维积分方程正演是航空电磁数值模拟的重要技术，它通过对异常区域进行剖分，建立异常区域内电场满足的第二类Fredholm积分方程，通过直接求解或利用迭代技术求解离散后的矩阵方程得到空间电磁场。在小规模异常体模拟中，系数矩阵占用存储小，求解速度快。然而，对于大型异常体，由于形成的线性方程组系数矩阵的密实性，计算耗时多、存储量大，限制了积分方程算法的应用。

为了避免大型线性方程组的求解，提高计算效率，学者提出许多近似方法。典型的方法包括扩展Born近似^[2]、准线性近似(quasi-linear or QL)^[3]、准线性级数(QLSE)^[4]、准解析(QA)近似和准解析级数(QASE)^[5]等。然而，上述近似方法只有在满足其特定的使用条件下才能获得理想的计算结果。通过对QL技术的改进，Zhdanov等^[6]实现了局部准线性近似(LQL)的频率域电磁正反演，并在对加拿大拉布拉多沃伊西湾镍-铜-钴矿航空电磁数据反演中取得了良好的效果。然而，该算法需要通过求解最小值问题来获得电反射系数，计算效率低。刘永亮等^[7]提出了快速准线性近似方法，并将其应用于三维频谱激电反演研究，取得了很好的效果。针对多源问题的海洋可控源电磁正演模拟，Ueda等^[8]结合积分方程(IE)和QL技术，提出了基于多重网格的准线性(MGQL)近似算法。MGQL算法的关键是假设背景场与异常场之间存在线性关系。基于MGQL技术，人们首先在粗网格上求解积分方程得到总场和异常场，进而求解粗网格上的电反射系数，并利用该系数在细网格上做线性插值得到离散后的电反射系数。在细网格上，利用电反射系数与背景场求解异常场，最后在细网格上获得接收点处的场值。与全积分方程法相比，对于大型异常体和网格剖分较多的情况，MGQL通过粗化网格可以大量节省计算时间。本文将该技术引入到多源航空电磁三维正演模拟中，由于反射系数计算在粗网格上进行，而细网格上的场值通过插值获得，因此可极大提高多源航空电磁的计算效率。

本文在MGQL算法的基础上，利用计算电磁学中的处理技术^[9]优化粗网格上积分方程的计算过程，将发射源首次发射时的格林函数存储为Toeplitz矩阵形式，避免了不同源格林函数的重复计算，进而利用快速傅里叶变换(FFT)技术实现矩阵与向量的乘积^[10]，提升MGQL的计算效率。MGQL算法网格粗化与计算精度有关，因此在精度允许的范围内，调节粗网格大小可有效降低存储量、加快计算速度。最后，通过将本文算法结果和已发表的典型模型计算结果对比验证精度，并通过与开源代码MarcoAir对比验证本文MGQL算法的有效性。

1 MGQL正演理论

根据积分方程理论，在水平层状各向同性介质中，对麦克斯韦方程组进行变换，可得电场E与磁场H的如下求解方程^[11-12]：

式中：D表示异常区域；为电并矢格林张量；为磁并矢格林张量；r和r'分别表示源和接收点位置；Δσ为异常电导率σ与背景电导率σ_b之差，即Δσ=σ-σ_b；E^b(r')和H^b(r')分别表示观测点处的背景电场和磁场。总电场可分为异常场E^a(r)与背景场E^b(r)两部分，即

根据Zhdanov等^[3]提出的准线性近似理论，假设背景场和异常场之间存在线性关系，即

式中，表示电反射系数。由式(1)、(3)、(4)可得异常场的准线性近似表达式为

式中：r_j表示观测位置；为单位张量。当r_j∈D时, 结合(4)和(5)式，可得电反射系数满足如下方程：

式中，E^B(r_j)表示异常场的Born近似^[13]，其表达式为

利用最小化求解方法，由式(6)可得电反射系数的求解表达式为

准线性近似求解方法的优势在于：可在一个剖分较粗的网格上求解电反射系数，进而通过对插值离散，实现细网格上电磁场计算^[14]。因此，该算法可以有效地应用于多重网格。MGQL算法使用两种网格(剖分较粗的网格和剖分较细的网格)，对电导率分布模型和电场进行离散。首先，我们利用积分方程在粗网格上求解出总场E(r_c)，则由式(3)可得异常场为

式中，r_c表示粗网格单元中心。则由式(4)可得粗网格上r_c处的电反射系数。需要指出的是，通常含有9个未知分量，不同学者基于不同假设分别给出相应的求解方法^{[3, 7]}。本文假设是对角张量，即

因此，式(4)可以写成独立的标量式，即

必须指出的是，依据公式(4)中的假设，反射系数计算虽然简单，但是如果背景场的某一个分量为零，将导致反射系数无法获取。为此，本文根据Gao等^[15]提出的针对各向异性介质三维电磁感应测井的模拟算法，假设异常场和层状半空间背景场绝对值存在如下线性关系：

式中，向量λ(r)=(λ_x, λ_y, λ_z)。式(12)表明，在粗网格上，由于|E^b(r_c)|≠0，则可以得到电反射系数的计算公式为

利用式(13)计算λ(r_c)之后，通过线性插值计算细网格上的λ(r_f)，则由式(12)可得单元r_f中心处的异常场为

而总场由式(3)给出，最后利用式(2)可计算观测点处的电磁场^[16-19]。

2 模型试验和结果分析 2.1 精度验证

针对Newman等^[20]给出的三维航空电磁模型，我们将本文算法的结果与其给出的全积分方程结果进行对比，以验证MGQL算法的精确性。模型如图 1所示，粗网格剖分数为2×10×10(x×y×z，下同)个，细网格剖分数为4×10×20个。测线沿着x方向，发射频率为900 Hz，系统飞行高度20 m。三维异常体的几何参数和埋深如图 1所示。

图 2a、b给出利用本文算法和Newman文中IE计算的电磁响应对比结果，可见MGQL的磁场H_x、H_z实虚分量曲线与Newman的结果吻合很好。图 2c、d给出的相对误差均保持在5%以内，说明本文算法具有很高的精度。

2.2 单个异常体模型

航空电磁勘探在地形条件比较复杂的区域应用较多，目标体上方的覆盖层会对异常响应产生影响^[21-22]。因此，本节采用的模型与图 1类似，只是在地表加一层电阻率为10 Ω·m、厚度为10 m覆盖层，异常体为30 m×180 m×70 m, 电阻率为1 Ω·m，埋深40 m，粗网格剖分网格数为2×10×10个，细网格剖分网格数为4×10×20个。测线沿x方向，飞行高度30 m，频率900 Hz。

图 3给出含覆盖层与不含覆盖层的单个异常体航空电磁响应对比。由3图可见：覆盖层不影响异常响应曲线形态，但是异常响应幅值发生明显变化，含覆盖层异常体电磁响应明显大于不含覆盖层的电磁响应，说明良导覆盖层增强了异常体的电磁响应(背景值得到提升)。

2.3 多个异常体模型

在三维航空电磁正演计算中，多个异常体之间存在电磁耦合。利用Zhdanov等^[14]提出的不均匀背景电导率法(IBC)，当存在多个异常体时，对每个异常体使用积分方程进行计算，然后通过耦合迭代求解电磁响应，但这种方法计算速度较慢。本文采用MGQL算法，由于积分是在粗糙网格上进行的，耦合迭代速度极大地提升。这种方法对网格数巨大的模型求解优势更为明显。为验证MGQL算法对多个异常体模拟的有效性，我们设计如图 4所示的模型，频率为900 Hz，粗网格剖分数为4×10×10个，细网格剖分数为8×10×20个，测线沿着x方向，飞行高度30 m。

为了方便对比，本文分别计算存在单个异常体和多个异常体时电磁响应。图 5a、b给出存在单个异常体时H_x和H_z响应，而图 5c、d给出存在两个异常体时H_x和H_z响应曲线。由图 5可以看出，由于两个异常体之间存在电磁耦合，单个异常体电磁响应的叠加明显大于两个异常体组合的电磁响应；说明异常体之间的二次耦合不可忽视。

3 计算效率分析

本文提出的MGQL算法的最大优点是在保证计算精度的前提下达到快速计算的目标。为了分析MGQL算法的计算效率，本文将其与开源软件MarcoAir (Version 2.3.1)及全积分方程算法结果作对比。本文设计与图 1类似的模型，异常体大小为30 m×200 m×100 m，电阻率1 Ω·m，埋深30 m，飞行高度30 m, 频率900 Hz，测线沿x方向。通常，全矩阵存储积分方程法需要占用大量的内存空间，用于存储格林系数矩阵。受内存限制，全矩阵存储积分方程法计算网格数有限。使用MGQL技术能够减少网格数、节约存储，可实现多网格数剖分的正演模拟。本文的计算在Intel(R) Xeon(R) CPU E5-2667 v3@3.20 Hz处理器上进行，内存128 GB，操作系统为64位Win7系统。

MarcoAir属于早期研发的积分方程算法，它利用格林函数的对称性和块迭代算法提升积分方程的计算速度^[22]。当异常体剖分单元较少时，直接解法优势较大；当剖分单元较多时，迭代法计算速度较快。我们对比时选用MarcoAir迭代法求解线性方程组。本文MGQL算法中将格林函数存储为Toeplitz矩阵形式，同时利用FFT实现矩阵与向量乘积，并采用稳定双共轭梯度法(BICGSTAB)求解线性方程组^[23-24]，以进一步加快计算速度。假设航空电磁系统发射源和接收机共移动60次。表 1给出在相同计算精度条件下，本文算法计算时间与MarcoAir的对比结果。由表 1可见：当网格数较少时，MarcoAir求解速度快于MGQL算法；随着模型剖分网格数增加，MGQL算法计算速度远远超过MarcoAir。例如，当网格数增加到10×20×20时，MGQL的计算速度是MarcoAir近54倍。另外，对于计算设备的有限存储空间和时间，MarcoAir可以计算的网格数不超过20×40×40个，而对于MGQL算法来说，可以计算的网格数为50×100×100。由于MGQL在粗网格上求解异常场，计算量极大减少，使得该算法在计算效率上拥有很大的优势。

对于表 1设计的模型，我们进一步在MGQL与全积分方程法(全矩阵存储和直接矩阵向量乘积算法)之间进行计算效率对比，结果如图 6所示。对比结果表明，MGQL算法明显具有计算速度快、占用内存少的优势。随着剖分网格增加，全积分方程计算时间和占用内存呈指数增加；但对MGQL算法，随网格数增加，其计算时间和占用内存增长缓慢，凸显出MGQL算法的高效性。

下面说明MGQL算法如何通过粗细网格结合，实现航空电磁响应高效三维数值模拟。采用的模型与表 1相同，航空电磁系统发射源和接收机共移动60次，异常体离散为2×10×10个粗网格和4×10×20个细网格。我们分别采用如下3种方法：1)在粗网格上应用积分方程法(BICGSTAB-FFT)；2)在精细网格上应用积分方程法；3)在粗细多重网格上应用MGQL算法。图 7是3种方法的计算的H_z响应对比。由7图可以看出：方法1)的结果与方法3)的计算结果差距较大，而方法2)和3)的响应曲线吻合度高。3种方法的计算时间分别为13.12、69.61和16.06 s。很显然，本文算法在保证计算精度的同时，大大提升了航空电磁数据正演模拟的计算效率。

4 结论

本文提出的MGQL算法有效地解决了积分方程计算时占用内存大、计算速度慢的问题。通过移动场源航空电磁三维模型试验，并与已发表的算法结果对比，表明响应误差均小于5%。另外本文算法模拟多个异常体也取得很好的结果，说明MGQL算法在满足精度要求的条件下能够求解大型复杂异常体模型，是解决多源航空电磁问题的有效算法。

在相同模型保证计算精度的前提下，本文算法相比于MarcoAir算法和全积分方程法，网格数越多，相对计算时间越少、占用存储越小，说明利用矩阵Toeplitz性质和BICGSTAB-FFT有效地提升了计算效率。探究精算精度与网格的关系发现，采用组合网格可以使用最小的内存和计算时间获得高的计算精度。因此，未来我们将尝试利用MGQL算法进行多源航空电磁数据的快速成像和反演。

致谢: 感谢开源软件MarcoAir的作者，感谢吉林大学“千人计划”电磁研究团队全体成员。