2016, Vol. 59
2. 纽芬兰纪念大学地球科学系, 圣约翰斯, 加拿大;
3. 中国矿业大学深部岩土力学与地下工程国家重点实验室, 江苏徐州 221116
2. Department of Earth Sciences, Memorial University of Newfoundland, St. John's, NL, A1B;
3. X5, Canada 3 State Key Laboratory for Geomechanics and Deep Underground Engineering, China University of Mining & Technology, Jiangsu Xuzhou 221116, China
The numerical accuracy of FE method depends on meshing quality to some extent, which is evaluated by the maximum allowable radius-edge ratio and the minimum allowable dihedral angle for the unstructured tetrahedral elements generated by the TetGen software. Taking the electric field excited by a HED within a 0.01 S·m-1 whole-space as examples, the results show that as the maximum allowable radius-edge ratio decreases or the minimum allowable dihedral angle increases, the numerical accuracy will be improved. By this model, we also validated the algorithm presented here.
For a conductive prism buried in a homogeneous half-space with a 100m-long grounded wire, the electric field calculated by our algorithm was compared with those calculated by the integral equation method based on secondary electric field, the finite volume method based on total electric field and the FE method based on magnetic vector potential. The results show that these four numerical solutions coincide well with each other, and the behavior of the electric field well indicates the conductive prism. Then, we applied the algorithm presented here to compute the amplitude and phase of the electric field for the model of disc-shape hydrocarbon buried in seabed. Through this model, the validity and ability of modelling the electromagnetic field for irregular bodies were tested simultaneously. Finally, an inclined fault, contact zone and metalliferous vein are always approximated by an inclined plate in three-dimensional modelling. For the inclined plate model with different conductivities, the electric field excited by a 1000m-long grounded wire source was calculated. The electric field for the inclined plate of 0.0333 S·m-1 has stronger anomaly responses than those for the inclined plate of 0.0167 S·m-1.
In the following study, this algorithm will be applied to three-dimensional modelling for a long grounded wire source laid on a rugged surface.
电磁法勘探中，电性源装置已成功应用于金属矿产勘探(Hu et al.，2013)、环境水文调查与监测(Fu et al.，2013；Grayver et al.，2014)、海底油气藏勘探等领域(Zhdanov et al.，2014).为了更加精细地反映地电结构、进一步改善勘探效果，开发精度高、速度快、能够模拟复杂地电模型的三维正演技术已成为必由之路.
电性源装置具体包括电偶源和接地长导线源.除一维解析法外，电性源的正演按照维度还可分为二维、2.5维和三维.二维正演问题是指地电模型为二维、接地长导线源沿地电模型走向无限延伸的情况，即二维模型、二维场源问题.阎述和陈明生(1999a，1999b)发展和完善了无限长接地导线的有限单元法二维正演；陈小斌和胡文宝(2002)利用有限元法二维正演分析了频率域电磁测深近区观测的可行性；王若等(2006)开展有限元法正演时，在计算区域外边界加载了一阶吸收边界条件，并将伪δ函数引入场源项，提高了近区数值解精度.
2.5维正演问题是指地电模型为二维、电性源为电偶源或接地长导线源的情况，即二维模型、三维场源问题.开展2.5维正演时，首先利用傅里叶变换将与地质构造走向一致的电磁场分量转换至波数域，再利用有限元法等数值方法求取一系列波数的电磁场响应；最后利用傅里叶逆变换将波数域电磁场响应转换至空间域，进而获得二维模型的三维电磁场响应.基于电磁场数值模拟的背景场/异常场算法，Unsworth等(1993)利用有限元法实现了水平电偶源的频率域电磁法2.5维正演问题，并有效克服了场源的奇异性问题；底青云等(2004)将该算法应用于复杂介质的正演计算，为反演解释提供了良好依据.Li和Key(2007)基于背景场/异常场算法，利用自适应有限元法实现了海洋可控源电磁法的2.5维正演问题，场源为水平电偶源.2.5维正演问题的另一种策略是采用电磁场数值模拟的总场算法，如何克服场源奇异性对数值解的不利影响是该策略的关键问题.Mitsuhata(2000)以伪δ函数代替δ函数处理场源项、克服了场源奇异性，实现了基于有限元法的水平电偶源频率域电磁法2.5维正演；Mitsuhata等(2002)延续上述思路，开展了接地长导线源频率域电磁法的2.5维正演和反演研究.
三维正演问题是指地电模型为三维、电性源为电偶源或接地长导线源的情况，即三维模型、三维场源问题.基于电磁场数值模拟的背景场/异常场算法，国内外许多学者开展了电性源频率域电磁法的三维正演问题，Farquharson和Oldenburg(2002)采用了积分方程法；Weiss和Constable(2006)、韩波等(2015)采用了有限体积法；Sasaki和Meju(2009)、Streich(2009)采用了交错网格有限差分法；Schwarzbach等(2011)、李勇等(2015)采用了矢量有限元法.采用电磁场数值模拟总场算法是电性源频率域电磁法三维正演的另一有效策略.阎述和陈明生(2000)对利用矢量有限元法开展电偶源频率域电磁法三维正演进行了探索研究.张继锋等(2009)基于电场矢量波动方程，利用节点有限元法实现了水平电偶源的频率域电磁法三维正演，以伪δ函数代替δ函数处理场源项，改善了近区数值解精度；并以水平电偶源的一维解析解作为计算区域外边界条件.徐志锋和吴小平(2010)基于磁矢量势Helmholtz方程，采用稳定化节点有限元法实现了水平电偶源的频率域电磁法三维正演，也以水平电偶源的一维解析解作为计算区域外边界条件.Ansari和Farquharson(2014)基于磁矢量势Helmholtz方程，采用基于非结构化四面体的矢量有限元法实现了水平电偶源和接地长导线源的频率域电磁法三维正演，将接地长导线看作有限多个水平电偶源的组合，消除了场源奇异性对近区数值解的不利影响.Jahandari和Farquharson(2014)基于电场总场Helmholtz方程，采用基于非结构化四面体和Voronoï网格的有限体积法实现了水平电偶源和接地长导线源的频率域电磁法三维正演，也将接地长导线看作有限多个水平电偶源的组合.
以接地长导线源作为发射装置的电磁法有广域电磁法(何继善，2010)、短偏移距瞬变电磁法(薛国强等，2014)、海洋可控源电磁法(韩波等，2015)，这些方法的共同点为测点不仅位于远区，而且还分布于过渡区和近区.如果仍将发射源看作水平电偶源开展数值模拟，将不利于获得真实的电磁场响应.本研究将接地长导线源看作有限多个水平电偶源的组合，基于电场总场Helmholtz方程，采用基于非结构化四面体网格剖分的矢量有限元法开展频率域电磁法三维正演研究.矢量有限元法将自由度定义在棱边上，自动满足电场切向分量连续、法向分量不连续的连续条件，因此可避免节点有限元法中强加法向分量连续而产生的伪解.目前，该方法已广泛应用于电磁法勘探的三维正演研究中(Schwarzbach et al.，2011；Li et al.，2011；Ansari and Farquharson，2014；李勇等，2015).非结构化四面体网格剖分不仅使得本文算法能够模拟更加复杂的地电模型，而且能够通过采取局部加密的方法，细化场源、测点和异常体所在区域网格剖分，进而提高数值解精度.
2 矢量有限单元法 2.1 控制方程在准静态、采用正谐时eiωt的条件下，麦克斯韦方程组中的电场E和磁场H旋度方程为：
以式(3)为控制方程，应用Galerkin法推导有限元方程.设余量r为
本研究中，将计算区域剖分为有限多个四面体单元e.每个四面体e中，点(x，y，z)处的电场矢量场可由各条棱边的电场和矢量基函数表示，有
 | 图 1 四面体单元中，节点与棱边编号 Fig. 1 The numbering diagram of the edges and nodes in a tetrahedral element |
当场源为接地长导线源时，为了消除场源对近区数值解的不利影响，需将接地长导线划分为有限多段，每段均可视为一个电偶源.假设沿x方向布置的接地长导线被划分为M段，且分布于M个单元中；第m个网格单元中的发射源x方向长度为dsmx，坐标为xm，而M段导线的y和z坐标分别为y0和z0，且保持不变，那么电流密度式(4)可改写为：
与迭代法相比，直接法求解大型方程组具有精度高、稳定性好的优点，但其内存需求大、求解效率低.近年来计算机硬件和并行计算的快速发展，直接法效率已有明显提升，并已大量应用于电磁法三维数值模拟(Streich，2009；Schwarzbach et al.，2011； Jahandari and Farquharson，2014).本研究将采用直接法开源软件MUMPS(Amestoy et al.，2006)求解大型线性方程组；由于矢量有限元法生成的系数矩阵具有对称性，因此只需存储和输入下三角矩阵.获得各棱边电场值后，即可通过(9)式计算单元内任意点的电场矢量场.
3 网格剖分与算法验证本研究采用的开源代码或程序如下：求解方程组采用MUMPS 5.0.0(Amestoy et al.，2006)、非结构化四面体网格剖分采用TetGen 1.5.1-beta1(Si，2015)，而非结构化四面体网格显示采用ParaView 4.3.1(Ayachit，2015).另外，计算环境为：采用Ubuntu 14.04系统和GNU Fortran编译器，硬件为Inter I5处理器、16 G内存台式机.
有限元法数值解精度一定程度上取决于网格剖分质量.在TetGen软件中，非结构化四面体形状和尺度由四面体外接圆半径与其最短棱边比值和四面体二面角(如图 1所示)角度控制(Si，2013，2015).Ansari和Farquharson(2014)采用外接圆半径与四面体最短棱边最大比值为1.414、二面角最小值为16°的网格剖分方案，使用矢量有限元法取得了理想的数值结果；而 Jahandari和Farquharson(2014)将四面体外接圆半径与其最短棱边的最大比值设置在1.1~1.2之间，并未限制四面体二面角角度，采用有限体积法也取得了高精度的数值解.
本研究将以置于0.01 S·m-1均匀全空间中的水平电偶源为例，分析四面体外接圆半径与最短棱边比值和二面角角度对有限元法数值解精度的影响.为了避免计算区域尺度不足引起的数值解误差，本例计算区域尺度为400 km ×400 km ×400 km.表 1为网格剖分方案技术参数，四面体外接圆半径与其最短棱边的最大比值(R)分别为1.2、1.3和1.4，二面角最小值(DA)分别为未限制、12°、14°、16°和18°.图 2为R值等于1.4、DA值未限制条件下的非结构化网格剖分示意图，其中图 2b为图 2a中白色线条所围区域的放大图.非结构化四面体网格剖分时，采用局部加密策略对电偶源和测点分布区域进行加密.加密方案具体如下：将电偶源或测点作为一个正四面体中心，该四面体的四个顶点即为加密点.本例中，正四面体棱边长度为1 m.
 | 表 1 基于非结构化四面体的网格剖分方案 Table 1 The meshing schemes of unstructured tetrahedrons |
 |
图 2 均匀全空间模型的非结构化网格示意图(y=0 m) 图(b)为图(a)中白色线条所围区域的放大图，其范围为16 km×16 km. Fig. 2 The section map (y=0 m) of tetrahedral meshing for the whole-space model Panel (b) is the enlarged view of the domain circled by the white lines in panel (a), and it is 16 km×16 km. |
在准静态、采用正谐时eiωt的条件下，沿x方向布置的电偶源在均匀全空间中激发的电场Ex分量解析解为(Ward and Hohmann，1988)：
图 3为电场Ex分量的有限元法数值解及其与解析解的相对误差曲线，图中黑色虚线表示3%的误差限.如图 3a所示：采用3种网格剖分方案，有限元法计算的数值解与解析解吻合良好；仅收发距为9000 m和10000 m，且R值为1.4的数值解精度略差.图 3b至3f展现了电场实部和虚部相对误差随二面角最小值变化的规律.不难发现实部相对误差均大于虚部相对误差，因此本文采用实部相对误差作为分析指标.综合分析图 3b、3c和3d，我们发现：随着收发距的增大，实部与虚部相对误差总体呈现增大趋势.当R值为1.4时，仅有收发距小于1000 m的少数测点实部相对误差小于3%；当R值为1.3时，收发距小于6000 m的大多数测点实部相对误差小于3%；当R值为1.2时，仅有收发距大于6000 m的少数测点实部相对误差大于3%.图 3e和3f中，仅有收发距大于7000 m的少数实部相对误差大于3%.如图 3b至3f所示：DA值固定时，电场实部与虚部相对误差随着R值减小而呈现出递减的变化规律.另外如表 1所示：DA值为18°的网格剖分过密，单元数与棱边数激增.比如R值为1.2、DA值为18°的网格棱边数高达74095813.这一网格剖分规模远大于本研究所采用硬件的允许计算规模(经测试约为1800000条棱边)，因此并未计算该网格剖分的有限元法数值解.
 |
图 3 置于0.01 S·m-1均匀全空间中的电偶源，频率为1 Hz的电场Ex分量数值解及其相对误差 紫色和绿色符号分别表示实部与虚部分量.图中R表示四面体外接圆半径与最短棱边的最大比值，而DA表示四面体二面角的最小值. Fig. 3 The numerical solutions and their relative errors for the Ex component excited by a horizontal electric dipole located in the 0.01 S·m-1 homogeneous whole-space for a frequency of 1 Hz The purple and green symbols represent the real and imaginary parts, respectively. The R and DA denote the maximum tetrahedral radius-edge ratio and the minimum dihedral angle. |
结合表 1和图 3综合分析，随着R值的减小或DA值的增大，非结构化网格剖分的四面体和棱边数目呈递增趋势、剖分网格逐渐精细，有限元法数值解精度也随之提高.综合考虑计算效率和精度，本文后续研究均采用R值为1.3和DA值为16°的网格剖分参数.
4 算例 4.1 块状高导体模型如图 4所示：一个电导率为0.2 S·m-1的块状体埋置于电导率为0.02 S·m-1的均匀半空间中，其三维尺度(依次为x、y、z方向)为120 m×200 m×400 m，中心位于(1000，0，300)m；空气电导率设置为10-8 S·m-1.测线布设于地表，起始位置为(400，0，0)m，终止位置为(1600，0，0)m.发射源是长度为100 m的接地长导线，其中心坐标为(50，0，0)m.发射电流为1 A、发射频率为3 Hz.本例计算区域为50 km×50 km×50 km，共剖分325207个四面体、378508条棱边.为了改善数值解精度，采用局部加密技术对接地长导线和测点所在区域网格进行了细分，并限制块状高导体所在区域的剖分单元最大体积为5000 m3(如图 5所示).另外，本例计算时间约为600 s.
 |
图 4 块状高导体模型示意图 (a) x-z平面断面图(y=0 m)； (b) x-y平面俯视图. Fig. 4 The diagrams for the conductive prism (a) The section view (y=0 m) in the x-z plane; (b) The top view in the x-y plane. |
 | 图 5 块状高导体模型的非结构化网格剖分断面图(y=0 m) Fig. 5 The section view (y=0 m) of the unstructured meshing for the conductive prism |
图 6为块状高导体的电场分量Ex响应曲线，其中实部为正值、虚部为负值.为了进一步验证本文算法的正确性，将本研究数值解与积分方程法(IE)(Farquharson and Oldenburg，2002)、有限体积法(FV)( Jahandari and Farquharson，2014)和基于磁矢量势Helmholtz方程的有限元法(FE)(Ansari and Farquharson，2014)数值解进行了对比.由图 6可看出：除积分方程法外，其他三种方法的数值解实部与虚部曲线形态一致、吻合程度良好；与其他三种方法相比，积分方程法数值解在块状体分布区域(x值为900至1100 m)偏大.积分方程法采用矩形块网格剖分方案，高导块状体沿x、y和z方向剖分的网格数目依次为8×8×8，被剖分为512个矩形块单元；其他三种方法均采用了非结构化四面体单元，比如本文算法采用TetGen软件将块状体剖分为7311个四面体单元.因此，我们推断积分方程法数值解与其他三种方法数值解存在上述差异的潜在原因为积分方程法求解时，块状高导体的剖分单元数目不足，影响了数值解精度.
 | 图 6 均匀半空间中块状高导体地电模型的电场总场Ex分量 Fig. 6 The Ex component of total electric field for the conductive prism in the homogeneous half-space |
另外，这些数值方法计算的电场响应对高导块状体反映显著，电场实部和虚部曲线在x为900至1100 m范围内均呈现向下凹陷现象.这是由于地下电流易被高导体“吸引”，流经地表的电流密度降低所致.
4.2 海底油气藏模型海底油气藏资源是近十年来电磁法勘探的热门领域，本例模型与典型油气藏模型(Weiss and Constable，2006)类似.如图 7所示，该地电模型包括电导率为3.2 S·m-1的海水层、1 S·m-1的海底以及0.01 S·m-1的油气藏.油气藏分布呈扁平圆柱体状，半径为1000 m，厚度为100 m，顶部埋深950 m.长度为100 m的接地长导线源分布于海水层，距离海底高度为100 m，其中心坐标为(-50，0，-100)m.发射电流为1 A，发射频率为1 Hz.计算区域为30 km×30 km×16 km，其中海水层厚度为1 km，该模型共剖分1342046个四面体以及1562622条棱边.如图 8a所示，我们依然对接地长导线、测点和电性异常体所在区域进行了局部加密，而图 8b说明了采用TetGen软件能够精细地剖分扁平圆柱体状油气藏资源，同时也避免了采用结构化网格易导致计算区域边界处出现扁平单元和狭长单元的现象.另外，本例计算时间约为3000 s.
 |
图 7 海底油气藏模型示意图 (a) x-z平面断面图(y=0 m)； (b) x-y平面俯视图. Fig. 7 The diagrams for the hydrocarbon buried in seabed (a) The section view (y=0 m) in the x-z plane; (b) The top view in the x-y plane. |
 |
图 8 海底油气藏模型的非结构化网格剖分结果 (a) x-z平面断面图(y=0 m)； (b) x-y平面切片图(z=1000 m). Fig. 8 The results of the unstructured meshing for the hydrocarbon buried in seabed (a) The section view (y=0 m) in the x-z plane; (b) The slice (z=1000 m) in the x-y plane. |
图 9为海底油气藏地电模型的电场分量Ex实部与虚部、振幅和相位曲线.由图 9可知，电场实部与虚部数量级从10-5降低到10-15，并存在多个极小值，因此对数值方法精度和稳定性要求极高.如图 9a所示，本文算法计算的电场实部和虚部与基于磁矢量势Helmholtz方程的有限元法(Ansari and Farquharson，2014)计算结果吻合程度良好.两种数值算法的主要区别为极小值的数值不同，即本研究计算结果偏大，这种现象由采样不足引起.本研究测点间距为50 m，而Ansari与Farquharson(2014)的测点间距为5 m，其更能反映真实场值.频率域海洋可控源电磁法中，通常结合航行数据，将实测信号转换为振幅-收发距曲线与相位-收发距曲线(李予国和段双敏，2014).因此，本研究将电场实部与虚部转换为了电场振幅与相位曲线.由图 9b和9c可知两种数值方法的电场振幅与相位完全吻合，可见采样间隔的这种差异对电场振幅与相位影响微小.
 |
图 9 海底油气藏模型的电场总场Ex分量 “FE”表示Ansari与Farquharson(2014)的有限元法数值解.(a) 实部与虚部；(b) 振幅；(c) 相位. Fig. 9 The Ex component of total electric field for the hydrocarbon buried in seabed The symbol of “FE” denotes the solution of the finite element method(Ansari and Farquharson,2014). (a) Real and imaginary parts; (b) Amplitude; (c) Phase. |
断层破碎带、接触带和金属矿脉往往可由倾斜板状体予以近似.如图 10所示：一顶部出露于地表的倾斜板状体分布于电导率为0.01 S·m-1的均匀半空间，其三维尺度为200 m×3000 m×1500 m，电导率分别为0.0333 S·m-1和0.0167 S·m-1；空气电导率设置为10-8 S·m-1.板状体沿走向方向(y方向)分布范围为[1500，4500] m；沿x方向分布范围： 顶部为[50，250] m，底部为[-250，-50] m.
 |
图 10 倾斜板状体模型示意图 (a) x-z平面断面图(y=3000 m)； (b) x-y平面俯视图. Fig. 10 The diagrams for the inclined plate (a) The section view (y=3000 m) in the x-z plane; (b) The top view in the x-y plane. |
本例在地表(z=0 m)布置了y坐标分别为1000、2000、3000、4000和5000 m的5条测线、共105个测点，其中每条测线长度为2000 m、含21个测点.布置于x方向的接地长导线源长度为1000 m，中心置于(0，0，0)m.发射电流为1 A、发射频率为1 Hz和32 Hz.本例计算区域为50 km×50 km×50 km，共剖分1518786个四面体、1766490条棱边.该地电模型的非结构化四面体网格剖分断面图(y=3000 m)如图 11所示.
 | 图 11 倾斜板状体的非结构化网格剖分断面图(y=3000 m) Fig. 11 The section view (y=3000 m) of the unstructured meshing for the inclined plate |
图 12为倾斜板状体模型的电场Ex分量振幅曲线.由图 12所示，两个频率的电场振幅曲线随测点x坐标变化的形态一致；相比于均匀半空间模型，电场振幅对倾斜高导板状体模型反映显著，其中0.0333 S·m-1模型的电场振幅比0.0167 S·m-1情形的数值更大，尤其是x坐标从-100 m和400 m的六个测点.两个频率电场振幅的主要区别在于远离测线中心的测点对高导板状体的反映能力.当x坐标为-1000~-200 m和500~1000 m时，倾斜板状体1 Hz的电场振幅略大于均匀半空间相应频率的电场振幅；而32 Hz时，两者电场振幅近似相等.
 | 图 12 倾斜板状体模型的电场总场Ex分量振幅(y=3 km) Fig. 12 The amplitude of the Ex component for the inclined plate (y=3 km) |
图 13和14为倾斜板状体模型的电场分量Ex振幅等值线图，其共同点是电场振幅随着y坐标的增大而逐渐衰减.图 13中1 Hz的电场振幅从1.2×10-5 V·m-1衰减至9.0×10-8 V·m-1，而图 14中32 Hz的电场振幅从1.9×10-5 V·m-1衰减至1.0×10-7 V·m-1.这一现象反映了高频电磁场随着传播距离的增大而衰减更快的客观规律.均匀半空间模型中，电场振幅等值线关于y轴对称；倾斜板状体模型中，电场振幅等值线并不对称.如图 13b、13c、14b和14c所示：由于受高导板状体影响，x方向[0，300]m、y方向[1800，4500]m区域内的等值线向下凹陷；与0.0167 S·m-1的板状体模型相比，0.0333 S·m-1板状体模型的等值线呈现出更加显著的向下凹陷现象.
 |
图 13 倾斜板状体模型的电场总场Ex分量振幅(1 Hz) (a) 均匀半空间； (b) 0.0167 S·m-1的板状体； (c) 0.0333 S·m-1的板状体. Fig. 13 The amplitude of the Ex component for the inclined plate (1 Hz) (a) The homogeneous half-space; (b) The plate of 0.0167 S·m-1; (c) The plate of 0.0333S·m-1. |
 |
图 14 倾斜板状体模型的电场总场Ex分量振幅(32 Hz) (a) 均匀半空间； (b) 0.0167 S·m-1的板状体； (c) 0.0333 S·m-1的板状体. Fig. 14 The amplitude of the Ex component for the inclined plate (32 Hz) (a) The homogeneous half-space; (b) The plate of 0.0167 S·m-1; (c) The plate of 0.0333 S·m-1. |
本研究将接地长导线源视为有限多个水平电偶源的组合，并以水平电偶源的形式实现了场源的分段加载，有效去除了电磁法三维正演中外加场源对近区数值解精度的不利影响.通过将本文算法数值解与均匀全空间中水平电偶源产生的电场解析解对比研究，分析得出非结构化网格剖分中，四面体外接圆半径与其最短棱边最大比值的减小或者四面体二面角最小值的增大都能提高有限元法数值解精度.然后，将本文算法数值解与块状高导体地电模型的积分方程法、有限体积法和基于磁矢量势Helmholtz方程的有限元法数值解对比分析，四种数值解吻合良好、对高导异常体反映显著，进一步验证了本文算法的正确性.本例中，基于矩形块网格剖分的积分方程法在异常体区域的数值解比基于非结构化四面体网格剖分的三种方法数值解偏大，说明了非结构化四面体网格能够更加精细地剖分电性异常体，有利于获得精确数值解.海底油气藏和倾斜板状体地电模型主要说明了基于非结构化四面体网格剖分的数值算法能够精确刻画复杂地质体，避免了结构化网格剖分复杂地质体易导致的扁平单元和狭长单元.后续研究中，我们将致力于起伏地表接地长导线源电磁法的三维正演算法.
致谢 感谢纽芬兰纪念大学地球科学系SeyedMasoud Ansari与Hormoz Jahandari博士在三维正演程序编写中的热情帮助.感谢国家留学基金委全额资助本文第一作者在纽芬兰纪念大学开展为期1年的博士后研究.特别感谢两位审稿专家的中肯建议，使得本文质量显著提升.| [1] | Amestoy P R, Guermouche A, L'Excellent J -Y, et al. 2006. Hybrid scheduling for the parallel solution of linear systems. Parallel Computing, 32(2): 136-156. |
| [2] | Ansari S, Farquharson C G. 2014. 3D finite-element forward modeling of electromagnetic data using vector and scalar potentials and unstructured grids. Geophysics, 79(4): E149-E165. |
| [3] | Ayachit U. 2015. The ParaView Guide: A Parallel Visualization Application. ParaView 4.3 ed. Carrboro, NC: Kitware, Inc. |
| [4] | Chen X B, Hu W B. 2002. Direct iterative finite element (DIFE) algorithm and its application to electromagnetic response modeling of line current source. Chinese J. Geophys. (in Chinese), 45(1): 119-130, doi: 10.3321/j.issn:0001-5733.2002.01.015. |
| [5] | Di Q Y, Unsworth M, Wang M Y. 2004. 2.5-D CSAMT modeling with the finite element method over 2-D complex earth media. Chinese J. Geophys. (in Chinese), 47(4): 723-730, doi: 10.3321/j.issn:0001-5733.2004.04.026. |
| [6] | Farquharson C G, Oldenburg D W. 2002. An integral equation solution to the geophysical electromagnetic forward-modelling problem. // Zhdanov M S, Wannamaker P E eds. Three-Dimensional Electromagnetics. Amsterdam: Elsevier, Inc., 3-19. |
| [7] | Fu C M, Di Q Y, An Z G. 2013. Application of the CSAMT method to groundwater exploration in a metropolitan environment. Geophysics, 78(5): B201-B209. |
| [8] | Grayver A V, Streich R, Ritter O. 2014. 3D inversion and resolution analysis of land-based CSEM data from the Ketzin CO2 storage formation. Geophysics, 79(2): E101-E114. |
| [9] | Han B, Hu X Y, Schultz A, et al. 2015. Three-dimensional forward modeling of the marine controlled-source electromagnetic field with complex source geometries. Chinese J. Geophys. (in Chinese), 58(3): 1059-1071, doi: 10.6038/cjg20150330. |
| [10] | He J S. 2010. Wide field electromagnetic sounding methods. Journal of Central South University (Science and Technology) (in Chinese), 41(3): 1065-1072. |
| [11] | Hu X Y, Peng R H, Wu G J, et al. 2013. Mineral Exploration using CSAMT data: Application to Longmen region metallogenic belt, Guangdong Province, China. Geophysics, 78(3): B111-B119. |
| [12] | Jahandari H, Farquharson C G. 2014. A finite-volume solution to the geophysical electromagnetic forward problem using unstructured grids. Geophysics, 79(6): E287-E302. |
| [13] | Jin J M. 2002. The Finite Element Method in Electromagnetics. 2nd ed. New York: John Wiley & Sons, Inc. |
| [14] | Li J H, Zhu Z Q, Liu S C, et al. 2011. 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, 8(4): 560-567. |
| [15] | Li Y, Wu X P, Lin P R. 2015. Three-dimensional controlled source electromagnetic finite element simulation using the secondary field for continuous variation of electrical conductivity within each block. Chinese J. Geophys. (in Chinese), 58(3): 1072-1087, doi: 10.3038/cjg20150331. |
| [16] | Li Y G, Duan S M. 2014. Data processing of marine controlled-source electromagnetic data. Periodical of Ocean University of China (in Chinese), 44(10): 106-112. |
| [17] | Li Y G, Key K. 2007. 2D marine controlled-source electromagnetic modeling: Part 1 — An adaptive finite-element algorithm. Geophysics, 72(2): WA51-WA62. |
| [18] | Mitsuhata Y. 2000. 2D electromagnetic modeling by finite-element method with a dipole source and topography. Geophysics, 65(2): 465-475. |
| [19] | Mitsuhata Y, Uchida T, Amano H. 2002. 2.5-D inversion of frequency-domain electromagnetic data generated by a grounded-wire source. Geophysics, 67(6): 1753-1768. |
| [20] | Sasaki Y, Meju M A. 2009. Useful characteristics of shallow and deep marine CSEM responses inferred from 3D finite-difference modeling. Geophysics, 74(5): F67-F76. |
| [21] | Schwarzbach C, Börner R U, Spitzer K. 2011. Three-dimensional adaptive higher order finite element simulation for geo-electromagnetics—a marine CSEM example. Geophysical Journal International, 187(1): 63-74. |
| [22] | Si H. 2013. TetGen, A quality tetrahedral mesh generator and 3D Delaunay triangulator, Version 1.5 User's Manual. Technical Report 13, Weierstrass Institute for Applied Analysis and Stochastics. |
| [23] | Si H. 2015. TetGen, a Delaunay-based quality tetrahedral mesh generator. ACM Transactions on Mathematical Software, 41(2), doi: 10.1145/2629697. |
| [24] | Streich R. 2009. 3D finite-difference frequency-domain modeling of controlled-source electromagnetic data: Direct solution and optimization for high accuracy. Geophysics, 74(5): F95-F105. |
| [25] | Unsworth M J, Travis B J, Chave A D. 1993. Electromagnetic induction by a finite electric dipole source over a 2-D earth. Geophysics, 58(2): 198-214. |
| [26] | Wang R, Wang M Y, Di Q Y. 2006. Electromagnetic modeling due to line source in frequency domain using finite element method. Chinese J. Geophys. (in Chinese), 49(6): 1858-1866, doi: 10.3321/j.issn:0001-5733.2006.06.035. |
| [27] | Ward S H, Hohmann G W. 1988. Electromagnetic theory for geophysical applications. // Nabighian M N ed. Electromagnetic Methods in Applied Geophysics: Volume 1, Theory.Tulsa, OK: Society of Exploration Geophysicists, 167-183. |
| [28] | Weiss C J, Constable S. 2006. Mapping thin resistors and hydrocarbons with marine EM methods, Part II — Modeling and analysis in 3D. Geophysics, 71(6): G321-G332. |
| [29] | Xu Z F, Wu X P. 2010. Controlled source electromagnetic 3-D modeling in frequency domain by finite element method. Chinese J. Geophys. (in Chinese), 53(8): 1931-1939, doi: 10.3969/j.issn.0001-5733.2010.08.019. |
| [30] | Xue G Q, Yan S, Chen W Y. 2014. Research prospect to grounded-wire TEM with short-offset. Progress in Geophysics (in Chinese), 29(1): 177-181, doi: 10.6038/pg20140124. |
| [31] | Yan S, Chen M S. 1999a. Line source frequency electromagnetic sounding two-dimensional forward modeling I. Coal Geology & Exploration (in Chinese), 27(5): 60-62. |
| [32] | Yan S, Chen M S. 1999b. Line source frequency electromagnetic sounding two-dimensional forward modeling II. Coal Geology & Exploration (in Chinese), 27(6): 56-59. |
| [33] | Yan S, Chen M S. 2000. Finite element solution of three-dimensional geoelectric models in frequency electromagnetic sounding excited by a horizontal electric dipole. Coal Geology & Exploration (in Chinese), 28(3): 50-56. |
| [34] | Zhang J F, Tang J T, Yu Y, et al. 2009. Three dimensional controlled source electromagnetic numerical simulation based on electric field vector wave equation using finite element method. Chinese J. Geophys. (in Chinese), 52(12): 3132-3141, doi: 10.3969/j.issn.0001-5733.2009.12.023. |
| [35] | Zhdanov M S, Endo M, Cox L H, et al. 2014. Three-dimensional inversion of towed streamer electromagnetic data. Geophysical Prospecting, 62(3): 552-572. |
| [36] | 陈小斌, 胡文宝. 2002. 有限元直接迭代算法及其在线源频率域电磁响应计算中的应用. 地球物理学报, 45(1): 119-130, doi: 10.3321/j.issn:0001-5733.2002.01.015. |
| [37] | 底青云, Unsworth M, 王妙月. 2004. 复杂介质有限元法2.5维可控源音频大地电磁法数值模拟. 地球物理学报, 47(4): 723-730, doi: 10.3321/j.issn:0001-5733.2004.04.026. |
| [38] | 韩波, 胡祥云, Schultz A等. 2015. 复杂场源形态的海洋可控源电磁三维正演. 地球物理学报, 58(3): 1059-1071, doi: 10.6038/cjg20150330. |
| [39] | 何继善. 2010. 广域电磁测深法研究. 中南大学学报(自然科学版), 41(3): 1065-1072. |
| [40] | 李勇, 吴小平, 林品荣. 2015. 基于二次场电导率分块连续变化的三维可控源电磁有限元数值模拟. 地球物理学报, 58(3): 1072-1087, doi: 10.3038/cjg20150331. |
| [41] | 李予国, 段双敏. 2014. 海洋可控源电磁数据预处理方法研究. 中国海洋大学学报, 44(10): 106-112. |
| [42] | 王若, 王妙月, 底青云. 2006. 频率域线源大地电磁法有限元正演模拟. 地球物理学报, 49(6): 1858-1866, doi: 10.3321/j.issn:0001-5733.2006.06.035. |
| [43] | 徐志锋, 吴小平. 2010. 可控源电磁三维频率域有限元模拟. 地球物理学报, 53(8): 1931-1939, doi: 10.3969/j.issn.0001-5733. 2010.08.019. |
| [44] | 薛国强, 闫述, 陈卫营. 2014. 接地源短偏移瞬变电磁法研究展望. 地球物理学进展, 29(1): 177-181, doi: 10.6038/pg20140124. |
| [45] | 阎述, 陈明生. 1999a. 线源频率电磁测深二维正演(一). 煤田地质与勘探, 27(5): 60-62. |
| [46] | 阎述, 陈明生. 1999b. 线源频率电磁测深二维正演(二). 煤田地质与勘探, 27(6): 56-59. |
| [47] | 阎述, 陈明生. 2000. 电偶源频率电磁测深三维地电模型有限元正演. 煤田地质与勘探, 28(3): 50-56. |
| [48] | 张继锋, 汤井田, 喻言等. 2009. 基于电场矢量波动方程的三维可控源电磁法有限单元法数值模拟. 地球物理学报, 52(12): 3132-3141, doi: 10.3969/j.issn.0001-5733.2009.12.023. |