2020, Vol. 63
2. 防灾科技学院, 地球科学学院, 三河 065201;
3. 河北省地震动力学重点实验室, 三河 065201
2. Schoolof Earth Sciences, Instituteof Disaster Prevention, Sanhe 065201, China;
3. Hebei Key Laboratory of Earthquake Dynamics, Sanhe 065201, China
大地电磁测深(Magnetotelluric sounding, MT)是资源勘查、能源勘探及深部构造探测的重要方法之一.MT法以天然的平面电磁波为场源, 依据不同频率的电磁波在导体中具有不同趋肤深度的原理，在地表测量由高频至低频的地球电磁响应序列，经过相关的数据处理和分析来获得大地由浅至深的电性结构(Tikhonov，1950；Cagniard, 1953).随着大地电磁测深二维正、反演算法研究的成熟和水平地形条件下大地电磁测深三维正演算法研究的进步, 三维反演算法也日渐趋于实际应用，近年来主要在地热、地震构造带、石油天然气的探测(王楠等，2016；胥博文等，2018)、地壳尺度的断层、俯冲带与缝合带的研究(王刚等，2017；崔腾发等，2020)以及岩石圈尺度的大范围电性结构反演解释(黄颖等，2019；杨文采等，2020)等领域得到应用, 并取得明显应用效果.但是，要真正实现大地电磁测深三维反演的实用化, 则必须解决复杂地形条件下大地电磁场的三维正演模拟和反演问题.
地形效应在感应电磁法中一般指由近地表地形引起的电流畸变现象, 这一现象对于感应电磁法的反演与解释有着较大的影响(Jiracek，1990；Vallianatos，2002).自从这一现象被发现以来, 地形效应一直是感应电磁法领域关注与研究重点之一.解决地形效应问题的方法大致可以分为两大类.一类是采用畸变校正的方式, 首先在观测数据中将地形效应消除, 得到校正后的数据, 再使用不带地形的传统数值方法来进行正反演(Jiracek et al., 1989; 晋光文等，1998; 胡祖志等，2008；Nam and Kim, 2010；薛国强等，2016).另一类是采用带地形反演方式, 即将地形纳入数值模型, 直接使用观测数据采用带地形数值方法进行正反演.畸变校正方法虽然可以不同程度的消除地形影响, 但也可能造成反演模型中异常体几何形态的变形, 难以完全解决地形效应问题(Gürer and Ìlkışık，1997).因此, 随着近年来计算机数值模拟与反演技术的发展, 越来越多的学者开始研究直接带地形正反演方法.
自20世纪80年代以来，国内外学者在起伏地形条件下大地电磁场数值模拟方面做了大量研究.Wannamaker等(1986)与Chouteau和Bouchard(1988)使用有限单元法对二维地形进行了模拟，并且引发了对于电磁法地形响应的广泛讨论；徐世浙等(1992)、陈小斌(2000)、刘云和王绪本(2010)分别实现了基于边界元和有限元法的大地电磁二维正演模拟算法，并利用数值方法探讨了地形对大地电磁测深法的影响；张翔和胡文宝(1999)则开发了二维带地形校正的联合反演算法；霍光谱等(2015)采用有限差分法实现了带地形大地电磁各向异性二维模拟.
随着计算条件和三维数值模拟方法的进步，对三维地形的MT数值模拟研究也逐年兴起，徐世浙等(1997)、阮百尧等(2007)采用边界元法实现了对三维地形的数值模拟；Chen等(1998)使用有限差分法对三维山峰模型进行了数值模拟；Baba和Seama(2002)将FS(FlatSeafloor)技术结合Mackie等(1994)的大地电磁三维有限差分模拟技术实现了海底起伏地形下的三维数值模拟；Sasaki(2003)采用交错网格有限差分法模拟了三维山丘地形的大地电磁响应；Shi等(2004)研究了结合散度校正的大地电磁矢量有限单元法三维正演；Nam等(2007, 2009)基于带散度校正的大地电磁矢量有限元法对三维地形和岛屿-海洋地形进行了模拟；董浩等(2014)研究了带散度校正的有限差分法带地形三维大地电磁正演，并在此基础上采用非线性共轭梯度法对三例不同复杂度的理论模型合成数据进行了带地形三维反演；顾观文等(2014)实现了大地电磁场三维地形影响的矢量有限元数值模拟.Ren等(2013)、殷长春等(2017)基于非结构有限元法模拟了三维地形条件下地电模型的大地电磁响应；Kordy等(2016a, b)研究了基于变形六面体的矢量有限元法带地形三维正演，在正演过程中引入散度校正技术以保证正演计算精度，进而在此基础上实现了基于数据空间的带地形三维反演；Yin等(2018)基于非结构有限元法三维正演，采用有限内存拟牛顿法实现了带地形MT三维反演，分别对水平地形和起伏地形下的三例理论模型合成数据进行反演，验证算法的有效性.
为提高带地形大地电磁测深资料的反演速度与反演结果的可靠性，本文首先在作者前期三维正演工作(顾观文等，2014)的基础上，将并行直接稀疏求解器PARDISO引入三维正演问题的求解，在正演过程中无需散度校正，实现快速三维正演计算.然后将快速三维正演方案应用于共轭梯度法三维反演中“拟正演”问题的求解，从而实现迭代反演中主要计算量：即雅可比偏导数矩阵与一个矢量的乘积、以及雅可比偏导数矩阵的转置与一个矢量的乘积的快速计算.最后，设计不同类型的五个理论模型进行正反演试算，验证本文开发的三维正反演算法的正确性和可靠性，进而利用本文实现的三维反演程序反演了某矿区大地电磁实测数据，进一步验证了本文反演算法的有效性.
1 带地形大地电磁三维矢量有限元数值模拟 1.1 大地电磁三维正演的边值问题本研究中使用的正演方法为基于六面体网格剖分的矢量有限元法(Shi et al., 2004; Nam et al., 2007)，在大地电磁研究的频率范围内，忽略位移电流的作用.取时谐场为e-iωt，麦克斯韦方程组的微分形式表示为
|
(1a) |
|
(1b) |
|
(1c) |
|
(1d) |
式中：E为电场强度矢量，H为磁场强度矢量，σ是介质的电导率，μ是介质的磁导率，ε是介质的介电常数，q为自由电荷密度.对(1a)式两边取旋度，再将(1b)式代入到(1a)式，则将电磁场满足的一阶微分方程变为二阶微分方程，可得电场E满足方程：
|
(2) |
求解偏微分方程组(2)得到电场分量Ex、Ey、Ez后，根据(1a)式求得磁场分量Hx、Hy、Hz.
由于空气中的电场和磁场受地形和不均匀体影响而不均匀，因此在数值模拟中必须包括空气层和地下介质，如图 1所示，模拟区域设为Ω，由空气层和地下介质两部分组成，Ω=Ω1∪Ω2.在模拟区域中，空气层的电导率一般在10-6~10-10 S·m-1之间，这时空气和地下介质的接触面变成内部边界.
|
图 1 带地形三维MT数值模拟区域剖面示意图 Fig. 1 Schematic section of numerical modeling domain for 3D MT with topography |
将(2)式写为
|
(3) |
取第一类边界条件：
|
(4) |
式(4)中g是边界上的矢量电场，可以用一维或者二维MT计算值(Mackie et al., 1993; Siripunvaraporn et al., 2002).这样，(3)式和(4)式构成了大地电磁三维正演的边值问题.
1.2 模型与地形剖分本项研究对正演模型区域采用矩形六面体剖分，矩形网格x、y和z方向的间距可以是不均匀的.地形剖分也采用矩形六面体，比如三维山峰地形的剖分如图 2所示.若在地形起伏界面(空气-地面界面)上加大网格剖分密度, 使其剖分足够精细时, 大地电磁场数值模拟响应将与理论场值相近, 甚至趋于一致(Chen et al., 1998).
|
图 2 三维地形及其网格剖分 (a)三维地形示意图；(b)地形网格剖分和MT测点分布示意图. Fig. 2 3D topography and grid (a) Sketch of 3D topography; (b) Topography meshing and distribution of MT measurement sites. |
有限元法求解上述区域(图 1)的电磁场问题，需要对研究区域离散化，即对研究区域进行六面体网格剖分(图 3a)，沿x、y和z轴方向分别剖分成Nx、Ny和Nz段，网格间距分别为Δx(i)(i=1, …, Nx)、Δy(j)(j=1, …, Ny)和Δz(k)(k=1, …, Nz).可以推导，每个六面体单元的内部电场分量用六面体的12条棱边的场值分量(图 3b)通过插值求取.公式为
|
(5) |
|
图 3 矢量有限元法的区域剖分示意图 (a)区域剖分示意图；(b)电场分量位置图. Fig. 3 Domain subdivision of the vector finite element method (a) Domain subdivision; (b) Location of electric field components. |
(5) 式写成矢量形式为
|
(6) |
其中：Nie=(Nxie, Nyie, Nzie)为矢量有限元中的形函数，Eie为单元棱边上的未知电场值，上标e表示第e个网格单元.即:
|
(7) |
|
(8) |
|
(9) |
式中，坐标转换函数为ξ=(x-xc)/a, η=(y-yc)/b, ζ=(z-zc)/c，(xc, yc, zc)是六面体的中心坐标，2a、2b、2c分别是六面体x、y、z方向的长度.(ξi, ηi, ζi)的取值与棱边编号有关.
由(3)式定义矢量余函数为
|
(10) |
把矢量基函数作为权函数，采用迦辽金方法使整个域内的积分矢量余函数为最小，即：
|
(11) |
式中：NE为计算域网格单元总数.
将(10)式代入到(11)式，对于第e个单元，得到：
|
(12) |
其中：n是外法线方向，Se是单元面积分，Ve是单元体积分.由于
|
(13) |
其中：Se表示场源项，Ee表示棱边上的电场，Ke为单元刚度矩阵，是一个12×12阶的复数矩阵，可按式(14)解析计算得出(金建铭，1998)：
|
(14) |
将每个单元电场满足的线性方程进行组合，可以得到整个计算域上电场满足的线性方程组：
|
(15) |
其中：K是系统刚度矩阵，E是整个计算域的网格单元棱边上的电场值向量，s是源向量，由计算域的上、下、左、右的边界场值与边界上的单元刚度矩阵计算得到.
根据大地电磁场线性方程组特有的稀疏性，开发有针对性的数据结构接口，将PARDISO快速求解器(Schenk and Gärtner, 2004；潘克家等，2016)应用到式(15)中的大型稀疏线性方程组的求解，得到计算域上网格单元棱边上的电场值，然后根据麦克斯韦方程组(1a式)微分求取磁场.
1.4 视电阻率及阻抗相位的计算根据(Newman and Alumbaugh, 2000；谭捍东等，2004)的论文，假设两种线性无关的场源激发的表面电场和磁场分别为Ex1，Ey1，Hx1，Hy1，Ex2，Ey2，Hx2，Hy2, 由此可以计算三维MT的张量阻抗：
|
(16) |
张量阻抗的每一个分量可表示为
|
(17a) |
|
(17b) |
|
(17c) |
|
(17d) |
其中：下标1、2表示极化模式，下标x、y表示x、y分量，E表示电场，H表示磁场，Z表示阻抗.式(17b)和(17c)定义的响应分别称为XY和YX模式响应，按照式(18)可求出三维介质的视电阻率和相位：
|
(18) |
其中i=X, Y, j=X, Y.
2 带地形大地电磁三维反演 2.1 目标函数根据Tikhonov正则化理论，定义三维带地形的MT反演问题的目标函数(Mackie and Madden, 1993；Newman and Alumbaugh, 2000)为
|
(19) |
目标函数的第一部分ϕd表征数据的拟合情况，第二部分ϕm表征模型的光滑度.λ是平衡数据拟合和模型光滑度的正则化因子，目前对于反演过程中λ的选取主要采用两种方法，一种是在反演过程中先采用比较大的正则化系数，随着迭代的进行逐步减小正则化系数，当λ值或拟合差小于阀值时终止迭代的自适应方法(刘云鹤和殷长春，2013)，另一种是通过实验的方法来确定最佳的λ(吴小平和徐果明，2000；吴小平和汪彤彤，2001)，本文反演过程中λ的选取采用后者的方法.dobs是观测数据向量，采用张量阻抗的四个元素的值或者归一化阻抗值，dfwd是当前模型正演数据向量，Wd为观测数据的方差，m是模型向量，取m=lnσ(σ为网格单元的电导率)，mref为先验模型向量，Wm为模型光滑度矩阵，一般可取为二阶微分算子例如拉普拉斯算子，上标H表示共轭转置，上标T表示转置.
采用高斯—牛顿法求取目标函数的极小值，第k次迭代模型的修改量等价于求方程(20)的解(Mackie and Madden, 1993)：
|
(20) |
式中：J为雅可比偏导数矩阵，
当式(20)左边部分(JHWd-1J+λWmTWm)·δmk括号中的矩阵为奇异矩阵时，采用Marquardt(1963)、陈小斌等(2005)、Shi等(2009)、Siripunvaraporn(2012)的方法增加一个阻尼因子，即将方程(20)式改为
|
(21) |
式中：ε为阻尼因子，I为单位矩阵，Δd=dobs-dfwd表示观测数据与预测数据的残差，Δm=mref-mk表示预测模型与先验模型的残差.
2.2 反演流程共轭梯度法(CG)的基本思想是把共轭性和梯度法结合，利用已知点处的梯度方向构造一组共轭方向，并沿这组方向进行搜索，求出目标函数的极小点.其实质是采用该方法求解高斯—牛顿反演方程，不需要直接计算雅可比偏导数矩阵的值，只需计算雅可比偏导数矩阵和一个矢量的乘积以及雅可比矩阵的转置和一个矢量的乘积.根据互换原理(Madden, 1972；Parasnis, 1988; Rodi, 1976；Mackie and Madden, 1993)，雅可比偏导数矩阵与一个矢量的乘积、以及雅可比偏导数矩阵的转置与一个矢量的乘积可以转换为两个不同虚场源的正演问题(Newman and Alumbaugh, 2000；Rodi and Mackie, 2001)，即“拟正演”问题.从而克服带地形三维MT反演问题中雅可比矩阵存储量大和计算时间长的两大困难.反演算法的流程(图 4)描述如下：
|
图 4 带地形三维MT共轭梯度反演流程图 Fig. 4 Flowchart of 3D MT inversion with topography using conjugate gradients |
(1) 输入初始模型、反演数据和模型剖分参数等，设置迭代次数k=1和最大迭代次数niters；
(2) 正演计算当前模型，并计算数据残差向量Δdk=dobs-dkfwd；
(3) 计算方差‖Δdk‖，如果方差小于或等于阀值，反演结束，否则，转下一步；
(4) 简记式(21)为A·x=b，则b=JHWd-1(Δdk)+λWmTWm(mref-mk)；
(5) 通过“拟正演”求取右端项b；
(6) 初始化共轭梯度的条件：x0=0, r0=b；
(7) 令共轭梯度的迭代次数i=1，最大迭代次数为ncgs；
(8) 计算共轭梯度修改量：
|
(9) 若达到最大迭代次数ncgs，结束CG迭代，求得x，即δmk，进入下一步反演；否则，回步骤(8)，继续共轭迭代；
(10) 计算模型修改量：mk+1=mk+δmk；
(11) 若达到反演最大迭代次数niters，结束反演；否则，回到步骤(2)，继续迭代反演.
从上述反演算法流程可以看出，在每一次反演迭代过程中，求数据的残差Δdk需要进行一次正演.由于雅可比矩阵的共轭转置与一个矢量的乘积为JHy=(JTy*)*，实际上就是求雅可比矩阵的转置与一个共轭向量的乘积，然后求其共轭，所以求反演方程的右端项b即为计算雅可比矩阵的转置和一个矢量的乘积，相当于一次“拟正演”，共轭梯度法计算模型修改量的时候，计算Api=(JHWd-1J+λWmTWm+εI)pi，这个过程相当于两次“拟正演”.因此，上述反演方法迭代一次所需要的的正演次数为(2+2×(ncgs-1))×nf，ncgs为共轭梯度的最大迭代次数，nf为频点数.由于三维MT的正演必须计算两种独立极化场源，上述算法实际的正演计算量需要乘以2.可以看出，共轭梯度法三维反演的主要计算量在于正演问题的求解，为了实现快速三维正演，将并行直接稀疏求解器PARDISO引入三维正演问题的求解，以提高三维反演计算效率.
2.3 雅可比矩阵的计算“拟正演”问题从共轭梯度反演算法可知，在求取反演方程(21)式右端项b和计算模型修改量过程中的Api均涉及到雅可比矩阵J，从b和Api表达式的形式上看，只需计算出J后直接代入两者相应的表达式即可，附录A给出了J的具体计算过程(Rodi, 1976；Mackie and Madden, 1993；Newman and Alumbaugh, 2000；Shi et al., 2009)，可以看出，计算J需要进行M次正演(M为模型个数)，其计算量相当巨大.为了减少共轭梯度法反演的计算量，有必要避免直接求取雅可比矩阵J，而只需要将J与向量的乘积和JT与向量的乘积转换为两次“拟正演”，下面介绍如何通过解“拟正演”问题来直接计算出J或JT与向量的乘积.
将雅可比矩阵与向量的乘积简记为J·x，采用向量(J的阶数为N×M，N为观测数据的个数，M为模型个数)表示，则有：
|
(22) |
将(A3)式代入(22)式，则J·x的任意行第i行可写为
|
(23) |
将式(A7a)和(A7b)代入(23)式，可得Ex、Ey极化模式的J·x计算公式为
|
(24a) |
|
(24b) |
令：
|
(25a) |
|
(25b) |
式中：
令：
|
(26) |
(26) 式两端同乘以系数矩阵K，则变为
|
(27) |
由(27)式可以看出，求u1和u2的过程就相当于把
基于上述三维正反演方法，我们采用C++语言开发了带地形大地电磁三维正反演程序.为了检验三维正演计算程序的正确性和三维反演程序的有效性，设计5个不同类型的地电模型进行正反演试算，并分析试算结果.为了进一步验证反演算法的有效性，反演了某矿区大地电磁实测数据，并分析三维反演结果.
本文所述的正反演算例均在曙光W560-G20工作站上完成，计算及程序编译环境如下：CPU为Intel E5-2643 V4(3.4 G)，内存为64GB, 操作系统为Windows 7(64位)；编译环境为Microsoft Visual Studio 2012(已集成Intel Parallel Studio XE 2013).
3.1 三维正演计算程序验证为了验证正演程序对MT地形影响数值模拟效果，采用(Wannamaker et al., 1986)的二维山峰模型.模型背景电阻率为100 Ωm，山峰地形如图 5所示.设二维山峰地形的走向为x方向，山峰倾向为y方向，z方向垂直向下.为了尽量避免对y方向的三维影响，将山峰地形沿x方向延伸3.43 km，整个模型区域(34300 m×34300 m×91993 m)沿x、y和z方向剖分为43×43×29(其中z方向山峰顶面以上7层为空气层)个网格单元.山峰地形采用9个纵向网格单元的划分，网格间距为50 m.
|
图 5 二维山峰地形示意图 Fig. 5 Sketch of 2D ridge |
对上述模型(图 5)进行三维正演模拟，计算频率f=2 Hz时TE极化模式和TM极化模式的视电阻率曲线如图 6所示.图 6中黑色虚线和实线分别是二维有限元计算的TE和TM极化模式的视电阻率曲线图；圆圈和方框是本次研究算法(PARDISO)计算的TE和TM极化模式的视电阻率曲线图.“+”号和“×”号是前期研究算法(BICG)计算的TE和TM极化模式的视电阻率曲线图.从图中可以看出，本次研究算法(PARDISO)计算的结果与二维有限元的计算结果一致，从而说明了本次研究算法计算结果准确可靠, 同时可以看出，在平地与山峰拐点处，PARDISO计算的TM结果比BICG更接近于二维有限元计算结果.计算时间方面，本文采用的PARDISO求解器正演一次上述模型仅耗时24 s，而BICG迭代算法耗时407 s，PARDISO与BICG法的计算速度比达17倍，有利于快速反演.
|
图 6 三维矢量有限元算法(PARDISO、BICG)计算的二维地形影响与二维有限元结果对比图 Fig. 6 Comparision between modeling results of 3DVFEM(PARDISO, BICG) and 2DFEM for 2D ridge |
设计的地电模型为山谷地形下的单个低阻棱柱体模型(图 7)，其地形的谷深相对于地面为2.4 km，山谷模型的宽度为6 km，地形倾角约为18°.低阻棱柱体的规模为10 km×10 km×2 km，顶界面距山谷底部1.9 km，底界面距山谷底部3.9 km，其电阻率为10 Ωm，围岩电阻率为100 Ωm.
|
图 7 三维山谷地形下的三维低阻体模型 (a)三维山谷地形示意图；(b) X-Y平面图；(c) X-Z方向剖面图；(d) Y-Z方向剖面图. Fig. 7 3D low-resistivity model with valley topography (a) Sketch of valley; (b) X-Y plane diagram; (C) X-Z profile; (d) Y-Z profile. |
网格剖分采用36×36×32(Z方向包括5层空气层，本文后面3个模型算例的空气层数设置与此相同，不再单独说明)，测点设为网格单元中间，观测数据采用算区中间的28×28个测点，在每个测点上大地电磁测深观测信号的频率范围为100~0.1 Hz，共取7个频点，反演输入的数据为XY和YX模式的阻抗，并施加5%的随机噪声.反演参数设置为：正则化因子λ=10-3，阻尼因子ε=10-5，共轭梯度内部迭代次数设为6次.
反演迭代20次耗时约88h(3.7天)，图 7是反演的拟合误差RMS曲线下降图，拟合误差计算表达式为：
|
图 8 反演拟合误差下降曲线图 Fig. 8 RMS variation with iteration number after 20 iterations |
|
图 9 山谷地形下低阻体模型三维反演结果 Fig. 9 3D inversion results of MT data on valley model |
设计的地电模型为山峰地形下的单个高阻棱柱体模型(图 10)，其地形的峰高相对于地面为2.4 km，山峰模型的宽度为4 km，地形倾角约为18°.高阻棱柱体的规模为10 km×10 km×1.8 km，顶界面距水平地表 0.9 km，底界面距水平地表 2.7 km，其电阻率为100 Ωm，围岩电阻率为10 Ωm.
|
图 10 三维山峰地形下的三维高阻体模型 (a)三维山峰地形示意图；(b) X-Y平面图；(c) X-Z方向剖面图；(d) Y-Z方向剖面图. Fig. 10 3D high-resistivity model with ridge (a) Sketch of 3D ridge; (b) X-Y plane; (c) X-Z profile; (d) Y-Z profile. |
网格剖分采用36×36×32，观测数据采用算区中间的28×28个测点，在每个测点上大地电磁测深观测信号的频率范围为100~0.1 Hz，共取7个频点，反演输入的数据为XY和YX模式的阻抗，并施加5%的随机噪声.
迭代反演10次耗时约34 h(1.4天)，图 11为迭代10次的反演电阻率模型图，可以看出，反演结果能较好地反映出高阻异常体模型的空间位置和形态，与理论模型电性结构(图 10)吻合.
|
图 11 山峰地形下高阻体模型三维反演结果 Fig. 11 3D inversion results of MT data on high-resistivity model with ridge |
设计的地电模型为峰谷组合地形下含两个低阻棱柱体模型(如图 12所示)，其地形的峰高和谷深相对于地面均为1.5 km，山峰和山谷模型的宽度分别为4 km和6 km，在山峰和山谷地形的下方各设计了两个规模均为10 km×10 km×1.5 km的低阻棱柱体，其电阻率为10 Ωm，围岩电阻率为100 Ωm.
|
图 12 峰谷组合地形下含低阻体模型 (a)三维峰谷组合地形示意图；(b)三维峰谷组合地形下含低阻体Y-Z方向剖面图. Fig. 12 Schematic of 3D low-resistivity model with ridge-valley combined terrain (a) Sketch of 3D ridge-valley topography; (b) Y-Z profile. |
网格剖分采用50×50×32，观测数据采用算区中间的42×42个测点，在每个测点上大地电磁测深观测信号的频率范围为100~0.01 Hz，共取10个频点，反演输入的数据为XY和YX模式的阻抗，并施加了5%的随机噪声.
迭代反演10次耗时约172 h(约7天)，图 13为迭代10次的反演电阻率模型图，不难看出，两个低阻异常体在横向和纵向都得到了很好的归位，反演的地电断面与理论模型(图 12b)吻合.
|
图 13 峰谷组合地形下低阻体模型三维反演结果 Fig. 13 3D inversion results of MT data on low-resistivity model with ridge-valley combined terrain |
设计一例起伏地形下较为复杂地电模型，并对该模型正演合成数据施加10%的高斯噪声，以实验本文反演方法对观测数据噪声的依赖程度，为该反演方法后续应用于野外实测数据的反演解释提供可行性依据.设计的地电模型为峰谷组合地形下含低阻体和高阻体、低阻和高阻处于接触状态(如图 14)，其地形起伏状况与以上峰谷组合下低阻体模型完全相同.该模型的围岩导电性为两层结构，上层电阻率为500 Ωm，下层产状为两边倾斜的高阻基底，其电阻率为2000 Ωm.在围岩上下接触层之间存在一个梯形台柱体，其电阻率为20 Ωm，该台柱体与下层高阻基岩顶界面接触.
|
图 14 三维峰谷组合地形下高低阻体接触模型Y-Z方向剖面图 Fig. 14 Y-Z profile of the combined conductive-resistive model with ridge-valley topography |
网格剖分采用50×50×38，观测数据采用算区中间的42×42个测点，在每个测点上大地电磁测深观测信号的频率范围为100~0.01 Hz，共取10个频点，反演输入的数据为XY和YX模式的阻抗，并对阻抗数据施加10%的高斯噪声.
迭代反演10次耗时约282 h(约12天)，图 15为含10%噪声的阻抗数据反演电阻率模型图，由图 15可以看出，围岩的二层结构特征、两边产状倾斜的高阻基底以及低阻异常体与高阻基岩的接触状态在反演断面中均得到了很好的重现；同时，反演模型中低阻异常体的规模和位置、高阻基底的深度与理论模型基本一致.反演的电性结构与理论模型(图 14)吻合，表明本文的电阻率三维反演结果可靠，反演方法对观测数据的噪声依赖性不大，具备应用于实测数据反演的可行性.
|
图 15 引入10%观测误差时三维峰谷组合地形下高低阻体接触模型三维反演结果 Fig. 15 3D Inversion results of MT data bearing 10% measurement error on combined conductive-resistive model with ridge-valley topography |
(1) 关于某矿区资料
该矿区分为一、二、三矿带，本次研究区域位于矿区北部的一矿带，据测区地质图(图 16)可以看出：测区北侧为震旦系洗肠井群第四岩组千枚岩，南侧为第三岩组大理岩、千枚岩.两者的接触部位(近东西错动断裂)为矿区的一矿带赋存位置.二矿带位于测区的东南角，位于三岩组第二岩性段千枚岩与第三岩性段大理岩夹千枚岩过渡部位.测区南端为三矿带，测点未完全覆盖该矿带.测区内沿近南北走向的F3构造与东西向错动断裂带相交叉，为后期的断裂带，也是岩体入侵通道，与矿体的走向与产状有重要的关系.
|
图 16 测区地质图与大地电磁测线布置 1震旦系洗肠井群四岩组千枚岩; 2震旦系洗肠井群三岩组三岩段千枚岩; 3印支期花岗岩; 4震旦系洗肠井群三岩组三岩段大理岩; 5震旦系洗肠井群三岩组二岩段大理岩; 6震旦系洗肠井群三岩组二岩段千枚岩; 7花岗斑岩脉; 8矿带; 9断裂; 10测线; 11板岩; 12钻孔. Fig. 16 Geological map and layout of MT survey lines in study area 1 Phyllite of Siyan formation from Sinian Xichangjing Group; 2 The third rock section phyllite of Sanyan formation from Sinian Xichangjing Group; 3 Indosinian granite; 4 The third rock section marble of Sanyan formation from Sinian Xichangjing Group; 5 The second rock section marble of Sanyan formation from Sinian Xichangjing Group; 6 The second rock section phyllite of Sanyan formation from Sinian Xichangjing Group; 7 Granite porphyry dike; 8 Metallogenic belt; 9 Fault; 10 Survey line; 11 Slate; 12 Borehole. |
相对地层与岩体而言，磁黄铁矿、黄铁矿与铅锌矿均具有低电阻率的特征，这可以作为本区找矿的地球物理标志.千枚岩本身与大理岩的电阻率差别不大，但是含炭千枚岩具有低电阻率的特征，连同测区北部的炭质板岩容易成为找矿的干扰因素.
为了进一步查明本区矿体赋存的有利部位，在研究区开展大地电磁探测工作，垂直于一矿带接触带构造方向布置9条大地电磁测线，测线方位为北偏东16°，测线从西至东依次编号为L240~L400(图 16), 测线间距为200 m.在每条测线上分布29个测点，点距为50 m.测点呈规则网格状分布，观测区共261个测点.
(2) 对反演所得三维电阻率模型初步分析
为了获得本研究区的三维精细电性结构，利用本文的三维反演程序对该区观测的大地电磁数据进行了三维反演.网格剖分采用41×41×31，观测数据采用观测区的9×29个测点，每个测点上大地电磁观测信号的频率范围为10000~5.86 Hz，共取33个频点，反演输入的数据为XY和YX模式的阻抗.根据观测数据的视电阻平均值，将初始模型电阻率设置为500 Ωm，迭代反演10次耗时约183 h(约7.6天).图 17为该测区的三维反演电阻率模型，从图中可以看出，测区的整体电性结构为北部呈现低阻特征、南部表现为高阻特征.结合测区地质图(图 16)可知，北部的低阻带为含碳质的千枚岩，南部的高阻带为大理岩，三维电性结构展示了大理岩与千枚岩的接触关系，并清晰地刻画出接触带由南西向至北东向延伸的展布特征.从分立于测区西北角和东北角两个低阻体的规模和形态可以看出，近南北走向的断裂F3对地层的破坏作用相当明显，成为岩体入侵的主要通道，与矿体的产状和走向有重要关系.图 17中异常①、②和③是电阻率等于30 Ωm的低阻体，异常①和②都处于大理岩与千枚岩两种岩体的接触带附近，钻孔ZK02(见图 16)位于异常②区内，该孔主要矿化蚀变类型为黄铁矿化、弱褐铁矿化，可认为低阻体②为矿致异常.根据一矿带的赋存条件可知，矿体位于大理岩与千枚岩的接触部位，而异常体①正处于岩体接触带附近，从而可推断低阻异常①也为矿致异常.异常③没有钻孔资料验证，由于临近三矿带位置，推测为隐伏硫化物矿床.
|
图 17 测区大地电磁观测数据三维反演结果 Fig. 17 3D inversion results of MT data in survey area |
依据本次反演获得的三维精细电性结构，揭示了研究区大理岩与千枚岩两大岩体的接触关系，并清晰刻画接触带的展布特征，查明已知断裂的走向及其与岩体的接触关系，推断矿致异常赋存的有利部位，在一定程度上有助于深化对研究区成矿作用的认识，为容矿有利层位的寻找、矿区潜在矿产资源的评价提供依据.
4 结论本项研究基于矢量有限元正演采用共轭梯度法实现了带地形大地电磁三维反演，并利用C++语言开发了三维正反演计算程序代码，设计不同类型的5个带地形三维理论模型开展了三维正反演试算，数值实验和算例分析表明所研发的带地形三维正、反演计算程序正确、可靠.同时利用本文实现的三维反演程序反演了某矿区大地电磁实测数据，进一步验证了本文反演算法的有效性.
为提高带地形大地电磁三维正演速度，采用基于并行直接稀疏求解器PARDISO且无需散度校正的正演方案，实现快速三维正演，对典型地形模型，在中等规模计算条件下，与前期工作的双共轭梯度法(BICG)正演计算结果比较，PARDISO与BICG法的计算速度比达17倍.
在共轭梯度法反演中，为了避免直接计算雅可比偏导数矩阵，基于互换原理将雅可比偏导数矩阵与一个矢量的乘积、以及雅可比偏导数矩阵的转置与一个矢量的乘积转换为两个不同虚场源的“拟正演”问题.在此基础上，进而将本文研究的快速三维正演方案应用于共轭梯度法三维反演中雅可比矩阵计算的“拟正演”问题求解，从而克服了带地形三维MT反演问题中雅可比矩阵存储量大和计算时间长的两大困难.
将本文开发的反演算法成功应用于某矿区大地电磁实测数据的三维反演，获得的三维精细电性结构清晰地反映了研究区的地电特征，为该区容矿有利层位的寻找、矿区潜在矿产资源的评价提供依据，表明该算法具备一定的实用性.
本文提出的三维快速正演方案仅针对于快速求解线性方程组这一层面，如在下一步的反演算法研究中，将MT的观测频点作为并行粒度，并将本文的快速正演方法融入其中开展并行算法研究，有望进一步提高三维反演效率.
附录A 基于矢量有限元的雅可比偏导数计算
对于∂Z/∂mk，具体到Z的某个分量(Zxx，Zxy，Zyx，Zyy)，根据式(17a—d)，可以得到其相对mk的导数，以Zxy为例，令：
|
(A1) |
设地表观测点的电场值为Ex1、Ex2、Ey1和Ey2，磁场值为Hx1、Hx2、Hy1、Hy2、Hz1和Hz2，据文献(Rodi，1976；吴小平，1999)，地面电磁场可表示为
|
(A2) |
式中：ax、ay、bx、by和bz为与测点位置有关的行向量，下标1和2分别表示Ex和Ey极化模式，E1、E2为地球内部主场矢量.
将式(A2)代入到(A1)式，对相同极化模式的同类项合并，可简化为
|
(A3) |
其中：
|
(A4a) |
|
(A4b) |
由式(15)可知，Ex和Ey极化模式的电场满足：
|
(A5a) |
|
(A5b) |
将(A5a)和(A5b)式两边对网格单元参数m求导：
|
(A6a) |
|
(A6b) |
将(A6a)和(A6b)式两边同乘以K-1，可得：
|
(A7a) |
|
(A7b) |
式中：K为矢量有限元的总体刚度矩阵，上标-1表示求逆运算，
令：
|
(A8a) |
|
(A8b) |
|
(A8c) |
|
(A8d) |
则(A7a)和(A7b)可变形为
|
(A9a) |
|
(A9b) |
通过解方程(A9a)和(A9b)可求取
Baba K, Seama N. 2002. A new Technique for the Incorporation of Seafloor Topography in Electromagnetic Modelling. Geophysical Journal International, 150(2): 392-402. DOI:10.1046/j.1365-246X.2002.01673.x |
Cagniard L. 1953. Basic theory of the magneto-telluric method of geophysical prospecting. Geophysics, 18(3): 605-635. DOI:10.1190/1.1437915 |
Cao X Y, Yin C C, Zhang B, et al. 2018. 3D magnetotelluric inversions with unstructured finite-element and limited-memory quasi-Newton methods. Applied Geophysics, 15(3-4): 556-565. DOI:10.1007/s11770-018-0703-8 |
Chen P F, Hou Z H, Fan G H. 1998. Three-dimensional topographic responses in MT USING finite difference method. Acta Seismologica Sinica, 11(5): 631-635. DOI:10.1007/s11589-998-0079-6 |
Chen X B. 2000. On the research of the influence of terrain to MT 2D forward compulation. Geophysical Prospecting for Petroleum (in Chinese), 39(3): 112-120. |
Chen X B, Zhao G Z, Tang J, et al. 2005. An adaptive regularized inversion algorithm for magnetotelluric data. Chinese Journal of Geophysics (in Chinese), 48(4): 937-946. |
Chouteau M, Bouchard K. 1988. Two-dimensional terrain correction in magnetotelluric surveys. Geophysics, 53(6): 854-862. DOI:10.1190/1.1442520 |
Cui T F, Chen X B, Zhan Y, et al. 2020. Characteristics of deep electrical structure and seismogenic structure beneath Anhui Huoshan earthquake area. Chinese Journal of Geophysics (in Chinese), 63(1): 256-269. DOI:10.6038/cjg2019M0458 |
Dong H, Wei W B, Ye G F, et al. 2014. Study of Three-dimensional magnetotelluric inversion including surface topography based on Finite-difference method. Chinese Journal of Geophysics (in Chinese), 57(3): 939-952. DOI:10.6038/cjg20140323 |
Gu G W, Wu W L, Li T L. 2014. Modeling for the effect of Magnetotelluric 3D topography based on the vector finite-element method. Journal of Jilin University (Earth Science Edition) (in Chinese), 44(5): 1678-1686. |
Gürer A, Ìlkışık O M. 1997. The importance of topographic corrections on magnetotelluric response data from rugged regions of Anatolia. Geophysical Prospecting, 45(1): 111-125. DOI:10.1046/j.1365-2478.1997.3160236.x |
Hu Z Z, He Z X, Sun W B, et al. 2008. An improved MT topographic correction method and discussion on relative issues. Oil Geophysical Prospecting (in Chinese), 43(3): 343-348. |
Huang Y, Wang X Q, Zhou Z K, et al. 2019. Application of clustering analysis algorithm in three-dimensional magnetotelluric interpretation. Progress in Geophysics (in Chinese), 34(2): 568-572. DOI:10.6038/pg2019CC0332 |
Huo G P, Hu X Y, Huang Y F, et al. 2015. MT modeling for two-dimensional anisotropic conductivity structure with topography and examples of comparative analyses. Chinese Journal of Geophysics (in Chinese), 58(12): 4696-4708. DOI:10.6038/cjg20151230 |
Jin G W, Zhao G Z, Xu C F, et al. 1998. The affection and correction on magnetotelluric response data for Inclination two Dimension terrain. Seismology and Geology (in Chinese), 20(4): 454-458. |
Jin J M.1998. The Finite Element Method in Electromagnecic Fields (in Chinese). Wang J G trans. Xi'an: Xidian University Press, 176-189.
|
Jiracek G R. 1990. Near surface and topographic distortions in electromagnetic induction. Surveys in Geophysics, 11(2-3): 163-203. DOI:10.1007/BF01901659 |
Jiracek G R, Reddig R P, Kojima R K. 1989. Application of the rayleigh-FFT technique to magnetotelluric modeling and correction. Physics of the Earth and Planetary Interiors, 53(3-4): 365-375. DOI:10.1016/0031-9201(89)90023-X |
Kordy M, Wannamaker K, Maris V, et al. 2016a. 3-D magnetotelluric inversion including topography using deformed hexahedral edge finite elements and direct solvers parallelized on SMP computers-Part Ⅰ:forward problem and parameter Jacobians. Geophysical Journal International, 204(1): 74-93. DOI:10.1093/gji/ggv410 |
Kordy M, Wannamaker K, Maris V, et al. 2016b. 3-dimensional magnetotelluric inversion including topography using deformed hexahedral edge finite elements and direct solvers parallelized on symmetric multiprocessor computers-Part Ⅱ:direct data-space inverse solution. Geophysical Journal International, 204(1): 94-110. DOI:10.1093/gji/ggv411 |
Liu Y, Wang X B. 2010. FEM using adaptive topography in 2D MT forward modeling. Seismology and Geology (in Chinese), 32(3): 382-391. |
Liu Y H, Yin CC. 2013. 3D inversion for frequency-domain HEM data. Chinese Journal of Geophysics (in Chinese), 56(12): 4278-4287. DOI:10.6038/cjg20131230 |
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. |
Mackie R L, Madden T R. 1993. Three-dimensional magnetotelluric inversion using conjugate gradients. Geophysical Journal International, 115(1): 215-229. DOI:10.1111/j.1365-246X.1993.tb05600.x |
Mackie R L, Smith J T, Madden T R. 1994. Three-dimensional electromagnetic modeling using finite difference equations:the magnetotelluric example. Radio Science, 29(4): 923-935. DOI:10.1029/94RS00326 |
Madden T R. 1972. Transmission systems and network analogies to geophysical forward and inverse problems, Technical Report NOW-14-67-A-0204-0045, M.Z.T.Report No.72-3, Dept. of Earth and Planetary Sciences, Massachusetts Institute of Technology.
|
Marquardt D W. 1963. Analgorithm for least-squares estimation of nonlinear parameters. Journal of the Society for Industrial and Applied Mathematics, 11(2): 431-441. DOI:10.1137/0111030 |
Nam M J, Kim H J. 2010. 3D MT inversion using an edge finite element modeling algorithm. Geosystem Engineering, 13(2): 43-52. DOI:10.1080/12269328.2010.10541308 |
Nam N J, Kim H J, Song Y, et al. 2007. 3D magnetotelluric modelling including surface topography. Geophysical Prospecting, 55(2): 277-287. DOI:10.1111/j.1365-2478.2007.00614.x |
Nam N J, Kim H J, Song Y, et al. 2009. Three-dimensional topographic and bathymetric effects on magnetotelluric responses in Jeju Island, Korea. Geophysical Journal International, 176(2): 457-466. DOI:10.1111/j.1365-246X.2008.03993.x |
Newman G A, Alumbaugh D L. 2000. Three-dimensional magnetotelluric inversion using non-linear conjugate gradients. Geophysical Journal International, 140(2): 410-424. DOI:10.1046/j.1365-246x.2000.00007.x |
Pan K J, Tang J T, Du H K, et al. 2016. Trust region inversion algorithm of high-resolution array lateral logging in axisymmetric formation. Chinese Journal of Geophysics (in Chinese), 59(8): 3110-3120. DOI:10.6038/cjg20160833 |
Parasnis D S. 1988. Reciprocity theorems in geoelectric and geoelectromagnetic work. Geoexploration, 25(3): 177-198. DOI:10.1016/0016-7142(88)90014-2 |
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 |
Rodi W L. 1976. A technique for improving the accuracy of finite element solutions for magnetotelluric data. Geophysical Journal International, 44(2): 483-506. DOI:10.1111/j.1365-246X.1976.tb03669.x |
Rodi W L, Mackie R L. 2001. Nonlinear conjugate gradients algorithmfor 2-D magnetotelluric inversion. Geophysics, 66(1): 174-187. |
Ruan B Y, Xu S Z, Xu Z F. 2007. Modeling the 3D terrain effect on MT by the boundary element method. Earth Science-Journal of China University of Geosciences (in Chinese), 32(1): 130-134. |
Sasaki Y. 2003. 3-D electromagnetic modelling and inversion incorporating topography. ASEG Extended Abstracts, 2003(1): 1-7. |
Schenk O, Gärtner K. 2004. Solving unsymmetric sparse systems of linear equations with PARDISO. Future Generation Computer Systems, 20(3): 475-487. DOI:10.1016/j.future.2003.07.011 |
Shi X, Utada H, Wang J, et al. 2004.Three dimensional magnetotelluric forward modelling using vector finite element method combined with divergence corrections(VFE++)//17th IAGA WG1.2 Workshop on electromagnetic Induction in the Earth. Hyderabad: [s.n.]: 465-473.
|
Shi X, Utada H, Wang J. 2009. 3-D magnetotelluric forward modeling and inversion incorporating topography by using vector finite-element method combined with divergence corrections based on the magnetic field (VFEH++).//American Geophysical Union, Fall Meeting 2009. AGU.
|
Siripunvaraporn W. 2012. Three-dimensional magnetotelluric inversion:an introductory guide for developers and users. Surveys in Geophysics, 33(1): 5-27. DOI:10.1007/s10712-011-9122-6 |
Siripunvaraporn W, Egbert G, Lenbury Y. 2002. Numerical accuracy of magnetotelluric modeling:A comparison of finite difference approximations. Earth, Planets and Space, 54(6): 721-725. DOI:10.1186/BF03351724 |
Tan H D, Wei W B, Deng M, et al. 2004. General use formula in MT tensor impedance. Oil Geophysical Prospecting (in Chinese), 39(1): 113-116. |
Tikhonov A N. 1950. On determining electrical characteristics of the deep layers of the earth's crust. Doklady Akademii Nauk, 73(2): 295-297. |
Vallianatos F. 2002. A note on the topographic distortion of magnetotelluric impedance. Annals of Geophysics, 45(2): 313-320. |
Wang G, Wei W B, Jin S, et al. 2017. A study on the electrical structure of eastern Gangdese metallogenic belt. Chinese Journal of Geophysics (in Chinese), 60(8): 2993-3003. DOI:10.6038/cjg20170808 |
Wang N, Zhao S S, Hui J, et al. 2016. Three-dimensional audio magnetotelluric sounding of coal-bed methane reservoirs in southern Qinshui basin. Progress in Geophysics (in Chinese), 31(6): 2664-2676. DOI:10.6038/pg20160642 |
Wannamaker P E, Stodt J A, Rijo L. 1986. Two-dimensional topographic responses in magnetotellurics modeled using finite elements. Geophysics, 51(11): 2121-2144. |
Wu X P, Wang T T. 2001. Study on some problems for 3D resistivity inversion using conjugate gradient. Seismology and Geology (in Chinese), 23(2): 321-327. |
Wu X P, Xu G M. 1999. Derivation and analysis of partial derivative matrix in resistivity 3-D inversion. Oil Geophysical Prospecting (in Chinese), 34(4): 363-372. |
Wu X P, Xu G M. 2000. Study on 3-D resistivity inversion using conjugate gradient method. Chinese Journal of Geophysics (in Chinese), 43(3): 420-427. |
Xu B W, Liu Z L, Ye G F, et al. 2018. Evaluation and analysis of geothermal resources in Ninghe uplift-evidence from magnetotelluric method. Progress in Geophysics (in Chinese), 33(6): 2278-2284. DOI:10.6038/pg2018BB0484 |
Xu S Z, Ruan B Y, Zhou H, et al. 1997. Numerical modeling of 3-D Terrain Effect on MT Field. Science in China Series D:Earth Sciences, 40(3): 269-275. DOI:10.1007/BF02877535 |
Xu S Z, Wang Q Y, Wang J. 1992. Modeling 2-D terrain effect on MT by the boundary element method. Acta Geophysica Sinica (in Chinese), 35(3): 380-388. |
Xue G Q, Yan S, Chen W Y. 2016. A fast topographic correction method for electromagnetic data. Chinese Journal of Geophysics (in Chinese), 59(12): 4408-4413. DOI:10.6038/cjg20161202 |
Yang W C, Jin S, Zhang L L, et al. 2020. The three-dimensional resistivity structures of the lithosphere beneath the Qinghai-Tibet Plateau. Chinese Journal of Geophysics (in Chinese), 63(3): 817-827. DOI:10.6038/cjg2020N0197 |
Yin C C, Zhang B, Liu Y H, et al. 2017. A goal-oriented adaptive algorithm for 3D magnetotelluric forward modeling. Chinese Journal of Geophysics (in Chinese), 60(1): 327-336. DOI:10.6038/cjg20170127 |
Zhang X, Hu W B. 1999. Joint 2-D inversion of magnetotelluric sounding data having landform effect. Oil Geophysical Prospecting (in Chinese), 34(2): 190-196. |
陈小斌. 2000. MT二维正演计算中地形影响的研究. 石油物探, 39(3): 112-120. DOI:10.3969/j.issn.1000-1441.2000.03.015 |
陈小斌, 赵国泽, 汤吉, 等. 2005. 大地电磁自适应正则化反演算法. 地球物理学报, 48(4): 937-946. DOI:10.3321/j.issn:0001-5733.2005.04.029 |
崔腾发, 陈小斌, 詹艳, 等. 2020. 安徽霍山地震区深部电性结构和发震构造特征. 地球物理学报, 63(1): 256-269. DOI:10.6038/cjg2019M0458 |
董浩, 魏文博, 叶高峰, 等. 2014. 基于有限差分正演的带地形三维大地电磁反演方法. 地球物理学报, 57(3): 939-952. DOI:10.6038/cjg20140323 |
顾观文, 吴文鹂, 李桐林. 2014. 大地电磁场三维地形影响的矢量有限元数值模拟. 吉林大学学报(地球科学版), 44(5): 1678-1686. |
胡祖志, 何展翔, 孙卫斌, 等. 2008. 一种改进的大地电磁地形校正方法及相关问题探讨. 石油地球物理勘探, 43(3): 343-348. DOI:10.3321/j.issn:1000-7210.2008.03.020 |
黄颖, 王雪秋, 周子坤, 等. 2019. 聚类分析算法在大地电磁三维解释中的应用. 地球物理学进展, 34(2): 568-572. DOI:10.6038/pg2019CC0332 |
霍光谱, 胡祥云, 黄一凡, 等. 2015. 带地形的大地电磁各向异性二维模拟及实例对比分析. 地球物理学报, 58(12): 4696-4708. DOI:10.6038/cjg20151230 |
晋光文, 赵国泽, 徐常芳, 等. 1998. 二维倾斜地形对大地电磁资料的影响与地形校正. 地震地质, 20(4): 454-458. |
金建铭. 1998.电磁场有限元方法.王建国译.西安: 西安电子科技大学出版社, 179-189.
|
刘云, 王绪本. 2010. 大地电磁二维自适应地形有限元正演模拟. 地震地质, 32(3): 382-391. DOI:10.3969/j.issn.0253-4967.2010.03.004 |
刘云鹤, 殷长春. 2013. 三维频率域航空电磁反演研究. 地球物理学报, 56(12): 4278-4287. DOI:10.6038/cjg20131230 |
潘克家, 汤井田, 杜华坤, 等. 2016. 轴对称地层中高分辨率阵列侧向测井信赖域反演法. 地球物理学报, 59(8): 3110-3120. DOI:10.6038/cjg20160833 |
阮百尧, 徐世浙, 徐志锋. 2007. 三维地形大地电磁场的边界元模拟方法. 地球科学-中国地质大学学报, 32(1): 130-134. DOI:10.3321/j.issn:1000-2383.2007.01.020 |
谭捍东, 魏文博, 邓明, 等. 2004. 大地电磁法张量阻抗通用计算公式. 石油地球物理勘探, 39(1): 113-116. DOI:10.3321/j.issn:1000-7210.2004.01.021 |
王刚, 魏文博, 金胜, 等. 2017. 冈底斯成矿带东段的电性结构特征研究. 地球物理学报, 60(8): 2993-3003. DOI:10.6038/cjg20170808 |
王楠, 赵姗姗, 惠健, 等. 2016. 沁水盆地南部煤层气藏三维音频大地电磁探测. 地球物理学进展, 31(6): 2664-2676. DOI:10.6038/pg20160642 |
吴小平, 汪彤彤. 2001. 利用共轭梯度方法的电阻率三维反演若干问题研究. 地震地质, 23(2): 321-327. DOI:10.3969/j.issn.0253-4967.2001.02.028 |
吴小平, 徐果明. 1999. 电阻率三维反演中偏导数矩阵的求取与分析. 石油地球物理勘探, 34(4): 363-372. |
吴小平, 徐果明. 2000. 利用共轭梯度法的电阻率三维反演研究. 地球物理学报, 43(3): 420-427. DOI:10.3321/j.issn:0001-5733.2000.03.016 |
胥博文, 刘志龙, 叶高峰, 等. 2018. 宁河凸起地热资源远景区评价分析——来自大地电磁测深法的证据. 地球物理学进展, 33(6): 2278-2284. DOI:10.6038/pg2018BB0484 |
徐世浙, 阮百尧, 周辉, 等. 1997. 大地电磁场三维地形影响的数值模拟. 中国科学(D辑), 27(1): 15-20. |
徐世浙, 王庆乙, 王军. 1992. 用边界单元法模拟二维地形对大地电磁场的影响. 地球物理学报, 35(3): 380-388. DOI:10.3321/j.issn:0001-5733.1992.03.012 |
薛国强, 闫述, 陈卫营. 2016. 电磁测深数据地形影响的快速校正. 地球物理学报, 59(12): 4408-4413. DOI:10.6038/cjg20161202 |
杨文采, 金胜, 张罗磊, 等. 2020. 青藏高原岩石圈三维电性结构. 地球物理学报, 63(3): 817-827. DOI:10.6038/cjg2020N0197 |
殷长春, 张博, 刘云鹤, 等. 2017. 面向目标自适应三维大地电磁正演模拟. 地球物理学报, 60(1): 327-336. DOI:10.6038/cjg20170127 |
张翔, 胡文宝. 1999. 带地形的大地电磁测深联合二维反演. 石油地球物理勘探, 34(2): 190-196. DOI:10.3321/j.issn:1000-7210.1999.02.008 |