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

doi:10.11720/wtyht.2020.0343

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

顾观文^,¹^,², 武晔¹^,², 石砚斌¹^,²

1.防灾科技学院地球科学学院,河北廊坊 065201

2.河北省地震动力学重点实验室,河北廊坊 065201

Research on fast three-dimensional forward algorithm of magnetotelluric sounding based on vector finite element

GU Guan-Wen^,¹^,², WU Ye¹^,², SHI Yan-Bin¹^,²

1. School of Earth Sciences, Institute of Disaster Prevention, Langfang 065201, China

2. Hebei Key Laboratory of Earthquake Dynamics, Langfang 065201, China

责任编辑: 沈效群

收稿日期: 2020-07-2 修回日期: 2020-09-10 网络出版日期: 2020-12-20

基金资助:

中央高校创新团队项目. ZY20180102
自然资源部“十五”重点科技项目. 20010211

Received: 2020-07-2 Revised: 2020-09-10 Online: 2020-12-20

作者简介 About authors

顾观文(1975-),男,教授级高级工程师,从事电磁法正反演研究及其软件研制工作。 Email:sun_ggw@163.com

摘要

本文采用直接求解器PARDISO且无需散度校正的正演方案,求解矢量有限元法对应的大型线性方程组,获得不同地形(水平和起伏)条件下地电模型的大地电磁响应,提高了大地电磁三维正演计算的速度。在中等规模计算条件下,通过本文的计算方法与带散度校正的迭代求解法对比,计算速度可提高十倍以上。

关键词： 大地电磁 ; 矢量有限元法 ; 三维正演 ; PARDISO

Abstract

The finite element method has the characteristics of strong adaptability in simulating the electromagnetic response of rugged topography and complex geological bodies. In recent years, it has been widely used in the three-dimensional (3D) forward modeling of magnetotelluric (MT) sounding. However, the finite element method also has some shortcomings in terms of computational efficiency. The large amount of calculation and long running time of the method are the main factors that lead to the lag of the practical process of the 3D MT inversion technology based on the finite element method compared with the 3D MT inversion technology based on the finite difference method. In order to improve the 3D forward speed of MT, the authors adopt the forward modeling scheme which uses the direct solver PARDISO and does not need divergence correction to solve the large-scale linear equations corresponding to the vector finite element method, and obtain the MT response of the geoelectric model under such different terrain conditions as flat and rugged topography. Under the conditions of medium-scale calculation, through the comparison between the direct solution method without divergence correction and the iterative solution method with divergence correction, the authors have detected that the direct solution method without divergence correction has advantages in calculation accuracy and calculation time, especially in the calculation. In terms of time, the ratio of the calculation speed of the direct solution and the iterative solution is raised by more than ten times.

Keywords： magnetotelluric ; vector finite element method ; 3D forward ; PARDISO

PDF (1541KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

顾观文, 武晔, 石砚斌. 基于矢量有限元的大地电磁快速三维正演研究. 物探与化探[J], 2020, 44(6): 1387-1398 doi:10.11720/wtyht.2020.0343

GU Guan-Wen, WU Ye, SHI Yan-Bin. Research on fast three-dimensional forward algorithm of magnetotelluric sounding based on vector finite element. Geophysical and Geochemical Exploration[J], 2020, 44(6): 1387-1398 doi:10.11720/wtyht.2020.0343

0 引言

大地电磁测深法(magnetotelluric sounding, MT)具有施工方便、勘探效率高、成本低 (相对于地震勘探)、勘探深度大等优点,目前已被广泛应用于资源勘查、能源勘探及深部构造探测等方面。对于大地电磁数据的反演和解释,由于受限于理论方法及计算设备的运算能力,早期以一维或二维反演技术为主。20 世纪90年代以来,随着计算设备运算能力的提高及数值计算方法的进步,大地电磁三维正反演技术开始逐步发展,不同的三维正演方法(积分方程、有限差分、有限元等)及其计算技术取得了巨大进展^{[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]}。目前在实际中得到应用的三维反演技术主要是基于有限差分法三维正演的反演方法,特别是国内实测大地电磁资料的三维解释基本上都采用基于有限差分法的三维反演技术^{[16,17,18,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)$\nabla \times E=i \omega \mu H$,

(1b)$\nabla \times H=\sigma E$,

(1c)$\nabla \times E=\frac{q}{\varepsilon}$,

(1d)$\nabla \times H=0$,

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

(2)$\nabla \times (\frac{1}{i \omega \mu} \nabla \times E)=\ sigma E$,

求解偏微分方程组(2)得到电场分量E_x、E_y、E_z后,根据式(1a)求得磁场分量H_x、H_y、H_z。

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

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 带地形三维MT数值模拟区域剖面示意^[11]

Fig.1 Section diagram of numerical modeling domain for 3D MT with topography^[11]

将式(2)写成

(3)$\nabla \times \nabla \times E- i \omega \mu \sigma E=0,(x,y,z) \in \Omega$,

取第一类边界条件:

(4)

E (x, y, z) |_{\partial Ω} = g (x, y, z) |_{\partial Ω},

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

1.2 矢量有限元分析

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

(5)

E_{x}^{e} = \overset{4}{\sum_{i = 1}} N_{xi}^{e} E_{xi}^{e}, E_{y}^{e} = \overset{4}{\sum_{i = 1}} N_{yi}^{e} E_{yi}^{e}, E_{z}^{e} = \overset{4}{\sum_{i = 1}} N_{zi}^{e} E_{zi}^{e},

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 矢量有限元法的区域剖分示意

a—区域剖分示意; b—电场分量位置示意

Fig.2 Domain subdivision of the vector finite element method

a—domain subdivision; b—location of electric field components

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

(6)

E^{e} = \overset{12}{\sum_{i = 1}} N_{i}^{e} E_{i}^{e},

式中: $N_{i}^{e}$ =( $N_{xi}^{e}$ , $N_{yi}^{e}$ , $N_{zi}^{e}$ )为矢量有限元中的基函数, $E_{i}^{e}$ 为单元棱边上的未知电场值,上标e表示第e个网格单元。即:

(7)

N_{i = [1,4]}^{e} = \frac{1}{4} (1 + η η_{i}) (1 + ζ ζ_{i}) x,

(8)

N_{i = [5,8]}^{e} = \frac{1}{4} (1 + ζ ζ_{i}) (1 + ξ ξ_{i}) y,

(9)

N_{i = [9,12]}^{e} = \frac{1}{4} (1 + ξ ξ_{i}) (1 + η η_{i}) z,

式中:坐标转换函数为ξ=(x-x_c)/a,η=(y-y_c)/b,ζ=(z-z_c)/c;(x_c,y_c,z_c)是六面体的中心坐标,2a、2b、2c分别是六面体x、y、z方向的长度;(ξ_i,η_i,ζ_i)的取值与棱边编号有关。

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

(10)$r^{e}=\nabla \times \nabla \times E^{e}- i \omega \mu \sigma E^{e}$,

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

(11)

R = \overset{NE}{\sum_{e = 1}} \underset{V_{e}}{\int \int \int} r^{e} \cdot N_{i}^{e} d v \to 0,

式中:NE为计算域网格单元总数。将式(10)代入式(11)中的r^e· $N_{i}^{e}$ ,则:

(12)$r^{e} · N^{e}_{i}=(\nabla \times \nabla \times E^{e}) · N^{e}_{i}- i \omega \mu \sigma E^{e} · N^{e}_{i}$

利用矢量恒等式:

(13)$(\nabla \times A)· B=\nabla (A \times B)+(\nabla \times B)· A$,

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

(14) $(\nabla \times \nabla \times E^{e}) · N^{e}_{i}=\nabla · (\nabla \times E^{e} \times N^{e}_{i})+\nabla \times N^{e}_{i} · \nabla \times \nabla \times E^{e}$

根据高斯散度定理,可知:

(15)

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

(16)

式中:n是外法线方向,S_e是单元面积分,V_e是单元体积分。由于$[n \times (\frac{1}{ i \omega \mu_{0}} \ nabla \times E^{e})]=[n \times H^{e}]$是连续的,因此式(16)右端第三项在单元的装配过程中互相抵消,第e单元的电场所满足的方程可以表示为

(17)

K^{e} E^{e} = S^{e},

式中:S^e表示场源项;E^e表示棱边上的电场;K^e为单元刚度矩阵,是一个12×12阶的复数矩阵,可按下式解析计算^[32]得出:

(18)

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

(19)

K \cdot E = s,

式中: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分解,融合了向左和向右算法的优点,采用超节点消去树进行填充元优化,降低算法的时空消耗。对于复线性方程组有

(20)

Ax = b,

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

① 矩阵重排与符号分解:针对稀疏矩阵A的结构,设计合适的行列交换矩阵,对原矩阵进行交换与重排,使得新矩阵 $\dot{A}$ 分解后含有尽量少的非零元素。

② 矩阵LU分解:对新矩阵 $\dot{A}$ 进行LU分解。

③ 方程求解与迭代:利用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)

[\begin{array}{l} 1 & 0 & 5 & 0 \\ 0 & 2 & 0 & 3 \\ 5 & 0 & 2 & 0 \\ 0 & 3 & 0 & 1 \end{array}] \begin{array}{l} a [6] = [152321] \\ ja [6] = [132434] \\ ia [5] = [13567] \end{array}

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

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

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

根据Newman等^[34]、谭捍东等^[35]的研究,假设两种线性无关的场源激发的表面电场和磁场分别为E_x₁、E_y₁、H_x₁、H_y₁、E_x₂、E_y₂、H_x₂、H_y₂,由此可以计算三维MT的张量阻抗:

(22)

Z = [\begin{array}{l} Z_{xx} & Z_{xy} \\ Z_{yx} & Z_{yy} \end{array}],

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

(23a)

$Z_{xx} = \frac{E_{x 1} H_{y 2} - E_{x 2} H_{y 1}}{H_{x 1} H_{y 2} - H_{x 2} H_{y 1}},$

(23b)

$Z_{xy} = \frac{E_{x 2} H_{x 1} - E_{x 1} H_{x 2}}{H_{x 1} H_{y 2} - H_{x 2} H_{y 1}},$

(23c)

$Z_{yx} = \frac{E_{y 1} H_{y 2} - E_{y 2} H_{y 1}}{H_{x 1} H_{y 2} - H_{x 2} H_{y 1}},$

(23d)

$Z_{yy} = \frac{E_{y 2} H_{x 1} - E_{y 1} H_{x 2}}{H_{x 1} H_{y 2} - H_{x 2} H_{y 1}},$

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

(24)

(ρ_{a})_{ij} = \frac{1}{ω μ_{0}} | Z_{ij} |^{2}, ϕ_{ij} = Arg (Z_{i, j}) 。

式中: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,需要模型的电阻率为有限值,空气层的电阻率设为一个较大的常数,一般在10⁶~10¹⁰ Ω·m之间,本文涉及的三维模型其空气层电阻率均取为10¹⁰ Ω·m,表1、图3为数值模拟结果。从表1可以看出,视电阻率误差最高为0.743 7%(频点为50 Hz),其余频点计算误差均未超过0.6%,尤其在0.1 Hz以后误差均在0.1%以下。相位计算误差最高为0.233 56%(频点为100 Hz),其余频点计算误差均在0.2%以下。

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

Table 1 Comparison of vector finite element solution and analytical solution of uniform half space model

频点 /Hz	视电阻率(ρ_xy,ρ_yx)/(Ω·m)			相位φ/(°)
频点 /Hz	矢量有限元解	解析解	误差/%	矢量有限元解	解析解	误差/%
10000	99.4934	100	0.5066	44.9824	45	0.039111
8000	99.5168	100	0.4832	44.9800	45	0.044444
5000	99.5580	100	0.442	44.9782	45	0.048444
2000	99.6149	100	0.3851	44.9802	45	0.044
1000	99.6266	100	0.3734	44.9895	45	0.023333
500	99.7002	100	0.2998	44.9601	45	0.088667
200	99.8467	100	0.1533	45.0661	45	0.14689
100	99.4868	100	0.5132	45.1051	45	0.23356
50	99.2563	100	0.7437	45.0452	45	0.10044
10	99.4329	100	0.5671	44.9565	45	0.096667
5	99.5295	100	0.4705	44.9596	45	0.089778
2	99.6083	100	0.3917	44.9717	45	0.062889
1	99.6484	100	0.3516	44.9776	45	0.049778
0.5	99.6599	100	0.3401	44.9970	45	0.006667
0.1	99.7267	100	0.2733	44.9301	45	0.155333
0.05	99.9013	100	0.0987	44.9393	45	0.134889
0.01	100.0090	100	0.009	44.9873	45	0.028222
0.005	100.0070	100	0.007	44.9949	45	0.011333
0.001	100.0010	100	0.001	44.9995	45	0.001111
0.0005	100.0000	100	0	44.9998	45	0.000444
0.0001	100.0000	100	0	45.0000	45	0

新窗口打开| 下载CSV

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 均匀半空间模型三维正演视电阻率(a)和相位(b)与解析解对比

Fig.3 Comparison of 3D forward apparent resistivity (a) and phase (b) of homogeneous half space model with analytical solution

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

3) H型层状模型

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

表2 H型层状模型参数

Table 2 Parameters of H-type layered model

层参数	第一层	第二层	第三层
电阻率/(Ω·m)	100	10	1000
层厚/m	370	268	∞

新窗口打开| 下载CSV

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

Table 3 Comparison of vector finite element solution and analytical solution of H-type layered model

频点 /Hz	视电阻率(ρ_xy,ρ_yx)/(Ω·m)			相位φ/(°)
频点 /Hz	矢量有限元解	解析解	误差/%	矢量有限元解	解析解	误差/%
10000	99.4944	100	0.505638	44.9825	45.00002	0.03894
8000	99.5156	99.99967	0.484074	44.9807	45.00007	0.00043
5000	99.564	100.0036	0.439562	44.9729	44.9985	0.000569
2000	99.1812	99.72298	0.543282	45.0089	45.02392	0.000334
1000	100.138	100.1248	0.01314	44.3103	44.43167	0.002732
500	108.223	107.9785	0.22646	44.9127	44.67756	0.00526
200	112.433	113.6883	1.104134	51.6513	51.46143	0.00369
100	94.4807	95.60484	1.175822	59.265	59.07526	0.00321
50	63.587	64.76429	1.817807	63.9834	63.85949	0.00194
10	28.0673	28.37355	1.079333	45.849	46.13936	0.006293
5	32.2082	32.33247	0.384363	32.3095	32.54673	0.007289
2	55.1794	55.21716	0.06839	21.4561	21.55953	0.004797
1	89.6415	89.67347	0.035652	18.7262	18.78408	0.003081
0.5	143.375	143.3732	0.00122	19.0537	19.09105	0.001956
0.1	348.002	347.851	0.04341	25.4995	25.50588	0.00025
0.05	457.716	457.5438	0.03763	29.036	29.03635	1.19E-05
0.01	692.088	691.9579	0.0188	36.1384	36.13554	7.9E-05
0.005	769.11	769.008	0.01327	38.3803	38.37769	6.8E-05
0.001	888.395	888.3428	0.00588	41.8079	41.80643	3.5E-05
0.0005	919.642	919.6043	0.0041	42.7015	42.70037	2.7E-05
0.0001	963.197	963.179	0.00187	43.9466	43.94614	1.1E-05

新窗口打开| 下载CSV

图4

新窗口打开| 下载原图ZIP| 生成PPT

图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型层状模型参数

Table 4 Parameters of K-type layered model

层参数	第一层	第二层	第三层
电阻率/(Ω·m)	100	1000	10
层厚/m	370	268	∞

新窗口打开| 下载CSV

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

Table 5 Comparison between vector finite element solution and analytical solution of K-type layered model

频点 /Hz	视电阻率(ρ_xy,ρ_yx)/(Ω·m)			相位φ/(°)
频点 /Hz	矢量有限元解	解析解	误差/%	矢量有限元解	解析解	误差/%
10000	99.4928	99.99995	0.507154	44.9823	44.99998	0.039291
8000	99.5179	100.0004	0.482464	44.9795	44.99994	0.04542
5000	99.5522	99.9958	0.443615	44.9822	45.00155	0.042987
2000	99.9938	100.3013	0.306591	44.9731	44.9962	0.05134
1000	98.5514	99.11594	0.569571	45.4799	45.46765	0.02695
500	94.8705	94.84442	0.0275	43.8641	44.11071	0.55908
200	111.446	109.5436	1.73666	42.1518	41.35312	1.93136
100	124.19	128.152	3.091671	47.226	45.89309	2.90438
50	115.34	122.0064	5.46397	54.8037	54.506	0.54618
10	55.4147	56.83843	2.504864	64.9559	65.42769	0.721085
5	38.6174	39.2693	1.660063	64.5842	64.96375	0.584241
2	25.6211	25.88921	1.03559	61.7858	62.03061	0.394665
1	20.052	20.20837	0.773778	58.9626	59.13016	0.28337
0.5	16.5801	16.68206	0.611172	56.1309	56.24364	0.200444
0.1	12.6035	12.65607	0.415342	50.8822	50.92886	0.091626
0.05	11.7748	11.82107	0.391403	49.346	49.36433	0.037132
0.01	10.7493	10.77999	0.284731	46.9822	47.06061	0.166611
0.005	10.5349	10.54576	0.102932	46.4109	46.47605	0.140175
0.001	10.2412	10.24057	0.00616	45.6579	45.67163	0.030065
0.0005	10.17	10.16953	0.00462	45.4711	45.47688	0.012716
0.0001	10.0755	10.07547	0.00035	45.2136	45.21445	0.001873

新窗口打开| 下载CSV

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 K型层状模型三维正演视电阻率(a)和相位(b)数值解与解析解对比

Fig.5 Comparison of 3D forward apparent resistivity (a) and phase (b) of K-type layered model with analytical solution

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方向。

图6

新窗口打开| 下载原图ZIP| 生成PPT

图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)计算结果,可以看出二种方法的计算结果基本一致。

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 本文矢量有限元正演算法的计算结果与IE方法的计算结果对比

a—Z_xy模式正演视电阻率;b—Z_yx模式正演视电阻率;c—Z_xy模式正演阻抗相位; d—Z_yx模式正演阻抗相位

Fig.7 Comparison between the calculation results of vector finite element forward algorithm and IE method

a—Z_xy mode forward apparent resistivity;b—Z_yx mode forward apparent resistivity;c—Z_xy mode forward impedance phase;d—Z_yxmode forward impedance phase

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

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

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8 三维地形及其网格剖分示意

a—三维地形示意; b—地形网格剖分和MT测点分布示意

Fig.8 3D topography and grid

a—sketch of 3D topography; b—topography meshing and distribution of MT measurement sites

为了验证正演程序对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

新窗口打开| 下载原图ZIP| 生成PPT

图9 二维山峰地形示意

Fig.9 Sketch of 2D ridge

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

图10

新窗口打开| 下载原图ZIP| 生成PPT

图10 三维矢量有限元算法(PARDISO)计算的二维地形影响与二维有限元结果对比

a—正演视电阻率; b—正演阻抗相位

Fig.10 Comparision between modeling results of 3DVFEM(PARDISO) and 2DFEM for 2D ridge

a—forward apparent resistivity; b—forward impedance phase

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倍。可见本文的无需散度校正直接求解法与带散度校正的迭代求解法相比,在计算精度和计算时间方面均有优势,特别是在计算时间方面表现出明显优势,适合应用于中等计算规模的三维反演算法中。

图11

新窗口打开| 下载原图ZIP| 生成PPT

图11 无需散度校正的直接解法与带散度校正的迭代解法计算结果对比曲线

a—正演视电阻率; b—正演阻抗相位

Fig.11 Comparison of calculation results of VFE-PARDISO without divergence correction and VFE-BICG with divergence correction

a—forward apparent resistivity; b—forward impedance phase

4 结论

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

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Wannmaker P

Advances in three-dimensional magnetotelluric modeling using integral equations

[J]. Geophysics, 1991,56:1716-1728.

[本文引用: 3]

[2]

徐凯军, 李桐林, 张辉, 等.

利用积分方程法的大地电磁三维正演

[J]. 西北地震学报, 2006(2):104-107.