2016, Vol. 59
2. 中国海洋大学海洋地球科学学院, 青岛 266100
2. College of Marine Geosciences, Ocean University of China, Qingdao 266100, China
海洋可控源电磁法(Marine Controlled Source Electromagnetic，简称MCSEM)是海底油气资源勘探的一种重要地球物理方法(Schwalenberg et al., 2005；Weitemeyer et al., 2006；Zach and Brauti, 2009).由于MCSEM对高阻油气储层具有很强的识别能力，近十几年来引起了勘探地球物理学家们极大的关注，其仪器设备、数据采集及反演解释技术均得到了迅猛发展(何展翔和余刚，2008).正演是反演的引擎，反演解释的准确性和高效性在很大程度上依赖于正演.由于海洋环境的特殊性，海底结构往往会非常复杂，在这种情况下，三维正演对于能否正确地解释海洋可控源电磁数据显得十分重要(da Silva et al., 2012).
在过去的几十年，海洋可控源电磁法三维数值模拟技术取得了很快发展和进步.早期主要采用有限差分方法(Newman and Alumbaugh, 1995；Streich，2009)和积分方程法(Xiong and Tripp, 1997；Zhdanov et al., 2006)；近年来有限元法(Badea et al., 2001；Schwarzbach et al., 2011；Puzyrev et al., 2013)和有限体积法(Weiss and Constable, 2006；韩波等，2015)受到了广泛关注.Avdeev (2005)、Börner (2010)及Everett (2012)对近年来地球电磁场数值模拟方法发展及应用做了较为全面的综述.
在众多的数值模拟方法中，有限元法在处理模型复杂度和计算精度方面更具优势，特别是该方法非常适合非结构网格剖分.由于非结构网格剖分能精确地模拟起伏地形和倾斜界面等复杂构造，基于非结构网格的有限元法非常适合复杂地电模型的正演模拟.Rücker等(2006)将非结构网格应用到起伏地形的直流电阻率法有限元正演中，并通过局部节点加密提高了有限元解的精度.Key和Weiss (2006)将非结构网格应用到了二维大地电磁的有限元正演中.Li和Key (2007)将非结构网格应用到海洋可控源电磁二维自适应有限元正演中，算法具有较高的计算精度.Li和Pek (2008)利用非结构网格研究了海底各向异性地层二维大地电磁响应.Liu等(2008)将非结构网格应用到了三维大地电磁矢量有限元正演中.任政勇和汤井田(2009)将局部加密的非结构网格应用到了电阻率法三维有限元正演中，并提出了一种局部节点加密方法.Ren和Tang (2010)应用自适应四面体网格研究了三维直流电阻率问题.Wang等(2013)基于非结构网格研究了电阻率三维各向异性正演问题.Ye等(2014)应用非结构化网格有限元法研究了激发极化法2.5维正演响应.蔡红柱等(2015)基于非结构网格研究了电导率各向异性介质海洋电磁有限元正演算法.杨军等(2015)研究了海洋可控源电磁法三维非结构网格矢量有限元数值模拟.
本文从Maxwell方程组出发推导出关于磁矢量势和电标量势的偏微分方程和有限元方程，然后采用非结构网格对模型区域剖分，并提出了一种非结构网格的局部加密方法，实现了基于局部加密非结构网格的海洋可控源电磁三维有限元正演算法.模型实例表明，基于局部加密的非结构网格可以有效提高有限元数值解的精度，要想获得更高精度的数值解，对异常区域体积约束加密是有必要的.
2 偏微分方程假定时间因子为e-iωt，在拟稳态情形下，电场(E)和磁场(H)满足如下的微分方程：
|
(1) |
|
(2) |
采用库伦规范势的定义，E和H可表示成如下形式：
|
(3) |
式中，A为磁矢量势，ϕ为电标量势.将(3)式代入到(1)和(2)式中，可得到如下关于磁矢量势和电标量势的双旋度方程：
|
(4) |
为了保证解的唯一性，采用库伦规范条件Δ·A=0(Badea et al., 2001)，则上式可变为如下Helmholtz方程：
|
(5) |
|
(6) |
采用合适的边界条件，方程(5)和(6)可以利用有限元法进行求解.为了消除电偶源引起的奇异性及数值模拟困难，将总的电导率分解为背景模型电导率(σp)和异常电导率(σs)，此时总势分解为一次势(Ap, ϕp)和二次势(As, ϕs)：
|
(7) |
对于背景电导率模型，一次势满足的偏微分方程为
|
(8a) |
|
(8b) |
利用式(5)、(6)、(7)、(8)，可得二次势满足的偏微分方程为
|
(9a) |
|
(9b) |
由于Ep=iω(Ap+Δϕp)，代入上式可得：
|
(10a) |
|
(10b) |
由方程(10)可知，在已知背景模型电场的情况下，采用有限元解法，选取合适的边界条件，可得各节点的二次势.
3 有限元方程将方程(10)的磁矢量势沿x、y、z三个方向展开，利用伽辽金法(Jin，2002)，同时结合矢量恒等式和散度定理，可得如下体积分方程组：
|
(11a) |
|
(11b) |
|
(11c) |
|
(11d) |
式中，η为测试函数，通常取插值函数为测试函数，Ω为体积分区域，
本文采用非结构四面体单元对求解区域进行离散，每个单元内电导率为常量，二次势呈线性变化，离散后得到如下形式的线性方程组：
|
(12) |
当采用节点有限元法时，每个节点的未知量为4，假设剖分区域内节点数为M，则未知数为4M，K矩阵为由以下4×4的分块矩阵构成的4M×4M的大型、稀疏、复对称系数矩阵：
|
(13) |
式中，Ni为与节点i对应的线性插值函数，I33为3×3的单位矩阵.右端源项可表示为b=(b1, …, bn)T，且
|
(14) |
在正演模拟中，模型的网格剖分方法按照生成单元类型可分为结构网格和非结构网格剖分.结构网格剖分生成的网格单元一般为正四边形或正六面体，而非结构网格剖分生成的网格单元一般为三角形或四面体.在众多完全非结构网格生成技术中(Shephard and Gcorges, 1991；Moller and Hansbo, 1995；Shewchuk，2002)，数学理论、单元质量和可靠性方面最为出众的为约束Delaunay三角算法(Constraint Delaunay Triangulation，CDT).CDT算法核心思想为：首先将给出的离散点进行三角化，基于给定的边线修订网格；然后参考给定的单元质量和大小约束挑选出需要加密的单元序列；选定单元e，做其外接圆C，添加其圆心点o，连接o和单元e的端点形成新的四面体，为了保证新四面体质量，其外接圆C包含的其他端点都要被移除；此过程一直循环，直到所有单元质量和大小都满足条件为止(Shewchuk，2002).
本文采用的Delaunay网格剖分由Tetgen (Si，2003)生成，它通过Q值(四面体的外接球半径R与四面体的最短边L之比，即Q=R/L)来控制所生成的四面体的质量.一般情况下，质量好的四面体Q值较小，质量差的四面体Q值较大，但是Q值太小，生成的四面体单元数会大大增加，这就会导致有限元计算量大大增加，降低了计算效率.因此，在实际应用Tetgen时，Q值取1~1.4为宜，本文非结构化网格剖分和加密采用的Q值为1.2.
4.2 网格局部加密网格剖分后单元的大小会影响有限元数值解的精度.单元太大，精度往往不高，单元太小，会造成节点太多，从而求解线性方程组的压力较大.使用规则网格时，网格的剖分和加密要靠手动调节.非结构网格的剖分和加密则更加灵活，非结构网格剖分为完全自动剖分，网格加密方法有：全局体积加密(Global Volume-Refinement，GVR)、局部体积加密(Local Volume-Refinement，LVR)和局部测点加密(Local Node-Refinement，LNR).GVR虽然可以提高有限元解的精度，但是其计算消耗非常大，过多的节点对计算效率提高不大，因此非结构网格的加密通常采用LNR和LVR.Rücker等(2006)指出在测点的正下方一定深度增加一个节点，作为生成单元时的控制顶点，额外增加的控制点会改善求解结果，但是它会破坏节点分布的均匀性，从而影响单元的质量，而且节点距离测线太近，会带来过多的节点，并且对精度没有显著提高.任政勇和汤井田(2009)提出通过在测点正下方加入一个小四边形，小四边形的四个顶点作为网格生成时的控制顶点，以达到改善求解精度的效果，并指出精度要求较高时或模型复杂时，LVR是高精度模拟的保证，但是对于起伏地形在测点下方加入四个顶点难以控制，而且这种人为的加入控制顶点，在网格生成时，为保证单元体的质量，有时会有顶点被剔除.为了克服以上弊端，本文非结构化网格的局部加密采用LNR+LVR的局部加密方式，LNR通过寻找测点所在的单元，然后对这些单元进行体积约束加密，由于本文采用二次势的方法，异常区域相当于源区，因此LVR采用对异常体区域进行体积约束加密，网格局部加密策略如图 1所示.
|
图 1 网格局部策略加密示意图 Fig. 1 Mesh locally refinement strategy |
采用图 2所示的海洋可控源一维储层模型对算法进行精度验证.该模型分五层，在深度2~2.4 km范围存在100 Ωm高阻储层.电偶极子源沿x方向激发，测线方向沿y方向，观测三个电磁场Ex、Hy和Hz分量.激发场源位于海底上方100 m，其坐标为(0 m，0 m，900 m)，发射频率为0.1 Hz.在海底布设51个接收站，分布范围为y=-5~5 km，接收站间距200 m.
|
图 2 一维储层模型 Fig. 2 1D reservoir model |
图 3为采用三种不同加密策略得到的测线附近非结构网格.图 3a为直接由TetGen产生的受约束Delaunay剖分网格(CDT)，节点个数为37893，单元个数为227607；图 3b为在测点附近施加体积约束的局部测点加密网格(LNR)，加密后节点个数为52904，单元个数为321698；图 3c为同时在测点附近和异常体区域施加体积约束加密后的网格(LNR+LVR)，由于对一维薄层异常区域进行体积加密，在离源点较远的区域会产生多余的节点，从而增加了计算负担，因此，本文选择了对半径为5 km的圆柱形区域进行局部体积加密，加密后节点个数为114966，单元个数为714748.从图 3可看出，通过LNR+LVR加密后，网格剖分在测点附近和离源点较近的异常区域均得到了加密，而在其他区域增加的节点很少.
|
图 3 一维储层模型测线附近非结构网格y-z切面 Fig. 3 y-z slice of the unstructured tetrahedral mesh along survey line of the 1D reservoir model |
图 4为三种不同加密策略得到的网格二次电场数值解与一维解析解比较，从图中可看出，LNR，LNR+LVR网格与CDT网格相比，精度有不同程度的提高，最大相对误差由CDT的22.3%下降到LNR的8.0%和LNR+LVR的2.3%，平均相对误差也由CDT的7.0%下降到LNR的2.5%和LNR+LVR的1.0%，表明经LNR+LVR网格局部加密后，有限元数值解的精度会有明显地提高.图 5为LNR+LVR网格的电磁场总场分量与一维解析解比较，从图中可看出，电磁场分量三维数值解与一维解析解吻合地很好，幅值最大相对误差不超过1%，相位最大误差不超过0.3°，表明该算法具有很高的计算精度.
|
图 4 一维模型不同加密网格的二次电场与一维解析解的比较 (a)实部；(b)虚部；(c)实部相对误差；(d)虚部相对误差. Fig. 4 A comparison of the secondary electric fields of different meshes and analytic solution for the 1D reservoir model Panel (a) and (b) are the real and imaginary parts, (c) and (d) are their relative errors |
|
图 5 一维模型电磁场总场分量与一维解析解的比较 (a)幅值；(b)相位；(c)幅值相对误差；(d)相位误差. Fig. 5 A comparison of the total electromagnetic fields and analytic solution for the 1D reservoir model Panel (a) and (b) are the amplitude and phase, (c) and (d) are the errors. |
设计的一个二维储层模型如图 6所示，为了对该算法的计算结果进行精度验证，与Li和Key (2007)的二维自适应有限元算法计算结果进行了比较.由图 6可见，二维异常体为梯形，走向沿x方向，上底y方向分布范围为2~4 km，下底y方向分布范围为1~5 km，二维异常体的厚度为400 m，电阻率为100 Ωm.发射场源沿x方向激发，位于海底上方50 m，其位置坐标为(0，-2000 m，950 m)，发射频率为0.1 Hz.剖面沿y方向布设，接收测站分布于海底y=-1.8~8 km之间，间距200 m，共50个接收测站.
|
图 6 二维储层模型 Fig. 6 2D reservoir model |
图 7为剖面附近三种不同加密方式的非结构网格及Li和Key (2007)算法得到的最终细化网格.从图 7可看出，LNR网格在剖面附近得到了加密，LNR+LVR网格在剖面附近和异常区域均得到了加密，而在其他区域增加的节点和单元数很少.图 7a为直接由TetGen产生的受约束Delaunay剖分网格(CDT)，节点个数为10579，单元个数为62023；图 7b为在测点附近施加体积约束的局部测点加密网格(LNR)，加密后节点个数为23321，单元个数为142149；图 7c为同时在测点附近和异常体区域施加体积约束加密后的网格(LNR+LVR)，加密后节点个数为104394，单元个数为660735.图 7d为二维自适应有限元法最终细化的网格，从图中可看出，剖面附近和异常体区域均得到了加密.
|
图 7 二维储层模型三维非结构网格剖分y-z切面及二维自适应有限元算法得到的最终细化网格 Fig. 7 y-z slice of the unstructured tetrahedral mesh along survey line of the 2D reservoir model and the 2D unstructured mesh refined by the adaptive FE method |
图 8为三种不同加密策略得到的网格二次电场数值解与二维有限元数值解的比较，从图 8中可看出，最大相对误差由CDT的47.0%下降到LNR的40.0%和LNR+LVR的9.0%，平均相对误差也由CDT的17.5%下降到LNR的15.0%和LNR+ LVR的3.5%，表明LNR法对算法精度的改善较小，而LNR+LVR法对算法精度的改善更明显，要想得到更高精度的数值解，对异常区域局部体积加密是有必要的.图 9为LNR+LVR网格的电磁场总场分量与二维数值解的比较，从图中可看出，电磁场三维数值解与二维数值解吻合得很好，电磁场幅值最大相对误差约2%，相位的最大误差约0.7°，进一步表明该算法具有较高的计算精度.
|
图 8 二维模型不同加密网格的二次电场与二维有限元解的比较 (a)实部；(b)虚部；(c)实部相对误差；(d)虚部相对误差. Fig. 8 A comparison of the secondary electric fields of different meshes and 2D FE solution for the 2D reservoir model Panel (a) and (b) are the real and imaginary parts, (c) and (d) are their relative errors. |
|
图 9 二维模型电磁场总场分量与二维有限元解的比较 (a)幅值；(b)相位；(c)幅值相对误差；(d)相位误差. Fig. 9 A comparison of the total electromagnetic fields and 2D FE solution for the 2D reservoir model. Panel (a) and (b) are the amplitude and phase, (c) and (d) are the errors. |
第三个模型实例为一个圆盘模型，该模型如图 10所示.Schwarzbach等(2011)和Puzyrev等(2013)曾对该模型进行过正演计算与分析，可作为本文计算结果的参考.如图 10所示，圆盘的半径为2 km，厚度为100 m，上底距海底埋深950 m，电阻率为100 Ωm.发射场源沿x方向激发，位于海底上方100 m，离圆盘中心的水平距离为3 km，其位置坐标为(0，0，900 m)，发射频率为1 Hz.剖面穿过圆盘中心，沿x方向布设，接收测站分布于海底，根据模型的对称性，Ey、Hx、Hz分量在该剖面上消失，测量Ex、Hy、Ez分量.模型的计算区域为66 km ×66 km ×66 km.
|
图 10 圆盘模型 Fig. 10 The Canonical disc model |
表 1列出了采用非结构网格(CDT，LNR，LNR+LVR)与结构网格的有限元计算参数统计.从表 1可看出，非结构网格经局部加密后，占用内存和消耗时间会相应增加，计算代价也会相应提高.但是采用非结构网格计算时，节点数较少，单元数是节点数的6倍左右，而采用结构网格计算时，节点数比单元数多.因此，在占用内存和消耗时间方面，采用非结构网格计算比采用结构网格计算要少的多，下降约一个数量级，表明采用非结构网格可以极大地提高节点有限元法的计算效率.
|
表 1 圆盘模型网格参数统计 Table 1 The mesh variants and their statistics of the Canonical disc model |
图 11为圆盘模型的总场分量(Ex、Hy、Ez)有限元计算结果，与Schwarzbach等(2011)和Puzyrev等(2013)中的计算结果相比，一致性非常好，表明该算法的计算结果稳定可靠.从图 10可看出，由于二次场与一次场相比幅值很小，当收发距小于2.5 km时，含有油气储层的响应与不含油气储层的背景场信息很难区分开来，当收发距大于2.5 km时，随着背景场的逐渐减小，油气储层响应能从背景场中分离出来，表现为曲线衰减变缓的特征，当收发距大于7 km时，由于远离异常体，此时二次场也逐渐变小，曲线逐渐与背景场重合.更多关于该模型的正演特征及空气波的影响，可见参考文献Schwarzbach等(2011)和Puzyrev等(2013).
|
图 11 圆盘模型总场分量(Ex、Hy、Ez)有限元解 Fig. 11 Non-vanishing total field components (Ex, Hy, Ez) for the Canonical disc model |
本算例对一复杂海底地形资料进行了数值模拟分析，海底地形如图 12a所示，海水深度在793~1933 m之间变化，设海水的电阻率为0.3 Ωm，海底的电阻率为1 Ωm.在海底沿y方向布置一条15 km长的剖面，从-5 km到10 km变化，发射源位置位于(0 m，-1000 m，1600m)，发射频率为0.1 Hz.分别计算broadside和inline装置的电磁响应.海底地形非结构网格剖分及对剖面附近区域进行局部加密处理，如图 12b所示.
|
图 12 海底地形(a)及其非结构网格剖分(b) Fig. 12 Complex seafloor topography (a) and unstructured tetrahedral mesh of the seafloor topography (b) |
为了分析海底地形对电磁场的影响，将正演计算结果与水平海底地形的响应作了比较，水平海底地形模型的海水深度设置为发射源位置的海水深度.图 13为broadside和inline装置的起伏海底地形与水平海底地形模型的正演计算结果比较，从图 13可看出，无论是inline装置，还是broadside装置，海底地形会对海底电磁场响应会产生较大的影响，因此，海底地形影响是海洋可控源电磁场反演解释时有必要考虑的一个因素.
|
图 13 起伏地形与水平地形模型电磁场响应的比较 (A)为broadside装置，(B)为inline装置. Fig. 13 The comparisons of the total electromagnetic fields between the complex seafloor and flat seafloor for the broadside (A) and inline (B) configuration |
本文实现了基于完全非结构网格的海洋可控源电磁场三维有限元正演算法.文中将电磁势分解为一次势和二次势，避免了总场算法的场源奇异性.针对非结构网格具有局部加密的特点，实现了几种不同非结构网格的局部加密，对比分析了不同局部加密网格的有限元数值解的精度.
一、二维数值模拟实例表明，直接采用CDT网格计算，精度较差，对网格进行LNR加密后，能够有限地提高数值解的精度，采用LNR+LVR加密后，可以大幅度地提高有限元数值解的精度，说明计算精度要求较高或模型较复杂时，采用LNR+LVR局部加密网格是有必要的.
三维数值模拟实例表明，采用非结构网格剖分能够以最少的节点数达到最优化的单元数，在保证精度的同时，提高计算效率；复杂海底地形模拟实例表明，该算法可以模拟复杂的地电模型，如复杂起伏海底地形，能够为实际的海洋电磁资料解释提供一种途径.
基于局部加密的非结构网格有限元法尽管可以得到高精度的数值解，但网格的剖分和加密仍具有一定的人为性，要想实现完全自动的网格加密，则需借助于自适应有限元法.自适应有限元法可根据每个单元的后验误差估计来调整单元的大小，从而得到高精度的数值解，这将是下一阶段的研究目标.
| Avdeev D B. 2005. Three-dimensional electromagnetic modelling and inversion from theory to application. Surveys in Geophysics, 26(6): 767-799. DOI:10.1007/s10712-005-1836-x | |
| Badea E A, Everett M E, Newman G A, et al. 2001. Finite-element analysis of controlled-source electromagnetic induction using Coulomb-gauged potentials. Geophysics, 66(3): 786-799. DOI:10.1190/1.1444968 | |
| Börner R U. 2010. Numerical modelling in geo-electromagnetics: advances and challenges. Surveys in Geophysics, 31(2): 225-245. DOI:10.1007/s10712-009-9087-x | |
| Cai H Z, Xiong B, Michael Zhdanov M. 2015. Three-dimensional marine controlled-source electromagnetic modelling in anisotropic medium using finite element method. Chinese J. Geophys. , 58(3): 2839-2850. DOI:10.6038/cjg20150818 | |
| da Silva N Y, Morgan J V, MacGregor L, et al. 2012. A finite element multifrontal method for 3D CSEM modeling in the frequency domain. Geophysics, 77(2): E101-E115. DOI:10.1190/geo2010-0398.1 | |
| Everett M E. 2012. Theoretical developments in electromagnetic induction geophysics with selected applications in the near surface. Surv. Geophys., 33(1): 29-63. DOI:10.1007/s10712-011-9138-y | |
| Han B, Hu X Y, Adam Schultz A, et al. 2015. Three-dimensional forward modeling of the marine controlled-source electromagnetic field with complex source geometries. Chinese J. Geophys. , 58(3): 1059-1071. DOI:10.6038/cjg20150330 | |
| He Z X, Yu G. 2008. Marine EM survey technology and its new advances. Progress in Exploration Geophysics , 31(1): 2-9. | |
| Jin J M. The Finite Element Method in Electromagnetics.(2nd ed).New Jersey: Wiley-IEEE Press, 2002. | |
| Key K, Weiss C. 2006. Adaptive finite element modeling using unstructured grids: The 2D magnetotelluric example. Geophysics, 71(6): G291-G299. DOI:10.1190/1.2348091 | |
| Li Y G, Key K. 2007. 2D marine controlled-source electromagnetic modelling, Part I-An adaptive finite-element algorithm. Geophysics, 72(2): WA51-WA62. DOI:10.1190/1.2432262 | |
| Li Y G, Pek J. 2008. Adaptive finite element modelling of two-dimensional magnetotelluric fields in general anisotropic media. Geophysical Journal International, 175(3): 942-954. DOI:10.1111/gji.2008.175.issue-3 | |
| Liu C S, Ren Z Y, Tang J T, et al. 2008. Three-dimensional magnetotellurics modeling using edge-based finite-element unstructured meshesd. Applied Ggeophysics, 5(3): 170-180. DOI:10.1007/s11770-008-0024-4 | |
| Moller P, Hansbo P. 1995. On advancing front mesh generation in three dimensions. Int. J. Numer. Mesth. Engrg., 38(21): 3551-3569. DOI:10.1002/(ISSN)1097-0207 | |
| Newman G A, Alumbaugh D L. 1995. Frequency-domain modelling of airborne electromagnetic responses using staggered finite differences. Geophysicsal Prospecting, 43(8): 1021-1042. DOI:10.1111/gpr.1995.43.issue-8 | |
| Puzyrev V J, Koldan J, de la Puente J, et al. 2013. A parallel finite-element method for three-dimenstional controlled-source electromagnetic forward modelling. Geophysical Journal International, 193(2): 678-693. DOI:10.1093/gji/ggt027 | |
| Ren Z Y, Tang J T. 2009. Finite element modeling of 3-D DC resistivity using locally refined unstructured meshes. Chinese J. Geophys. , 52(10): 2627-2634. DOI:10.3969/j.issn.0001-5733.2009.10.023 | |
| Ren Z Y, Tang J T. 2010. 3D direct current resistivity modeling with unstructured mesh by adaptive finite-element method. Geophysics, 75(1): H7-H17. DOI:10.1190/1.3298690 | |
| Rücker C, Günther T, Spitzer K. 2006. Three-dimensional3D modelling and inversion of DC resistivity data incorporating topography-I. Modelling. Geophysical Journal International, 166(2): 495-505. DOI:10.1111/gji.2006.166.issue-2 | |
| Schwalenberg K, Willoughby E, Mir R, et al. 2005. Marine gas hydrate electromagnetic signatures in Cascadia and their correlation with seismic blank zones. First Break, 23(4): 57-63. | |
| 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. DOI:10.1111/gji.2011.187.issue-1 | |
| Shephard M S, Gcorges M K. 1991. Automatic three dimensional mesh generation by the finite Octree technique. Int. J. Numer. Mesth. Engrg., 32(4): 709-749. DOI:10.1002/(ISSN)1097-0207 | |
| Shewchuk J R. 2002. Delaunay refinement algorithms for triangular mesh generation. Computational Geometry: Theory and Applications, 22(1-3): 21-74. DOI:10.1016/S0925-7721(01)00047-5 | |
| Si H. 2003. TETGEN: A 3D Delaunay Tetrahedral Mesh Generator. http://tetgen.berlios.de. | |
| 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. DOI:10.1190/1.3196241 | |
| Wang W, Wu X P, Spitzer K. 2013. Three-dimensional DC anisotropic resistivity modelling using finite elements on unstructured grids. Geophysical Journal International, 193(2): 734-746. DOI:10.1093/gji/ggs124 | |
| 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. DOI:10.1190/1.2356908 | |
| Weitemeyer K, Constable S, Key K. 2006. Marine EM techniques for gas-hydrate detection and hazard mitigation. The Leading Edge, 25(5): 629-632. DOI:10.1190/1.2202668 | |
| Xiong Z H, Tripp A C. 1997. 3-D electromagnetic modeling for near-surface targets using integral equations. Geophysics, 62(4): 1097-1106. DOI:10.1190/1.1444210 | |
| Yang J, Liu Y, Wu X P. 2015. 3D simulation of marine CSEM using vector finite element method on unstructured grids. Chinese J. Geophys. , 58(8): 2827-2838. DOI:10.6038/cjg20150817 | |
| Ye Y X, Li Y G, Deng J Z, et al. 2014. 2. 5D induced polarization forward modeling using the adaptive finite-element method. Applied Geophysics, 11(4): 500-507. | |
| Zach J J, Brauti K. 2009. Methane hydrates in controlled-source electromagnetic surveys-analysis of a recent data example. Geophysical Prospecting, 57(4): 601-614. DOI:10.1111/gpr.2009.57.issue-4 | |
| Zhdanov M S, Lee S K, Yoshioka K. 2006. Integral equation method for 3D modeling of electromagnetic fields in complex structures with inhomogeneous background conductivity. Geophysics, 71(6): G333-G345. DOI:10.1190/1.2358403 | |
| 蔡红柱, 熊彬, Michael ZhdanovM. 2015. 电导率各向异性的海洋电磁三维有限单元法正演. 地球物理学报, 58(8): 2839–2850. DOI:10.6038/cjg20150818 | |
| 韩波, 胡祥云, Adam SchultzA, 等. 2015. 复杂场源形态的海洋可控源电磁三维正演. 地球物理学报, 58(3): 1059–1071. DOI:10.6038/cjg20150330 | |
| 何展翔, 余刚. 2008. 海洋电磁勘探技术及新进展. 勘探地球物理进展, 31(1): 2–9. | |
| 任政勇, 汤井田. 2009. 基于局部加密非结构化网格的三维电阻率法有限元数值模拟. 地球物理学报, 52(10): 2627–2634. DOI:10.3969/j.issn.0001-5733.2009.10.023 | |
| 杨军, 刘颖, 吴小平. 2015. 海洋可控源电磁三维非结构矢量有限元数值模拟. 地球物理学报, 58(8): 2827–2838. DOI:10.6038/cjg20150817 | |