2017, Vol. 60
2. 中国科学院大学, 北京 100049
2. Graduate University, Chinese Academy of Sciences, Beijing 100049, China
大地电磁法(Magnetotelluric，MT)作为研究地球深部电性结构的地球物理方法，在对我国的岩石圈研究中发挥了重要作用(Wei et al.，2001;马晓冰等，2005; 汤吉等，2005; Bai et al.，2010;闫永利等，2012).近年来，随着计算机技术的发展和数据解释的需要等，三维大地电磁正反演问题成为了大地电磁法研究的热点(Siripunvaraporn et al.，2002; Ledo，2006; Börner，2010).为了快速精确的实现三维大地电磁正演模拟，国内外学者进行了大量的研究(Avdeev，2005)，采用的数值方法有积分方程法(Hohmann，1975)、有限单元法(徐世浙，1994)和有限差分法(Lines and Jones，1973)等，其中有限差分法凭借计算快速、适应性强等优点，被大量的应用于三维大地电磁正反演算法中(Siripunvaraporn et al.，2002; Miensopust et al.，2013).Mackie等(1993，1994)在传统有限差分法的基础上，通过改进电场和磁场的分布，提出了三维大地电磁交错网格有限差分方法，提高了计算精度.同时为了加快求解速度，Mackie等(1993) 和Smith(1996) 引入了电磁场散度校正方法，研究表明该方法可有效的提高三维电磁场正演模拟的计算速度和精度(陈辉等，2011; Farquharson and Miensopust，2011).谭捍东等(2003) 较系统的论述了三维大地电磁有限差分算法的实现过程，并详细的介绍了电磁场的边界条件，但为了满足截断边界条件，需要将外边界设置在足够远处，以此来避免截断边界的影响，从而势必增加了正演的计算量和计算时间.因此，在保证计算精度的前提下，改进大地电磁场正演模拟的边界条件，压缩数值模拟求解空间，是提高三维大地电磁正演效率的有效途径.
针对电磁场正演模拟中的截断边界问题，一些学者研究了多种吸收边界条件理论，其中最佳匹配层(Perfectly Matched Layer，PML)吸收边界条件应用比较广泛(李展辉和黄清华，2014).1994年，Berenger提出了 PML吸收边界条件理论，Chew和Weedon(1994) 推导了Berenger 的PML吸收边界条件下的Maxwell修正方程，随后，Katz等(1994) 将其应用于三维散射电磁场数值计算中，并对比了Mur吸收边界和PML吸收边界，显示PML吸收边界具有更好的吸收效果.Alumbaugh等(1996) 将Berenger的PML吸收边界条件引入到偶极子场源下的电磁场模拟中，取得了较好的效果，但PML吸收边界条件的加入增大了系数矩阵的条件数，进而大幅度的增加了迭代次数和求解时间.Sacks等(1995) 研究Berenger的PML吸收边界条件后，提出了各向异性PML吸收边界条件，各向异性PML吸收边界条件不需要对Maxwell方程进行修改，且更适用于频率域电磁场的计算.Gedney(1996) 对比分析Berenger的PML和各向异性PML吸收边界条件，认为各向异性PML吸收边界条件具有更高的计算效率和稳定性，但是Wu等(1997) 在计算各向异性PML吸收边界条件下的电磁场时，结果显示各向异性PML边界条件仍较大幅度的增加了迭代次数.所以在应用PML吸收边界条件时，需要克服其对系数矩阵条件数的影响.
为了进一步提高三维大地电磁正演的计算效率，本文研究将各向异性PML吸收边界条件耦合到基于二次场的三维大地电磁有限差分正演模拟方法中.文章首先分析了三维大地电磁二次场法理论及其有限差分格式，然后讨论了PML吸收边界条件在大地电磁正演模拟中的设置和散度校正等，最后应用二次场算法正演模拟两个三维模型，分析计算结果，总结二次场法和PML吸收边界条件的适用性和计算经验等.
2 大地电磁二次场有限差分法考虑到大地电磁法(MT)研究的频率范围，忽略位移电流，设时谐因子为eiωt，则频率域Maxwell方程组形式为
|
(1) |
|
(2) |
式中ω为角频率，μ和σ分别为磁导率和电导率，E和H分别为电场强度和磁场强度.
根据叠加原理，可将电磁总场分解为一次场(Primary field)和二次场(Secondary field)之和，公式为
|
(2) |
其中，大地电磁一次场可以通过解析式精确求得，二次场则需要利用数值方法进行求解.
将式(3) 代入方程(1) ，得到大地电磁二次场控制方程为
|
(4) |
|
(5) |
式中，${\tilde{\mu }}$和${\tilde{\sigma }}$分别为有效磁导率和电导率，在研究区域内部，${\tilde{\mu }}$和${\tilde{\sigma }}$为内部介质磁导率μ和电导率σ，但在截断边界区域，${\tilde{\mu }}$和${\tilde{\sigma }}$则需根据边界条件的不同而进行相应的设置.Δσ=1/ρ－1/ρp，为异常体与背景介质电导率差值，ρ为介质电阻率，ρp是均匀半空间或层状介质背景电阻率.与总场方程(Mackie et al.，1993; 谭捍东等，2003)相比，二次场方程多了电场场源项，电磁场变为由内部产生并向四周传播的散射电磁场.
本文采用交错网格有限差分法求解二次场方程和，二次电磁场采样位置如图 1，以单元(i，j，k)为例，电场采样点位于网格单元表面的中心，磁场采样点位于网格棱边的中心上.首先对方程按电流密度Js三个分量x，y，z进行离散，公式为
|
(6) |
|
(7) |
|
(8) |
|
图 1 交错采样网格单元示意图 Fig. 1 Sketch of staggered grid |
其中电流密度Js=σEs，右端场源项ΔJs与一次场Ep离散关系为
|
其中，ρx，ρy，ρz和Δρx，Δρy，Δρz分别为x，y，z方向的电阻率和异常电阻率差.
同样的，对方程(4) 进行有限差分离散，公式为
|
(9) |
|
(10) |
|
(11) |
将式(6) —(11) 相互联立便可建立总数为Nx×(Ny－1) ×(Nz－1) +Ny×(Nx－1) ×(Nz－1) +Nz×(Ny－1) ×(Nx－1) 的二次磁场线性方程组.二次场有限差分线性方程中的各项系数与总场法完全相同(谭捍东等，2003)，只是在方程右端多了二次场源项.
3 PML吸收边界条件在有限的区域数值模拟三维电磁场时，必须在计算区域的截断边界处给出边界条件，边界条件的处理是影响电磁场计算精度和效率的重要因素(Mackie et al.，1993).如图 2，在三维大地电磁有限差分数值模拟中，存在6个外边界区域，顶部为空气边界，下部是地下底边界，和四个方向的侧边界.在总场法中，6个外边界的边界场值各不相同，在数值计算中需要对6个边界进行不同的处理，如在四个方向侧边界上加载对应的电磁场值(谭捍东等，2003).而在二次场法中，当6个边界离二次场源足够远时，大地电磁二次场在边界处场值衰减至零，因此二次场法的边界条件比较简洁.但是不管总场法还是二次场法，为了满足边界条件，需要将边界设置在“足够远”处，以此避免截断边界的影响，都需要比较大的计算区域.
|
图 2 三维PML边界分布 Fig. 2 3D boundaries of the PML region |
为了有效减小三维大地电磁数值模拟的求解空间，本文引入各向异性PML吸收边界条件作为二次场边界条件.根据各向异性PML吸收边界条件理论(Sacks et al.，1995; Gedney，1996)，在三维计算区域的6个截断边界处设置PML吸收边界层(图 2)，当二次场传播到截断边界处时，如图 3，为了使得二次场无反射的进入PML吸收边界层，边界层介质的电导率${\tilde{\sigma }}$和磁导率${\tilde{\mu }}$需满足公式为(Wu et al.，1997)
|
(12) |
其中Λ为匹配矩阵，si(i=x，y，z)为不同方向的PML吸收因子.吸收因子s的取值可以为复数，但本文考虑到大地电磁的研究频率和大量试验计算，建议取吸收因子s为实数.
PML吸收边界条件通过对不同方向的吸收因子s进行取值，使得进入PML边界的电磁波快速衰减，最终减小或避免截断边界的影响(Wu et al.，1997).三维电磁数值模拟中存在三个方向的截断边界，每个方向有两个 PML平面区Λside，不同方向的PML平面区Λside吸收因子s的取值见表 1.如图 3，当二次场进入垂直于y轴的截断边界区域时，为使得二次场无反射的进入PML层并快速衰减，区域内的有效电导率和磁导率须为
|
对于两个或三个PML边界层相交的区域Λedge和Λcorner，只需对应的平面区匹配矩阵相乘即可(葛德彪，闫玉波，2002).
|
表 1 不同PML平面区的参数取值 Table 1 Absorbing factor values in different areas of PML |
|
图 3 二次场入射PML边界示意图 Fig. 3 Sketch of secondary wave incident on a boundary interface |
在实际计算中，一般需要设置多层PML边界，对于每层PML边界的吸收因子设置存在两种方式(Alumbaugh et al.，1996)，一种是每一层PML吸收因子都固定不变，一种是逐步增大吸收因子值.Berenger(1994) 和Gedney(1996) 研究发现，在每一层PML中设置相同的吸收因子s容易产生伪解，建议逐步增大s值，因此本文采用逐步增大方式，在PML边界层中逐步增大吸收因子s值，如沿着x轴方向，每层吸收因子s(x)值为
|
(13) |
式中，x0为交界面，d是PML厚度，smax是最大吸收因子.理论上，不同方向的吸收因子s值可以各不相同，但为了保持线性方程矩阵的对称和防止条件数增大(Alumbaugh et al.，1996)，smax一般作为全局的最大吸收因子.
4 求解器和散度校正加入边界条件后，便可以形成一个完整的大型线性方程组为
|
(14) |
其中，A是大型稀疏系数矩阵，f为待求解的二次磁场三分量，b是二次场场源.
对于方程(14) 的求解，目前存在两种应用较广泛的求解器：双共轭梯度求解器(BICG)和准最小残差求解器(QMR)(Alumbaugh et al.，1996).Alumbaugh等(1996) 、沈金松(2003) 和Puzyrev等(2013) 对比了这两种求解器，发现两种算法都可以得到较好的计算结果，其中BICG的迭代收敛较快，但是存在残差振荡现象，而QMR算法则更加稳定，更加适合求解复数方程，同时对内存的需求也小.本文由于采用了二次场法，方程与Alumbaugh等(1996) 和Puzyrev等(2013) 文章中的方程结构更加相似，而且引入的吸收边界条件增加了系数矩阵A的条件数和复杂性，因此采用QMR求解器求解方程，并加入不完全LU分解预处理器.
在二次场中加入PML吸收边界后，会增大系数矩阵的条件数，从而大量的增加迭代次数和求解时间(Alumbaugh et al.，1996; Wu et al.，1997).为了改善PML边界条件下方程系数矩阵A的条件数，提高迭代收敛速度，本文在QMR算法迭代求解中加入散度校正技术(Mackie et al.，1994; Smith，1996).由于引入了PML吸收边界条件，边界区域介质磁导率值变为张量，因此不能只对磁场强度Hs散度进行校正，而需要对磁感应矢量Bs进行散度校正，即：
|
(15) |
求解方程组(14) 获得大地电磁二次场，与一次场相加后得到三维大地电磁场，继而计算出地表电磁场值和阻抗张量响应.
5 模型计算分析为了验证基于二次场有限差分算法(Secondary Field Method，SFM)的正确性，以及耦合各向异性PML吸收边界条件的二次场算法SFM-PML有效性和计算效率，本文计算了两个三维地电模型：低阻异常模型TM1和都柏林测试模型TM2，分析二次场法和PML吸收边界条件的有效性和适用性等.本文算法的计算环境是至强四核3.40 GHz处理器、8G内存、64位Windows7系统计算机.
5.1 低阻异常模型低阻异常模型TM1是在100 Ωm的均匀半空间中赋存一个1 Ωm低阻异常体(图 4)，低阻异常体尺寸为5 km×5 km×10 km，埋深为5 km.分别应用基于总场的有限差分算法mackie(Mackie et al.，1994)和二次场算法SFM进行计算，采用相同的网格剖分70×60×34，截断边界均远离异常体.图 5展示了沿着y=0测线、周期100 s下的两种算法的计算结果，可以看出，两种算法的结果吻合的较好，视电阻率均正确的反映了地下的电性变化.
|
图 4 模型TM1 Fig. 4 3D model TM1 |
|
图 5 模型TM1在100 s周期下沿测线y=0的视电阻率 (a)XY模式；(b)YX模式. Fig. 5 Apparent resistivity along y=0 under 100 s period of model TM1 (a)XY mode；(b)YX mode. |
在以上二次场算法正演模拟网格模型的基础上，大幅减小三个方向外边界距离，压缩模型TM1的计算区域，剖分网格数由70×60×34减为42×32×25，此时，x和y方向的计算范围缩小为－6.5~6.5 km.然后在三个方向截断边界附近设置两层网格为PML吸收边界层，应用算法SFM-PML进行计算.图 5中展示了小边界情况下，吸收因子smax=15时的计算结果和不加PML边界条件(no PML)的计算结果，可以看出，在小边界情况下即截断边界距离异常体较近时，正演响应中存在较严重的截断边界效应.而加载了PML吸收边界条件后，则较好的减小了截断边界的影响，说明了PML吸收边界条件的有效性.
图 6是不同吸收因子smax下沿着测线y=0的视电阻率误差分布情况，随着吸收因子的变大，两种模式的视电阻率误差快速降低，当吸收因子smax=15时，误差全部降至2%以下，说明调节吸收因子smax的取值，可以大幅度的减小截断边界的影响.但当吸收因子继续增大至20，视电阻率Rxy误差部分开始出现增大现象，分析认为过大的吸收因子引起了离散误差的增大，因此吸收因子取值不宜过大.
|
图 6 不同吸收因子下沿测线y=0周期100 s的视电阻率误差的变化 (a)XY模式；(b)YX模式. Fig. 6 Error variation with different absorbing factors along y=0 (a)XY mode；(b)YX mode. |
都柏林测试模型TM2是一个在电阻率100 Ωm的均匀半空间中赋存三个异常体的三维地电模型，各异常体的电性与空间分布见图 7.图 8是坐标中心原点x=y=0处大地电磁测深曲线，包括基于总场的有限单元算法Han/Lee、有限差分算法Miensopust计算结果(Miensopust et al.，2013)和二次场法SFM计算结果.其中有限单元算法Han/Lee的网格剖分数为48×47×31，有限差分算法Miensopust网格剖分数85×62×44，二次场算法SFM网格剖分数60×57×40，所有算法计算0.1 s～10000 s范围内21个周期的电磁响应，外边界设置在足够远处.从图 7可以看出，二次场法SFM计算结果与其他两种算法的计算结果吻合的较好，进一步说明了二次场算法SFM的可靠性.
|
图 7 都柏林测试模型TM2 Fig. 7 Three dimensional model TM2 |
|
图 8 模型TM2中心点x=y=0视电阻率和相位曲线 (a)XY模式；(b)YX模式. Fig. 8 Apparent resistivity and phase curves of model TM2 at XY0 site x=y=0 km (a)XY mode；(b)YX mode. |
同样的，大幅收缩二次场算法SFM正演模型的外边界距离，减小三维地电模型TM2的计算空间，x和y方向的计算范围分别压缩为-24～24 km和-26.5～26.5 km.然后在截断边界附近设置厚度为3.5 km的PML吸收边界层，剖分网格数由二次场算法SFM计算时的60×57×40减少为40×40×30，应用算法SFM-PML计算小边界情况下的模型TM2响应，以二次场算法SFM的计算结果作为参考结果进行对比分析.
图 9和图 10分别是中心测点与边界附近测点的测深曲线，从图中可以看出，随着周期的变大，截断边界开始影响计算结果，尤其随着测点越靠近边界，截断边界的影响越大(图中no PML)，而PML吸收边界条件则有效的减小了截断边界的影响，当吸收因子smax=5时，视电阻率和相位有了较大幅度的改善，继续增大吸收因子至smax=15，计算结果则进一步逼近参考结果，上述的计算比较进一步说明了耦合PML吸收边界的二次场算法SFM-PML的有效性.
|
图 9 不同吸收因子下模型TM2中心点x=y=0视电阻率和相位曲线 (a)XY模式；(b)YX模式. Fig. 9 Apparent resistivity and phase curves of model TM2 at x=y=0 km (a)XY mode；(b)YX mode. |
|
图 10 不同吸收因子下模型TM2边界附近测点的视电阻率和相位曲线 (a)和(b)分别为边界处测点x=0，y=-22.5 km的XY和YX模式计算结果，(c)和(d)分别为边界处测点x=20 km，y=0的XY和YX模式计算结果. Fig. 10 Apparent resistivity and phase curves of model TM2 at boundary sites (a)and(b)XY and YX mode results of boundary site x=0，y=-22.5 km, (c)and(d)XY and YX mode results of boundary site x=20 km, y=0. |
Alumbaugh等(1996) 研究发现，PML吸收边界条件的加入，会引起线性方程系数矩阵条件数的增大，从而可能大幅度的增加迭代求解的难度.图 11为不同吸收因子下的求解器QMR的迭代次数，可以看出，虽然加入PML吸收边界条件后，总体趋势上，随着吸收因子的增大，迭代次数也逐步增加，但迭代求解次数并没有大幅度的上升，而且迭代次数总体上仍小于参考结果SFM计算的迭代次数.表 2是不同吸收因子下模型TM2正演计算21个周期所用的计算时间，也可以看出PML吸收边界条件的加入增加了计算时间，但是增加幅度有限，仍远小于大边界情况下参考结果SFM的正演计算时间.PML吸收边界条件的加入并未引起算法SFM-PML正演模拟迭代次数和计算时间的大幅增加，分析认为考虑PML参数的磁感应矢量散度校正技术有效的改善了系数矩阵A的条件数，避免了PML吸收边界条件下的方程求解迭代的大量增大.
|
图 11 不同PML吸收因子下模型TM2的方程迭代次数 Fig. 11 Iterations numbers for different frequencies and PML absorbing factors of model TM2 |
|
|
表 2 不同边界条件下21个周期的模型TM2正演计算时间 Table 2 Calculation time of model TM2 with different boundary conditions |
图 12是吸收因子smax=15下不同频率的模型TM2视电阻率Rxy剖面以及对应的参考结果SFM，其中图中每个剖面中的方框内部是研究区，外部为截断边界区域.从图 12可见，在PML吸收边界区域，电阻率变化剧烈，而在研究区域内部，从高频到低频，PML吸收边界条件下的视电阻率Rxy剖面和参考结果REF基本一致.综合以上的计算分析，本文认为耦合PML吸收边界条件的三维大地电磁二次场有限差分算法具有较高的计算效率和精度，可以有效的提高了三维大地电磁求解的效率.
|
图 12 不同周期的模型TM2视电阻率Rxy剖面对比 (a)PML吸收边界条件下计算结果;(b)参考结果SFM. Fig. 12 Comparison of resistivity profiles of model TM2 with PML boundary conditions and REF profiles (a)PML boundary conditions;(b)SFM. |
除了上述的计算分析，本文又进行了大量的不同背景电阻率、不同异常体、不同PML边界厚度和不同吸收因子等的正演模拟，正演结果显示耦合PML吸收边界条件的三维大地电磁二次场有限差分算法在计算精度和效率上具有普遍适用性.一般，异常体与背景电阻率差异越小、异常体越小、截断边界放置的越远，则所需的吸收因子越小.对于PML吸收因子在电磁场正演模拟中合适的取值，需要综合考虑多种因素，除了需要考虑异常体与背景电阻率差异、截断边界放置的远近等，还需要考虑边界网格的大小、边界网格的变化幅度和正演误差容忍度等.
6 结论本文研究分析了三维大地电磁二次场法理论及其有限差分法离散格式，引入各向异性最佳匹配层(PML)吸收边界条件作为二次场数值计算的边界条件，讨论了PML吸收边界的设置，然后应用QMR求解器和磁感应矢量散度校正技术对线性方程进行迭代求解，实现了耦合PML吸收边界条件的基于二次场的三维大地电磁有限差分算法.最后对模型的正演模拟结果进行分析，得出以下结论：
(1) 三维大地电磁二次场有限差分算法计算结果与其他算法均吻合的较好，具有较高的精确度和可靠性.
(2) 耦合PML吸收边界条件的二次场算法可以有效的减小截断边界的影响，大幅压缩模型计算的求解空间，提高了三维大地电磁正演模拟的效率.
(3) PML吸收因子对正演计算的误差和线性方程的迭代求解具有重要影响，PML吸收因子的取值需要综合考量，不宜过大.
致谢非常感谢审稿专家和编辑提出的宝贵意见和建议，对文章的修改起到了很大的帮助.
| Alumbaugh D L, Newman G A, Prevost L, et al. 1996. Three-dimensional wideband electromagnetic modeling on massively parallel computers. Radio Science, 31(1): 1-23. DOI:10.1029/95RS02815 | |
| Avdeev D B. 2005. Three-dimensional electromagnetic modelling and inversion from theory to application. Surv. Geophys., 26(6): 767-799. DOI:10.1007/s10712-005-1836-x | |
| Bai D H, Unsworth M J, Meju M A, et al. 2010. Crustal deformation of the eastern Tibetan plateau revealed by magnetotelluric imaging. Nat. Geosci., 3(5): 358-362. DOI:10.1038/ngeo830 | |
| Berenger J P. 1994. A perfectly matched layer for the absorption of electromagnetic waves. Journal of Computational Physics, 114(2): 185-200. DOI:10.1006/jcph.1994.1159 | |
| Börner R U. 2010. Numerical modelling in geo-electromagnetics:advances and challenges. Surv. Geophys., 31(2): 225-245. DOI:10.1007/s10712-009-9087-x | |
| Chen H, Deng J Z, Tan H D, et al. 2011. Study on divergence correction method in three-dimensional magnetotelluric modeling with staggered-grid finite difference method. Chinese J. Geophys. (in Chinese), 54(6): 1649-1659. | |
| Chew W C, Weedon W H. 1994. A 3D perfectly matched medium from modified maxwell's equations with stretched coordinates. Microwave and Optical Technology Letters, 7(13): 599-604. DOI:10.1002/(ISSN)1098-2760 | |
| Farquharson C G, Miensopust M P. 2011. Three-dimensional finite-element modelling of magnetotelluric data with a divergence correction. Journal of Applied Geophysics, 75(4): 699-710. DOI:10.1016/j.jappgeo.2011.09.025 | |
| Ge D B, Yan Y B. Electromagnetic Wave and FDTD (in Chinese).Xi'an: Xidian University Press, 2002. | |
| Gedney S D. 1996. An anisotropic perfectly matched layer-absorbing medium for the truncation of FDTD Lattices. IEEE Transactions on Antennas and Propagation, 44(12): 1630-1639. DOI:10.1109/8.546249 | |
| Hohmann G W. 1975. Three-dimensional induced polarization and electromagnetic modeling. Geophysics, 40(2): 309-324. DOI:10.1190/1.1440527 | |
| Katz D S, Thiele E T, Taflove A. 1994. Validation and extension to three dimensions of the Berenger PML absorbing boundary condition for FD-TD meshes. IEEE Microwave and Guided Wave Letters, 4(8): 268-270. DOI:10.1109/75.311494 | |
| Ledo J. 2006. 2-D versus 3-D magnetotelluric data interpretation. Surv. Geophys., 27(1): 111-148. DOI:10.1007/s10712-006-0002-4 | |
| Li Z H, Huang Q H. 2014. Application of the complex frequency shifted perfectly matched layer absorbing boundary conditions in transient electromagnetic method modeling. Chinese J. Geophys. (in Chinese), 57(4): 1292-1299. DOI:10.6038/cjg20140426 | |
| Lines L R, Jones F W. 1973. The perturbation of alternating geomagnetic fields by an island near a coastline. Canadian Journal of Earth Sciences, 10(4): 510-518. DOI:10.1139/e73-050 | |
| Ma X B, Kong X R, Liu H B, et al. 2005. The electrical structure of northeastern Qinghai-Tibet Plateau. Chinese J. Geophys. (in Chinese), 48(3): 689-697. | |
| Mackie R L, Madden T R, Wannamaker P E. 1993. Three-dimensional magnetotelluric modeling using difference equations-Theory and comparisons to integral equation solutions. Geophysics, 58(2): 215-226. DOI:10.1190/1.1443407 | |
| Mackie R L, Smith J T, Madden T R. 1994. Three-dimensional electromagnetic modeling using finite difference equations:the magnetotelluric example. Radio Sci., 29(4): 923-935. DOI:10.1029/94RS00326 | |
| Miensopust M P, Queralt P, Jones A G, et al. 2013. Magnetotelluric 3-D inversion-a review of two successful workshops on forward and inversion code testing and comparison. Geophys. J. Int., 193(3): 1216-1238. DOI:10.1093/gji/ggt066 | |
| Puzyrev V, Koldan J, de la Puente J, et al. 2013. A parallel finite-element method for three-dimensional controlled-source electromagnetic forward modelling. Geophys. J. Int., 193(2): 678-693. DOI:10.1093/gji/ggt027 | |
| Sacks Z S, Kingsland D M, Lee R, et al. 1995. A perfectly matched anisotropic absorber for use as an absorbing boundary condition. IEEE Transactions on Antennas and Propagation, 43(12): 1460-1463. DOI:10.1109/8.477075 | |
| Shen J S. 2003. Modeling of 3-D electromagnetic responses in frequency domain by using staggered-grid finite difference method. Chinese J. Geophys. (in Chinese), 46(2): 281-288. | |
| Siripunvaraporn W, Egbert G, Lenbury Y. 2002. Numerical accuracy of magnetotelluric modeling:A comparison of finite difference approximations. Earth Planets Space, 54(6): 721-725. DOI:10.1186/BF03351724 | |
| Smith J T. 1996. Conservative modeling of 3-D electromagnetic fields, Part II:biconjugate gradient solution and an accelerator. Geophysics, 61(5): 1319-1324. DOI:10.1190/1.1444055 | |
| Tan H D, Yu Q F, Booker J, et al. 2003. Magnetotelluric three-dimensional modeling using the staggered-grid finite difference method. Chinese J. Geophys. (in Chinese), 46(5): 705-711. | |
| Tang J, Zhan Y, Zhao G Z, et al. 2005. Electrical conductivity structure of the crust and upper mantle in the northeastern margin of the Qinghai-Tibet plateau along the profile Maqên-Lanzhou-Jingbian. Chinese J. Geophys. (in Chinese), 48(5): 1205-1216. | |
| Wei W B, Unsworth M, Jones A, et al. 2001. Detection of widespread fluids in the Tibetan crust by magnetotelluric studies. Science, 292(5517): 716-718. DOI:10.1126/science.1010580 | |
| Wu J Y, Kingsland D M, Lee J F, et al. 1997. A comparison of anisotropic PML to Berenger's PML and its application to the finite-element method for EM scattering. IEEE Transactions on Antennas and Propagation, 45(1): 40-50. DOI:10.1109/8.554239 | |
| Xu S Z. The Finite Element Method in Geophysics (in Chinese).Beijing: Science Press, 1994. | |
| Yan Y L, Ma X B, Chen Y, et al. 2012. The study of magnetotelluric sounding on Coqên-Xainza profile in Tibet. Chinese J. Geophys. (in Chinese), 55(8): 2636-2642. DOI:10.6038/j.issn.0001-5733.2012.08.015 | |
| 陈辉, 邓居智, 谭捍东, 等. 2011. 大地电磁三维交错网格有限差分数值模拟中的散度校正方法研究. 地球物理学报, 54(6): 1649–1659. | |
| 葛德彪, 闫玉波. 电磁波时域有限差分方法.西安: 西安电子科技大学出版社, 2002. | |
| 李展辉, 黄清华. 2014. 复频率参数完全匹配层吸收边界在瞬变电磁法正演中的应用. 地球物理学报, 57(4): 1292–1299. DOI:10.6038/cjg20140426 | |
| 马晓冰, 孔祥儒, 刘宏兵, 等. 2005. 青藏高原东北部地区地壳电性结构特征. 地球物理学报, 48(3): 689–697. | |
| 沈金松. 2003. 用交错网格有限差分法计算三维频率域电磁响应. 地球物理学报, 46(2): 281–288. | |
| 谭捍东, 余钦范, BookerJ, 等. 2003. 大地电磁法三维交错采样有限差分数值模拟. 地球物理学报, 46(5): 705–711. | |
| 汤吉, 詹艳, 赵国泽, 等. 2005. 青藏高原东北缘玛沁-兰州-靖边剖面地壳上地幔电性结构研究. 地球物理学报, 48(5): 1205–1216. | |
| 徐世浙. 地球物理中的有限单元法.北京: 科学出版社, 1994. | |
| 闫永利, 马晓冰, 陈赟, 等. 2012. 西藏错勤-申扎剖面大地电磁测深研究. 地球物理学报, 55(8): 2636–2642. DOI:10.6038/j.issn.0001-5733.2012.08.015 | |