复杂形态回线源激发电磁场的矢量有限元解

2017-12-18 10:27李建慧胡祥云刘亚军郭士明

石油地球物理勘探 2017年6期

关键词：回线有限元法脉冲响应

李建慧胡祥云陈斌刘亚军郭士明

(①中国地质大学(武汉)地球物理与空间信息学院，湖北武汉 430074；②中国矿业大学深部岩土力学与地下工程国家重点实验室，江苏徐州 221116)

·非地震·

复杂形态回线源激发电磁场的矢量有限元解

李建慧①②胡祥云*①陈斌①刘亚军①郭士明①

(①中国地质大学(武汉)地球物理与空间信息学院，湖北武汉 430074；②中国矿业大学深部岩土力学与地下工程国家重点实验室，江苏徐州 221116)

回线源瞬变电磁法实际工作中，发射回线易发生形变，这将影响资料处理解释的精度。为此，基于电磁场数值模拟的总场算法，利用矢量有限元法实现了铺设于水平地表的复杂回线源电磁场的三维正演。复杂回线源可视为由一系列首尾相接的电偶源组成，这些电偶源可进一步分解为与笛卡尔坐标系中水平坐标轴平行的一对正交电偶源。通过对比铺设于均匀半空间表面的圆形回线圆心处磁场脉冲响应的一维解析解和矢量有限元法解，验证了三维正演算法的正确性。对于复杂形态高导体地电模型，计算了矩形回线、圆形回线和普通四边形回线激发的磁场脉冲响应。结果表明，发射回线形状主要影响早期电磁响应。

复杂形态回线源矢量有限元法非结构化四面体

1 引言

目前，回线源瞬变电磁法已广泛应用于金属矿产勘探[1,2]、环境水文调查与监测[3-8]等领域。回线源瞬变电磁法，尤其是大定源装置，野外工作易受起伏地形、河流、公路等影响，导致回线源并非按照既定方案铺设，而是存在一定畸变，比如发射回线在垂直方向起伏或水平方向倾斜等。如果三维正演计算不考虑发射回线的这种畸变，将对资料处理解释造成不利影响。

瞬变电磁法的三维正演问题大致有两种策略：时步法和频谱法。时步法是以某一时刻电磁场三维分布状态作为初始值，用显式差分法或隐式差分法离散时间域电磁场微分方程，并逐步迭代直接获取时间域电磁场的一种数值方法。其中，每一时间步均需采用数值方法，比如有限差分法[9-11]、有限元法[12]等，求解一次空间域电磁场。频谱法通常由傅里叶逆变换或拉普拉斯逆变换将频率域或拉普拉斯域电磁场数值解转换至时间域，这些电磁场数值解可由积分方程法[13]、有限元法[14,15]等获得。上述研究中，发射源形态规则，比如铺设于水平地表的矩形回线或接地长导线等，并没有考虑形态复杂的发射源，比如圆形回线、普通四边形回线等。目前，有人研究了复杂形态发射源的正演计算。Yin等[16]提出了频率域航空电磁法，对发射—接收装置姿态变化引起电磁场畸变进行校正；嵇艳鞠等[17]将该校正技术应用于航空瞬变电磁法中心回线装置，为航空电磁测量数据的精确处理和解释奠定了理论基础；刘云鹤等[18]采用一维正演分析了海洋可控源电磁法中电偶源姿态变化引起的电磁场误差分布特征；韩波等[19]和严波等[20]基于电磁场数值模拟的背景场/异常场算法，分别利用有限体积法和有限差分法实现了长导线源海洋可控源电磁法的三维正演，并考虑了长导线垂向弯曲引起的电磁场畸变。

基于电磁场数值模拟的总场算法，利用矢量有限元法开展水平铺设的复杂回线源的电磁场三维正演，比如圆形回线和普通四边形回线。目前，总场算法已广泛应用于电磁场数值模拟，其不仅可避免计算异常体区域的背景电磁场，还有利于高效地模拟复杂地电模型和处理复杂发射源[21-23]。将复杂回线源视为由一系列首尾相接的水平电偶源组成，但这些电偶源的空间展布并非全都平行于笛卡尔坐标系的水平坐标轴。采用矩形块单元或规则四面体单元无法精确处理这些电偶源，因此需要采用非结构化四面体网格结合局部加密技术精确描述这些复杂场源。

2 矢量有限元法

矢量有限元法将自由度定义在棱边上，自动满足电场切向分量连续、法向分量不连续的条件，避免了节点有限元法中强加法向分量连续的条件而产生的伪解。目前，基于非结构化四面体网格的矢量有限元法已广泛应用于可控源电磁法三维正演[23-25]。

2.1 矢量有限元法方程

在准静态采用正谐时eiω t的条件下，麦克斯韦方程组中的电场E和磁场H的旋度方程为

(1)

(2)

式中：B为磁感应强度，B=μH，μ为磁导率，不考虑磁导率变化对电磁场的影响，故μ取真空中的磁导率μ0；ω为角频率，其值等于2πf，f为频率；σ为地电模型电导率；Js为场源电流密度。将式(2)带入式(1)，得到频率域电场总场双旋度方程

(3)

采用矢量有限元法离散式(3)，最终可得

(4)

2.2 加载场源

(5)

式中：δ表示狄拉克函数;I表示电流;X表示x方向的单位矢量。

图1 回线源划分为多个电偶源的示意图

(6)

式中Y表示y方向的单位矢量。

图2 电偶源ds分解示意图

3 算法验证

从矢量有限元法加载场源的角度分析，圆形回线可看作复杂回线的一个特例。本文以电导率为0.01S/m的均匀半空间为例，通过计算圆形回线圆心处磁场的一维解析解(见附录A)来验证矢量有限元解的正确性。计算环境为Ubuntu 14.04系统和GNU Fortran编译器，硬件为Intel Xeon CPU E5-2640 v3处理器和64GB内存台式机。

非结构化四面体网格剖分和显示分别采用TetGen 1.5.1-beta1[30]和ParaView 4.3.1[31]。采用TetGen网格剖分时，四面体外接圆半径与其最短棱边的最大比值设置为1.4，二面角最小值设置为16°，并局部加密测点和圆形回线附近网格。网格局部加密时，将测点和每条线段的两个端点分别作为一个正四面体中心，该四面体的四个顶点即为加密点。本例中正四面体棱边长度设置为1m。

计算区域为20km×20km×20km，其中空气层和均匀半空间厚度各为10km，空气电导率设置为10-8S/m。矢量有限元法求解时，首先将圆形回线近似等分为首尾相接的48、100或200段。局部加密每条直线段两个端点附近的网格，因此每条直线段会分布于多个四面体之中，用以近似圆形回线的电偶源数目也远多于48、100或200。表1给出了网格剖分的详细参数及相应的矢量有限元法计算时间。图3是100段近似情况下，半径100m圆形回线的非结构化网格剖分示意图。

表1 均匀半空间模型的网格剖分方案及相应的

图3 100条线段近似情况下，半径为100m的圆形回线的非结构化网格剖分示意图

图4是半径为50m的圆形回线在均匀半空间表面激发磁感应强度脉冲响应的矢量有限元解，可见矢量有限元解与一维解析解吻合良好，相对误差在4.5%以内，相对误差并没有随着直线段数目的增多而下降，并且各个时刻的数值起伏较大。图5为三个尺度的圆形回线在均匀半空间表面激发脉冲响应的矢量有限元解，这三个圆形回线均被近似为100段。结果表明：三个尺度圆形回线的矢量有限元解均具有较高精度，与一维解析解的相对误差均小于4%。

图4 对于0.01S/m的均匀半空间模型，半径a=50m的圆形回线在其圆心处激发的磁感应强度脉冲响应及相对误差

图5 半径a=50m、100m和200m的圆形回线在0.01S/m均匀半空间表面激发的磁感应强度脉冲响应及相对误差

4 算例

以复杂形态高导体模型为例，研究矩形回线、圆形回线和普通四边形回线激发的磁场脉冲响应特征。该高导体三维形态与加拿大纽芬兰和拉布拉多省的Ovoid大规模硫化物矿体一致[22]，如图6a、图6b和图6c所示；与发射回线的位置关系如图6c、图6d和图6e所示。其中，矩形回线与圆心回线面积相同，矩形回线边长为500m，圆形回线半径为282.1m；普通四边形回线的发射面积比矩形回线略小。121个测点布置于三类发射回线内外，详见图6f。地下半空间和高导体电导率分别设置为0.01S/m和1S/m，计算区域为50km×50km×50km，不同发射回线的网格剖分情况详见表2。

表2 复杂形状高导体模型的网格剖分方案

图7为测点(555800m,6243165m,110m)和(555500m,6242965m,110m)处的磁感应强度脉冲响应。前一个测点位于矩形回线中心和圆形回线圆心，后一测点位于发射回线外部(图6f)。对于这两个测点，三类发射回线激发的脉冲响应总体趋势一致，但在早期存在明显差异。对于前一个测点，普通四边形回线激发的早期脉冲响应幅值比其他两个回线更大，这是由于普通四边形回线的两条倾斜边距离测点更近。而对于后一个测点，圆形回线激发的脉冲响应由负到正的变号时刻略晚，这是由于圆形回线距离该测点更远。另外，普通四边形回线激发的脉冲响应比其他两类回线在晚期时刻的数值略小。这是由于不同形状的发射回线在发射周期的末期均可视为一个磁偶源，电磁场分布状态和幅值大小与发射回线形状几无关系，仅与发射磁矩有关。普通四边形回线的磁矩略小，因此其激发的晚期脉冲响应也略小。

图8～图10分别为矩形回线、圆形回线和普通四边形回线激发磁感应强度脉冲响应的不同时刻等值线图。在0.053ms时刻，三类发射回线在其内部的脉冲响应等值线形态与其回线形态基本一致，受高导体影响微弱；在0.356ms时刻，等值线形态受红色、白色和绿色线框分别表示矩形回线、圆形回线和普通四边形回线；两个红色实心点分别表示测点(555800m,6243165m,110m)和(555500m,6242965m,110m)，白色实心点表示其他测点回线形态影响减弱，高导体的影响凸显；在2.395ms和32.9ms时刻，等值线形态主要受高导体影响。这些现象说明：脉冲响应等值线形态在早期与发射回线形态密切相关，受高导体影响较小；而在晚期与发射回线形态几乎无关，主要受高导体影响。

图6 复杂形态高导体模型、发射回线以及测点布置示意图

(a)高导体在z-x平面的投影； (b)高导体在z-y平面的投影； (c)矩形回线与矿体在海拔110m平面的位置关系示意图； (d)圆形回线与矿体在海拔110m平面的位置关系示意图； (e)普通四边形回线与矿体在海拔110m平面的位置关系示意图； (f)三类发射回线与测点的位置关系示意图

图7 复杂形状高导体模型在不同发射回线激发时的磁感应强度脉冲响应垂直分量

图8 对于复杂形态高导体模型，矩形回线激发的不同时刻磁感应强度脉冲响应

图9 对于复杂形态高导体模型，圆形回线激发的不同时刻磁感应强度脉冲响应

图10 对于复杂形态高导体模型，普通四边形回线激发的不同时刻磁感应强度脉冲响应

5 结束语

水平铺设的复杂回线源可分解为一系列首尾相接的电偶源，并可将这些电偶源进一步分解为与笛卡尔坐标系中水平坐标轴平行的一对正交电偶源。本文以这种方式在矢量有限元法分析中加载了复杂回线源。在铺设于均匀半空间表面的圆形回线模型中，解析法和矢量有限元法分别计算的圆心处磁感应强度脉冲响应吻合良好，说明了上述场源加载方式正确可行。以复杂形态高导体模型为例，采用矢量有限元法计算了矩形回线、圆形回线和普通四边形回线激发的磁感应强度脉冲响应，并发现早期脉冲响应主要受发射回线形态的影响，而晚期脉冲响应几乎不受回线形态影响。

[1] Xue G Q,Qin K Z,Li X et al.Discovery of a large-scale porphyry molybdenum deposit in Tibet through a modified TEM exploration method.Journal of Environmental and Engineering Geophysics,2012,17(1):19-25.

[2] Yang D,Oldenburg D W.Survey decomposition:A scalable framework for 3D controlled-source electromagnetic inversion.Geophysics,2016,81(2):E31-E49.

[3] 刘树才,刘志新,姜志海.瞬变电磁法在煤矿采区水文勘探中的应用.中国矿业大学学报,2005,34(4):414-417.

Liu Shucai,Liu Zhixin,Jiang Zhihai.Application of TEM in hydrogeological prospecting of mining district.Journal of China University of Mining & Technology,2005,34(4):414-417.

[4] 薛国强,李貅,郭文波等.大回线源瞬变电磁场响应特性.石油地球物理勘探,2007,42(5):586-590.

Xue Guoqiang,Li Xiu,Guo Wenbo et al.Characters of response of large-loop transient electro-magnetic field.OGP,2007,42(5):586-590.

[5] Xue G Q,Bai C Y,Yan S et al.Deep sounding TEM investigation method based on a modified fixed central-loop system.Journal of Environmental and Engineering Geophysics,2012,76(17):23-32.

[7] 李建慧,刘树才,焦险峰等.地—井瞬变电磁法三维正演研究.石油地球物理勘探,2015,50(3):556-564.

Li Jianhui,Liu Shucai,Jiao Xianfeng et al.Three-dimensional forward modeling for surface-borehole transient electromagnetic method.OGP,2015,50(3):556-564.

[8] 胡祖志,石艳玲,何展翔等.瞬变电磁法在西部黄土层勘探中的应用.石油地球物理勘探,2016,51(增刊):131-136.

Hu Zuzhi,Shi Yanling,He Zhanxiang et al.Application of transient electromagnetic method to detect loess layer in Western China.OGP,2016,51(S):131-136.

[9] Wang T,Hohmann G W.A finite-difference,time-domain solution for three-dimensional electromagnetic modeling.Geophysics,1993,58(6):797-809.

[10] Commer M,Newman G.A parallel finite-difference approach for 3D transient electromagnetic modeling with galvanic sources.Geophysics,2004,69(5):1192-1202.

[11] 邱稚鹏,李展辉,李墩柱等.基于非正交网格的带地形三维瞬变电磁场模拟.地球物理学报,2013,56(12):4245-4255.

Qiu Zhipeng,Li Zhanhui,Li Dunzhu et al.Non-orthogonal-grid-based three dimensional modeling of transient electromagnetic field with topography.Chinese Journal of Geophysics,2013,56(12):4245-4255.

[12] Um E S,Harris J M,Alumbaugh D L.3D time-domain simulation of electromagnetic diffusion phenomena:A finite-element electric-field approach.Geophysics,2010,75(4):F115-F126.

[13] Newman G A,Hohmann G W,Anderson W L.Transient electromagnetic response of a three-dimensional body in a layered earth.Geophysics,1986,51(8):1608-1627.

[14] Sugeng F.Modelling the 3D TDEM response using the 3D full-domain finite-element method based on the hexahedral edge-element technique.Exploration Geophysics,1998,29(4):615-619.

[15] Li J H,Zhu Z Q,Liu S C et al.3D numerical simulation for the transient electromagnetic field excited by the central loop based on the vector finite-element method.Journal of Geophysics and Engineering,2011,8(4):560-567.

[16] Yin C C,Fraser D C.Attitude corrections of helicopter EM data using a superposed dipole model.Geophysics,2004,69(2):431-439.

[17] 嵇艳鞠,林君,关珊珊等.直升机航空TEM中心回线线圈姿态校正的理论研究.地球物理学报,2010,53(1):171-176.

Ji Yanju,Lin Jun,Guan Shanshan et al.Theoretical study of concentric loop coils attitude correction in helicopter-borne TEM.Chinese Journal of Geophysics,2010,53(1):171-176.

[18] 刘云鹤,殷长春,翁爱华等.海洋可控源电磁法发射源姿态影响研究.地球物理学报,2012,55(8):2757-2768.

Liu Yunhe,Yin Changchun,Weng Aihua et al.Attitude effect for marine CSEM system.Chinese Journal of Geophysics,2012,55(8):2757-2768.

[19] 韩波,胡祥云,Schultz A等.复杂场源形态的海洋可控源电磁三维正演.地球物理学报,2015,58(3):1059-1071.

Han Bo,Hu Xiangyun,Schultz A et al.Three-dimensional forward modeling of the marine controlled-source electromagnetic field with complex source geo-metries.Chinese Journal of Geophysics,2015,58(3):1059-1071.

[20] 严波,李予国,韩波等.任意方位电偶源的MCSEM电磁场三维正演.石油地球物理勘探,2017,52(4):859-868.

Yan Bo,Li Yuguo,Han Bo et al.3D marine controlled-source electromagnetic forward modeling with arbitrarily oriented dipole source.OGP,2017,52(4):859-868.

[21] Weiss C J.Project AphiD:A Lorenz-gaugedA-Φde-composition for parallelized computation of ultra-broadband electromagnetic induction in a fully heterogeneous earth.Computer & Geosciences,2013,58:40-52.

[22] Jahandari H,Farquharson C G.A finite-volume solution to the geophysical electromagnetic forward problem using unstructured grids.Geophysics,2014,79(6):E287-E302.

[23] Ansari S,Farquharson C G.3D finite-element forward modeling of electromagnetic data using vector and scalar potentials and unstructured grids.Geophysics,2014,79(4):E149-E165.

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

Yang Jun,Liu Ying,Wu Xiaoping.3D simulation of marine CSEM using vector finite element method on unstructured grids.Chinese Journal of Geophysics,2015,58(8):2827-2838.

[25] 李建慧,Colin G Farquharson,胡祥云等.基于电场总场矢量有限元法的接地长导线源三维正演.地球物理学报,2016,59(4):1521-1534.

Li Jianhui,Farquharson C G,Hu Xiangyun et al.A vector finite element solver of three-dimensional mode-lling for a long grounded wire source based on total electric field.Chinese Journal of Geophysics,2016,59(4):1521-1534.

[26] Jin J M.The Finite Element Method in Electromagnetic II.John Wiley & Sons Inc,New York，2002.

[27] Amestoy P R,Guermouche A,L’Excellent J Y et al.Hybrid scheduling for the parallel solution of linear systems.Parallel Computing,2006,32(2):136-156.

[28] Anderson W L.Fourier cosine and sine transforms using lagged convolutions in double-precision (subprograms DLAGF0/DLAGF1).Open-File Report 83-320,U.S.Geological Survey,1983.

[29] Li J,Farquharson C G,Hu X.Three effective inverse Laplace transform algorithms for computing time-domain electromagnetic responses.Geophysics,2016,81(2):E113-E128.

[30] Si H.TetGen, a delaunay-based quality tetrahedral mesh generator.ACM Transactions on Mathematical Software,2015,41(2)：11.

[31] Ayachit U.The ParaView Guide:A Parallel Visualization Application.Kitware Inc,2015.

[32] Ward S H,Hohmann G W.Electromagnetic theory for geophysical applications // Electromagnetic Me-thods in Applied Geophysics:Volume 1,Theory.SEG,1988,167-183.

[33] Guptasarma D,Singh B.New digital linear filters for HankelJ0andJ1transforms.Geophysical Prospecting,1997,45(5):745-762.

附录A圆形回线在圆心处激发的解析解

采用正谐时eiω μ条件下，铺设于水平地表的圆形发射回线在其圆心处激发的磁场垂直分量为[32]