基于矢量有限元的大地电磁快速三维正演研究

顾观文，武晔，石砚斌

(1.防灾科技学院 地球科学学院，河北 廊坊 065201; 2.河北省地震动力学重点实验室，河北 廊坊065201)

0 引言

大地电磁测深法(magnetotelluric sounding, MT)具有施工方便、勘探效率高、成本低 (相对于地震勘探)、勘探深度大等优点，目前已被广泛应用于资源勘查、能源勘探及深部构造探测等方面。对于大地电磁数据的反演和解释，由于受限于理论方法及计算设备的运算能力，早期以一维或二维反演技术为主。20 世纪90年代以来，随着计算设备运算能力的提高及数值计算方法的进步，大地电磁三维正反演技术开始逐步发展，不同的三维正演方法(积分方程、有限差分、有限元等)及其计算技术取得了巨大进展[1-15]。目前在实际中得到应用的三维反演技术主要是基于有限差分法三维正演的反演方法，特别是国内实测大地电磁资料的三维解释基本上都采用基于有限差分法的三维反演技术[16-19]。不同于有限差分法，有限单元法在模拟起伏地形以及复杂地质体的电磁响应方面具有明显优势，特别是近些年发展迅速的矢量有限元法，由于其能有效地解决传统节点有限元法存在的伪解问题，目前已成为复杂地形和复杂地质体三维电磁响应模拟的主要方法。但有限单元法也存在一些不足，运算量大、计算时间长是导致基于有限元法的大地电磁三维反演技术实用化进程相对滞后(相对于基于有限差分法的三维反演技术)的主要因素。为此，开展基于有限元的MT快速三维正演算法研究，提高三维正演计算效率，对于大地电磁三维反演技术的实用性具有重要的意义。

MT三维正演模拟最终需要求解复数大型线性方程组，如何快速、准确地求解此线性方程成为大地电磁三维正演实现中的关键工作，此线性方程的求解时间约占整个正演计算时间的80%[20]。目前求解大型线性方程组主要有两种方法，一类为Krylov子空间迭代法，另一类为直接解法。迭代算法从一个初始解出发，通过逐步修正来逼近真实解，在每一步的计算过程中只需要计算矩阵向量乘积或矩阵的转置与向量的乘积，因此其内存需求较小。但是，迭代算法可能存在不收敛的问题：一方面，方程组阶数的增加导致矩阵条件数增大；另一方面，求解低频电磁场问题时，方程的定解条件变弱，也加大了矩阵的病态性，使迭代算法收敛变慢，导致过多的迭代求解步骤，从而增加了计算时间。采用迭代法求解线性方程组时，需要引入散度校正技术以改善在低频范围(尤其在0.01 Hz及以下)求解的收敛性[5,11,21-22]。直接解法对计算机的内存容量要求较高，但直接法求解稳定，对于给定的方程组，直接解法总能获得其精确解，尤其是对于大地电磁的低频正演模拟问题，无需做散度校正也可获得高精度的数值解。值得关注的是，近些年直接求解器获得了非常大的突破，产生了有效的直接法求解器，如基于LU矩阵分解法和OpenMP的PARDISO(parallel direct solver)求解器[23]，基于波前法和MPI的MUMPS(multifrontal massive parallel solver)求解器[24-25]。基于直接解法在解决大规模电磁模拟问题方面取得了明显效果[13,26-28]。直接求解法得到了越来越多的重视，有逐步取代迭代求解器的发展趋势[29]。基于PARDISO直接求解器的直接求解法具有计算精度高和无需散度校正等优势，同时PARDISO求解器支持复数双精度数据类型、64位整数索引、并可搭配共享内存环境下的OpenMP并行方式，可以方便地应用在电磁有限元分析后的大型稀疏矩阵的求解中。

根据大地电磁场满足的麦克斯韦方程组推导了大地电磁三维正演的边值问题，利用矢量有限元法基于六面体网格单元对计算区域进行离散，采用无需散度校正的直接求解方法(PARDISO)求解矢量有限元法对应的大型线性方程组。在对三维地电模型进行数值模拟的过程中，通过数值解与解析解对比以及本文模拟结果与国际公认检验模型的计算结果对比，验证三维矢量有限元方法及程序的正确性。通过无需散度校正直接求解法与带散度校正的迭代求解法对比，证明本文快速三维正演算法的有效性。

1 大地电磁三维正演算法

1.1 大地电磁三维正演的边值问题

正演方法基于六面体网格剖分的矢量有限元法[11-12]，在大地电磁研究的频率范围内，忽略位移电流的作用。取时谐场为e-iωt，麦克斯韦方程组的微分形式表示如下：

(1a)

(1b)

(1c)

(1d)

式中：E为电场强度矢量，H为磁场强度矢量，σ是介质的电导率，μ是介质的磁导率，ε是介质的介电常数，q为自由电荷密度。对式(1a)两边取旋度，再将式(1b)代入到式(1a)，则将电磁场满足的一阶微分方程变为二阶微分方程，可得电场E满足下面方程：

(2)

求解偏微分方程组(2)得到电场分量Ex、Ey、Ez后，根据式(1a)求得磁场分量Hx、Hy、Hz。

由于空气中的电场和磁场受地形和不均匀体影响而不均匀，因此在数值模拟中必须包括空气层和地下介质。如图1所示，模拟区域设为Ω，由空气层和地下介质两部分组成，Ω=Ω1∪Ω2。在模拟区域中，空气层的电导率一般在10-6～10-10S·m-1之间，这时空气和地下介质的接触面变成内部边界。

图1 带地形三维MT数值模拟区域剖面示意[11]Fig.1 Section diagram of numerical modeling domain for 3D MT with topography[11]

将式(2)写成

(3)

取第一类边界条件：

E(x,y,z)|∂Ω=g(x,y,z)|∂Ω,

(4)

式(4)中g是边界上的矢量电场，可以用一维或者二维MT计算值[4,30]。这样，式(3)和式(4)构成了大地电磁三维正演的边值问题。

1.2 矢量有限元分析

用有限元法求解上述区域(图1)的电磁场问题需要对研究区域离散化，即对研究区域进行六面体网格剖分。如图2a所示，沿x、y和z轴方向分别剖分成Nx、Ny和Nz段，网格间距分别为Δx(i)(i=1,…,Nx)、Δy(j)(j=1,…,Ny)和Δz(k)(k=1,…,Nz)。可以推导，每个六面体单元的内部电场分量用六面体的12条棱边的场值分量(如图2b所示)通过插值求取，公式为

(5)

a—区域剖分示意; b—电场分量位置示意a—domain subdivision; b—location of electric field components图2 矢量有限元法的区域剖分示意Fig.2 Domain subdivision of the vector finite element method

式(5)写成矢量形式如下：

(6)

(7)

(8)

(9)

式中：坐标转换函数为ξ=(x-xc)/a,η=(y-yc)/b,ζ=(z-zc)/c；(xc,yc,zc)是六面体的中心坐标，2a、2b、2c分别是六面体x、y、z方向的长度；(ξi,ηi,ζi)的取值与棱边编号有关。

由式(3)定义矢量余函数：

(10)

把矢量基函数作为权函数，采用迦辽金方法[31-32]使整个域内的积分矢量余函数为最小，即：

(11)

(12)

利用矢量恒等式：

(13)

则式(12)右边第一项分为两项：

(14)

根据高斯散度定理，可知：

(15)

对于第e个单元，式(11)可化为

(16)

KeEe=Se,

(17)

式中：Se表示场源项；Ee表示棱边上的电场；Ke为单元刚度矩阵，是一个12×12阶的复数矩阵，可按下式解析计算[32]得出：

(18)

将每个单元电场满足的线性方程进行组合，可以得到整个计算域上电场满足的线性方程组：

K·E=s,

(19)

式中：K是系统刚度矩阵；E是整个计算域的网格单元棱边上的电场值向量；s是源向量，由计算域的上、下、左、右的边界场值与边界上的单元刚度矩阵计算得到。

1.3 线性方程组求解

方程(19)为复数大型线性方程组，PARDISO是针对大规模稀疏线性方程组开发的高效并行直接求解器，采用PARDISO直接求解器求解方程(19)。该求解器采用BLAS软件包，实现算法中线性代数操作的并行计算。大量数值测试表明，PARDISO是目前最快的线性稀疏矩阵求解方法之一[33]。Intel公司获得授权后，Intel®Math Kernel Library(Intel MKL)提供了PARDISO的优化版本，计算效率高，稳定性好。

1.3.1 PARDISO求解过程

PARDISO求解器基于LU分解，融合了向左和向右算法的优点，采用超节点消去树进行填充元优化，降低算法的时空消耗。对于复线性方程组有

Ax=b，

(20)

式中：A矩阵是一个稀疏矩阵，矩阵中的元索值多数为零，只有与当前节点直接相连的节点处为非零元素。利用PARDISO求解方程(20)的具体步骤如下：

③ 方程求解与迭代：利用LU分解结果，求解方程。如对结果精度有进一步要求，使用迭代法进一步提高精度。

④ 迭代结束，释放计算过程所占内存。

1.3.2 压缩稀疏矩阵存储

PARDISO求解器采用目前广泛流行的压缩稀疏行格式(compressed sparse row, CSR)对稀疏矩阵进行压缩存储，大大降低了矩阵元素的访问时间和空间存储成本。该方法以行为单位存储每个非零数据。对于一个对称、稀疏矩阵A，PARDISO对矩阵的存储包括三个数组：

① 双精度数组a—矩阵A上三角部分的非零元素。A的非零元素通过下面的ja与ia映射到a数组中。

② 整型数组ja—数组a中每个元素在矩阵A中的列号。

③ 整型数组ia—矩阵A中上三角部分每行第一个非零元素在数组a中的位置，其最后一个元素为上三角矩阵的非零元素总数加一。

以式(21)为例说明CSR存储格式：

(21)

大地电磁三维正演问题(式(3)、式(4))经有限元离散后，得到的刚度矩阵(式(19) 中的K)为稀疏、对称矩阵，可采用对称CSR格式存储。

根据大地电磁场线性方程组特有的稀疏性，开发有针对性的数据结构接口，将PARDISO快速求解器应用到方程组(19)的求解，得到计算域上网格单元棱边上的电场值，然后根据麦克斯韦方程组(式(1a))微分求取磁场。

1.4 视电阻率及阻抗相位的计算

根据Newman等[34]、谭捍东等[35]的研究，假设两种线性无关的场源激发的表面电场和磁场分别为Ex1、Ey1、Hx1、Hy1、Ex2、Ey2、Hx2、Hy2，由此可以计算三维MT的张量阻抗：

(22)

张量阻抗的每一个分量可表示为

(23a)

(23b)

(23c)

(23d)

式中：下标1、2表示极化模式，下标x、y表示x、y分量，E表示电场，H表示磁场，Z表示阻抗；式23b和23c定义的响应分别称为XY和YX模式响应，按照式24可求出三维介质的视电阻率和相位：

(24)

式中:i=X,Y;j=X,Y。

2 三维正演结果的正确性验证

为了验证本文三维正演算法的正确性，分别对水平地形地电模型和起伏地形地电模型进行正确性验证。对于水平地形条件下的验证，采用均匀半空间、典型层状模型和国际电磁学术讨论会COMMEMI研究小组Zhdanov等设计的COMMEMI 3D-2模型[36]，对于起伏地形条件下的验证，采用Wannamaker等设计的目前国际上普遍认可的二维山峰模型[37]。

2.1 水平地形条件下的模型验证

2.1.1 均匀半空间模型和层状模型

1) 模型网格剖分及观测频率

以下的均匀半空间模型、H型和K型层状模型均采用相同网格剖分，并选取相同的观测频率。沿x、y和z方向剖分为41×41×35(其中z方向地面以上10层为空气层)个网格。

观测频率范围为10 000～0.000 1 Hz，共取21个频点，分别为(单位：Hz)：10 000，8 000，5 000，2 000，1 000，500，200，100，50，10，5，2，1，0.5，0.1，0.05，0.01，0.005，0.001，0.000 5，0.000 1。

2) 均匀半空间模型

大地的电阻率设定为100 Ω·m，空气通常认为是绝缘体，但为了求解方程2，需要模型的电阻率为有限值，空气层的电阻率设为一个较大的常数，一般在106～1010Ω·m之间，本文涉及的三维模型其空气层电阻率均取为1010Ω·m，表1、图3为数值模拟结果。从表1可以看出，视电阻率误差最高为0.743 7%(频点为50 Hz)，其余频点计算误差均未超过0.6%，尤其在0.1 Hz以后误差均在0.1%以下。相位计算误差最高为0.233 56%(频点为100 Hz)，其余频点计算误差均在0.2%以下。

从图3给出的视电阻率、相位的对比曲线可以看出，在10 000～0.000 1 Hz频段范围，三维矢量有限元正演(vector finite element with PARDISO solver(VFE-PARDISO)曲线与解析解曲线重合。

表1 均匀半空间模型矢量有限元解与解析解对比

图3 均匀半空间模型三维正演视电阻率(a)和相位(b)与解析解对比Fig.3 Comparison of 3D forward apparent resistivity (a) and phase (b) of homogeneous half space model with analytical solution

3) H型层状模型

H型模型参数设置如表2所示。表3、图4为H型层状模型三维数值模拟结果，可以看出，矢量有限元数值解与解析解的最大误差为1.8%(观测频率50 Hz)，其他频点误差均在2%以下，表明模拟精度符合要求。

表2 H型层状模型参数

表3 H型层状模型矢量有限元解与解析解对比

图4 H型层状模型三维正演视电阻率(a)和相位(b)数值解与解析解对比Fig.4 Comparison of apparent resistivity (a) and phase (b) of 3D forward modeling of H-type layered model with analytical solutions

(4) K型层状模型

K型模型参数设置如表4所示。表5、图5为K型层状模型三维数值模拟结果，可以看出，矢量有限元数值解与解析解的最大误差为5.5%(观测频率50 Hz)，其他频点误差均在3.5%以下，表明模拟精度符合要求。

表4 K型层状模型参数

2.1.2 COMMEMI 3D-2模型

COMMEMI 3D-2 模型是国际电磁学术讨论会COMMEMI研究小组设计的一个三维地电模型[37]，该模型包含了多个电性差异巨大的分界面，比较复杂，已成为国内外学者普遍认可用于MT三维正演算法验证的三维模型[1,4,9,11-12,38]。图6为COMMEMI 3D-2模型示意，背景是一个三层K型地电断面，第一层电阻率10 Ω·m，厚度10 km；第二层电阻率100 Ω·m，厚度20 km；第三层电阻率为0.1 Ω·m。在剖面的第一层中，镶嵌有电阻率分别为10 Ω·m的低阻棱柱体和100 Ω·m的高阻棱柱体，对称地分布在x轴的两侧，棱柱体的规模为20 km×40 km×10 km，长轴为y方向，短轴为x方向。

表5 K型层状模型矢量有限元解与解析解对比

图5 K型层状模型三维正演视电阻率(a)和相位(b)数值解与解析解对比Fig.5 Comparison of 3D forward apparent resistivity (a) and phase (b) of K-type layered model with analytical solution

图6 COMMEMI3D-2 模型示意Fig.6 Schematic diagram of COMMEMI3D-2

将该模型沿x、y和z方向剖分为28×21×19(其中z方向地面以上8层为空气层)个网格，采用矢量有限元法对该模型进行三维正演模拟，并将矢量有限元模拟结果与Wannamaker等的积分方程法(integral equation，IE)三维模拟结果[1]进行比较。图7为观测频率f=0.001 Hz的模拟结果对比，图中实心黑色圆圈是本文的无需散度校正基于PARDISO直接求解的矢量有限元法(VFE-PARDISO)计算结果，方框是积分方程法(IE)计算结果，可以看出二种方法的计算结果基本一致。

2.2 起伏地形条件下模型验证

对带地形模型采用矩形六面体剖分，矩形网格x，y和z方向的间距可以是不均匀的，比如图8所示的三维山峰地形剖分。若在地形起伏界面(空气—地面)上加大网格剖分密度,使其剖分足够精细时,大地电磁场数值模拟响应将与理论场值相近,甚至趋于一致[39]。

为了验证正演程序对MT地形影响数值模拟效果，采用Wannamaker等设计的二维山峰模型[37]。模型背景电阻率为100 Ω·m，山峰地形如图9所示。设二维山峰地形的走向为x方向，山峰倾向为y方向，z方向垂直向下。为了尽量避免对y方向的三维影响，将山峰地形沿x方向延伸3.43 km，整个模型区域(34 300 m×34 300 m×91 993 m)沿x、y和z方向剖分为43×43×29(其中z方向山峰顶面以上7层为空气层)个网格单元。山峰地形采用9个纵向网格单元的划分，网格间距为50 m。

对上述模型(图9)进行三维正演模拟，计算频率f=2 Hz时TE极化模式和TM极化模式的视电阻率和相位曲线。图10中黑色虚线和实线分别是二维有限元(2D FE)法的TE、TM极化模式的计算结果，圆圈和方框是本文的矢量有限元法(VFE-PARDISO)的xy模式、yx模式的计算结果。从图中可以看出，本研究算法(VFE-PARDISO)与二维有限元的视电阻率、相位的计算结果一致，从而说明了本文研究的三维正演算法对起伏地形地电模型的计算结果准确可靠。

a—Zxy模式正演视电阻率；b—Zyx模式正演视电阻率；c—Zxy模式正演阻抗相位； d—Zyx模式正演阻抗相位a—Zxy mode forward apparent resistivity；b—Zyx mode forward apparent resistivity；c—Zxy mode forward impedance phase； d—Zyxmode forward impedance phase图7 本文矢量有限元正演算法的计算结果与IE方法的计算结果对比Fig.7 Comparison between the calculation results of vector finite element forward algorithm and IE method

a—三维地形示意； b—地形网格剖分和MT测点分布示意a—sketch of 3D topography； b—topography meshing and distribution of MT measurement sites图8 三维地形及其网格剖分示意Fig.8 3D topography and grid

图9 二维山峰地形示意Fig.9 Sketch of 2D ridge

a—正演视电阻率； b—正演阻抗相位a—forward apparent resistivity； b—forward impedance phase图10 三维矢量有限元算法(PARDISO)计算的二维地形影响与二维有限元结果对比Fig.10 Comparision between modeling results of 3DVFEM(PARDISO) and 2DFEM for 2D ridge

3 无需散度校正的直接解法与带散度校正的迭代解法计算对比

为了对比无需散度校正的PARDISO直接解法(VFE-PARDISO without divergence correction)和带散度校正的BICG迭代解法(VFE-BICG with divergence correction)的计算精度和计算时间，分别采用这两种求解方法对如图9所示的二维山峰地形模型进行三维正演模拟，并对比模拟结果。两种求解方法的三维正演模拟均在曙光W560-G20工作站上完成，计算及程序编译环境如下：CPU为Intel E5-2643 (3.4G)，内存为64GB，操作系统为Windows 7(64位)；编译环境为Microsoft Visual Studio 2012(已集成Intel Parallel Studio XE 2013)。整个模型区域沿x、y和z方向剖分为43×43×29(其中z方向山峰顶面以上7层为空气层)个网格单元。

两种算法的计算结果(观测频率为2 Hz)对比如图11所示，图11中黑色虚线和实线分别是二维有限元计算的TE和TM极化模式的计算结果(视电阻率和相位)，圆圈和方框是本文的无需散度校正直接求解的矢量有限元法(VFE-PARDISO)计算的xy模式(视电阻率、相位)和yx模式(视电阻率、相位)曲线图，菱形和三角形是带散度校正迭代求解的矢量有限元法(VFE-BICG)计算的XY模式(视电阻率、相位)和YX模式(视电阻率、相位)曲线图。从图11中可以看出，本文的算法(VFE-PARDISO)和迭代求解算法(VFE-BICG)计算的结果均与二维有限元的计算结果一致，但在平地与山峰拐点处，PARDISO直接求解的TM结果比BICG迭代求解的结果更接近于二维有限元计算结果。在计算时间方面，无需散度校正直接求解的矢量有限元法(VFE-PARDISO)正演一次上述模型仅耗时24 s，而带散度校正迭代求解的矢量有限元法(VFE-BICG)耗时407 s；本文的直接解法(VFE-PARDISO)与迭代解法(VFE-BICG)的计算速度比达17倍。可见本文的无需散度校正直接求解法与带散度校正的迭代求解法相比，在计算精度和计算时间方面均有优势，特别是在计算时间方面表现出明显优势，适合应用于中等计算规模的三维反演算法中。

a—正演视电阻率； b—正演阻抗相位a—forward apparent resistivity； b—forward impedance phase图11 无需散度校正的直接解法与带散度校正的迭代解法计算结果对比曲线Fig.11 Comparison of calculation results of VFE-PARDISO without divergence correction and VFE-BICG with divergence correction

4 结论

采用基于并行直接稀疏求解器PARDISO且无需散度校正的正演方案，并利用C++语言编制正演计算程序，实现MT三维快速正演。为了检验本文的快速正演算法及计算程序的正确性，通过数值解与解析解对比以及本文模拟结果与国际公认检验模型[36-37]的计算结果对比，结果表明该算法及程序在水平地形和起伏地形条件下均满足三维正演计算的精度要求。在中等规模计算条件下，通过本文的无需散度校正直接求解法与带散度校正的迭代求解法对比，本文的无需散度校正直接求解法在计算精度和计算时间方面均有优势，特别是在计算时间方面表现出明显优势，直接解法与迭代解法的计算速度比达17倍。本文实现的快速三维正演算法对于促进MT三维反演技术的实用性具有现实意义。