电阻率各向异性介质大地电磁二维非结构有限元数值模拟

0 引言

介质的电性各向异性现象普遍存在[1]，在地电模型和地球动力学解释中发挥着重要作用。早在20世纪70年代，文献[2]和文献[3]就已经开始了对各向异性大地电磁(MT)的研究。在各向异性MT数值模拟领域， Pek等[4]公布了二维各向异性MT有限差分公式；Li[5]实现了二维各向异性介质的MT有限元数值模拟，其采用线性插值基函数，数值精度在电磁场变化剧烈的地区受限，且其数值模拟结果没有得到各向异性模型解析解的验证；QIN L[6]给出了二维各向异性断层模型的MT解析解，为各向异性二维MT的数值模拟提供了比对结果。依托MT有限差分正演模拟，胡祥云和霍光谱等 [7-8]在实际大地电磁数据解释中引入了各向异性。

在MT的应用方面，20世纪80年代，Cox[9]和Ranganayaki等 [10]应用MT展开对海岸效应的研究，并给出了该模型同性介质TM模式观测的特征频率和特征距离的经验公式；Constable等 [11]在的TE模式观测中也发现了海岸效应特征异常。

但是这些研究均基于各向同性介质模型，与各向异性普遍存在的地球模型有一定出入，因此有必要研究各向异性介质下的海岸效应现象。笔者从麦克斯韦方程组出发，推导了任意各向异性介质二维大地电磁的非结构二阶有限元公式，实现了二维大地电磁有限元正演模拟。相比较一阶插值，二阶插值有限元模拟结果更精确，而非结构有限元的优势在于，可以模拟任意复杂地形的MT观测结果。

1 有限元理论分析

任意各向异性介质下的麦克斯韦方程组如式(1)所示。

(1)

其中:是电导率张量，为一个3×3的矩阵：

电导率参数 与电阻率参数 互逆，即将麦克斯韦方程组(1)的旋度展开，并运用矢量公式变换进行转化，得到了关于Ex和Hx的方程组(2)和方程组(3)，以及辅助场Ex和Hx关于Ex和Hx的表达式(4)。

(2)

▽

(3)

(4)

其中：

C=σ11+σ12B+σ13A

D=σ22σ33-σ23σ32

求解公式(2)和公式(3)，可以得到水平电磁场量Ex和Hx。

1.1 插值函数

图1是二维非结构网格剖分图，采用Gmsh软件进行剖分。图2是非结构网格中的某个三角形单元，编号为该单元内插值节点的内部编号。图中的坐标系y轴为走向方向，x轴垂直于走向方向，z轴正方向垂直向下。我们在三角形单元内进行二次差值，即假定电场Ex和磁场Hx是y和z的二次函数，可以近似为：

i=1,…,6

(6)

其中:Ei和Hi分别是全局坐标系yoz中三角形单元的第i个节点的电场和磁场;二次差值的插值函数Ni采用文献 [12]中的定义，Ni公式为式(7)。

Ni(y,z)=2(Li-1)Li i=1,2,3

N4(y,z)=4L1L2

N5(y,z)=4L2L3

N6(y,z)=4L3L1

(7)

Li公式为式(8)。

c1=y2z3-z2y3;a1=z2-z3;b1=y3-y2

c2=y3z1-z3y1;a2=z3-z1;b2=y1-y3

(8)

1.2 边界条件

笔者利用有限元方法求解的是一个关于电磁场的边值问题，内部边界的电场切向分量Et和磁场切向分量Ht连续。外部边界条件选用Dirichlet边界条件。

左右两侧边界加载时，将模型的边界设置离异常体足够远，大概异常体规模10倍左右的距离，采用一维程序计算边界场值并加载。上边界(下边界)场值的加载，则将左右两侧的上边界(下边界)值读取出来，以y轴的距离为插值变量进行线性插值，从而计算出上边界(下边界)的场值。

1.3 单元分析

将公式(2)乘以电场的任意变分δ Ex，公式(3)乘以电场的任意变分δ Hx，并在模拟区域Ω上求积分。求解区域Ω由ne个三角形单元组成，单元编号记为e=1、2、…、ne，电阻率参数赋值到每一个三角形单元。利用矢量计算公式化简可得到如下积分公式:

▽Ex▽

(9)

▽

▽

(10)

公式(9)和公式(10)中的Γe为单元e的边界，ey和ez是方向单位向量。在内边界Γe上，电场的切向分量Et和磁场的切向分量Ht连续。在求积分的过程中，沿着每一个边界都进行了两次积分，积分的方向相反，因而，所有的内部边界积分之和为零。又因为在外边界上我们采用Dirichlet边界条件，故电场的变分δ Ex和磁场的变分δ Hx为零，故所有线性积分之和等于零，公式(9)和公式(10)可以进一步化简为式(11)与式(12)。

▽Ex▽

(11)

▽

▽δHxdΩ=0

(12)

图1 二维非结构单元
Fig.1 2D unstructured element

图2 三角形单元
Fig.2 Triangle element

1.4 刚度矩阵

将插值基函数公式(6)代入式(11)和式(12)，并将单元矢量Exe、Hxe、δ Exe和δ Hxe移到积分外面(因为它们不再依赖于y和z)，可得到式(13)和式(14)。

(13)

(14)

式中：在进行单元矩阵求和之前，先将各个单元的矢量Exe、Hxe、δ Exe和δ Hxe扩展成nd×1的系统矢量Ex、Hx、δ Ex和δ Hx，把6×6的单元矩阵K1e、K2e、K3e、K4e、K5e和K6e扩展成nd阶矩阵和为节点总数。然后，将所有单元上的扩展矩阵相加，考虑到δ Ex和δ Hx的任意性，则由式(13)和式(14)，可得线性方程组式(15)。

K U=0

(15)

其中，

由于单元矩阵K1e、K2e、K4e和K5e是对称矩阵，故矩阵K11和K22也是对称阵。由得所以系统矩阵K为对称矩阵。又因为K11、K22、K21和K22是nd阶的具有带状结构且含有大量零元素的稀疏矩阵，故系统矩阵K为2nd阶的复数的、对称的稀疏矩阵。图3中剖分结构由八个三角形单元23个节点构成，其对应的稀疏矩阵如图4所示。

本文中调用Fortran中的pardiso库函数，运用直接法求解系数矩阵。

图3 网格剖分实例
Fig.3 Grid example

图4 矩阵稀疏
Fig.4 Sparse matrix

1.5 视电阻率的计算

相对于各向同性介质而言，各向异性介质中的电场和磁场一般不能解耦，则阻抗张量Zxx和Zyy不为零，所以需要在两种不同的极化方式下求解这4个阻抗，其公式如下：

(16)

E和H下标中的1和2代表两种不同的极化方式。公式(16)中的辅助场Ey和Hy依据公式求得。根据公式(16)求得的阻抗值，就可以进一步求解视电阻率和相位，二者的计算公式为式(17)。

(17)

2 算法验证

为了验证本文任意各向异性介质二维大地电磁非结构有限元正演模拟算法(MT2DFEANI)的正确性，模拟了一维和二维测试模型的大地电磁响应，并依据公式(18)，计算模拟数据Dc与对比数据Dr的相对误差r。

(18)

2.1 一维模型验证

一维层状各向异性模型的大地电磁响应，存在解析解，因此设了一个三层的电阻率参数为“高低高”的模型和一个三层的电阻率参数为“低高低”的模型。两个模型的第一层和第二层深度均为3 km，其中“高低高”模型(D模型)的第一层和第三层的电阻率为1 000 Ω·m，第二层为电阻率各向异性层，其主轴电阻率一次为参数30 Ω·m、100 Ω·m、30 Ω·m，水平旋转角度αs为30°。低高低模型(H模型)的第一层和第三层的电阻率为100 Ω·m，第二层为电阻率各向异性层，其主轴电阻率一次为参数300 Ω·m、1 000 Ω·m、300 Ω·m，水平旋转角度αs为30°。

对比图5和图 6，可以得到三点结论：①模拟结果对比中，高频段的模拟结果出现振荡现象，而在低频段模拟结果很好，最大相对误差在2%以内；②在周期T=10 s之后，本文算法的模拟结果相当准确，相对误差接近零；③相较一次插值有限元模拟结果，二次插值有限元模拟的结果能够在低频段出现较小的震荡现象，且相对较快的收敛到解析解结果，证明了本文二次插值程序的优势性。

2.2 二维模型验证

图7是d'Erceville等[13]提出的断层模型，介质1的主轴电阻率依次为20 Ω·m、20 Ω·m、100 Ω·m，介质2的主轴电阻率依次为1 000 Ω·m、1 000 Ω·m、20 Ω·m。该断层模型，TE模式观测结果只与主轴电阻率σ11有关，即该各向异性断层模型的TE模式观测结果与以σ11为电导率的各向同性介质断层模型观测结果相同，应用MT2DFEANI计算了该模型在周期T=10 s时的ρyx和φyx。断层模型的边界设置为180 km，上下边界顶点的剖分尺寸为80 km，测线两侧顶点处和异常体剖分尺寸为500 m，网格剖分结果如图8所示。本文模拟结果与文献[6]的解析解结果对比如图9所示。

图9的结果可以看出，MT2DFEANI的结果与解析解的结果吻合一致。经计算，解析解与MT2DFEANI模拟的ρyx和φyx的最大相对误差分别为2.12%、2.07%，这验证了本文MT2DFEANI模拟的准确性。

为了进一步检测本文有限元程序的正确性，笔者参考Pe等[4]设计的一个较为复杂的模型。模型如图10所示，具有电阻率水平各向异性的水平地层位于一个二维异常块体下方，该二维体出露至地表，且其电阻率亦是水平各向异性的，主轴电阻率依次为30 Ω·m、100 Ω·m、30 Ω·m，其各向异性主轴x′与走向方向x轴之间的夹角为30°，而水平地层的各向异性主轴x′与走向方向x轴之间的夹角为120°，主轴

图5 高低高模型模拟结果
Fig.5 Results for 1D D-model

图6 低高低模型模拟结果
Fig.6 Results for 1D H-model

电阻率依次为10 Ω·m、100 Ω·m、10 Ω·m，即二维异常体的各向异性水平主轴与水平地层的各向异性主轴相互垂直。该模型的外边界设置为230 km，上下边界顶点的剖分尺寸为80 km，测线两侧顶点剖分尺寸为2 km，异常体剖分尺寸为200 m，网格剖分结果如图11所示，异常体剖分放大图如图图12所示。MT2DFEANI模拟的结果和文献[4]的有限差分模拟结果对比如图13所示，模拟周期为T=30 s。从图13的对比结果可以看出，在复杂各向异性的介质中，本文的算法模拟结果和有限差分的模拟结果吻合一致。视电阻率ρxx、ρxy、ρyx、ρyy和φxx、φxy、φyx、φyy最大相对误差均控制在2.5%以内，验证了MT2DFEANI模拟结果的正确性。

图7 各向异性半空间内无限延伸的断层模型
Fig.7 Anisotropic semi-infinite fault model

图8 断层模型的非结构网格剖分
Fig.8 Grid for the fault model

图9 断层模型的解析解和本文有限元结果的对比图
Fig.9 Comparison between analytic results and the
MT2DFEANI results for the fault model

图10 Pek and Verner设计的二维各向异性模型
Fig.10 Two dimensional anisotropic model
designed by Pek and Verner

图11 Pek and Verner设计模型的网格剖分图
Fig.11 Grid for model by Pek and Vernerr

图12 二维各向民性模型
Fig.12 Two dimensional anisotropic model

图13 模型的有限差分(FD)和本文
有限元(FE)模拟结果
Fig.13 Results comparison between finite-
difference and finite
element algorithm

3 各向异性影响研究

模型为含有一个二维各向异性异常体的半无限空间介质，异常体埋深2 km，半无限空间介质的电阻率为1 000 Ω·m。

3.1 倾斜各向异性

倾斜各向异性介质是y,z轴在yoz平面内绕x轴旋转αd角度得到，其电阻率张量与主轴电阻率张量和αd的关系式表达如下：

(19)

将表达式(19)代入式(2)和式(3)，电场和磁场耦合在一起的方程组可以化简为：

▽2Ex+iωμ0σ11Ex=0

(20)

▽▽Hx)+iωμ0Hx=0

(21)

由式(20)可以看出，在倾斜各向异性情况下，TE模式已经解耦，且其求解方式与电阻率为σ11的各向同性介质求解方式完全相同。式(21)表达的TM模式也已经解耦，但是其仍然较为复杂，不能用各向同性介质求解方式求解。倾斜各向异性的模型如图14所示，各向异性异常体的主轴电阻率ρx'、ρy'、ρz'，依次为300 Ω·m、10 Ω·m、300 Ω·m。为讨论αd对大地电磁响应的影响程度，我们设置αd依次为0°、30°、45°、60°、90°，该模型在T=10 s下的大地电磁响应如图15所示。从图15我们可以得到几个结论：①ρxy和φxy不受αd的影响，仅与σ11有关，这与偏微分方程的表达一致；②ρyx和φyx会随着αd变化，且αd越大，ρyx和φyx的变化也越大，其曲线也不再是关于中心对称；③当αd=90°时，介质转变成垂直各向异性介质，ρyx和φyx与垂直各向异性介质响应一致，曲线又恢复了中心对称形式。

图14 倾斜各向异性模型
Fig.14 Dipping anisotropic model

3.2 主轴各向异性

主轴各向异性的电阻率张量中，只有对角线存在非零元：

将代入和中，可以得到：

▽2Ex+iωμ0σ11Ex=0

(22)

图15 倾斜各向异性模型的MT响应
Fig.15 Modeling results for dipping anisotropic model

图16 TIH和TIV介质的视电阻率ρyx和相位φyx
Fig.16 Apparent and phase results for the TIH and TIV media

(23)

由式(22)可以看出，TE模式解耦，可以用各向同性的求解方式求解电场。由式(23)可以看出，TM模式也已经解耦，但是由于σ22≠σ33，故此时的TM模式无法按照各向同性的MT响应求解。垂直各向异性(TIV)和方位各向异性(TIH)是其中两种特殊的情况，TIV介质下介质下根据公式(21)和倾斜各向异性的讨论，得到主轴各向异性下的ρxy和φxy仅与σ11有关。图16是TIH介质和TIV介质的ρyx和φyx，模型中ρ1为300 Ω·m，ρ2依次取为10 Ω·m、100 Ω·m、300 Ω·m、500 Ω·m、1 000 Ω·m。

对于TIH介质，随着ρ2的增大，ρyx和φyx逐渐趋于平缓，直到当ρ2=1 000 Ω·m时，此时ρyx和φyx变成一条直线，即TIH介质中当ρ2与围岩介质电阻率相同时，此时无法检测到异常电阻率体。TIV介质的ρyx和φyx也随着ρ2的增大趋于平缓，但是相对于TIH介质的曲线变化微乎其微。根据以上分析得出以下几个结论：①在主轴各向异性情况下，相较于σ33而言，ρyx和φyx受σ22的影响更大；②ρyx和φyx不受σ11的影响，这一点由式(23)可知；③ρxy和φxy只与σ11有关。

3.3 水平各向异性

水平各向异性是x轴和y轴在xoy平面内绕z轴旋转一定角度αs，其电阻率张量与主轴电阻率张量和αs的关系式表达如下：

将其电阻率张量的表达式代入式(2)和式(3)，方程组可以化简为：

(24)

(25)

由式(23)和式(24)可以看到Ex和Hx依然耦合在一起，此时的电场Ex不仅和σ11、σ12、σ22有关，且阻抗张量中的Zxx和Zyy为非零元。水平各向异性的模型如图17所示，各向异性异常体的主轴电阻率依次为300 Ω·m、10 Ω·m、300 Ω·m，为讨论αs对大地电磁响应的影响程度，设置αs依次为0°、30°、45°、60°、90°。该模型在T=10 s下的大地电磁响应如图18所示，图19是其阻抗张量图。从图18可以看出，随着αs的变化，ρxy、φxy、ρyx和φyx均发生变化。从图19的阻抗张量图中可以得出三个结论，一是随着观测点逐渐远离异常体，阻抗Zxx和Zyy趋于零，在异常体中心位置出现最大值，二是在αs从0°增长到90°的过程中，Zxx和Zyy先增大后减小，且在αs=0°和αs=90°是均为零，三是αs从0°增长到90°的过程中，Zxy和Zyx在异常体附近减小，而在αs=90°时，Zxy和Zyx曲线出现较大变化，因为此时σ11和σ22发生互换。

图17 水平各向异性模型
Fig.17 Horizontal anisotropic model

图18 水平各向异性模型的大地电磁响应
Fig.18 Modeling results for horizontal anisotropic model

4 结论

通过本文算法的推导、验证及各向异性参数的分析研究，可以得到几点结论。

1)通过两个一维及两个二维各向异性模型算例，验证了本文的二阶插值非结构有限元算法的正确性。通过一维模型算例中，与一次插值非结构有限元结果的对比，验证了本文二次插值非结构有限元算法的优势性。

图19 水平各向异性模型的阻抗张量图
Fig.19 Impedance results for horizontal anisotropic model

2)对于较为特殊的倾斜各向异性介质和主轴各向异性介质而言，其电场和磁场已经完全解耦，其阻抗分量中的Zxx和Zyy仍然为零，这与各向同性介质相同，但是其Zxy≠Zyx。其TE模式只依赖于σ11。而其TM模式虽已经解耦，但是由于各向异性电阻率张量的存在，其磁场的传播与多个电阻率分量有关，使得其TM模式求解不再像各向同性介质一样简单，也是其视电阻率和相位有着复杂的变化。

电阻率各向异性现象非常复杂，但却符合真实的地质模型，因此在后续工作中，将进一步讨论电阻率各向异性对大地电磁观测的影响。