大地电磁三维矢量有限元正演模拟

石 明, 冯德山, 李开鹏，王 珣

(1.中南大学 地球科学与信息物理学院，湖南 长沙，410083；2. 有色资源与地质灾害探查湖南省重点实验室, 湖南 长沙，410083；3. 贵州省有色金属和核工业地质勘查局物化探总队，贵州 都匀，558004)



大地电磁三维矢量有限元正演模拟

石 明1,3, 冯德山1,2*, 李开鹏3，王 珣1,2

(1.中南大学 地球科学与信息物理学院，湖南 长沙，410083；2. 有色资源与地质灾害探查湖南省重点实验室, 湖南 长沙，410083；3. 贵州省有色金属和核工业地质勘查局物化探总队，贵州 都匀，558004)

从Maxwell方程出发，开展了三维大地电磁场所满足的边值问题研究，利用加权余量法导出了三维大地电磁有限元方程.介绍了三维矢量有限元六面体网格剖分方式、插值基函数选取，推导了三维大地电磁矢量有限元正演的单元刚度系数矩阵及离散格式.编制了三维矢量有限元大地电磁正演的Matlab程序.三维COMMEMI 3D-1模型的视电阻率曲线与国际通用的标准测试数据能很好地拟合，验证了作者编写的矢量有限元正演程序的正确性.通过对高、低阻异常体的阻抗张量形态分析，说明张量阻抗等值线图能用以大致判断异常体特性，丰富了大地电磁响应特征的表达方式.

矢量有限元；大地电磁；正演模拟；张量阻抗

大地电磁(MT)是以电离层激发的天然交变电磁场为场源，在地表观测相互正交的电场、磁场分量来获取地电构造信息的一种重要地球物理勘探方法[1].MT不需要庞大的发射源设备，只需采用比较轻便的接收设备，野外工作方便、成本低，被广泛应用于地壳和上地幔电性结构的研究,在石油天然气勘探、矿产资源勘探、工程与环境普查等领域,发挥着举足轻重的作用[2-9].可以预见，三维MT勘探技术是地球物理中深层领域的研究热点及今后MT的发展趋势，而三维MT正演是理解MT勘探物理现象并认识地质体电磁响应规律的有效手段，显然尤其重要.

尽管矢量FEM拥有诸多优点，但在地球物理的电磁法正演领域中，其应用并不多见，尚需要进一步完善.目前的研究主要包括：Yoshimura 等[10]开展了矢量FEM的MT响应数值模拟，并将矢量FEM的计算结果与交错网格FDM的计算结果进行了对比；Mitsuhata等[11]利用矢量FEM和节点FEM耦合的方法对三维MT数值模拟；Nam[12]采用不规则六面体矢量FEM直接计算电场，研究了起伏地形下MT的电阻率和相位的变化规律；刘长生等[13]将完全非结构化四面体单元引入到矢量有限元中，实现了三维大地电磁h-型自适应矢量有限元正演；王烨[14]开展了高频率大地电磁法矢量有限元正演，并采用改进的威尔金森方法求解大型病态方程组，提高了迭代速度；顾观文等[15]开展了矢量有限元法MT三维地形数值模拟，研究了地形起伏下三维阻抗张量的变化规律；杨军等[16]采用非结构四面体单元的三维矢量FEM实现了海洋可控源电磁数值模拟；苏晓波等[17]采用规则六面体单元的三维矢量FEM实现了大地电磁数值模拟，并对网格剖分的重要性进行了研究.

在前人基础上，作者推导了三维大地电磁矢量FEM正演的离散形式，应用矢量FEM算法计算了三维COMMEMI 3D-1国际模型[18]的MT视电阻率及模型张量阻抗，研究了高低阻异常体的电磁响应特性，有效地指导了MT的资料解释.

1 三维大地电磁边值问题

以e-iwt表示谐变场的时间因子，大地电磁满足的Maxwell方程组可表示为：

(1)

式中：ω为角频率；μ为介质的磁导率；σ为电导率；ε为介电常数.联立两式，消去H，并假定岩石和空气中的μ为常数得：

(2)

式(2)为电场矢量E所应满足的微分方程,式中k为波数,k2=iωμσ+ω2εμ.将电场E写成分量形式为：

E=Exex+Eyey+Ezez．

(3)

在三维情况下，Ex，Ey和Ez是相互联系的，无法将其中一个分量单独分离出来.根据广义变分原理，将微分方程和边界条件都考虑在内，利用加权余量法[9]推导得到的三维大地电磁边值问题相关的加权余量方程为：

(4)

2 矢量FEM网格离散及单元分析

矢量FEM与常规节点FEM求解过程非常类似，仅仅在区域剖分、插值基函数选取及单元分析方面略有不同，主要步骤包括：(ⅰ) 区域剖分；(ⅱ) 选用矢量插值基函数；(ⅲ) 单元分析；(ⅳ) 总体合成；(ⅴ) 多元函数求偏导数；(ⅵ) 解线性代数方程组.

图1 矢量有限元三维区域剖分图

图2 矢量六面体单元中的节点与棱边编号

通过分配常切向场分量给单元的每一条边，单元内x,y,z3个方向上场分量可分别表示为：

(5)

其中：

(6)

矢量FEM定义的基函数具有零散度和非零旋度的特性.而且，构成离散单元各个小平面上的切向场仅由组成小平面的棱边上的切向场决定，不仅保证了穿越棱边的切向场连续，而且保证了穿越离散单元表面时切向场的连续性.用六面体单元对整个区域进行剖分.单元的坐标变换关系为：

x=xc+ξ·a/2,

y=yc+η·b/2,

z=zc+ζ·c/2．

(7)

式中：xc，yc，zc是子单元中心的坐标，微分关系为：

(8)

将式(4)中的区域积分分解为各单元积分之和：

(9)

由插值基函数形式，可得到

(10)

δEx,δEy,δEz分别为电场Ex,Ey,Ez的变分.将式(10)带入，可得到式(9)中的第一项积分为：

∫e×δE·×EdΩ=∫eδET·(×NT)×

(11)

(12)

将形函数代入并进行积分求解，可以得到式(9)第一积分项的三维单元刚度矩阵系数：

其中N1，N2，N3均为经过三重积分后得到的4×4的单元矩阵，结果如下：

将形函数代入式(9)的第二项单元积分中：

(13)

式中：

K11=K22=K33=∫ek2NNTdΩ，

(14)

求解式(14)，可以得到式(9)第二项的单元刚度矩阵系数为：

(15)

当单元的边界面5678落在区域底面边界EFGH上时，边界积分

(16)

式中：

K11=∫χNNTdΓ，

(17)

这里Ni，Nj为四边形单元的矢量基函数，不是六面体单元的基函数.求解式(17)，可以得到式(9)第3积分项的三维矢量FEM单元刚度矩阵系数：

(18)

在一个单元内，将三项积分得出的矩阵K1e，K2e，K3e相加，再将单元的系数矩阵扩展成由全体棱边组成的系数矩阵，再将各单元相加，对(9)式进行离散、扩展，将顶界面ABCD上的边界条件代入线性方程组后，解方程组，就可以得到各个棱边的Ex,Ey,Ez的值.根据式(1)中第一式，进行有限元求解即可得到相应的Hx，Hy，Hz的值.

3 三维矢量FEM大地电磁正演模拟

为了验证三维矢量FEM程序的正确性和精度，选取国际通用COMMEMI 3D-1 标准测试模型，模型示意图如图3所示.

图3中异常体大小为1 000 m×2 000 m×2 000 m，异常体上界面据地面250 m，电阻率为0.5·m，背景电阻率为100·m，MT测线分别布置在x，y坐标轴上，范围为0 km～2.5 km.其中x，y，z方向网格数为 41×41×39(8 层空气层)，计算区域为30 000 m×25 000 m×35 000 m， 空气层厚度为 5 000 m，空气层电阻率设为 1015·m .

图3 COMMEMI 3D-1标准测试模型

图4为f=10 Hz时XY模式与YX模式下矢量FEM正演模拟所得到的视电阻率曲线.VFEM3D表示三维矢量FEM计算结果，黑色竖线为COMMEMI提供的误差值.由图中可见，f=10 Hz时矢量FEM的计算结果与COMMEMI提供的结果尽管存在细小误差，但相差不大，且在误差允许范围之内，表明在该频率下矢量FEM的计算结果可信.在图4(a)XY模式下矢量FEM计算结果较COMMEMI偏小，特别是第一个点已经超出误差棒的范围.这一现象普遍存在，推断为矢量FEM后处理精度不足引起的.而图4(b)YX模式下VFEM3D结果与COMMEMI结果吻合相当好.

x/mm

x/mm

图5(a)，(b)分别为f=0.1 Hz时XY模式与YX模式下矢量FEM正演视电阻率曲线.分析图5(a),(b)可知，两幅图中的矢量FEM曲线与COMMEMI所提供的数据都能够很好地吻合，说明无论是在低频还是高频部分，应用矢量FEM开展三维大地电磁正演，都具有较高的精度，同时也验证了矢量FEM算法及程序的正确性.

x/mm

x/mm

图6为应用矢量FEM正演计算COMMEMI3D-1模型得到的张量阻抗.由图可见，10 Hz与0.1 Hz两个频率下的张量阻抗形态基本一致，10 Hz的数值较0.1 Hz要大.对比图中4个不同的张量阻抗，可以发现，图6(a),(d),(e)和(h)中两个频率下的Zxx与Zyy分为四瓣，且阻抗值较小，四瓣的中心反映了异常体的边界，而图6(b),(c),(f)和(g)中的Zxy与Zyx阻抗值较大，反映了入射场的特性.根据张量阻抗理论可知，当构造为二维构造时，Zxx和Zyy为零，即当异常体走向方向越长，Zxy与Zyx越小，Zxy与Zyx差异也越大.由此，张量阻抗分解后，无需做反演即可以判断出异常体的简单特性.

图6 COMMEMI 3D-1A模型单个低阻异常体张量阻抗等值线图

为了进一步认识大地电磁的响应特性，对比高、低阻异常体张量阻抗的不同，在图3中COMMEMI3D-1测试模型的基础上，仅将低阻异常体改为1 000 Ω·m高阻异常体.其他参数均与国际模型相同.应用三维矢量FEM开展三维高阻异常体模型的张量阻抗研究.

图7为应用三维矢量FEM正演的10 Hz大地电磁张量阻抗图.分析图7(a)与图7(d)可知，高阻异常体张量阻抗中的Zxx与Zyy同样分为四瓣，且阻抗值较小，其四瓣的中心反映了异常体的边界.由于异常体x方向与y方向的比值为1∶2，图7(b)中的Zxy与图7(c)中的Zyx差异较大.对比高低阻异常10 Hz时的阻抗相位Zxx，虽然两者都为四瓣，但是阻抗值正负值的分布正好相反，低阻异常体四瓣的中心向外辐射，幅值变小趋于0；而高阻异常体四瓣的中心向外辐射，幅值变小趋于0之后会发生反转之后再次趋于0.Zyy具有相同的规律.对比Zxy和Zyx，高阻异常体中心仅出现一个闭合异常形态，而低阻异常体则形态更为复杂.

图7 COMMEMI 3D-1A模型单个高阻异常体张量阻抗等值线图

4 结 论

1) 介绍了三维矢量有限元区域剖分方式，对矢量FEM插值基函数以及单元插值方式进行了阐述，应用Galerkin算法，推导了三维矢量FEM大地电磁方程离散格式，编制了矢量FEM三维MT的Matlab模拟程序.

2) 设置三维COMMEMI 3D-1模型进行矢量FEM的计算，模拟结果与COMMEMI提供的数据拟合效果很好，验证了矢量有限元程序的正确性.通过对比高低阻异常体的张量阻抗，分析了不同异常下张量阻抗的特点，进一步认识了MT的响应特性.

[1] 柳建新，童孝忠，郭荣文,等. 大地电磁测深勘探：资料处理反演与解释[M]. 北京：科学出版社，2012：1-12.

LIU jian-xin, TONG Xiao-zhong, GUO Rong-wen,etal. Magnetotelluric sounding exploration: data processing inversion and interpretation [M].Beijing: Science Press, 2012:1-12.(In Chinese)

[2] 底青云，王若. 可控源音频大地电磁数据正反演及方法应用[M]. 北京：科学出版社，2008：1-8.

DI Qing-yun, WANG Ruo. Controlled source audio magnetotelluric data inversion and methods application [M].Beijing: Science Press, 2008: 1-8.(In Chinese)

[3] 谭捍东，余钦范，JOHN B，等. 大地电磁法三维交错采样有限差分数值模拟[J]. 地球物理学报，2003，46(5)：705-711.

TAN Han-dong, YU Qin-fan, JOHN B,etal. Magnetotelluric three-dimensional modeling using the staggered-grid finite difference method[J]. Chinese Journal of Geophysics, 2003, 46(5):705-711.(In Chinese)

[4] 陈辉，邓居智，谭捍东,等. 大地电磁三维交错网格有限差分数值模拟中的散度校正方法研究[J]. 地球物理学报，2011，54(6)：1649-1659.

CHEN Hui, DENG Ju-zhi, TAN Han-dong,et,al. Study on divergence correction method in three-dimensionalmagnetotelluric modeling with staggered-grid finite difference method [J]. Chinese Journal of Geophysics, 2011, 54(6): 1649-1659. (In Chinese)

[5] 李焱，胡祥云，杨文采,等. 大地电磁三维交错网格有限差分数值模拟的并行计算研究[J]. 地球物理学报，2012，55(12)：4036-4043.

LI Yan, HU Xiang-yun, YANG Wen-cai,etal. A study on parallel computation for 3D magnetotelluric modeling using the staggered-grid finite difference method [J]. Chinese Journal of Geophysics, 2012, 55(12): 4036-4043. (In Chinese)

[6] 徐凌华，童孝忠，柳建新，等. 基于有限单元法的二维/三维大地电磁正演模拟策略[J]. 物探化探计算技术，2009，31(5)：421-425.

XU Ling-hua, TONG Xiao-zhong, LIU Jian-xin,etal. Solution strategies for 2D and 3D magnetotelluric forward modeling based on the finite element method [J]. Computing Techniques for Geophysical and Geochemical Exploration, 2009, 31 (5): 421-425. (In Chinese)

[7] MOGI T. Three-dimensional modeling of magnetotelluric data using finite- element method [J]. Journal of Applied Geophysics, 1996, 35(2): 185-189.

[8] 梁生贤，张胜业，吾守艾力,等. 复杂三维介质的大地电磁正演模拟[J]. 地球物理学进展，2012，27(5)：1981-1988.

LIANG Sheng-xian, ZHANG Sheng-ye, WU Shouaili,etal. Magnetotelluric forward modeling in complex three- dimensional media [J]. Progress in Geophys, 2012, 27(5): 1981-1988. (In Chinese)

[9] 金建铭.电磁场有限元方法[M]. 西安：西安电子科技大学出版社, 1998:164-177.

JIN Jian-ming. Finite element method of electromagnetic field [M]. Xi’an: Xi’an Electronic University Press, 1998:164-177.(In Chinese)

[10]YOSHIMURA R，OSHIMAN N. Edge-based finite element approach to the simulation of geoelectromagnetic induction in a 3-D sphere [J]. Geophysical Research Letters, 2002, 29(3): 1039-1042.

[11]MITSUHATA Y, UCHIDA T. 3D magnetotelluric modeling using theT-Ωfinite-element method [J].Geophysics, 2004, 69(1):108-119.

[12]NAM M J, KIM H J, SONG Y,etal. 3D magnetotelluric modelling including surface topography [J]. Geophysical Prospecting, 2007, 55(2): 277-287.

[13]刘长生，汤井田，任政勇,等. 基于非结构化网格的三维大地电磁自适应矢量有限元模拟[J]. 中南大学学报:自然科学版，2010，41(5)：1855-1860.

LIU Chang-sheng, TANG Jing-tian, REN Zheng-yong,etal. Three-dimension magnetotellurics modeling by adaptive edgefinite-element using unstructured meshes [J]. Journal of Central South University: Science and Technology,2010, 41(5): 1855-1860. (In Chinese)

[14]王烨. 基于矢量有限元的高频大地电磁法三维数值模拟[D]. 长沙：中南大学地球科学与信息物理学院,2008：35-62.

WANG Ye. A study of 3D high frequency magnetotelluric modeling by edge-based finite element method [D]. Changsha: Central South University. School of Geosciences and Info-Physics, 2008:35-62. (In Chinese)

[15]顾观文，吴文鹂，李桐林. 大地电磁场三维地形影响的矢量有限元数值模拟[J]. 吉林大学学报:地球科学版,2014,44(5):1678-1686.

GU Guan-wen, WU Wen-li, LI Tong-lin. Modeling for the effect of magnetotelluric 3D topography based on the vector finite-element method [J]. Journal of Jilin University: Earth Science Edition, 2014,44(5):1678-1686. (In Chinese)

[16]杨军，刘颖，吴小平. 海洋可控源电磁三维非结构矢量有限元数值模拟[J]. 地球物理学报，2015，58(8)：2827-2838.

YANG Jun, LIU Ying, WU Xiao-ping. 3D simulation of marine CSEM using vector finite element method on unstructured grids [J]. Chinese Journal of Geophysics, 2015, 58(8): 2827-2838. (In Chinese)

[17]苏晓波，李桐林，朱成,等．大地电磁三维矢量有限元正演研究[J].地球物理学进展，2015，30(4)：1772-1778.

SU Xiao-bo, LI Tong-lin, ZHU Cheng,etal. Study of three-dimensional MT forward modeling using vector finite element method [J].Progress in Geophysics, 2015,30(4):1772-1778. (In Chinese)

[18]ZHDANOV M S, VARENTSOV I M, WEAVER J T,etal. Methods for modelling electromagnetic fields Results from COMMEMI—the international project on the comparison of modelling methods for electromagnetic induction [J]. Journal of Applied Geophysics, 1997, 37(3/4):133-271.

Three-dimensional Magnetotelluric Forward Modeling Using Vector Finite Element Method

SHI Ming1,3，FENG De-shan1,2†，LI Kai-peng3，WANG Xun1,2

(1．School of Geosciences and Info-Physics，Central South Univ，Changsha, Hunan 410083, China； 2. Key Laboratory of Non-ferrous Resources and Geological Detection of Hunan Province, Changsha, Hunan 410083, China； 3. Geophysical and Geochemical Prospecting Team, Non-ferrous Metals and Nuclear Industry Geological Exploration Bureau of Guizhou Province, Duyun, Guizhou 558004, China)

Starting from the Maxwell equations, this article studied the boundary conditions of 3D MT. By using the weighted residual method, we derived the three-dimensional MT finite element equation. The three-dimensional vector finite element hexahedral meshing mode was introduced and the basis functions were selected. Then we derived the three-dimensional magnetotelluric vector finite element stiffness coefficient matrix and discrete format. A three-dimensional vector finite element magnetotelluric forward Matlab program was done. The apparent resistivity curve of the dimensional COMMEMI 3D-1 model matches the international standard test data, which proves the correctness of 3D magnetotelluric forward program. With the analysis of high and low resistivity anomalies, it shows that tensor impedance map can roughly determine the anomaly characteristics, which enriches the magnetotelluric response characteristics of expression.

vector finite element method；magnetotelluric; forward modeling; impedance tensor

1674-2974(2016)10-0119-07

2016-03-28

国家自然科学基金资助项目(41574116), National Natural Science Foundation of China(41574116)；中南大学创新驱动项目(2015CX008)；中南大学教师研究基金资助项目(2014JSJJ001); 中南大学升华育英人才计划(2012)；湖湘青年创新创业平台培养对象共同资助项目(2013)

石 明(1969-)，男，湖南益阳人，中南大学博士研究生，高级工程师

†通讯联系人，E-mail: fengdeshan@126.com

P631

A