2019, Vol. 62
海洋可控源电磁法是一种重要的海底油气勘探的地球物理方法.海洋可控源电磁法可以获得海底地层的电性信息，因此可对海底地质情况和高阻油气储层进行评价.Chave和Cox(1982)对水平电流源和垂直电流源在一维海洋模型中产生的频率域响应进行了研究.Edwards和Chave(1986)采用频率域到时间域转换方法获得水平电偶极子在海水和基底双半空间模型中产生的阶跃和脉冲响应.Cheesman等(1987)对不同瞬变电磁系统的海底电导率成像理论进行了研究.Everett和Edwards(1992)用有限元法对2.5D海洋电磁模型进行了时间域正演模拟.20世纪90年代以来，三维正演模拟技术得到了快速发展.Newman和Alumbaugh(1995)和Commer和Newman(2004)分别提出了三维频率域有限差分算法和直接求解时间域响应的有限差分算法；Avdeeva等(2007)应用这两种算法进行对比研究，探究了频率域和时间域海洋可控源方法对油气储层的探测能力，发现时间域海洋可控源方法对于浅海油气勘探效果更好；Zhdanov等(2006)将积分方程方法应用到频率域海洋电磁中取得了良好的应用效果；Swidinsky和Edwards(2010)用积分方程方法模拟了时间域三维海洋背斜高阻油气藏的电磁响应；Um等(2010)提出了基于非结构有限元的时间域海洋电磁模拟算法，并讨论了使用垂直电流源的响应特点(Um et al., 2012).
随着地电模型复杂程度的增加，非结构有限元由于可对任意复杂模型进行模拟，优势变得越来越明显.然而，有限元模拟精度非常依赖网格质量.为快速获取最优化网格，许多学者近年来研究了自适应网格加密算法.Ren和Tang(2010)使用基于梯度恢复的自适应方法对直流问题进行了研究；Schwarzbach(2011)使用全局自适应高阶有限元法进行了频率域海洋可控源电磁模拟；Ren等(2013)使用面向目标自适应有限元对3D大地电磁问题进行了研究.Yin等(2016)使用面向目标的自适应有限元方法对频率域航空电磁问题进行了正演模拟；Zhang等(2018)使用全局自适应有限元方法对时间域航空电磁问题进行了研究.本文将Zhang等(2018)提出的时间域自适应有限元方法应用到时间域海洋电磁正演模拟中.
本文使用基于自适应非结构网格的有限元算法对时间域海洋电磁响应进行模拟, 并在自适应过程中采用基于法向电流连续的后验误差评估方法.为进行时间离散，我们采用后推欧拉方法.该方法对大时间步长也可稳定求解(Um, 2011; Oldenburg et al., 2013)，使得我们在求解晚期道电磁响应时可以使用较大的时间步长，提高计算效率.为了控制自适应过程的网格数同时保证正演模拟的稳定性，我们对所有时间道记录的需要加密的网格进行融合，从而得到总的需要加密的网格，并按照某一比例对需要加密的网格进行随机挑选并加密，生成适用于多时间道海洋电磁正演模拟的有效网格.我们使用大型稀疏矩阵直接求解器MUMPS进行方程组的求解.当发射电流为下阶跃波时，初始电场用节点有限元方法求解泊松方程进行计算.
为验证本文算法精度，我们计算了均匀海底半空间的时间域电磁响应并与均匀半空间的半解析解进行对比.进而，对三维海底起伏地形及起伏海底下存在高阻异常体的情况计算了半拖曳式和双船拖曳式时间域海洋电磁响应，并对海底起伏地形及高阻层的识别特征进行了分析和研究.
1 正演算法 1.1 时间域非结构矢量有限元方法在一个给定的计算域，时间域电场扩散方程为
|
(1) |
式中，e(r, t)为t时刻位置r上的电场，μ为介质的磁导率，σ为电导率，js为源电流密度.
将计算域剖分为与之匹配的四面体单元.在每一个四面体单元内，电场可展开为
|
(2) |
式中，nie(r)和eie(t)分别为第e个小单元内的矢量基函数和第i个棱边上的电场(金建铭, 1998).
定义一个残差矢量
|
(3) |
并且确保每个四面体的加权残差积分为零
|
(4) |
其中，Re为四面体单元e的残差矢量，Ve为四面体单元e的体积.将(2)式代入(4)式可得
|
(5) |
矩阵Ae和Be的第(i, j)个元素、矩阵Se的第i个元素可表示为
|
(6) |
|
(7) |
|
(8) |
对于时间离散我们采用后推欧拉法.通过Taylor展开可近似得到二阶形式(Um et al., 2010)
|
(9) |
将式(9)代入式(5)并组装所有单元得
|
(10) |
而边界条件使用Dirichlet边界条件(Um et al., 2010)
|
(11) |
将公式(10)简写为如下矩阵方程：
|
(12) |
式中，K为系数矩阵，e为棱边电场，而b是与源项和前两个时刻棱边电场相关的已知项.本文使用直接求解器MUMPS对(12)式进行直接求解，得到四面体棱边上的电场，并通过(2)式插值得到任意点的电场.如果Δt保持不变，(12)式求解时系数矩阵只需分解一次.
1.2 初始电场计算为求解方程(12)，我们需要给出初始时刻棱边上的电场.当采用电性源发射上阶跃波或脉冲波时，初始电场为零；当采用电性源发射下阶跃波时，初始电场为一个直流电阻率问题的解(Um et al., 2010).
点源产生的电位满足泊松方程.由于海洋电磁的求解区域非常大，在边界处的电位满足Dirichlet边界条件，则电位满足的边值问题如下：
|
(13) |
|
(14) |
式中φ(r)为在位置r处的电位.我们采用非结构节点有限元进行求解.为此，在每一个四面体单元内将电势展开为
|
(15) |
其中，ωie和φie分别为第e个四面体单元内的标量插值基函数和第i个点上的电位.
定义如下标量残差：
|
(16) |
并且使每一个四面体的加权残差积分为零
|
(17) |
其中pe为四面体单元e的残差.将(15)式代入(17)式有
|
(18) |
矩阵Me第(i, j)个元素、矩阵qe的第i个元素为
|
(19) |
|
(20) |
将每个四面体满足的方程(18)组合起来可得如下矩阵方程组：
|
(21) |
式中，F为系数矩阵，Φ为四面体节点处的电位，g为源项.我们同样使用直接求解器MUMPS求解式(21).在得到四面体节点处的电位后，可以计算四面体棱边上的电场.
1.3 后验误差评估首先定义f为四面体Ti和Tj的交界面，四面体Ti和Tj中的电场分别为ei和ej，则根据电流密度法向连续条件有
|
(22) |
式中，nf为交界面f的单位法向量，Ji和Jj为穿过交界面f的电流密度.根据公式J=σe，(22)式可写为
|
(23) |
由于有限元的解存在误差，导致方程(23)会有一个残差，我们使用这个残差对后验误差进行评估(Ren et al., 2013)
|
(24) |
式中，m表示四面体单元i的面的局部编号，j表示与四面体i交界面为m的四面体.
1.4 自适应过程
时间域海洋可控源电磁正演模拟的非结构网格采用如下自适应方式生成：(1)读入初始网格信息，并给出最大迭代次数Ndmax和最大网格数Nsmax.为了保持计算的稳定性，我们给出最小的单元体积Vmin；(2)当迭代次数Nd小于Ndmax且网格数Ns小于Nsmax时，计算每个单元每一时间道的相对后验误差
为了更直观的展示自适应网格加密过程，图 1给出了自适应过程的流程图.
|
图 1 自适应过程流程图 Fig. 1 Flow chart for the adaptive process |
为了验证本文自适应算法的精度，我们将本文计算结果与海底半空间模型的半解析解进行对比.为此，首先建立如下直角坐标系：x轴沿测线方向，y轴垂直测线方向，x, y-平面位于海水表面，z轴向下为正.计算模型如图 2所示，海水厚度为400 m，电阻率为0.3 Ωm；海底均匀半空间电阻率为1 Ωm.发射源是一个沿x方向长度为300 m的电偶源，发射电流为1 A，发射信号为下阶跃脉冲.发射源位于距海底50 m的海水中，接收点位于距海底30 m的海水中，16个接收点的坐标为：同线x=1000~8000 m，y=0 m，旁线x=0 m，y=1000~8000 m，接收点距均为1000 m.计算使用电脑的处理器为3.20 GHz，运行内存为128 GB.计算过程中的相关信息如表 1所示.图 3给出了本文自适应有限元结果与一维半解析解的对比，其中(a)和(b)、(c)和(d)及(e)和(f)分别给出同线Ex、同线Ez和旁线Ex场值对比及误差曲线.对于同线Ex最大相对误差小于1.5%，同线Ez最大相对误差小于3.5%，而旁线Ex除了场值变号处，其他测点相对误差均小于5%，说明本文算法具有较高精度.另外，从图中给出的结果可以看出，对于旁线装置，观测点位于发射偶极的赤道方向，电磁场随时间衰减过程中发生变号现象.
|
图 2 海底均匀半空间模型 Fig. 2 A homogeneous half-space under the ocean |
|
表 1 计算过程中相应参数 Table 1 The parameters in the process of calculation |
|
图 3 本文计算结果和一维半解析解对比(模型参见图 2) (a)(b)同线Ex和相对误差；(c)(d)同线Ez和相对误差；(e)(f)旁线Ex和相对误差.图中虚线表示场值变号. Fig. 3 Comparison of FE solution from this paper and 1D semi-analytical solution for the homogeneous half-space model in Fig. 2. (a), (b) In-line Ex and relative error; (c), (d) In-line Ez and relative error; (e), (f) Broadside Ex and relative error. The dashed line represents the sign reversal of the field. |
为了展示本文网格加密的合理性，图 4中我们展示了不同时间道的网格加密和后验误差.由于本文使用的是全局加密方法，为了提高接收点处响应的计算精度，我们在自适应网格细化过程中，通过设定一个阈值对接收点处单元棱边长度进行限制.从图 4中可以看出，由早期到晚期时间道，加密范围逐渐增大；在早期时间道，近发射源处的后验误差较大，该区域网格得到加密；而在晚期时间道，离发射源较远的位置后验误差较大，因此晚期时间道远离发射源的位置网格也得到有效加密.这与电磁场的传播特征完全一致，早期道电磁场主要集中在发射源附近，因此源附近区域需要得到加密，而在晚期道随着电磁场向远处扩散，影响范围逐渐增大，远处区域影响增大，必须进行相应的网格加密.
|
图 4 不同时间道的自适应网格加密 Fig. 4 Adaptive grids for different time channels |
为了研究海洋可控源电磁法对三维高阻层的探测能力，我们建立如图 5所示的地电模型.高阻异常体大小为3000 m×3000 m×200 m，顶部埋深为1000 m，异常体中心坐标为(4000 m, 0 m, 1500 m)，电阻率为100 Ωm.其余计算参数与图 1相同.
|
图 5 海底高阻异常体模型 Fig. 5 A resistive body buried in a half-space under the ocean |
我们对模型采用同线装置进行观测，测线范围为y=0 m，x=1000~8000 m，点距为100 m.计算结果如图 6所示，其中(a)(b)为存在高阻异常体时总的电场响应，而(c)(d)为相对于均匀半空间模型的相对百分异常.由图可见，仅从电场响应中我们难以识别高阻异常体的存在和几何位置.然而，在背景场去除后从相对百分异常中可以清晰地看出高阻异常体的存在，Ex最大相对百分异常为19.6%，而Ez的相对百分异常成倍增加.最大异常位置与异常体边界对应关系良好，因此从相对异常的分布特征可以很好地确定异常体位置和边界.必须指出的是，相比于Ex，Ez晚期道信号发生变号现象(图 6b)，但电场响应和相对百分异常均非常光滑，反映电磁场在海底的缓慢扩散过程.
|
图 6 高阻异常体模型电场响应(a)(b)及与背景模型的相对百分比异常(c)(d).图中虚线表示场值变号 Fig. 6 Electric fields response for a resistive target (a) (b) and the relative anomalies with background model (c) (d). The dashed line represents the sign reversal of the field |
实际观测中，由于船上安装有声纳系统，无论海底是平坦还是起伏，海底地形均可以有效测定，因此我们可以事先计算出海底半空间响应，并从实际观测数据中减除，获得异常体相对百分异常，从而实现对异常体的有效识别.
2.3 起伏海底三维高阻异常体模型为了研究海底地形对高阻异常体响应的影响特征，我们设计了如图 7所示的隆起和凹陷海底地形.该隆起和凹陷的高度均为200 m，x方向长2000 m，y方向长1000 m.在该隆起和凹陷地形下存在一个高阻异常体(如图中白色区域所示)，其参数与图 5中的高阻异常体一致.对于隆起和凹陷地形，我们的接收机到海底表面的距离始终保持30 m.图 8给出了隆起和凹陷海底地形下无高阻体和存在高阻体时电场响应及相对百分异常.其中(a)和(b)分别为隆起和凹陷纯地形电场响应，(c)和(d)分别为隆起和凹陷地形下埋有高阻异常体的电场响应，(e)和(f)分别为隆起和凹陷地形下埋有高阻异常体时电场相对百分异常.从图 8a和8b中可以看出：当海底地形存在起伏时，电场曲线的形态随地形发生改变，且曲线形态与地形形态一致.对比图 8a和8c，两者的差异不是十分明显.然而，将电场响应转换成相对百分异常发现其最大值可达20%，且百分异常极值点位置与远离发射源的异常体边界相对应.在近源处的异常体边界处，出现负异常的最大值.因此，在海底存在地形起伏的条件下，通过计算相对百分异常，并根据最大正负异常的位置可确定异常体的几何位置，这与前述水平海底的结论一致，说明地形的存在不会影响高阻异常体的识别特征.同样，通过对比图 8b、8d和8f我们也可得出相似的结论.
|
图 7 隆起和凹陷海底地形下有一高阻异常体 Fig. 7 A hill or valley ocean bottom with an embedded resistive body |
|
图 8 起伏海底地形下存在高阻异常体的电场响应 (a)和(b)分别为隆起和凹陷地形的纯地形响应；(c)和(d)分别为隆起和凹陷地形下埋有高阻异常体的响应；(e)和(f)分别为隆起和凹陷地形下存在高阻异常体的相对百分异常. Fig. 8 Electric field response for a resistive target buried under a topographic ocean floor (a) and (b) Pure topographic response for the hill and the valley; (c) and (d) Response for a resistive target embedded under the hill and the valley; (e) and (f) Relative anomaly of a resistivity target embedded under the hill and the valley. |
为了进一步对时间域海洋电磁方法进行研究，我们设计了如图 9所示的双船拖曳式海洋电磁测量模式.首先考虑图 5中给出的水平海底中存在高阻体的模型.同线和旁线测量收发距均为6000 m，测量过程中双船保持同步移动，观测点距为200 m.其中，同线装置测线设定为y=0 m，发射源从x=-8000~9000 m变化，相应地接收机范围为x=-2000~15000 m.对于旁线装置，发射源位于y=0 m处，接收机设定位于y=6000 m处，发射源和接收机变化范围为x=-1000~9000 m.图 10给出了双船拖曳海洋电磁响应计算结果，其中(a)(c)(e)和(b)(d)(f)分别为同线和旁线的总电场响应、异常响应和相对百分异常.对于同线测量，从图 10a可以看出，当发射源驶过异常体和接收机驶过异常体时均出现明显异常.通过图 10c和10e可以更为直观地对异常进行分析.事实上，当发射源逐渐靠近高阻体时，在接收机处出现负异常，当发射源到达异常体右边界处时负异常达到最大，当发射源驶过异常体上方时，异常由负变到正，当发射源达到异常体左边界处，正异常达到最大，随着发射源逐渐远离异常体，异常逐渐减小.但随着接收机逐渐接近异常体，异常又开始增大，当接收机到达异常体右边界附近时，异常再次达到最大，随着接收机驶过异常体，正异常逐渐减小，且由正异常变到负异常，当接收机到达异常体的左边界处负异常达到最大.对比发射源越过异常体和接收点越过异常体时的异常曲线可以发现，异常信号极大值与异常体边界之间存在明确的对应关系，因此异常体位置可以很好地被确定.对于旁线测量，通过图 10b可以看出，电磁异常反映高阻体不明显，并且信号在某一时间道发生变号与图 3e得出的结论相似.然而，从图 10d和10f给出的异常响应和相对异常可以看出，当发射源驶过异常体时，在异常体边界处异常出现极大值，很好的界定了异常体大小和位置信息.
|
图 9 双船拖曳测量模式示意图 Fig. 9 Sketch map for a dual-ship towed MCSEM |
|
图 10 双船拖曳海洋电磁水平海底下存在高阻体电磁响应 (a)(c)(e)为同线测量响应，而(b)(d)(f)为旁线测量电磁响应.图中虚线表示场值变号. Fig. 10 EM responses for the resistive target model under a level ocean floor for dual-ship towed configuration (a), (c), and (e) are for in-line mode, while (b), (d), and (f) are for the broadside mode. The dashed lines represent the sign reversal of the field. |
为进一步讨论地形对海洋电磁响应的影响及高阻体识别特征，我们针对图 9给出的双船拖曳作业方式对图 7给出的模型进行正演模拟.图 11和12分别给出同线和旁线装置的计算结果，其中(a)(c)(e)(g)针对隆起海底地形模型，而(b)(d)(f)(h)针对凹陷海底地形模型.首先，从图 11给出的纯地形响应(a)(b)和带地形异常体响应(c)(d)可以看出，无论对于隆起还是凹陷地形，同线装置海洋电磁响应与地形之间存在一定的对应关系，异常极值点对应地形拐点.然而, 由于地形效应比异常体响应强得多，异常体响应被掩盖.相比之下，由图 11e—11h给出的异常响应和相对百分异常可以清晰看出异常响应和异常体位置之间的对应关系.另外，对比海底不存在(图 10c，10e)和存在地形(图 11e—11h)两种情况发现，虽然地形对总响应影响较大，但同线装置异常响应及百分异常特征基本不受地形影响.这为双船拖曳海洋电磁勘探海底高阻目标体提供良好的物理前提.
|
图 11 双船拖曳同线观测电磁响应 图中虚线表示场值变号. Fig. 11 MCSEM responses for dual-ship towed in-line measurement The dashed line represents the sign reversal of the field. |
|
图 12 双船拖曳旁线测量电磁响应图中虚线表示场值变号. Fig. 12 MCSEM responses for dual-ship towed broadside measurementThe dashed line represents the sign reversal of the field. |
图 12a—12d给出双船拖曳海洋电磁旁线观测电磁响应, 从中可以看出与图 11相似的特征：由于海底地形的存在导致响应变得复杂，异常体响应被掩盖.另外，由于电磁信号从某一时间道开始发生变号，导致电场响应曲线与地形起伏的对应关系变得模糊.然而，从图 12e—12h给出的异常响应和百分异常可以看出，在异常体的边部异常响应幅值较大，异常体的中部异常幅值相对较小，两个负向异常峰值分别对应高阻体的左右边界.对比图 12e—12h与图 10e, 10f发现, 由于高阻体引起的异常受到了地形的影响，中部区域实质上为地形响应与异常体响应的叠加.这种特征从百分异常表现更为明显，特别是在场值发生变号的时间道附近，相对百分异常曲线与无地形时存在很大区别.对于隆起地形，相对异常在地形中心呈现正向极大值异常，而对于凹陷地形，相对异常在地形中心呈现负向极大值异常.因此，在进行全拖曳式海洋电磁勘探时，同线装置较之于旁线装置异常形态简单、异常幅值较大，对地下高阻体识别能力较强.
3 结论本文将时间域自适应有限元方法应用于海底复杂构造的时间域海洋电磁正演模拟取得了很好的应用效果.该方法不仅有效地改善了网格质量提高了计算精度，同时所采用的随机网格优选方法能很好地控制自适应过程中网格细化的数量，因而提高计算效率.后推欧拉方法适合于大时间步长，并且时间步长不变时，使用直接求解器只需对系数矩阵进行一次分解，大大提高了求解速度.需要指出的是，当发射波形为下阶跃时，正演模拟过程中还需要计算初始直流电场.时间域海洋电磁对海底存在高阻体具有较强的识别能力，特别是异常响应和百分异常响应的形态和峰值可以有效确定异常体边界.当海底存在起伏地形时，虽然异常响应和百分异常响应仍能在一定程度上反映海底高阻异常体的分布特征，但海底地形造成电磁响应和异常体几何位置之间的对应关系复杂化.全拖曳式海洋电磁同线装置较之于旁线装置对海底高阻目标体具有更强的勘探能力.
致谢 作者向各位审稿人及《地球物理学报》编辑表示衷心感谢.另外，也向殷长春“千人计划”电磁研究团队的其他成员表示感谢.
Avdeeva A, Commer M, Newman G A. 2007. Hydrocarbon reservoir detectability study for marine CSEM methods: time-domain versus frequency-domain.//2007 SEG Annual Meeting. San Antonio: SEG.
|
Chave A D, Cox C S. 1982. Controlled electromagnetic Sources for measuring electrical conductivity beneath the oceans: 1. Forward problem and model study. Journal of Geophysical Research: Solid Earth, 87(B7): 5327-5338. DOI:10.1029/JB087iB07p05327 |
Cheesman S J, Edwards R N, Chave A D. 1987. On the theory of sea-floor conductivity mapping using transient electromagnetic systems. Geophysics, 52(2): 204-217. DOI:10.1190/1.1442296 |
Commer M, Newman G A. 2004. A parallel finite-difference approach for 3D transient electromagnetic modeling with galvanic sources. Geophysics, 69(5): 1192-1202. DOI:10.1190/1.1801936 |
Edwards R N, Chave A D. 1986. A transient electric dipole-dipole method for mapping the conductivity of the sea floor. Geophysics, 51(4): 984-987. DOI:10.1190/1.1442156 |
Everett M E, Edwards R N. 1992. Transient marine electromagnetics: the 2.5-D forward problem. Geophysical Journal International, 113(3): 545-561. DOI:10.1111/j.1365-246X.1993.tb04651.x |
Jin J M. 1998. Electromagnetic Finite-Element Method (in Chinese). Wang J G, trans. Xi′an: Xidian University Press.
|
Newman G A, Alumbaugh D L. 1995. Frequency-domain modelling of airborne electromagnetic responses using staggered finite differences. Geophysical Prospecting, 43(8): 1021-1042. DOI:10.1111/j.1365-2478.1995.tb00294.x |
Oldenburg D W, Haber E, Shekhtman R. 2013. Three dimensional inversion of multisource time domain electromagnetic data. Geophysics, 78(1): E47-E57. DOI:10.1190/geo2012-0131.1 |
Ren Z Y, Kalscheuer T, Greenhalgh S, et al. 2013. A goal-oriented adaptive finite-element approach for plane wave 3-D electromagnetic modelling. Geophysical Journal International, 194(2): 700-718. DOI:10.1093/gji/ggt154 |
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 |
Schwarzbach C, Börner R, Spitzer K. 2011. Three-dimensional adaptive higher order finite element simulation for geo-electromagnetic—a marine CSEM example. Geophysical Journal International, 187(1): 63-74. DOI:10.1111/j.1365-246X.2011.05127.x |
Swidinsky A, Edwards R N. 2010. The transient electromagnetic response of a resistive sheet: an extension to three dimensions. Geophysical Journal International, 182(2): 663-674. DOI:10.1111/j.1365-246X.2010.04649.x |
Um E S, Harris J M, Alumbaugh D L. 2010. 3D time-domain simulation of electromagnetic diffusion phenomena: A finite-element electric-field approach. Geophysics, 75(4): F115-F126. DOI:10.1190/1.3473694 |
Um E S. 2011. Three-dimensional finite-element time-domain modeling of the marine controlled-source electromagnetic method [Ph. D. thesis]. San Francisco: Stanford University.
|
Um E S, Alumbaugh D L, Harris J M, et al. 2012. Numerical modeling analysis of short-offset electric-field measurements with a vertical electric dipole source in complex offshore environments. Geophysics, 77(5): E329-E341. DOI:10.1190/GEO2011-0442.1 |
Yin C C, Zhang B, Liu Y H, et al. 2016. A goal-oriented adaptive finite-element method for 3D scattered airborne electromagnetic method modeling. Geophysics, 81(5): E337-E346. DOI:10.1190/GEO2015-0580.1 |
Zhang B, Yin C C, Ren X Y, et al. 2018. Adaptive finite element for 3D time-domain airborne electromagnetic modeling based on hybrid posterior error estimation. Geophysics, 83(2): WB71-WB79. DOI:10.1190/geo2017-0544.1 |
Zhdanov M S, Lee S K, Yoshioka Ken. 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 |
金建铭. 1998.电磁场有限元方法.王建国译.西安: 西安电子科技大学出版社.
|