基于频域有限差分的三维大地电磁的低感应数预条件求解方法研究

蒋涛 （油气资源与勘探技术教育部重点实验室 （长江大学），湖北 武汉430100 长江大学工程技术学院，湖北 荆州434020）

胡文宝 （油气资源与勘探技术教育部重点实验室 （长江大学），湖北 武汉430100）

正演数值计算的算法研究是地球物理勘探的基础之一，只有充分了解和认识了已知场源所产生的地球物理异常，才可能从获得的数据中反演引起异常的异常体电磁学性质。常用的数值计算方法主要有基于积分方程的积分方程法、边界单元法和基于微分方程的有限差分法、有限元方法。随着计算机技术的发展和新算法技术的引进，近十年来大地电磁三维正反演取得了丰硕的成果［1，2］。有限差分法因其物理意义明确、简单直观、计算量小、计算速度快被很多学者选用［3～5］①Newman G A，Alumbaugh D L．Three-dimensional electromagnetic modeling on massively parallel computers．Sandia National Laboratories Report，SAND96-0582，1996．。谭捍东等［6］采用双共轭梯度稳定解法求解三维交错采样有限差分离散后得到的大型系数矩阵方程组，从迭代收敛稳定、计算精度高、速度快等方面，通过2个理论模型检验了该算法的正确性和计算精度。肖骑彬等［7］从网格精度、形成系统方程的方式、边界条件以及预条件线性算子等，对大地电磁有限差分数值解作了对比，结果表明，合适的预条件再配合好的线性算子，不仅可以加速收敛，而且可以降低迭代次数。李焱等［8］为了更有效地提高大地电磁三维正演的计算速度，根据各频率之间求取电磁场值的过程是相互独立的特点，引入并行处理技术，对2个理论模型的大地电磁响应采用三维交错网格有限差分数值模拟，高效地实现了正演计算。在上述研究中，如何快速、准确地求解经过有限网格离散后的大型稀疏线性方程组成为关键。

该次研究采用Yee［9］提出的交错网格有限差分法导出了频域Maxwell方程的离散化关系式，将三维电磁响应的模拟问题转化成求解大型线性代数方程组的问题。笔者将总场分离成一次场 （背景场）和二次场 （散射场），背景场用解析方法求取，二次场用有限差分实现。导出的线性代数方程组采用预条件方法的Krylov子空间类迭代方法求解。

为了解决频域低于1MHz时采用频域有限差分方法得到的线性方程组在迭代求解过程中收敛缓慢的问题，Newman等［3，4］提出了一种能显著减少三维频域有限差分数值计算时间的方法，即基于电场的Helmholtz分解定理的低感应数预条件方法。Maxwell微分方程经过频域有限差分网格离散后得到的线性方程组的系数矩阵具有大型复数稀疏对称的性质，广泛采用QMR （quasi-minimal residual，准最小残差法）进行迭代求解。而低感应数预条件能显著加快线性方程组迭代求解过程中的收敛速度，比简单的Jacobi预条件加速时间快1／10，可以将迭代次数从未做预条件处理的上千次降低到几十次。笔者在此基础上，研究了在利用QMR迭代求解线性方程组的过程中，为构建低感应数预条件矩阵，2种不同算法 （共轭梯度法和最小残差法）对正演计算迭代次数和时间的差别，分析了原因并给出了结论。

1 频域有限差分法

对于谐变电磁场，采用时谐因子eiωt（其中，i为虚数；ω为角频域，rad／s；t为时间，s）在各向同性介质中散射电场Es分布的控制方程可以从Maxwell方程组获得：

式中：△为哈密顿算子，后面接 “×”表示旋度；ES为散射电场，V／m；μ为介质的磁导率，H／m；σ为介质的电导率，S／m；σ0为背景介质的电导率，S／m；EP为背景介质中的一次场，V／m。按照Newman等［3］的方法，式 （1）在交错笛卡尔网格中离散，ES位于元胞边的中心，利用中心差分代替微分，可得到大型病态线性方程组：

式中：A为经过频域有限差分离散后得到的大型复数对称稀疏矩阵，每一行最多含有13个非零元素；E为散射电场矢量；b为与源或者边界条件相关的矢量。线性方程组的左边第m行可以写成二矢量之积αTmum（其中，αm为系数矩阵A的第m行13个非零元素组成的向量；u为待求电场，V／m）。在磁导率变化不大的条件下 （μ＝μb）（其中，μb为背景介质的磁导率，H／m），X方向的散射电场的线性方程可写为αTmum：

式中：i、j、k分别为直角坐标X、Y、Z方向上的网格节点序号；Δ为系数矩阵A的对角元素；δ为相邻网格节点之间的距离，m；δ＊为相邻网格节点中心之间的距离，m；Ex、Ey、Ez分别为待求的散射电场三分量，V／m；x （i）、y （j）、z （k）分别为X、Y、Z方向上网格节点坐标，m。

系数矩阵A的对角元素为：

式中：σ（i＋1／2，j，k）为沿着元胞棱边的4个相邻单元按体积权重所取的电导率，S／m。

式 （4）～ （6）部分地给出了式 （2）的复数对称系数矩阵，容易通过无矩阵内积方法实现。在利用QMR的每一步迭代过程中，需要计算一次矩阵－矢量积。近似解xi收敛到实际解的速度依赖于条件数（其中，下标P为LP范数），这里系数矩阵条件数很大，所以收敛过程缓慢。

2 预条件技术

对中高频率的大地电磁波，采用简单的Jacobi预处理的QMR能够有效地求解式 （2）。但是，当场源频率低到接近似稳场时，用频域有限差分离散得到的方程组条件数很大，若用Krylov子空间类迭代方法求解式 （2），则因为方程的病态程度严重，使得迭代过程的收敛速度很缓慢。Druskin等［10］利用与频域有限差分近似的SLDM （spectral Lanczos decomposition method）报道了求解过程中的问题。为了加快迭代收敛的速度，必须要使用预条件技术 （又称预优技术）。

Newman等［4］和Weiss等［12］根据Helmholtz原理把电场 （E）分解成无散场 （F）与无旋场 （△f）之和，即：

将式 （7）代入式 （1），“△×△×”算符可化为矢性拉普拉斯算符：

式中：▽后不接符号表示梯度。

而式 （9）的解并不满足电流的连续性方程，为了满足该条件，对式 （1）取散度然后代入式 （7）可得：

式中：▽后接 “·”表示散度。

这样，当频率很低即低感应数条件满足时，经过交错网格频域有限差分离散，可采用共轭梯度类算法先后求解式 （9）、（10），可以使式 （8）的近似解达到限制的误差范围内。

改善收敛缓慢的方法是，定义一预条件矩阵M，把式 （2）的求解可以有效转换为求解（AM－1）x＝b（其中，x＝ME为解矢量）。若M为A的主对角线，则M称为Jacobi预条件矩阵。如果A＝A1＋A2，当A1比A2占优势时，可取M ＝A1为预条件矩阵。在这里经过预条件后的线性方程组变为：

其中：A1E＝x。

这种预条件可以被应用到频率较低的准静态场的计算中，因为此时电场的旋度部分比起散度部分占优势。Newman等［4］提出的低感应数预条件方法正是基于此，将无散场部分F的扩散方程的解与无旋场▽f的泊松方程的解之和作为预条件矩阵。低感应数法对QMR迭代收敛速度的加速效果可以使迭代次数从未做预条件处理的上千次降低到几十次，优于不完全LU分解和简单Jacobi预条件方法。

3 低感应数预条件矩阵的求取方法

由于线性方程组的系数矩阵A是一个大型稀疏复对称的非埃尔米特矩阵，Krylov子空间类迭代方法是求解该类线性方程最为有效的方法。常用的求解稀疏复线性方程组的Krylov子空间类迭代方法主要是QMR，因为其稳定性好且收敛速度也较快。在利用QMR求解线性方程组时所需的低感应数预条件矩阵的求取方法主要有2种：共轭梯度法 （conjugate gradient）和最小残差法 （minimum residual）。该次研究对采用2种方法得到的低感应数预条件矩阵对大地电磁三维正演结果的影响进行了对比研究（2种方法的详细内容见文献 ［5］）。

共轭梯度法和最小残差法均属于共轭方向算法的一部分，等同于使误差函数在每次迭代过程的子空间中最小。考虑二次函数f（x）＝xTHx－hTx（其中，H为对称正定矩阵；h是线性方程组Hx＝h的右端向量）的最小值出现在＝H－1h （其中，是函数f（x）取最小值时的向量）并且值为f（）＝－hTH－1h （其中，f（＾x）是函数f（x）的最小值），因此找到二次函数的最小值相当于找到方程Hx＝h的解。

3.1 共轭梯度法

若A为对称正定，对于共轭梯度法有H＝A和h＝b。因此在最小化f （x）的每一步骤中，解矢量x和残差矢量r＝b－Ax根据下列方程更新：

式中：αi是沿着pi方向进行第i次迭代搜索时的步长；pi为第i次搜索方向。

选择搜索方向pi，使得函数f （x）对于沿着搜索方向的距离保持不变。当rTi＋1pi＝0，给定，可以得到搜索方向的更新方程：

式中：βi是共轭梯度方向更新因子。

选择βi使得Api＝0（搜索方向是线性独立且与A是共轭的），同时也得到：

通过选择p0＝r0，最小值出现在由 ｛r0，Ar0，A2r0，… ｝分隔的Krylov子空间中。

3.2 最小残差法

而最小残差法是选择H＝ATA和h＝ATb来最小化函数f （x）。在该方法中，对于x、r、p的更新与共轭梯度法一致，α和β需要重新定义，以使得Api＝0和 （Api＋1）TApi＝0，这也就给出了：

最小残差法与共轭梯度法一样，也是选择p0＝r0，使最小值出现在由 ｛r0，Ar0，A2r0，… ｝ 分隔的Krylov子空间中。

4 模型研究

为了比较共轭梯度法和最小残差法2种求解算法在构建低感应数预条件矩阵后对正演结果的影响，笔者设计了一种理论模型进行试算。设计的三维棱柱体模型如图1所示，几何尺寸为4km×4km×0．1km，电阻率为1Ω·m，顶面埋深为1000m，围岩电阻率为100Ω·m。

图1 三维棱拄体模型

研究了大地电磁频率为10Hz时 （因大地电磁的TE（正交极化或水平极化）、TM （平行极化或垂直极化）模式对该次研究没有影响，故未区分），在QMR迭代计算地表电场的过程中，分别采用最小残差法和共轭梯度法获取的低感应数预条件矩阵作为预条件，研究相对误差和时间随每一次迭代过程的变化关系。从图2、3可以看出，在每一次线性方程组的迭代求解过程中，采用最小残差法获得的相对误差及需要的时间都要小于共轭梯度法，并且线性方程组的迭代次数比共轭梯度法所需要的62次少1次。究其原因在于每一次迭代过程中需要通过先后计算矢量势 （无散场部分F）和标量势 （无旋场部分▽f），从而来构建低感应数预条件矩阵，而利用最小残差法需要的迭代次数均要小于共轭梯度法，如图4、5所示。

图2 2种方法求低感应数预条件矩阵时相对误差随QMR迭代次数变化的曲线对比

图3 2种方法求低感应数预条件矩阵时求解时间随QMR迭代次数变化的曲线对比

图4 2种方法分别求标量势时的迭代次数与QMR迭代次数的关系

图5 2种方法分别求矢量势时的迭代次数与QMR迭代次数的关系

5 结语

在对研究区域经过长方体网格的频域有限差分离散后，三维大地电磁的正演计算，时间主要耗费在线性方程组的迭代求解过程中，以此来达到所要求的计算精度。为了加快迭代收敛的速度，必须使用预条件技术。当场源频率低到接近似稳场时，低感应数预条件矩阵具有明显的加速效果。采用最小残差法构建低感应数预条件矩阵后，线性方程组的迭代次数和和求解时间都较同一类型的共轭梯度法要少，具有一定的速度优势；同样，在构建低感应数预条件矩阵的过程中，最小残差法需要先行计算的标量势和矢量势的迭代次数均比采用共轭梯度法所需的迭代次数要少。

［1］Siripunvaraporn W，Egbert G，Lencury Y．Three-dimensional magnetotelluric inversion：data-space method ［J］ ．Physics of the Earth and Planetary Interiors，2005，150：3～14．

［2］Siripunvaraporn W，Sarakorn W．An efficient data space conjugate gradient Occam’s method for three-dimensional magnetotelluric inversion ［J］．Geophys J Int，2011，186：567～579．

［3］Newman G A，Alumbaugh D L．Frequency-domain modeling of airborne electromagnetic responces using staggered finite differences［J］．Geophys Prosp，1995，43：1021～1042．

［4］Newman G A，Alumbaugh D L．Three-dimensional induction logging problems，Part 2：A finite-difference solution ［J］ ．Geophysics，2002，67：484～491．

［5］Mackie R L，Smith J T，Madden T R．Three-dimensional electromagnetic modeling using finite difference equations：The magnetotelluric example［J］．Radio Sci，1994，29：923～935．

［6］谭捍东，余钦犯，Booker John，等 ．大地电磁法三维交错采样有限差分数值模拟 ［J］．地球物理学报，2003，46（5）：705～711．

［7］肖骑彬，赵国泽 ．大地电磁有限差分数值解对比 ［J］．地球物理学报，2010，53（3）：622～630．

［8］李焱，胡祥云，杨文采，等 ．大地电磁三维交错网格有限差分数值模拟的并行计算研究 ［J］．地球物理学报，2012，55（12）：4036～4043．

［9］Yee K S．Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media ［J］．IEEE Trans Ant Prop，1996，14：302～307．

［10］Druskin V L，Knizhnerman L A，Lee P．New spectral Lanczos decompostion method for induction modeling in arbitrary 3-D geometry［J］．Geophysics，1999，64：701～706．

［11］Weiss C J，Newman G A．Electromagnetic induction in a generalized 3Danisotropic earth，Part 2：The LIN preconditioner ［J］．Geophysics，2003，68：922～930．