三维大地电磁响应各向异性特征研究

王 磊, 岳明鑫,,3, 杨晓冬,, 李 勇, 吴小平

(1.中国科学技术大学 地球和空间科学学院，合肥 230026；2.中国地质科学院 地球物理地球化学勘查研究所,自然资源部地球物理电磁法探测技术重点实验,廊坊 065000；3.安徽国科骄辉科技有限公司,合肥 230040)

0 引言

地球介质的各向异性现象是普遍存在的[1-4]，微观上，电导率各向异性与晶体的优势取向以及地球动力学过程密切相关,能提供有关岩体内部介质裂隙,矿物晶体排列和应力场、形变带特征等信息[5-6]；宏观上，成矿作用时地下岩浆活动会引起地壳应力变化，进而影响含流体微裂隙以及孔隙的定向排列从而引发强烈的电各向异性特征[7]。因此,研究地下结构的电导率各向异性特征，有助于揭示岩石圈壳幔物质的分布特征以及变形机制，对于理解构造运动和岩浆活动等球动力学过程具有重要的科学意义。

大地电磁测深法(MT)以天然交变电磁场为场源，工作效率高、探测深度大，可以获取地下介质的电性参数,是深部地球物理探测的一种主要方法技术，近年来被广泛应用于寻找地热资源[8]、深部地壳结构研究[9]以及上地幔结构等研究中[10]。现阶段大地电磁数据的处理主要基于各向同性理论[11-13]，各向同性假设会对大地电磁实测数据的处理与解释带来较大的偏差，无法反应真实的地质结构，甚至会导致错误的结果。现阶段野外观测数据中的某些特征，例如相位超象限现象已经被认为对电各向异性结构起到一定的指示意义，然而更多的因电各向异性产生的数据特征，需要电各向异性理论模型数值实验的验证，因此迫切需要研究准确高效的电各向异性三维正演模拟算法，模拟各向异性地质体对观测大地电磁场的影响。

数值模拟是反演方法研究和数据解释的基础，关于电各向异性介质的数值模拟研究最早可以追溯到上世纪六十年代[14]。Pek等[15]实现了各向异性介质中二维MT有限差分法，解释了相位超象限现象； Li[16]实现了各向异性介质中二维MT有限元法，能更好的模拟复杂地下介质条件。在三维MT各向异性的数值模拟方面，目前主要包括积分方程法[17]、有限差分法[18]、有限体积法和有限元法[19]。其中基于有限元方法的MT正演能更好模拟地下不规则异常体，矢量有限元法因其完美模拟单元边界上电磁场的不连续性，并严格满足散度为零的条件，故在MT数值模拟领域得到越来越广泛的应用。近几年，曹晓月等[20]、Liu 等[21]先后提出基于总场的各向异性介质中三维MT非结构矢量有限元模拟方法；秦策等[22]、李志旋等[23]先后提出基于二次场的各向同性介质三维正演模拟方法，而基于二次场任意各向异性介质三维非结构矢量有限元大地电磁正演模拟尚少见报道。

笔者从麦克斯韦方程出发，推导了任意各向异性介质中三维大地电磁的二次场变分方程，使用非结构矢量有限元法进行求解，对于大型稀疏线性方程组的求解，选用直接求解器PARDISO。使用本算法计算COMMEMI模型[24]和二维各向异性模型，与前人已发表结果进行对比，验证了算法的准确性，设计板状体模型计算倾斜各向异性条件下的MT响应，为突显非结构网格优势，设计球型异常体模型并计算其在主轴各向异性条件下的MT响应，分析异常体主轴电导率对视电阻率的影响，详细讨论了水平各向异性条件下，主轴电导率、欧拉角对视电阻率的影响。

1 二次场的三维MT正演算法

假设时谐因子为e-iωt，频率域的麦克斯韦方程为：

(1)

(2)

电导率各向异性可以由各向异性参数σx、σy、σz，即αs、αd、αl三个主轴分量和三个欧拉角来表示。以电场为例，式(1)可以写成双旋度方程如下：

(3)

(4)

这里背景模型采用一维层状介质模型。联合式(3)、式(4)，不考虑介电常数的变化，可得二次场方程：

(5)

(6)

模型空间可以离散为四面体单元，如图1所示。四面体单元内任一点的电场可以近似表达为式(7)。

图1 四面体单元及棱边编号图Fig.1 Diagram for tetrahedron element and its edge numbers

(7)

其中：ei表示第i条边上的电场值；Ni为型矢量形函数[26]表达式为式(8)。

Ni=(Li1Li2-Li2Li1)li

(8)

将式(7)、式(8)带入式(6)，消去变分项可得：

(9)

其中：

(10)

计算四种单元积分，将结果带入式(9)。一次场值Ep由解析公式计算，得到每个节点的一次场值，使用投影法将节点上一次场值转化到每条棱边上。将单元刚度矩阵合成为总体刚度矩阵，表达式为式(11)。

KEs=S

(11)

式中:K为大型对称稀疏矩阵，采用CSR格式进行存储以节约内存空间,为获得快速稳定的计算结果，选用PARDISO直接求解器进行方程组求解。

2 MT响应计算

由式(11)可以得到计算域内每条棱边上的二次场值，带入式(1)可得二次电场值。进而通过式(7)可以得到任意测点的二次场值，将其与此点的一次场相加，即得总场值。另外，各向异性介质中的电磁场一般不能解耦，故阻抗张量Zxx与Zyy一般不为零，这里在两个正交源下分别模拟得到两组数据，求解阻抗，公式如下：

(12)

式中：下角标1和2分别表示x、y方向的源。最后可以通过式(13)计算视电阻率与相位：

(13)

图2 一维层状模型示意图Fig.2 Diagram of one-dimensional layered model

图3 一维层状解析解与模拟结果对比Fig.3 1D layered analytical solution compared with simulation

3 算法验证

3.1 一维层状解析解

使用一维层状解析解来验证算法的准确性。一维层状模型示意图如图2所示，对比结果如图3所示，可以看到对于不同频点，本算法模拟结果与解析结果相吻合，验证了算法的准确性。

3.2 COMMEMI-3D1模型

一维层状介质和简单二维模型具有解析解，一般的三维模型很难得到其解析解，因此对于三维模型，当不同的数值模拟方法得到不同的结果时，需要一个客观标准COMMEMI(Comparison of Modeling Methods for Electroma -gnetic Induction)是一个众多学者参与的国际合作项目[24]，各个机构采用不同的数值模拟方法计算标准模型的MT响应，为新的算法提供一个对比标准。

图4 COMMEMI-3D1模型Fig.4 COMMEMI-3D1 model(a)YOX剖面;(b)XOZ剖面

图5 模拟结果与前人结果对比Fig.5 Comparisons of modelling results and previous results(a)x测线结果；(b)x测线结果;(c)y测线结果；(d)y测线结果

笔者采用COMMEMI-3D1模型验证三维正演算法的正确性。模型如图4所示，背景电阻率为100 Ω·m，异常体电阻率为0.5 Ω·m，测线布置为过原点的沿x、y轴两条测线，模型生成186 627个非结构四面体单元，设置频点为0.1 Hz，图5显示本文算法计算结果与前人结果相吻合，曲线平滑且落在误差区间内，验证了本算法对于三维模型模拟的正确性。

3.3 水平各向异性模型

为了验证算法对各向异性介质模拟计算的准确性，参考Pek等[15]设计了一个较为复杂的二维各向异性模型(图6)，各向异性块体出露地表，下方与各向异性水平层相接，上下两个异常体的各向异性水平主轴相互垂直异常体1的参数为：ρ1=30,ρ2=100,ρ3=30,αs=30°；异常体2的参数为：ρ1=10,ρ2=100,ρ3=10,αs=120°。本程序模拟的视电阻率模拟对比如图7所示。可以看到模拟结果与有限差分结果吻合，验证了算法的正确性。

图6 二维水平各向异性模型Fig.6 Two-dimensional horizontal anisotropy model(a)电性介质分布图;(b)水平各向异性示意图

图7 水平各向异性模型视电阻率的有限差分结果与本文矢量有限元结果Fig.7 Comparison of apparent resistivity between finite difference method and the finite element method proposed by this paper(a)ρxx;(b)ρxy;(c)ρyx;(d)ρyy

图8 三维倾斜各向异性模型Fig.8 Three-dimensional dipping anisotropy mode(a)XOY剖面;(b)倾斜各向异性示意图

4 算例

4.1 倾斜各向异性

倾斜各向异性介质是y、z轴在yoz平面中绕x轴旋转αd得到，电导率张量与主轴电导率张量以及旋转角αd的关系式为式(14)。

(14)

其中：

(15)

如图7所示，其主轴电导率为σ1=0.01 S/m,σ2=0.1 S/m,σ3=0.01 S/m。模型空间被分为186 627个非结构四面体单元，依次设置倾角大小为αd=0°、30°、60°、90°，频点为1 Hz，讨论倾角αd对测点MT响应的影响。二维情况下，视电阻率ρxy与φxy相位不受倾角αd影响，仅与主轴电导率σ1有关[18]。而三维情况下，不同于二维模拟的结论，视电阻率ρxy与相位φxy曲线随着倾角αd会发生较小变化。视电阻率ρxy与相位φxy受到较大影响，随着倾角变大，其异常体范围(-1 000 m,1 000 m)处的视电阻率值逐渐变小，当倾角为0°、90°时，曲线保持其对称性，当倾角为30°、90°时，曲线关于y=0非对称。两种模式视电阻率的差异以及非对称性与倾角αd有显著关联，因此可以作为判断异常体各向异性特征的依据。

4.2 主轴各向异性

主轴各向异性介质的电导率张量，只有对角线存在非零元素：

(16)

设计球状异常体模型，讨论主轴各向异性以及水平各向异性，背景电导率为0.01 S/m，球体半径为500 m，模型空间被分为183 835个非结构四面体单元，网格剖分见图10。VTI介质与VHI介质是主轴各向异性介质中两种常见情况，分别对两种介质进行模拟计算MT响应。

图9 倾斜各向异性模型视电阻率与相位随倾角αd变化曲线Fig.9 Dipping anisotropy model apparent resistivity and phase changing with different αd(a)ρxy;(b)ρyx;(c)φxy;(d)φyx

图10 球状异常体模型网格剖分Fig.10 Spheroidal anomaly model and discretation using unstructured girds

4.2.1 VTI介质

设计球型异常体主轴电阻率σ1=σ2=0.1 S/m保持不变，垂直主轴电导率σ3=1、0.1、0.01、0.001 S/m相相对背景电导率，经历几个量级的转变，频点设为1 Hz其对测点MT响应的影响见图11，可以看到：对于VTI介质，虽然垂直主轴电导率σ3经历了几个量级的变化，但测点处的视电阻率与相位曲线仅有微小变化，说明主轴各向异性条件下，σ3对地面测点的MT响应影响微弱。

4.2.2 VHI介质

设计异常体主轴电导率为：σ1=σ3=0.1 S/m保持不变，水平y方向主轴电导率σ2=1、0.1、0.01、0.001 S/m从从大到小变化，频点设为1 Hz，观察其对测点MT响应的影响。如图12所示，σ2经历几个数量级的变化，视电阻率ρxy与相位φxy几乎不受影响，然而视电阻率ρxy与相位φxy则发生了明显的变化，随着主轴电导率分量σ2从大到小变化，视电阻率曲线的异常体位置范围(-1 000 m,1 000 m)内也表现由低阻到高阻的变化特征，这与文献[16]的结论是一致的。

图11 VTI介质视电阻率与相位随σ3变化曲线Fig.11 Apparent resistivity and phase of VTI media with different σ3(a)ρxy;(b)ρyx;(c)φxy;(d)φyx

4.3 水平各向异性

水平各向异性是x轴与y轴在xoy平面内绕z轴旋转一定的角度，其余角度为零，其电导率张量表达式为式(18)。

(17)

(18)

沿用球模型，均匀半空间电阻率与球体位置大小均不变，研究水平旋转角αs变化对测点MT响应的影响。设异常体的电导率参数为：σ1=0.1 S/m,σ2=0.001 S/m,σ3=0.01 S/m，水平旋转角由零逐渐变大。过模型空间中心点的x方向测线在1 Hz下的视电阻率与相位曲线(图13)可以看到，视电阻率与相位均受到αs的显著影响，当αs=0°时，此时介质为主轴各向异性，ρxy和ρyx对应图13(a)绿色曲线；随着αs变大，ρxy曲线中心由低阻向高阻转变，相反ρyx曲线由高阻逐渐转变为低阻，原因在于，σxx的值逐渐向σ2过渡，而σyy逐渐向σ1变化所导致，直到αs=90°，此时σxx=σ2，σyy=σ1。

图12 VHI介质视电阻率与相位随σ2变化曲线Fig.12 Apparent resistivity and phase of VHI media with different σ2(a)ρxy;(b)ρyx;(c)φxy;(d)φyx

图13 过中心点测线的视电阻率与相位随αs变化曲线Fig.13 Central survey line apparent resistivity and phase with different αs(a)ρxy;(b)ρyx;(c)φxy;(d)φyx

图14 视电阻率在1 Hz平面内随αs变化Fig.14 Apparent resistivity in the 1 Hz plane with different αs

图16 αs=45°、γ不同情况下、ρxy旋转程度大小Fig.16 In condition of αs=45°，γ changing , the degree of ρxy rotation(a)γ=1.4;(b)γ=3.2;(c)γ=7.7;(d)γ=14.1

5 结论

实现了电导率任意各向异性大地电磁三维非结构矢量有限元法正演数值模拟，一维层状解析解、COMMEMI模型以及各向异性模型的计算结果与前人模拟结果很好吻合，验证了算法的准确性。在此基础上，进行了不同电导率各向异性模型的数值模拟，结果表明，倾斜各向异性条件下，视电阻率ρxy曲线随着倾角αs会发生较小变化，而视电阻率ρyx受到较大影响而产生非对称特征。主轴各向异性条件下，α3对视电阻率产生较小影响，而ρxy主要反映α1的信息，ρyx主要反映α2的信息。特别在水平各向异性条件下，视电阻率形态随着走向角αs变化而发生旋转，只有当水平各向异性系数γ非常大时，视电阻率形态的旋转角能够比较准确地反映αs的大小；当水平各向异性系数γ相对较小时，视电阻率形态的旋转角与αs差别较大，必须同时考虑水平各向异性系数γ的影响。实际地球深部介质的电导率各向异性比较复杂、大小不一，笔者大地电磁各向异性数值模拟获得的一些新的认识，对深部岩石圈电性结构探测有重要意义。