0 引言

大地电磁法以天然平面电磁波为场源，通过在地表观测相互正交的电磁场分量来获取地下的电性结构信息(Tikhonov, 1950; Cagniard, 1953;汤井田等2017, 2018)，是目前研究岩石圈及上地幔电性结构及构造特征的一种重要方法(魏文博等, 2006; Candansayar and Tezkan, 2010; Tang et al., 2013;强建科等, 2014; Li G et al., 2017; Hu et al., 2018; Mo et al., 2017; Li J et al., 2019).由于地球构造应力场、地球形变带、岩石裂隙、空隙水及地质沉积等因素的综合影响，造成地球电导率的分布呈各向异性现象.而基于各向同性假设的模型，无法合理解释含各向异性介质的大地电磁勘探数据(Li, 2002;Wannamaker, 2005).随着计算机技术和各向异性理论的发展，大地电磁各向异性数值模拟已成为研究的前沿(Li and Pek, 2008; Qin et al., 2013a, b ;霍光谱等, 2015;嵇艳鞠等，2016；邱长凯等，2018).

对于简单的二维各向异性地电模型，电磁响应存在解析计算表达式(Qin et al., 2013a, b )；而对于复杂的各向异性地电模型，须借助高精度的数值模拟方法来求解.目前，常用于各向异性电磁法数值模拟的方法有：积分方程法、有限差分法和有限单元法.积分方程法具有半解析解的精度，只需要对异常体进行网格剖分，具有未知数少的特点，但存在求解耗时和不适合处理带地形问题等缺陷(陈桂波等, 2009, 2010; Sun, 2015;任政勇等, 2017)；有限差分法具有简单的数学表达形式，且生成的是定带宽线性方程组，便于求解，但是该方法不适合模拟带地形的复杂模型(Saraf et al., 1986; Pek and Verner, 1997;霍光谱等, 2015)；有限单元法可以采用非结构化网格，适合模拟带任意地形的复杂模型(Reddy and Rankin, 1975; Li, 2002; Li and Pek, 2008; Ren et al., 2013; Fu and Gao, 2017).

有限体积法是基于有限差分和有限元两种方法发展而来的一种数值模拟方法，它兼顾有限差分法和有限元法的优点：1)线性方程组表达形式简单；2)能够采用非结构化网格剖分，适合处理带任意地形的复杂模型(Caughey and Jameson, 1981; Trompert and Hansen, 1996;谢德馨和杨仕友, 2009; Du et al., 2016).近几年来，该方法在地球物理勘探中的应用逐渐发展起来.众多学者采用基于正交网格的有限体积法，研究了电磁法的正演模拟问题(周建美等, 2014;彭荣华等, 2016)，但由于这种网格本身的缺陷，对于带地形的复杂模型难以进行精确的模拟.而基于非结构化网格的有限体积法，虽然能够很好的处理带任意地形的复杂模型(Jahandari and Farquharson, 2014, 2015; Du et al., 2016)，但是对于带任意地形的复杂各向异性模型的电磁响应，利用有限体积法仍未得到有效的模拟.

为了解决复杂地形情况下二维各向异性大地电磁正演问题，本文开发了基于非结构化三角形网格的有限体积正演算法.并采用本文算法、解析算法、有限元算法分别对层状各向异性模型和复杂各向异性模型进行计算，结果表明本文算法能够稳定、可靠地求解带任意地形的复杂电导率各向异性模型的大地电磁响应.

1 正演理论 1.1 边值问题

大地电磁场满足如下麦克斯韦方程(取时间因子e^－iωt)：

(1)

(2)

式中E和H分别为电场和磁场，ε为介电常数，ω=2πf为角频率，f为频率，为虚数单位，μ为磁导率，σ为电导率张量.

如图 1所示，假设电性主轴i与直角坐标系x方向(构造走向)一致，考虑电性主轴k, m分别与y, z坐标轴成一角度α(徐世浙和赵生凯, 1985; Li, 2002).在直角坐标系中，电导率张量σ可表示为：

(3)

式中A为坐标变换张量：

σ′为对角矩阵：

其中，σ_i、σ_k和σ_m分别为电性主轴i、k和m方向上的电导率.

对于二维介质，电磁场沿构造走向x方向保持不变(即)，Maxwell方程可解耦成两个独立的极化方程(Li, 2002)：

TE极化：

(4)

TM极化：

(5)

方程组(4)和方程组(5)可以写成统一的表达形式：

(6)

式中，∇为梯度算子.对于TE极化模式，u=E_x，.对于TM极化模式，u=H_x，其中：

1.2 节点中心有限体积法

有限体积法存在两种方式：单元中心方式和节点中心方式(Caughey and Jameson, 1981;谢德馨和杨仕友, 2009).单元中心方式将单一的网格单元作为控制体单元，并假设单元内电磁场为常数，该方式存在难以求解辅助电磁场的缺陷(Trompert and Hansen, 1996; Du et al., 2016).而节点中心方式以网格节点为中心形成控制体单元，并将未知场置于节点处，借助于线性形函数，可以容易地计算辅助电磁场.因此，本文采用节点中心方式的有限体积法进行模拟.

为了处理任意地形和复杂电导率结构模型，本文采用开源剖分程序Triangle(Shewchuk, 1996)对求解区域Ω进行三角单元网格剖分.在Ω上选取任意节点v_i(如图 2所示)，依次相连以v_i为顶点的三角单元的形心和三角单元各边的中点，形成一个封闭多边形c_i，称为包含节点v_i的控制体单元.

方程(6)在控制体单元c_i上的积分可以表示为：

(7)

利用散度定理(Jin, 2002)，方程(7)可化为如下积分：

(8)

其中，为多边形c_i的边，n是c_i的单位外法向量.假设求解区域Ω共有N_n个节点，那么控制体单元的个数必是N_n，方程(8)的个数也为N_n，从而形成定解的有限体积问题.方程(8)和方程(6)的场u通过各节点的未知场u_k和线性形函数S_k展开：

(9)

将(9)式代入方程(8)可得：

(10)

(10) 式写成矩阵的形式是：

(11)

其中，U为N_n个节点上的未知场向量，A为大型稀疏系数矩阵，其系数表达式为：

(12)

其中，i, k=1, 2, …, N_n.

事实上，控制体单元c_i是由一系列多边形P_t(t=1, 2, …, n_p，n_p是以节点v_i为顶点的三角形的总数)组成.因此，A矩阵的系数又可写为：

(13)

其中，是多边形P_t与控制体单元c_i所重合的边，A_ij在多边形P_t上的求解过程见附录A.

求解方程(11)中的U需要添加适当的边界条件.当电导率各向异性异常体距离截断边界足够远(5~10倍趋肤深度)时，二次场在截断边界逐渐衰减为零.因此，对于方程(6)和方程(11)的总场公式，将截断边界扩展到距离地下异常体5~10倍趋肤深度处，取1D各向异性模型的解析解作为边界条件.结合方程(11)，从而得到最后的线性方程组.利用并行直接求解器Pardiso对线性方程组进行快速求解，得到未知场向量U(E_x, H_x).通过方程(4)和(5)进一步求得电磁场的其他分量：

(14)

(15)

从而计算出视电阻率和相位值：

(16)

(17)

其中

2 正演分析 2.1 层状各向异性模型

首先，采用具有解析解的层状各向异性模型来验证本文算法的正确性.模型参数如图 3所示，其中第一层和背景层均为各向同性层，电阻率分别为ρ₁=10 Ωm和ρ₃=1000 Ωm，第二层为ρ_2i=200 Ωm，ρ_2k=2 Ωm，ρ_2m=40 Ωm，α₂=60°，h₂=4 km的各向异性层，顶部埋深1 km.采用三角形单元进行网格剖分(如图 4)，单元个数为550338，节点个数为275408，22个计算频点按对数均匀分布在10^－5~10²Hz.如图 5所示，将本文计算结果与解析解(Li, 2002)进行对比，可以看出两者的计算结果十分吻合，视电阻率相对误差小于1.3%，相位残差小于0.3°，证明了本文算法的正确性.

2.2 各向异性矩形棱柱体模型

为了进一步验证有限体积法求解大地电磁各向异性问题的能力，选取较为复杂的二维各向异性异常体模型进行计算，将有限体积法的结果与已发表的有限单元法的结果(Li, 2002; Ren, 2014)进行对比.模型参数如图 6所示，背景是ρ₁=1000 Ωm的各向同性均匀半空间，各向异性矩形棱柱体电阻率为ρ_2i=500 Ωm，ρ_2k=10 Ωm，ρ_2m=500 Ωm(其中电性主轴i与走向x轴平行)，α₂=60°.柱体截面宽为2 km，厚为8 km，顶部埋深1 km，沿x轴方向无限延伸.模型计算区域为y∈[－500 km, 500 km]，z∈[－500 km, 500 km]，计算频率是0.1 Hz，41个测点均匀分布在地表y∈[－4 km, 4 km]上.

采用相同的网格剖分单元(如图 7所示)，图 8对比了本文计算结果、Li(2002)的有限元法结果与Ren(2014)的有限元法结果.从图中可以看出，三者的计算结果十分吻合，本文计算的视电阻率与Li的相对误差小于1.5%，与Ren的相对误差小于0.01%.图 9对比了Y_A=0 km、Y_B=0.8 km、Y_C=-0.8 km处三种算法在整个计算频域中的结果，从图中可以看出，有限体积法的计算结果与Li的相对误差均小于1.4%，与Ren的相对误差小于0.38%.该测试结果进一步验证了本文算法的稳定性及准确性.

表 1为有限体积法和有限单元法(Ren, 2014)的计算参数.从表中可以看出，在相同的计算平台和网格条件下，对于二维大地电磁各向异性模型的正演，有限体积法与有限单元法的消耗相当.

2.3 复杂背斜模型

选用背斜山谷模型进一步测试本文算法处理复杂模型的能力.背斜顶部山谷宽为4 km，谷底深为0.5 km，背景电阻率为ρ₁=100 Ωm的各向同性结构，各向异性背斜的电阻率为ρ_2i=10 Ωm，ρ_2k=100 Ωm，ρ_2m=500 Ωm，夹角α₂=30°，如图 10所示.计算区域为y∈[－1000 km, 1000 km]，z∈[－1000 km, 1000 km].19个计算频率按对数均匀分布于10^－3~10³ Hz，测点位于地表.

采用相同的网格剖分单元(如图 11所示)，分别使用有限体积法和有限元法(Ren, 2014)计算了该模型的大地电磁响应，图 12为有限体积法计算的结果，图 13是两种方法计算结果的相对误差.从图中可以看出，二者的计算结果非常吻合.TE模型的视电阻率相对误差小于0.95%，相位残差小于0.22°；TM模型的视电阻率相对误差小于0.52%，相位残差小于0.075°.表明本文开发的有限体积算法，能够高精度的模拟带地形的复杂各向异性模型.

由图 12可知，TE模式视电阻率和相位拟断面图清楚地显示了背斜的位置和形态.但TM模式中，由于各向异性的影响，视电阻率和相位在背斜核部发生严重畸变，使得TM模式的视电阻率和相位拟断面图无法显示背斜的形态.

3 结论

本文开发了一种新的二维大地电磁正演模拟算法.从麦克斯韦方程出发，推导了大地电磁场各向异性介质的边值问题，在控制体单元上，基于散度定理，采用线性插值函数，建立了节点型有限体积法的计算流程，得到电磁场的响应值.通过三种不同复杂度的模型算例对比，得到以下结论：

(1) 本文提出的基于非结构化网格的大地电磁各向异性有限体积算法，是一种能够适合模拟带任意复杂地形模型的高精度、稳定的算法.

(2) 有限体积法基本思路易于理解，基于守恒原理，利用散度定理可直接将复杂的赫姆霍斯方程转换成线性方程组进行求解，原理简单，易于编程实现.

(3) 对于同一模型，在相同的计算平台和网格条件下，有限体积法与广泛应用的有限单元法的计算精度相似，且消耗相当.因此，有限体积法是处理电磁法各向异性问题的另一种有效方法.

致谢

感谢Shewchuk提供的非结构三角形网格剖分开源代码Triangle，感谢Schenk提供的方程组求解器Pardiso，感谢网格剖分可视化软件PARAVIEW的作者.

附录A

首先考虑式(13)中的第一部分线积分.对于TE模式，积分可以表示为：

(A1)

其中，为多边形p_t(i)的边(如附图A1)，n为的外法向量，n₂和n₃分别表示边L₂₃和L₃₄的外法向量，l₂₃和l₃₄分别表示边L₂₃和L₃₄的长度，e_y和e_z分别表示y和z方向的单位向量.

令，有：

(A2)

其中：

是三角形单元的面积.

对于TM模式，，式(13)中的第一部分线积分可以表示为：

(A3)

其具体表达式为：

(A4)

然后，计算式(13)中的第二部分在多边形p_t(i)上的面积分.建立二维自然坐标系L_i, L_j, L_m进行坐标转换(如附图A2)：

在多边形p_t(i)上1点处：L_j=0, L_m=0；2点处：L_j=1/2, L_m=0；3点处：L_j=1/3, L_m=1/3；4点处：L_j=0, L_m=1/2，所以yz平面上的多边形p_t(i)变换成L_j, L_m平面上的四边形.

那么，积分可以表示为：

(A5)

令，有：

(A6)