2012, Vol. 55
本文所研究的地面可控源频率测深是沿用朴化荣[1]以及Kaufman 和Keller等[2]关于地面可控源电磁测深的描述，源采用有限长电流线，频率域发射，观测源周围电磁场分布.在一般意义上，反演的数据不是以单个测点为单位，而是以整个工区内在所有测点上对所有观测分量做的观测为输入数据.因此，这里的反演区别于目前流行的可控源音频大地电磁测深法反演[3]，具体表现在:观测可以在靠近源的测区进行，此时，源的存在和形状以及源与测点间几何关系必须加以考虑；反演的数据可以是直接观测得到的电磁场各个分量，也可以是由这些分量导出的其它量，如视电阻率等，但不必是大地电磁测深意义上的阻抗，因为在那种情况下源的部分信息可能丢失[4].因此，关于地面可控源频率测深数据三维反演，与可控源音频大地电磁测深法采用平面波源的大地电磁测深三维反演有本质的区别，即这里必须采用考虑场源的电磁测深三维反演技术.
近年来，电磁勘探三维数据反演技术研究取得了长足进步，Avdeev[5]，汤井田等[6]，赵国泽等[7]以及Brner等[8]系统回顾了国内外在这个研究领域所取得的成果，并概括体现在高精度三维数值模拟、灵敏度元素计算以及反演方法技术等三个核心方面，里程碑是Alumbaugh 和Newman等[9]在Mackie等[10]工作基础上提出的三维共轭梯度并行电磁反演技术，其理论框架是用交错网格有限差分完成正演模拟，采用线性共轭梯度技术实现反演，借助并行计算技术提高反演速度.结合地球物理模型非线性特点，Newman等[11]进一步提出采用非线性共轭梯度反演方法，并以大地电磁测深三维反演为例进行了方法验证.由于非线性共轭梯度反演过程中将灵敏度计算融合到梯度计算中，不仅总的数值计算量成比下降，同时极大降低了对计算资源的苛刻要求，从而为在相对小型服务器上进行三维电磁反演研究奠定了理论基础[5].正是由于非线性共轭梯度反演的实现，三维电磁反演研究进入崭新的时代[12-15]，并应用在大地电磁测深技术中[10, 11, 16-18].正如Commer等[19]所指出，非线性共轭梯度方法是目前三维电磁反演最有效的方法，因此该方法最近被应用到海洋可控源电磁法的反演当中.
由于电磁理论的普适性以及反演理论数学上的通用性，本文将目前非线性共轭梯度反演所取得的研究成果应用到地面可控源电磁测深数据反演中.在反演过程中，基于异常电场满足的微分方程，借助交错网格有限差分技术实现三维正演模拟.非线性共轭梯度反演针对电场的Ex分量进行.这些数据由理论模型响应叠加噪音后得到.反演的模型包括:(1)层状背景模型中含有高阻和低阻两个异常体模型；(2)包含阴影效应的模型.反演结果验证了非线性共轭梯度反演在地面可控源电磁法资料三维反演中的有效性.
2 正演数值模拟 2.1 数学物理模型在地面可控源频率电磁勘探中，通常采用有限长导线作为发射源，在源周围观测电磁场分布.图 1是该观测方式的数学物理描述.三维异常体位于层状均匀各向同性背景电阻率中.异常电阻率不必是均匀的，可以表示为空间位置的函数ρ(r).N层层状背景模型各层电阻率和厚度分别为(ρi，hi，i=1，…，N)，各层磁导率和介电常数与自由空气的相等.数值计算中采用直角坐标系，坐标系原点在地面.设地面导线源沿x轴方向布置，中心在(xs，ys，0)处，长度为2L，电流I，频率ω，电偶极矩PE=2LI.电流的时谐因子取e-iωt.观测在地面上P(x，y，0)点进行，到源中心点距离矢量为r′.这样，任意一点电场E满足如下的Helmholtz方程[20]:
|
(1) |
|
图 1 可控源频率电磁法三维反演数学物理模型 层状背景模型中包含电阻率三维不均匀性，有限长电流源激发，地面观测电磁场分量. Fig. 1 Schematical model of three dimensional controlled source electromagnetic method A current line of finite length is located at the surface of a layered background，where there exists 3D resistivity abnormal. The electric or magnetic fields will be received at the earth surface. |
式中k2 = (iωμ0)/ρ，μ0 是空气磁导率，i=
为了克服地面可控源电磁法在数值模拟过程中源描述的困难以及源奇异性对数值稳定性的影响，正演数值模拟针对异常二次电场进行[19].为此，将总电场E分解为由层状背景介质产生的基本场Ep和三维异常体产生的异常(二次)场Es 之和，即
|
(2) |
则总电场满足的Helmholtz方程可以转化为如下的异常电场满足的微分方程:
|
(3) |
其中kp2 为背景介质复波数，k2 是含异常后介质复波数.利用(3)式计算三维电磁场分布的优势是，二次场只和基本场有关，而与产生基本场的源无关.因此，该方程在形式上提供了通用的三维电磁模拟框架，即对于不同的源，只需要提供基本场，通过(3)式就可以得到其三维电磁场分布.对于地面可控源频率电磁法，地面导线源在层状介质中任意点产生的基本场计算采用虚界面方法[21-22]，并采用离散复镜像法完成数值计算[23].
2.3 交错网格有限差分法利用交错网格有限差分法求解(3)式时，以测区为中心把地下空间(实际还包括上面的空气)剖分成长方体网格，并用节点编号(i，j，k)标注网格中心位置.这样，每一个长方体单元的电阻率值标记为ρ(i，j，k).这里，i的范围是1到Nx，代表x方向的网格序号.j和k代表y和z方向的网格序号，范围分别从1 到Ny和1 到Nz.当数值模拟电场分布时，电场位于网格单元边界中点，磁场分布在长方体面中心[24-25].(3)式是矢量微分方程，按照Alumbaugh等[25]的思路，可以分解为对应于3 个坐标矢量的3个微分方程.合并分别离散化后得到的矩阵方程，并进行对称化处理，最终形成如下的矩阵方程:
|
(4) |
K是复系数大型稀疏对称矩阵；s包括由背景场以及电阻率异常确定的(3)式右手源矢量和边界条件.假设反演模拟区域远大于三维异常分布范围，可以假设二次场在网格外边界为0，因此s可不必包含虑边界条件.Es 为各长方体单元边界中点处异常电场矢量(包含直角坐标三个方向分量).图 2 给出了离散化(3)式在第(i，j，k)单元上需要的电场分量在网格空间分布的规律.对于节点(i，j，k)(图中粗黑点)，其对应的电场分量为(Ex(i+1/2，j，k)，Ey(i，j+1/2，k)，Ez(i，j，k+1/2))，节点编号取±1/2，是由于电场位于单元边界中点.图中粗箭头对应的电场矢量用于构成x方向矢量满足的微分方程.
|
图 2 基于电场Helmholtz方程的有限差分交错网格 中心处的黑点代表节点（i，j，k），该处电场有三个对应的分量. 13个大的黑色箭头代表构成（i，j，k）处x方向方程的未知电场，其它的小的黑色箭头为附加的电场用来构成y和z方向的方程，据文献[11]. Fig. 2 Staggered grid used to modeling Helmholtz equation based on electric field The bold black point at the grid center represents node (i,j,k),associated with three electrical components Ex,Ey，and Ez. The big thirteen arrows are used to construct the finite differential equation in x direction,while the small arrows are used to form the equations in y and z direction,respectively. (After references [11]). |
在数值模拟过程中，先后引入了预条件技术[26-27]、Krylov 子空间迭代方法[16, 28]以及散度校正[29-30]等处理手段，极大提高了(4)式的计算效率和精度.
2.4 观测数据插值算子一旦求出网格节点上的总电场，利用插值算子可以获得空间任意一点处的电场，并借助法拉第定律微分形式计算出磁场分量，进而可以由电、磁场分量组合得到观测点处理论响应:
|
(5) |
式中f是ζi的函数；而ζi为插值算子，利用它由三维模拟得到的总电场E(参见(2)式)插值计算观测点处理论数据所需要的某个电、磁场分量；Kd为获得观测数据需要使用的电磁场分量个数.m为离散的三维模型数据矢量.在试验中，反演直接针对电场的Ex分量进行.此时，(5)式可以简化成
|
(6) |
式中，Ψ 为由总电场插值得到Ex分量的单个线性插值算子.
3 非线性共轭梯度反演 3.1 基本原理将非线性共轭梯度用于地面可控源频率观测数据的反演，遵从非线性共轭梯度反演数学形式的通用性.通常定义如下目标函数[9]:
|
(7) |
式中，m为待反演模型电阻率参数矢量，长度M；m0为初始模型(先验模型)矢量；mref为参考模型矢量；Cm为模型约束矩阵；Cd为数据方差矩阵；d是长度为N的观测数据矢量，这里是电场的Ex分量.长度N取决于发射频点个数、观测点数量以及观测参数个数.
非线性共轭梯度法通过迭代获得该目标函数的极点.首先通过计算该目标函数某参考模型的负梯度来得到共轭梯度，并在这个方向上而不是简单的梯度方向上搜索一个最优步长使目标函数达到极小值.这样得到的解就应该是目标函数在目前共轭梯度方向的极点.由于观测响应是模型参数的非线性函数，通常利用这个极点作为新的模型进行迭代，重复上述步骤直到收敛.具体计算流程如下[11]:
(1) 令i=1，选择初始模型m(i)=mref并计算r(i) =-ΔV(m(i))；
(2) 令u(i) = M(i)-1r(i)，M为预处理算子；
(3) 搜索最优步长α(i)来最小化V(m(i)+αu(i))；
(4) 令m(i+1)=m(i)+α(i)u(i)，r(i+1)=-ΔV(m(i+1))；
(5)如果‖r(i+1)‖ 满足要求则停止计算，否则跳到(6)；
(6) 令β(i+1)=
(7) 令u(i+1) = M(i+1)-1r(i+1)+β(i+1)Tu(i)；
(8) 令i=i+1，返回到步骤(3).
从上面的计算过程可以看出，应用非线性共轭梯度法时主要计算有两个，一个是目标函数梯度的计算，另一个是对于某个特定的模型参数m和指定的共轭搜索方向u，要找出一个最优步长α(i)使V(m+αu)最小化.最优步长α(i)的确定是一维优化问题，可以采用割线法进行[15].因此，梯度计算是非线性共轭梯度反演的关键.
3.2 梯度计算由式(7)，目标函数梯度g表示为(8)
式中J ik=
|
(9) |
则式(8)可写成
|
(10) |
从表面上看，为计算梯度g，只需要计算出J 后直接带入(10)式即可.附录A 给出了J 的具体计算公式[31].但如果显式计算J 并存储矩阵元素，需要庞大的内存和正演计算[32-33].下面从J 或者J T 与向量乘积的角度来分析说明为了计算梯度，只需要相当于2 次正演的计算量并且不需要显式存储J [9].
如果简记(10)式中右手第二项括号中的第二项灵敏度矩阵与模型向量积为y=J x，将(A7)式代入有
|
(11) |
式中xP =Px.继续令乘积K-1xP =v，两边乘以K，有
|
(12) |
比较(4)式可见，计算v事实上就是利用新的源项进行一次正演.最终有y=Lv.同样，令z=J Ty，由(A8)式，有z=PT (KT)-1LTy.若记yL=LTy，并令 (KT)-1yL=u，则
|
(13) |
上式是梯度计算的第二次正演.这样，通过(12)式和(13)式的2 次正演以及4 次矩阵和向量的乘积运算，完成了梯度计算.
综上所述，在1次非线性共轭梯度反演迭代中，为得到共轭梯度需要相当于2 次单独正演的计算量，外加在第(3)步中，3 次正演确定最优步长以及一次模型响应正演计算，完成一次迭代最少总共需要6次正演计算.在发射源个数少于观测数据个数时，整体计算时间大大小于直接计算J 的方法[5].这也是非线性共轭梯度在三维电磁反演中得以推崇的根本原因之一.
4 反演结果分析 4.1 理论模型一图 3模型说明本文反演技术的效果.该物理模型包含一个高阻和一个低阻异常体，它们被放置在三层层状背景模型的第一层中.高阻异常体靠近地表，长和宽都是800m，厚度为200m，电阻率为300Ωm；低阻异常体埋深200 m，长和宽都是400 m，厚度200m，电阻率为50Ωm.线电流源长度1000 m，距离测区最近测线2700m，发射频率分别为5.12 Hz、32、128Hz和8192 Hz.测量在图 3a 中2000 m×2000m 的灰色区域进行，观测参数为Ex分量的实部和虚部，点距100m，线距100m，共441个物理测点.因此，反演数据总个数为3528个.反演区域与测量范围一致，并在x，y和z三个方向被剖分为41×41×21共35301个网格，网格大小为100m×100m×100m.
|
图 3 包含两个异常体的反演理论模型 (a)异常及几何装置顶视图；(b)模型过x = 4600 m的纵断面. Fig. 3 Theoretical model used to check the validation of the inversion method Two abnormals are embedded into the first layer of a stratified background. The shallower one is a resistive slab，and the deeper smaller body is a relative conductor. Source AB is 2700m away from the survey zone (grayed rectangle) of size 2000m×2000 m. The abnormals are in the central area of the zone. Spacing between observing stations is 100 m. (a) The top view； (b) The section view from x=4600 m |
理论模拟数据的实部和虚部分别被施以5%的0均值高斯噪音后作为观测数据参与反演.反演在4CPU 的兼容工作站上进行串行运算，IntelXeonE5620处理器，主频2.4GHz，内存16G.反演从层状背景电阻率模型开始，迭代了120次，用时24.2h.图 4给出了迭代过程中均方相对误差(RMS)的变化情况.从图可以看出，随着迭代的进行，数据拟合误差逐渐减小，在一定迭代次数后，误差减小变慢.因此，对于非线性共轭梯度反演，为了将拟合差控制到理论上的1附近，可以预见，还需要较多的迭代计算.这也反映了非线性共轭梯度反演收敛速度较慢的缺点.
|
图 4 120次迭代过程中RMS变化曲线 Fig. 4 RMS variation with iteration number after 120 iterations |
图 5为反演结果模型在两个深度上的切片及相关的纵剖面.图 5a切片深度在100m，该平面横切浅部的高阻异常中心.在反演电阻率的平面上，高阻异常分布在(3900，3900)×(4600，4300)范围内，其平面位置与图 3 浅部高阻异常体对应.同时，从图 5b与图 5c中过高阻异常体的纵剖面可以看出，反演的高阻异常主要分布在地面近200m 深范围内，在垂向上也基本与理论模型一致.
|
图 5 三维反演电阻率分布图 (a) z=100 m处平面图；(b) y = 4100 m处z方向剖面；(c)z= 4100 m处y方向剖面； (d) z=300 m处平面图；(e) y=4900 m处z方向剖面；（f）z= 4900 m处y方向剖面. Fig. 5 3D resistivity distribution at the end of inversion (a) is the plan view at z=100 m，while (b) and (c) are the cross sections at y=4100 m and x=4100 m，respectively. (c) is the resistivity isotropic contour inner the earth at depth z=300 m，with sections crossed at y=4900 m in (e) and x=4900 m in (f). |
深度300m 的水平切片上的电阻率分布如图 5d.图中低阻异常平面上大致分布在(4600，4400)×(5000，5200)范围内，和设计的异常基本吻合.过x=4900m和y=4900m 的垂直正交纵剖面(图 5e和图 5f)中反演得到的低阻体埋深在200~500 m之间，也反映了理论模型的垂直分布.但综合图 5的结果，需要指出的是，在x方向上，异常体的位置和大小反演控制较好；而在y方向上，高阻体被压短，低阻体被拉长.这可能是由于在反演过程中，只采用Ex进行反演所致.
4.2 理论模型二图 6所示的模型用来说明源在可控源电磁数据三维反演中的作用.在虚线灰色测区范围之外，在源方向离测区100m 位置上，在三层背景模型的第一层中放置一个电阻率为10Ωm，顶部埋深100m 的低阻异常体，其在x、y和z三个方向的维度为500m×200m×400m.测区无异常电阻率存在.
|
图 6 测区外存在异常体的理论模型 (a)异常及几何装置顶视图；(b)模型过x = 4600 m的纵断面. Fig. 6 A model with an anormal located between the source and survey area and outside of the area (a) Shows the position of the anormal, along with that of transmitter and receivers on the ground; (b) The section crossing perpendicular to the abnormal in x,y,z plan through x= 4600 m |
图 7 给出了测区外的异常体对测量结果的影响.从结果可以看出，由于异常体的存在，测区电场受到了屏蔽而变形，靠近异常体的影响大，随测点距异常体距离增大屏蔽影响逐渐减小.显然，如果不考虑场源影响，直接基于音频大地电磁反演技术进行反演，畸变的电场必然会在测区造成虚假电阻率异常.
|
图 7 测区外异常体对测区内电场Ex分量的影响 The effects of the anormal on the real part of the electrical fields in the survey zone Fig. 7 图中参数为Re[Ex一Ex0]×1010，Re表示实部，Ex为存在局部异常体时电场，Ex0为背景电场.(a)f-2 Hz;(b)f=512 Hz. Two frequencies are transmitted. (a) is 2 Hz, and (b) is 512 Hz. Parameter in the figures is Re[Ex-Ex0] , where is the electrical field of background, while Ex represents the electrical field from the layered model with an abnormal； and Re denotes real part of a complex. |
反演中为了考虑局部异常体分布，将反演范围扩到到(3500，3000)×(5500，6000)，所用其它参数同图 3.图 8为由5%噪音扰动后的理论数据三维反演得到的结果，图中过y=4000 m 的横线表示测量区域源侧边界线.从图 8a和图 8b两个不同深度电阻率切片可以看出低阻异常向实际位置聚焦，低阻异常只分布在测区源侧边界200m 范围内，更远处电阻率基本不受影响.比较图 7可以看出，虽然异常体对测区内电场响应影响很大，但是在反演结果中该影响基本得到有效控制.此外，在图 8中的不同位置y向纵断面(c、d、e)和(f)中，异常在深度上主要分布在测区边界靠源一侧，对测区下部基本没有影响.这说明采用考虑场源的真正可控源三维反演基本不存在所谓的阴影效应[34].
|
图 8 测区外存在异常体的反演结果 (a) z=100m处平面图；(b) z=300 m处平面图；(c) x = 4300 m处y方向剖面；(d) x=4500 m处 y方向剖面；(e) x = 4600 m处y方向剖面；(f) x = 4700 m处y方向剖面. Fig. 8 Inversed resisitivity distribution of the model in Fig. 6 (a) and (b) are the plan views at y=100 m and 300 m,while (c)，(d)，(e) and (f) are the cross sections at x=4300 m，4500 m，4600 m and 4700 m，respectively. |
本论文将非线性共轭梯度反演算法用于地面可控源频率域电磁法数据的三维反演中，通过数值模拟试验，获得了如下认识:
(1) 理论模型反演结果表明论文中提出的反演技术线路是可行的.
(2) 非线性共轭梯度反演技术收敛较慢，完成反演需要较多的迭代次数.
(3) 对于地面可控源频率勘探资料，必须在反演中考虑场源几何信息，才能压制测区外异常对目标异常的影响，得到可靠的反演结果.
在本文工作的基础上，今后还要开展多分量观测数据的三维反演，并讨论影响反演结果的各种因素，如初始模型对反演结果的影响等问题.
附录A 灵敏度表达式灵敏度矩阵表达式推导遵从Egbert等[31]的框架.由观测数据响应表达式(6)式，对于第i个观测，有
|
(A1) |
则利用连续求导法则有
|
(A2) |
若令
|
可以写出
|
(A3) |
由于L可以解析计算，而观测装置参数通常是独立模型参数，故Q=0.因此，为计算J ，关键是计算F.如果源被看作是与模型参数无关的(实际情况往往如此)，由于E满足(4)式，两边对模型参数求导，有
|
(A4) |
上式两边同乘以K-1有
|
(A5) |
进一步令
|
(A6) |
从而最终有
|
(A7) |
及
|
(A8) |
| [1] | 朴化荣. 电磁测深法原理. 北京: 地质出版社, 1990 . Piao H R. Electromagnetic Soundings (in Chinese). Beijing: Geological Publishing House, 1990 . |
| [2] | Kaufman A A, Keller G V. Frequency and Transient Soundings. New York: Elsevier, 1983 . |
| [3] | 汤井田, 何继善. 可控源音频大地电磁法及其应用. 长沙: 中南工业大学出版社, 2005 . Tang J T, He J S. Controlled Source Audio Frequency Electromagnetic Method and Its Applications (in Chinese). Changsha: Central South University Press, 2005 . |
| [4] | 万扬Л Л著, 孔祥儒, 曾治权 译. 电磁测深. 北京: 海洋出版社, 2001. Ваньян Л Л. Electromagnetic Soundings (in Chinese). Translated by Kong X R, Zeng Z Q. Beijing: Ocean Press, 2001. |
| [5] | Avdeev D B. Three-dimensional electromagnetic modelling and inversion form theory to application. Surveys in Geophysics , 2005, 26(6): 767-799. DOI:10.1007/s10712-005-1836-x |
| [6] | 汤井田, 任政勇, 化希瑞. 地球物理学中的电磁场正演与反演. 地球物理学进展 , 2007, 22(4): 1181–1194. Tang J T, Ren Z Y, Hua X R. The forward modeling and inversion in geophysical electromagnetic field. Progress in Geophysics (in Chinese) , 2007, 22(4): 1181-1194. |
| [7] | 赵国泽, 陈小斌, 汤吉. 中国地球电磁法新进展和发展趋势. 地球物理学进展 , 2007, 22(4): 1171–1180. Zhao G Z, Chen X B, Tang J. Advanced geo-electromagnetic methods in China. Progress in Geophysics (in Chinese) , 2007, 22(4): 1171-1180. |
| [8] | Börner R U. Numerical modelling in geo-electrcmagnetics: Advances and Challenges. Surveys in Geophysics , 2010, 31(2): 225-245. DOI:10.1007/s10712-009-9087-x |
| [9] | Alumbaugh D L, Newman G A. Three-dimensional massively parallel electromagnetic inversion-II. Analysis of a crosswell electromagnetic experiment. Geophysical Journal International , 1997, 128(2): 355-363. DOI:10.1111/gji.1997.128.issue-2 |
| [10] | Mackie R L, Madden T R. Three-dimensional magnetotelluric inversion using conjugate gradients. Geophysical Journal International , 1993, 115(1): 215-229. DOI:10.1111/gji.1993.115.issue-1 |
| [11] | Newman G A, Alumbaugh D L. Three-dimensional magnetotelluric inversion using non-linear conjugate gradients. Geophysical Journal International , 2000, 140(2): 410-424. DOI:10.1046/j.1365-246x.2000.00007.x |
| [12] | Zhdanov M S. Geophysical Inverse Theory and Regularization Problems. New York: Elsevier, 2002 . |
| [13] | Newman G A, Boggs P T. Solution accelerators for large-scale there-dimensional elect romagnetic inverse problem. Inverse Problems , 2004, 20(6): 151-170. DOI:10.1088/0266-5611/20/6/S10 |
| [14] | Siripunvaraporn W, Egbert G. Data space conjugate gradient inversion for 2-D magnetotelluric data. Geophysical Journal International , 2007, 170(3): 986-994. DOI:10.1111/gji.2007.170.issue-3 |
| [15] | Kelbert A, Egbert G D, Schultz A. Non-linear conjugate gradient inversion for global EM induction: resolution studies. Geophysical Journal International , 2008, 173(2): 365-381. DOI:10.1111/gji.2008.173.issue-2 |
| [16] | Rodi W, Mackie R L. Nonlinear conjugate gradients algorithm for 2-D magnetotelluric inversion. Geophysics , 2000, 66(1): 174-187. |
| [17] | 谭捍东, 余钦范, BookerJ, 等. 大地电磁法三维快速松弛反演. 地球物理学报 , 2003, 46(6): 850–855. Tan H D, Yu Q F, Booker J, et al. Three-dimensional rapid relaxation inversion for the magnetotelluric method. Chinese Journal of Geophysics (in Chinese) , 2003, 46(6): 850-855. |
| [18] | 胡祖志, 胡祥云, 何展翔. 大地电磁非线性共轭梯度拟三维反演. 地球物理学报 , 2006, 49(4): 1226–1234. Hu Z Z, Hu X Y, He Z X. Pseudo-Three dimensional magnetotelluric inversion using nonlinear conjugate gradients. Chinese J. Geophys. (in Chinese) , 2006, 49(4): 1226-1234. |
| [19] | Commer M, Newman G A. New advances in three-dimensional controlled-source electromagnetic inversion. Geophysical Journal International , 2008, 172(2): 513-535. DOI:10.1111/gji.2008.172.issue-2 |
| [20] | Nabighian M N. Electromagnetic Methods in Applied Geophysics. Vol.1, SEG, Tulsa, Oklahoma, 1987. |
| [21] | Da U C. A reformalism for computing frequency- and time-domain EM responses of a buried finite-loop. 65th Ann. Internat. Mtg., Soc. Expi. Geophys. Expanded Abstracts, 1995: 811-814. |
| [22] | 刘云鹤, 刘国兴, 翁爱华, 等. 基于二级近似离散复镜像法的低频电磁格林函数计算. 地球物理学报 , 2011, 54(4): 1114–1121. Liu Y H, Liu G X, Weng A H, et al. Calculation of low-frequency electromagnetic Green's functions using two-level approximate discrete complex image method. Chinese J. Geophys. (in Chinese) , 2011, 54(4): 1114-1121. |
| [23] | 刘云鹤. 三维可控源电磁法非线性共轭梯度反演研究. 长春: 吉林大学, 2011 . |
| [24] | Yee K S. Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media. IEEE Trans. Ant. Prop. , 1966, 14(3): 302-307. |
| [25] | Alumbaugh D L, Newman G A, Prevost L, et al. Three-dimensional wide band electromagnetic modeling on massively parallel computers. Radio Science , 1996, 31(1): 1-23. DOI:10.1029/95RS02815 |
| [26] | Avdeev D B, Kuvshinov A V, Pankratov O V, et al. Three-dimensional induction logging problems. Part I: An integral equation solution and model comparisons. Geophysics , 2002, 67(2): 484-491. DOI:10.1190/1.1468608 |
| [27] | Siripunvaraporn W, Egbert G, Lenbury Y. Numerical accuracy of magnetotelluric modeling: A comparison of finite difference approximations. Earth Planets Space , 2002, 54: 721-725. DOI:10.1186/BF03351724 |
| [28] | Mackie R L, Madden T R. Conjugate direction relaxation solutions for 3-D magnetotelluric modeling. Geophysics , 1993, 58(7): 1052-1057. DOI:10.1190/1.1443481 |
| [29] | Mackie R L, Smith J T, Madden T R. Three-dimensional electromagnetic modeling using finite difference equations: The magnetotelluric example. Radio Science , 1994, 29(4): 923-935. DOI:10.1029/94RS00326 |
| [30] | Smith J T. Conservative modeling of 3-D electromagnetic fields, Part II: Biconjugate gradient solution and an accelerator. Geophysics , 1996, 61(5): 1319-1324. DOI:10.1190/1.1444055 |
| [31] | Egbert G D, Kelbert A. Computational recipes for electromagnetic inverse problems. Geophysical Journal International , 2012, 189(1): 251-267. DOI:10.1111/gji.2012.189.issue-1 |
| [32] | 吴小平, 徐果明. 电阻率三维反演中偏导数矩阵的求取与分析. 石油地球物理勘探 , 1999, 34(4): 363–372. Wu X P, Xu G M. Derivation and analysis of partial derivative matrix in resistivity 3-D inversion. Oil Gas Prospecting (in Chinese) , 1999, 34(4): 363-372. |
| [33] | 吴小平, 徐果明. 利用共轭梯度法的电阻率三维反演研究. 地球物理学报 , 2000, 43(3): 420–427. Wu X P, Xu G M. Study on 3-D Resistivity inversion using conjugate gradient method. Chinese J. Geophys. (in Chinese) , 2000, 43(3): 420-427. |
| [34] | Zonge K L, Hughes L J. Controlled source audio-frequency magnetotellurics //Nabighian M N, ed. Electromagnetic Methods in Applied Geophysics. Society of Exploration Geophysicists, 1988, 2: 713-809. |