地球物理学报  2018, Vol. 61 Issue (8): 3488-3498   PDF    
任意各向异性介质中三维可控源音频大地电磁正演模拟
邱长凯1, 殷长春1, 刘云鹤1, 陈辉1,2, 刘玲1, 蔡晶1     
1. 吉林大学地球探测科学与技术学院, 长春 130026;
2. 东华理工大学放射性地质与勘探技术国防重点学科实验室, 南昌 330013
摘要:在一些地层层理发育的地区,地下介质存在显著的电各向异性,此时基于各向同性模型解释含各向异性效应的可控源音频大地电磁(CSAMT)测深观测数据会导致错误的结果.本文通过引入3×3的对称正定张量表征电导率各向异性,采用非结构四面体网格和矢量有限元方法离散电场满足的矢量Helmholtz方程,并将电磁场源等效为系列电偶极子,实现任意各向异性介质中CSAMT高效数值模拟.本文首先通过层状各向异性模型检验三维有限元算法的精度和有效性,进一步建立三维地电模型研究异常体各向异性和围岩各向异性对CSAMT响应的影响,最后使用视电阻率极性图来识别各向异性电导率主轴方向.数值模拟结果表明,各向异性电导率对CSAMT视电阻率幅值及分布规律都有很大影响,视电阻率极性图能够很好地识别各向异性主轴方向.
关键词: 可控源音频大地电磁      各向异性      矢量有限元      三维正演     
3D forward modeling of controlled-source audio-frequency magnetotellurics in arbitrarily anisotropic media
QIU ChangKai1, YIN ChangChun1, LIU YunHe1, CHEN Hui1,2, LIU Ling1, CAI Jing1     
1. College of Geo-exploration Science and Technology, Jilin University, Changchun 130026, China;
2. Key Laboratory of Radioactive Geology and Exploration Technology Fundamental Science for National Defense, East China Institute of Technology, Nanchang 330013, China
Abstract: Controlled-source audio-frequency magnetotellurics (CSAMT) is playing an important role in shallow subsurface exploration. Most current three-dimensional (3D) forward modeling for CSAMT are based on an isotropic earth model which neglects the electrical anisotropy in the earth. Here we implement a 3D vector finite-element modeling algorithm to solve the frequency-domain Helmholtz equation for the total electrical field and analyze the effect of earth anisotropy on tensor CSAMT apparent resistivity responses. First we introduce a 3×3 positive-definite conductivity tensor to represent the arbitrarily electrical anisotropy. Then we discretize the weak-form of vector Helmholtz equation with an unstructured tetrahedral grid and lowest order of Nédélec vector basis functions. Besides, we approximate the horizontal transmitting sources for CSAMT by a series of electrical dipoles and solve the final system of equations with the direct solver MUMPS that can speed up the forward modeling while guaranteeing the numerical precision. We compare our 3D numerical results with one-dimensional solutions for a layered anisotropic earth to verify the accuracy and reliability of the proposed 3D anisotropic finite element algorithm. To further explore the anisotropic effect of a 3D anomalous body and host rock, we design a conductive brick model embedded in a homogeneous half-space. Then, we rotate the principal conductivity tensor of the anomalous body and host rock around x-axis and z-axis separately, and simulate the apparent resistivity for tensor CSAMT. Numerical results show that both the amplitude and distribution pattern of the CSAMT apparent resistivity vary with the anisotropic conductivity tensors. To recognize the electrical anisotropy of the earth, we show the CSAMT apparent resistivity in polar plots, thus one can clearly identify the principal axis of the anisotropic conductivity.
Key words: CSAMT    Anisotropy    Vector finite-element    3D forward modeling    
0 引言

可控源音频大地电磁测深法(Controlled-source audio-frequency magnetotellurics,缩写为CSAMT)是在大地电磁法和音频大地电磁法基础上发展起来的一种人工场源频率域测深方法.该方法具有抗干扰能力强、勘探深度范围大、观测效率高等特点(汤井田和何继善,2005),在石油、天然气、地热、固体矿产以及水文和环境工程勘查等领域得到了越来越广泛的应用.

由于CSAMT张量场源复杂的三维特性,为提高数据处理质量和解释水平,迫切需要开展CSAMT多维正反演研究.底青云等(2004)结合直立异常、倾斜异常和断陷模型研究了复杂介质2.5维CSAMT有限元法正演模拟.雷达(2010)使用三角网格剖分建立了起伏地形条件下基于Occam方法的CSAMT二维反演,有效地消除了地形对CSAMT的影响.邓居智等(2011)采用三维交错网格有限差分求解二次电场方程,实现CSAMT三维正演快速迭代求解.韩波等(2015)基于并行化直接求解器MUMPS和交错网格有限体积法实现频率域可控源三维正演,证明了直接求解对于模拟多场源可控源问题的有效性.Hu等(2015)基于六面体网格的矢量有限元法研究了带激电效应的三维张量CSAMT正演,并分析了典型模型的阻抗和视电阻率响应特征.王亚璐等(2017)采用非结构网格离散异常电场方程,引入散度条件约束实现CSAMT有限元数值模拟.

传统的CSAMT正演模拟主要基于各向同性模型,各向异性电阻率的CSAMT问题研究相对较少,大多局限于一维正反演.Li和Pedersen(1992)使用矢量势方法计算层状各向异性介质中的CSAMT响应,数值模拟结果表明即使是适中的各向异性也会产生和各向同性明显不同的结果.Yin和Maurer(2001)研究了CSAMT层状介质任意各向异性正演,分析了各向异性对CSAMT电场和视电阻率的影响.Li等(2000)使用非线性最小二乘法实现层状各向异性地层中CSAMT反演,实测数据反演表明各向异性模型比各向同性模型能更好地拟合数据.Zhou等(2014)研究了层状横向各向同性介质中的CSAMT数据正则化反演,发现CSAMT数据对于横向电阻率更敏感.曹强等(2016)运用Occam方法实现直立各向异性层状介质中CSAMT数据反演.三维正演方面,陈桂波(2009)利用积分方程算法研究了各向异性对CSAMT的场源效应与静态效应的影响.

考虑到在地层层理发育地区各向异性现象非常普遍,人们逐渐认识到电各向异性对大地电磁、海洋电磁和航空电磁等存在很大的影响并积极开展相关领域各向异性模拟研究(Yin and Fraser, 2004; Yin, 2006; Liu and Yin, 2014).相比之下,三维CSAMT各向异性正反演研究相对比较薄弱,各向异性对三维CSAMT的影响尚没有明确的认识.为此,本文研究各向异性介质中三维CSAMT正演模拟算法,通过引入3×3的对称正定矩阵表征各向异性电导率张量,使用基于非结构四面体网格的矢量有限元法离散电场Helmholtz方程,结合矩阵直接求解技术实现三维任意各向异性介质中CSAMT正演模拟.本文首先使用层状各向异性模型验证有限元算法的计算精度,然后分析各向异性异常体和各向异性围岩的CSAMT响应特征,最后提出使用视电阻率极性图识别地下介质各向异性电阻率的主轴方向.

1 控制方程和有限元离散

从微分形式的麦克斯韦方程组出发,取时谐因子为eiωt,在有源区域电磁场满足的方程为(Jin, 2014)

(1)

(2)

(3)

式(1)—(3)中,E表示电场强度,H表示磁场强度,Js表示外加电流源项,J表示总电流密度,i是虚数单位,ω=2πf是角频率,f是频率,σ是各向异性的电导率张量,μ是磁导率,ε是介电常数.为了将电导率张量从参考坐标系变换到计算坐标系,对主电导率张量σr进行三次欧拉旋转得到任意各向异性介质的电导率张量(Yin, 2000),即

(4)

其中主电导率张量σr

(5)

式中σxσyσz称为主电导率,ρxρyρz称为主电阻率.(4)式中D=DxDyDzDxDyDz分别为绕不同坐标轴的旋转矩阵(Yin, 2000).

由(1)和(2)式,可得在整个计算区域Ω上电场强度E满足的方程为

(6)

式中,上式即为电场总场Helmholtz方程.

为了保证电磁场存在唯一解,假设计算区域Ω的外边界Γ距离源足够远,电场强度满足Dirichlet边界条件:

(7)

其中n是外边界面上单位法向量.

定义为矢量Hilbert空间,L2是空间内平方可积的函数空间,空间内的内积(Ren et al., 2013)定义为

(8)

选取矢量测试函数VH(curl, Ω),对(6)式的第一项运用第一矢量格林定理并进行分部积分,其弱形式(Ren et al., 2013, 2014)表述为

(9)

其中,

(10)

(11)

本文运用矢量有限元法求取电场强度E.为此,将整个计算区域Ω剖分成四面体单元组成的网格NT是四面体单元的个数.选用一阶Nédélec矢量形函数(Nédélec, 1980)近似每个单元内的电场,即

(12)

其中,是单元内第j条棱边的切向电场,Nj(r)是第j个棱边的矢量形函数.

选取矢量测试函数V为矢量形函数N,将(12)式代入(10)和(11)式得到单元矩阵:

(13)

其中,

(14)

将(13)和(14)式中的单元矩阵组装可得到如下大型对称稀疏复线性代数方程组:

(15)

由于张量CSAMT使用两个分离或重叠场源,需要求解两个极化方向的电磁场,直接求解仅需对系数矩阵K进行一次矩阵分解,然后回代两个极化方向的源项即可.因此,本文采用多波前直接求解器MUMPS(Amestoy et al., 2006)对方程组(15)进行求解,得到棱边上的电场,进而利用(12)式进行插值求出观测点处的,并利用(1)式计算观测点处的磁场.

(14)式中的单元系数矩阵Ke的计算可参考(Jin, 2014).为计算使用总场方法时的单元源项矩阵Se,首先将水平导线发射源置于四面体网格中,计算每段导线与每个四面体单元的交点,并将各单元中所包含的线源部分近似为一个电偶极子.假设电偶极子的中点位置位于(x0, y0, z0),则通过单元的电流密度Jse可表示为(齐彦福等,2017)

(16)

式中Ie表示穿过该单元的源电流,dl表示该段偶极子长度,δ为Dirac-δ函数.由此,(14)式中的源项矩阵可简化为

(17)

对于任意形状、姿态和位置的电性源,均可利用电偶极子近似,进而由(17)式计算得到源项,(17)式中是矢量点乘,需要计算每个矢量形函数和外加电流源的点乘积.

为使单元电流源项满足电偶极子假设,剖分四面体网格时,在源所在的线段上嵌入一系列等距离点,通过施加网格棱边长度约束实现源附近网格的细化,如图 1a所示.为保证测点附近电磁场计算精度,同时尽可能减少网格单元个数,对测点附近的单元施加棱边长度约束实现网格的局部细化,如图 1b所示.

图 1 源(a)和测点(b)处四面体网格示意图 Fig. 1 Local view of the tetrahedral mesh near the sources (a) and receivers (b)

张量CSAMT通常在远区观测,测区场源可近似为垂直入射的平面波.可根据两次不同极化方向的电场和磁场进行转换得到阻抗张量,阻抗与电场和磁场分量关系式(Yin and Maurer, 2001)为

(18)

(18)式中上标表示极化模式.由(18)式可得,阻抗张量为

(19)

式中,ζ=Hy(2)Hx(1)Hy(1)Hx(2).进一步根据阻抗张量可定义视电阻率和相位:

(20)

式中ij可分别表示xy.

2 算例分析 2.1 层状各向异性模型精度验证

首先,设计一个两层各向异性模型验证有限元算法的精度.如图 2所示,第一层为各向同性地层,电阻率为10 Ωm,层厚度为100 m;第二层为各向异性地层,xyz方向电阻率分别为200 Ωm、50 Ωm和50 Ωm.以坐标原点为中心沿x轴和y轴方向布设两个互相垂直的源,发射频率为0.1 Hz到10000 Hz,按对数等间隔取21个频点,测点位于(3000, 3000, 0)m处.

图 2 两层各向异性模型 第一层为各向同性地层; 第二层为方位各向异性地层. Fig. 2 A two-layered anisotropic earth model The top layer is isotropic, while the second layer is azimuthally anisotropic.

将计算结果与Yin和Maurer(2001)的一维层状各向异性半解析解比较,如图 3所示,本文有限元算法结果与Yin和Maurer(2001)一维结果吻合很好,对于等效视电阻率和相位(Yin和Maurer, 2001),两种算法的相对差均在4%以内,验证了本文提出的数值算法的可靠性.

图 3 两层各向异性模型本文算法与半解析解对比 (a)等效视电阻率和相位;(b)相对误差. Fig. 3 Comparison of 3D results from this paper with semi-analytical solutions for a two-layered anisotropic earth model (a) For effective apparent resistivity and phase; (b) For relative error.
2.2 各向异性异常体的CSAMT响应特征

下面研究各向同性半空间中存在各向异性异常体的张量CSAMT响应特征.均匀半空间电阻率为100 Ωm,其中含有一个大小为600 m×600 m×600 m的异常体,异常体中心位于(6500, 13500, 700)m处,其主轴电阻率ρxρyρz分别为1 Ωm、50 Ωm和1 Ωm.将异常体主电导率张量绕xz轴分别旋转0°、45°和90°,计算视电阻率ρxyaρyxa,并分析异常体各向异性对CSAMT视电阻率的影响.

作为对比,首先在图 4中给出电阻率分别为1 Ωm和50 Ωm的各向同性异常体对应的视电阻率响应.由于异常体电阻率较低,根据电流密度在介质分界面的连续性,异常体内电场幅值较小,导致图 4中视电阻率在低阻异常体上方出现极小值,清晰地指示了异常体的水平位置;且异常体电阻率值越小,视电阻率值越低.根据(19)式,ρxya主要由Ex决定,反映x方向电阻率变化.对于图 4a的情况,沿x方向的电流穿过电性分界面时,电流密度连续,而电导率存在很大差异,导致法向电场发生严重的不连续性,因此ρxya等值线沿x方向非常密集. ρyxa的分布特征可作出相似的解释,只是由于此时Ey起决定作用,ρyxa等值线分布特征与ρxya正交.

图 4 各向同性异常体的张量CSAMT视电阻率响应 (a)和(b)为ρxya; (c)和(d)为ρyxa. Fig. 4 Apparent resistivities of tensor CSAMT for an isotropic anomalous body embedded in an isotropic half-space (a) and (b) are for ρxya, while (c) and (d) are for ρyxa.
2.2.1 各向异性主轴绕x轴旋转

图 5给出各向异性异常体主电导率张量绕x轴旋转时的张量CSAMT视电阻率响应.从图 5(a-c)可以看出,随着异常体主电导率张量绕x轴旋转角α的增大,ρxya基本不变,且ρxya图 4a中电阻率为1 Ωm的各向同性模型结果基本一致.这是由于ρxya主要由x方向电场决定,异常体电阻率变化对电场的影响较大而对磁场的影响较小,当异常体主电导率张量绕x轴旋转时,主轴电阻率ρx为1 Ωm保持不变,不影响x方向电场分布,因此ρxya对异常体主电导率张量绕x轴旋转不敏感.从图 5(d-f)可以看出,异常体各向异性对ρyxa有着很大影响,随着旋转角α增大,ρyxa逐渐减小.当α=0°时,ρyxa的幅值和分布规律和图 4d中电阻率为50 Ωm的各向同性模型结果相似;当α=90°时,ρyxa的幅值和分布规律和图 4c中电阻率为1 Ωm的各向同性模型结果相似.这是由于异常体主电导率绕x轴旋转时,主轴电阻率ρy由50 Ωm减小到1 Ωm,进而ρyxa分布特征由50 Ωm的各向同性异常体向1 Ωm的各向同性异常体的视电阻率分布特征转变.

图 5 各向异性异常体主电导率张量绕x轴旋转时张量CSAMT视电阻率响应 (a)—(c) ρxya; (d)—(f) ρyxa. Fig. 5 Apparent resistivities of tensor CSAMT for the anisotropic conductivity tensor of an anomalous body rotated around x-axis. (a)—(c) are for ρxya, while (d)—(f) are for ρyxa.
2.2.2 各向异性主轴绕z轴旋转

图 6可以看出,异常体主电导率张量绕z轴旋转对ρxyaρyxa的幅值和分布规律都产生了很大影响.当旋转角χ=0°时,ρxya图 4a中电阻率为1 Ωm的各向同性模型结果基本一致,ρyxa图 4d中电阻率为50 Ωm的各向同性模型结果基本一致.当旋转角χ=90°时,ρxyaρyxa分别和图 4b中电阻率为50 Ωm和图 4c中电阻率为1 Ωm的各向同性模型结果基本一致.当旋转角在0°到90°(比如图中χ=45°)变化时,ρxyaρyxa分布特征均随旋转角度发生改变,然而对不同旋转角度的极性图研究发现,ρxyaρyxa极轴的平均方向指示了异常体各向异性主轴方向.

图 6 异常体主电导率张量绕z轴旋转时张量CSAMT视电阻率响应 (a)—(c) ρxya; (d)—(f) ρyxa. Fig. 6 Apparent resistivities of tensor CSAMT for the anisotropic conductivity tensor of an anomalous body rotated around z-axis (a)—(c) are for ρxya, while (d)—(f) are for ρyxa.
2.3 各向异性围岩的CSAMT响应特征

在层理明显发育地区,围岩也可能呈现各向异性.本节研究各向异性围岩对CSAMT响应的影响特征.模型几何参数与2.2节相同.假设异常体为各向同性、电阻率为1 Ωm;围岩为各向异性,其主轴电阻率ρxρyρz分别为50 Ωm、100 Ωm和50 Ωm.将围岩主电导率张量分别绕xz轴旋转0°、45°和90°,计算视电阻率ρxyaρyxa,并分析围岩不同各向异性时的CSAMT视电阻率响应.

2.3.1 各向异性主轴绕x轴旋转

图 7可以看出,低阻异常体正上方,ρxyaρyxa的形态和各向同性情况相似,分别在垂直极化电场方向(其中,ρxya在y方向,而ρyxax方向)被拉伸.从图 7(a-c)可进一步看出,随着围岩主电导率绕x轴旋转角α的增大,异常体上方视电阻率ρxya增大,围岩上方视电阻率ρxya为50 Ωm不变.这是因为当围岩主电导率张量绕x轴旋转时,围岩主轴电阻率ρx保持不变(50 Ωm),x方向电场分布受影响很小,因此ρxya对围岩主电导率张量绕x轴旋转不敏感;而视电阻率是真电阻率的综合反应,围岩主轴电阻率ρx由50 Ωm变为100 Ωm,导致视电阻率ρxya增大.相比之下,由图 7(d-f)可以看出,随着旋转角α的增大,围岩上方和异常体上方视电阻率ρyxa均减小.其中的物理原因与ρxya相似,这里不做赘述.

图 7 围岩主电导率张量绕x轴旋转时张量CSAMT视电阻率响应 (a)—(c) ρxya; (d)—(f) ρyxa. Fig. 7 Apparent resistivities of tensor CSAMT for the anisotropic conductivity tensor of the host rock rotated around x-axis (a)—(c) are for ρxya, while (d)—(f) are for ρyxa.
2.3.2 各向异性主轴绕z轴旋转

图 8(a-c)可以看出,随着围岩主电导率绕z轴旋转角χ的增大,异常体和围岩上方视电阻率ρxya均增大.这是因为当围岩主电导率张量绕z轴旋转时,围岩主轴电阻率ρx由50 Ωm变为100 Ωm,围岩上方电场增大,视电阻率ρxya增大;视电阻率是异常体和围岩真电阻率的综合反映,围岩主轴电阻率ρx增大,导致异常体上方视电阻率ρxya也增大.同理,从图 8(d-f)可以看出,随着旋转角χ的增大,围岩上方和异常体上方视电阻率ρyxa均减小.当旋转角为γ=45°时,ρxyaρyxa均沿异常中心发生了偏转,且偏转方向和绕z轴旋转方向相反.和前面讨论得出类似的结论,ρxyaρyxa极性轴的平均方向指示围岩主电导率的旋转方向.

图 8 围岩主电导率张量绕z轴旋转时张量CSAMT视电阻率响应 (a)—(c) ρxya; (d)—(f) ρyxa. Fig. 8 Apparent resistivities of tensor CSAMT for the anisotropic conductivity tensor of the host rock rotated around z-axis (a)—(c) are for ρxya; (d)—(f) are for ρyxa.
2.4 围岩和异常体各向异性识别

参考Yin(2006)关于海洋大地电磁的研究,我们将阻抗张量和视电阻率从直角坐标系转换到极坐标系,使用极坐标系下的视电阻率分布(极性图)来识别地下介质的各向异性电阻率主轴方向.设极轴φ=0°时与直角坐标系x轴方向重合,在φ=270°时与直角坐标系y轴方向重合,则极坐标系中的张量阻抗定义为(Yin, 2006)

(21)

式中,ZZφr等价于使用极坐标下径向电场和切向磁场或者切向电场和径向磁场计算的阻抗分量.考虑到ZZφr的正交性,这里仅给出ρa的计算结果.

图 9展示针对上节三维各向异性异常体模型计算得到的视电阻率ρa.从图 9a可以看出,主电导率张量绕x轴由0°旋转到90°时,主轴电阻率ρx保持不变(1 Ωm),因此ρa短轴位于x方向,且沿该方向ρa保持不变.相比之下,随着旋转角α增大,主轴电阻率ρy由50 Ωm变为1 Ωm,视电阻率长轴随之减小.从图 9b可以看出,当异常体各向异性主轴绕z轴旋转时,ρa发生旋转,其长轴方向指示各向异性异常体电阻率较大的主轴方向.

图 9 异常体主电导率张量绕坐标轴旋转时ρa极性图 (a)绕x轴旋转;(b)绕z轴旋转. Fig. 9 Polar plots of ρa for the conductivity tensor of 3D anisotropic anomalous body rotated around x- or z-axes. (a) Rotation around x-axis; (b) Rotation around z-axis.

图 10展示针对上节各向异性围岩模型计算的视电阻率响应ρa.从图 10a可以看出,当各向异性主轴绕x轴旋转时,主轴电阻率ρx保持不变(50 Ωm),由此ρa短轴位于x方向,且该方向上ρa保持不变;相比之下,主轴电阻率ρy则由100 Ωm变到50 Ωm,导致视电阻率ρa长轴随着旋转角度增大明显变小.进一步从图 10b可以看出,当各向异性主轴绕z轴旋转时,ρa发生旋转,但与各向异性异常体的情况不同,此时ρa长轴方向垂直于电阻率较大的主轴方向.这可以给出如下解释:当围岩为各向异性时,沿层理方向的电流被低阻异常体吸引产生通道效应,导致垂直层理方向电流向异常体集中,因此围岩中电流密度减小,体现高阻特征.

图 10 围岩主电导率张量绕坐标轴旋转时ρa极性图 (a)绕x轴旋转; (b)绕z轴旋转. Fig. 10 Polar plots of ρa for the conductivity tensor of the anisotropic host rock rotated around x- or z- axes (a) Rotation around x-axis; (b) Rotation around z-axis.
3 结论与展望

本文基于非结构四面体网格和矢量有限元方法成功实现任意各向异性介质中张量CSAMT三维高效数值模拟.和层状各向异性模型的半解析解对比验证了本文提出的算法有效性和精确性.各向异性异常体和各向异性围岩模型的数值模拟结果表明,各向异性对视电阻率幅值及其地表分布特征均产生很大影响.极坐标下视电阻率ρa极性图能很好地识别各向异性主轴方向.因此,在各向异性效应较为显著的地区开展CSAMT工作必须考虑各向异性的影响.通过采用张量数据采集并结合相应的基于极坐标的视电阻率转换技术可有效识别地下目标体的各向异性特征.本文的研究结果可为认识CSAMT各向异性效应提供依据.如何利用反演技术从张量CSAMT数据中求解出各向异性参数将是我们未来的重要研究方向.

致谢

作者向齐彦福博士和吉林大学电磁“千人计划”研究团队其他成员在文章准备过程中提供的帮助表示感谢

References
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. DOI:10.1016/j.parco.2005.07.004
Cao Q, Lin C H, Tan H D, et al. 2016. Occam inversion of CSAMT data in layered media with azimuthal anisotropy. Geophysical and Geochemical Exploration (in Chinese), 40(3): 594-602. DOI:10.11720/wtyht.2016.3.23
Chen G B. 2009. Integral equation method for 3-D electromagnetic modeling in layered anisotropic earth and its applications[Ph. D. thesis] (in Chinese). Changchun: Jilin University.
Deng J Z, Tan H D, Chen H, et al. 2011. CSAMT 3D modeling using staggered-grid finite difference method. Progress in Geophysics (in Chinese), 26(6): 2026-2032. DOI:10.3969/j.issn.1004-2903.2011.06.017
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 Journal of Geophysics (in Chinese), 47(4): 723-730.
Han B, Hu X Y, Huang Y F, et al. 2015. 3-D frequency-domain CSEM modeling using a parallel direct solver. Chinese Journal of Geophysics (in Chinese), 58(8): 2812-2826. DOI:10.6038/cjg20150816
Hu Y C, Li T L, Fan C S, et al. 2015. Three-dimensional tensor controlled-source electromagnetic modeling based on the vector finite-element method. Applied Geophysics, 12(1): 35-46. DOI:10.1007/s11770-014-0469-1
Jin J M. 2014. The Finite Element Method in Electromagnetics. 3rd ed. Hoboken, NJ: Wiley-IEEE Pres.
Lei D. 2010. Studies and applications of 2-D CSAMT modeling and inversion with a dipole source and topography. Chinese Journal of Geophysics (in Chinese), 53(4): 982-993. DOI:10.3969/j.issn.0001-5733.2010.04.023
Li X B, Pedersen L B. 1992. Controlled-source tensor magnetotelluric responses of a layered earth with azimuthal anisotropy. Geophysical Journal International, 111(1): 91-103. DOI:10.1111/gji.1992.111.issue-1
Li X B, Oskooi B, Pedersen L B. 2000. Inversion of controlled-source tensor magnetotelluric data for a layered earth with azimuthal anisotropy. Geophysics, 65(2): 452-464. DOI:10.1190/1.1444739
Liu Y H, Yin C C. 2014. 3D anisotropic modeling for airborne EM systems using finite-difference method. Journal of Applied Geophysics, 109: 186-194. DOI:10.1016/j.jappgeo.2014.07.003
Nédélec J C. 1980. Mixed finite elements in R3. Numerische Mathematik, 35(3): 315-341. DOI:10.1007/BF01396415
Qi Y F, Yin C C, Liu Y H, et al. 2017. 3D time-domain airborne EM full-wave forward modeling based on instantaneous current pulse. Chinese Journal of Geophysics (in Chinese), 60(1): 369-382. DOI:10.6038/cjg20170131
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, Kalscheuer T, Greenhalgh S, et al. 2014. A finite-element-based domain-decomposition approach for plane wave 3D electromagnetic modeling. Geophysics, 79(6): E255-E268. DOI:10.1190/GEO2013-0376.1
Wang Y L, Di Q Y, Wang R. 2017. Three-dimensional modeling of controlled-source audio-frequency magnetotellurics using the finite element method on an unstructured grid. Chinese Journal of Geophysics (in Chinese), 60(3): 1158-1167. DOI:10.6038/cjg20170326
Yin C, Maurer H M. 2001. Electromagnetic induction in a layered earth with arbitrary anisotropy. Geophysics, 66(5): 1405-1416. DOI:10.1190/1.1487086
Yin C C. 2000. Geoelectrical inversion for a one-dimensional anisotropic model and inherent non-uniqueness. Geophysical Journal International, 140(1): 11-23. DOI:10.1046/j.1365-246x.2000.00974.x
Yin C C, Fraser D C. 2004. The effect of the electrical anisotropy on the response of helicopter-borne frequency-domain electromagnetic systems. Geophysical Prospecting, 52(5): 399-416. DOI:10.1111/gpr.2004.52.issue-5
Yin C C. 2006. MMT forward modeling for a layered earth with arbitrary anisotropy. Geophysics, 71(3): G115-G128. DOI:10.1190/1.2197492
Zhou J M, Wang J X, Shang Q L, et al. 2014. Regularized inversion of controlled source audio-frequency magnetotelluric data in horizontally layered transversely isotropic media. Journal of Geophysics and Engineering, 11(2): 025003. DOI:10.1088/1742-2132/11/2/025003
曹强, 林昌洪, 谭捍东, 等. 2016. 可控源音频大地电磁法方位非各向同性层状介质Occam反演. 物探与化探, 40(3): 594-602. DOI:10.11720/wtyht.2016.3.23
陈桂波. 2009. 各向异性地层中电磁场三维数值模拟的积分方程算法及其应用[博士论文]. 长春: 吉林大学.
邓居智, 谭捍东, 陈辉, 等. 2011. CSAMT三维交错采样有限差分数值模拟. 地球物理学进展, 26(6): 2026-2032. DOI:10.3969/j.issn.1004-2903.2011.06.017
底青云, UnsworthM, 王妙月. 2004. 复杂介质有限元法2.5维可控源音频大地电磁法数值模拟. 地球物理学报, 47(4): 723-730.
韩波, 胡祥云, 黄一凡, 等. 2015. 基于并行化直接解法的频率域可控源电磁三维正演. 地球物理学报, 58(8): 2812-2826. DOI:10.6038/cjg20150816
雷达. 2010. 起伏地形下CSAMT二维正反演研究与应用. 地球物理学报, 53(4): 982-993. DOI:10.3969/j.issn.0001-5733.2010.04.023
齐彦福, 殷长春, 刘云鹤, 等. 2017. 基于瞬时电流脉冲的三维时间域航空电磁全波形正演模拟. 地球物理学报, 60(1): 369-382. DOI:10.6038/cjg20170131
汤井田, 何继善. 2005. 可控源音频大地电磁法及其应用. 长沙: 中南大学出版社.
王亚璐, 底青云, 王若. 2017. 三维CSAMT法非结构化网格有限元数值模拟. 地球物理学报, 60(3): 1158-1167. DOI:10.6038/cjg20170326