地球物理学进展  2015, Vol. 30 Issue (2): 829-835   PDF    
引用本文
杨悦, 李亚彬, 翁爱华, 李大俊, 李斯睿, 杨大方. 2015. 频率选择对三维可控源反演的影响[J]. 地球物理学进展, 30(2): 829-835, doi: 10.6038/pg20150246  
YANG Yue, LI Ya-bin, WENG Ai-hua, LI Da-jun, LI Si-rui, YANG Da-fang. 2015. Effect of the choice of frequency on three-dimensional controlled source electromagnetic inversion. Progress in Geophysics, 30(2): 829-835, doi: 10.6038/pg20150246   
频率选择对三维可控源反演的影响
杨悦, 李亚彬, 翁爱华 , 李大俊, 李斯睿, 杨大方    
吉林大学地球探测科学与技术学院, 长春 130026
收稿日期: 2014-09-08 修回日期: 2015-01-06 .
基金项目：国家重大科研仪器专项(2011YQ05006010)资助.
作者简介：杨悦,女,1989年生,硕士,主要从事电磁法勘探理论模拟与应用研究.(E-mail:jlyangyue@126.com)
通讯作者：翁爱华,男,1969年生,教授,博导,从事电磁法勘探理论问题、三维电磁反演方法技术与地面核磁共振的数据分析与应用等研究.(E-mail:wengah@jlu.edu.cn)
摘要：本文研究了频率选择对三维可控源方法勘探效果的影响.研究的三维模型抽象自泥河铁矿.该铁矿规模为1000 m×1000 m×400 m, 埋深900 m.反演采用有限内存拟牛顿法, 反演过程计算采用有限差分法.观测数据为与发射源平行的Ex实部和虚部.反演初始模型为两层背景模型时, 选取不同频率组合进行反演, 迭代120次.数值试验结果表明, 不同频率组合反演效果不同, 存在最佳的频率组合.在反演初始模型合适时, 存在关键频率组合.对于本次模型, 关键频率组合为16 Hz+64 Hz, 即反演频率个数最少, 反演结果与理论模型最接近的频率组合.同时从反演结果可看出, 对于相对低阻盖层, 可控源的探测深度可以达到 1000 m.因此, 在野外观测时, 如果能获得合适的初始模型, 可以先根据该地区的地质特点设计模型, 进行数值模拟, 选出最佳频率组合, 指导三维可控源电磁测深方法野外数据采集, 提高工作的效率和准确性.
关键词： 可控源电磁法     三维     有限内存拟牛顿反演法     频率组合    
Effect of the choice of frequency on three-dimensional controlled source electromagnetic inversion
YANG Yue, LI Ya-bin, WENG Ai-hua , LI Da-jun, LI Si-rui, YANG Da-fang    
College for Geoexploration of Science and Technology, Jilin University, Changchun 130026, China
Abstract: Research on the effect of the choice of frequency on three-dimensional controlled source electromagnetic inversion.Three-dimensional model is designed based on the Nihe icon. The icon scale is 1000 m×1000 m×400 m,900 m depth.During inversion,Helmholtz equation based on secondary electric field governed by primary field as source is first computed by using staggered grid finite difference scheme.Inversion uses limited memory quasi-Newton method.The x-component of electric fields Ex with constant noises are put into inversion. When the background model is the two-layers model,choose combinations of different frequencies into inversion and terminate after 120 iterations.Numerical results show that it is different between inversions use combinations of different frequencies,and there is a combination of optimal frequencies. If the background model is appropriate, a combination of critical frequencies exits.For this model,the combination of critical frequencies is 16 Hz+64 Hz,which the number of frequencies used in inversion is minimum,and three-dimensional resistivity distributions from inversion coincide with the theoretical models. At the same time,the results also show that the exploration depth of controlled source electromagnetic can reach 1000 m. Therefore in field surveys,to guide the collection of field data,and improve the efficiency and accuracy,we can select the combination of critical frequencies using numerical algorithm of the three-dimensional electromagnetic model designed by the geologicalcharacteristics of the area.
Key words: controlled source electromagnetic method     three dimension     limited memory quasi-Newton method     combination of frequencies    
0 引 言

可控源电磁测深方法(CSEM)克服了天然场源信号的随机性和强度弱的缺点(汤井田等，2005)，已成为解决深部资源勘探问题的一种有效手段，因此对于可控源电磁测深资料的分析与认识变得越来越重要.近年来，可控源资料一维、二维的正反演已趋于成熟(底青云等，2008)，并取得一定效果，但考虑地下地质体的复杂性，采用三维(3D)反演将是一个更加有效的手段(Mackie et al., 1993；Newman et al., 1997；Newman et al., 2000).3D反演计算的基础是3D正演计算.目前常用的正演方法主要有有限差分法(Yee，1966；Newman et al., 1995；Commer et al., 2004；谭捍东等，2004；孙怀凤等，2013)、有限元法(Reddy et al., 1977；王绪本等，2009；徐志峰等，2010)、积分方程法(Ting et al., 1981；Wannamaker et al., 1984；Avdeev，2002a；Avdeev，2002b)，反演方法主要有牛顿法、拟牛顿法(Avdeeva et al., 2006；Avdeev et al., 2009)、非线性共轭梯度法(刘云鹤，2011；翁爱华等，2012)、快速松弛反演法(Smith et al., 1991)、Occam反演法(Siripunvaraporn et al., 2011)等.随着3D正反演方法与技术的不断发展，这些技术已被应用到航空电磁(刘云鹤等，2013)、海洋电磁(He et al., 2008；Newman et al., 2010)等领域，并取得了突破性的进展(胡祖志等，2005；汤井田等，2007).但由于三维模型模拟的计算量非常大，计算时间长，限制了其在三维可控源电磁勘探中的实际应用.因此，如何在反演过程中提高效率和质量具有重要的研究和现实意义.

本文通过数值试验研究频率选择对可控源电磁测深三维反演的影响.期望通过反演讨论不同频率组合对应的观测数据的反演效果，提出满足节约计算量、提高反演速度和不牺牲反演效果的条件下，给出最紧频率组合，为可控源实际应用提出一种解决方案.本文拟采用有限差分法正演，有限内存拟牛顿法反演.文中理论模拟的模型抽象于泥河地区铁矿的地球物理模型(吴明安，2011)，矿区地层由上到下分别为砖桥组、双庙组和浮山组，岩性主要以泥质粉砂岩和凝灰质粉砂岩为主，厚度约为900 m，这些地层电阻率较低，总体平均大约为66 Ω·m.基底侵入岩主要以辉石闪长玢岩为主，电阻率高阻特征约为800 Ω·m.泥河铁矿产在基岩和盖层的接触带上，规模大约为1000 m×1000 m×400 m，属特大型深部隐伏铁矿，其总体电阻率约为30 Ω·m.观测与源平行的Ex实部和虚部，组合不同频率进行反演，研究不同频率组合反演结果.针对该模型的可控源可探测性分析表明，该深部矿产生的异常相对较弱，但反演结果显示出较为可靠的异常.本模型的计算结果同时证明了可控源的勘探深度可以达到1000 m，为可控源深部探矿的应用提供了直接证据.

1 数值理论 1.1 数值模拟

当发射源为有限长电偶源时，假设场随时间变化的因子为e－iωt，三维频率域可控源电磁正演计算基于下面的Helmholtz方程(Nabighian，1987)：

式中 k2=(iωμ0)/ρ，μ0为空气磁导率，i=√ －1 ，Jp为发射导线中的空间域电流.

正演计算主要针对感应场，通常求解如下的微分方程：

式中k2p为背景场介质复波数，k2为含异常后介质复波数，E p为背景场，E s为异常体产生的感应场(Commer et al.， 2008).

背景场Ep计算采用虚界面法(Da，1995；刘云鹤，2011)，数值上采用直接积分法计算(翁爱华等，2003).采用交错网格有限差分法求解感应场Es即(2)式，将三维介质剖分成正六面体单元，即沿x、y、z三个坐标方向分别用若干平行的平面以不同的间距将研究区域划分成若干个小的长方体网格单元.电场取在各网格单元边界的中心，定义Ex、Ey、Ez分别为采样点所在边缘电场的平均值.磁场取在各网格单元表面的中点，Hx、Hy、Hz分别为采样点所在的平面磁场的平均值(如图 1所示).(2)式为矢量微分方程，可以分解为3个坐标矢量的3个微分方程，将离散化后的矩阵方程合并，并且进行对称化处理，得到如下的矩阵方程：

K是复系数大型稀疏对称矩阵，s 包括由背景场以及电阻率异常确定的(2)式右手源矢量和边界条件.求得每个边中点的电场值，再根据法拉第定律微分形式计算磁场分量.
图 1 交错采样网格示意图Fig. 1 Schematic diagram of Finite Difference
1.2 有限内存拟牛顿反演法

设正则化目标函数为

其中d为观测数据和理论模型响应的拟合差，m为模型范数，λ为正则化因子.

设目标函数二次可微，牛顿法是取目标函数在某已知点 x k处的二阶泰勒展开式

其中：k=(x k)；gk= g(x k)为目标函数的一阶导数； p k表示模型搜索方向；Hk为目标函数的二阶导数，即Hessian矩阵，为了寻找极小化目标函数的搜索方向，令，可得

即牛顿法模型更新方式为

因为利用了目标函数二阶导数的信息，能快速达到局部收敛，但直接求解(7)式，需要计算和存储Hessian矩阵 H并求其逆，需要很长的计算时间和较大的存储空间.拟牛顿法改进了这一缺点，它用给定的对称正定矩阵B k作为 Hessian矩阵Hk的近似，即

为了使Bk合理的逼近Hk，需要符合条件：

其中，s k= x k+1－ x k，y k=gk+1－gk.

矩阵Bk是由前一次的近似矩阵Bk-1校正得出.本文采用的BFGS校正公式，即给定一个对称正定初始矩阵Bo∈ R n×n，通过(10)给出的校正公式计算为

本文采用有限内存拟牛顿法就是基于BFGS校正公式，但它给定一个正整数m，只使用最近m次迭代的信息，即使用向量对{ s i，y i}i=k－mk计算Bk+1，丢掉早期迭代的信息，从而极大的减少了内存，提高了计算速度(袁亚湘等，1997).


2 数值结果分析

2.1 装置描述

根据泥河铁矿地质特点，设计一个两层介质的三维地质模型.其中盖层介质电阻率为66 Ω·m，厚900 m，基底电阻率为800 Ω·m.在埋深900 m处有一个1000 m×1000 m×400 m的低阻异常体，电阻率为30 Ω·m(如图 2).

图 2 装置示意图
(a)模型及装置顶视图；(b)模型过y=12000 m的纵断面.Fig. 2 Geometry of a CSEM survey
(a)The top view;(b)The section view from y=12000 m.

发射源使用1000 m的长导线源，距离测区最近测线10000 m，发射频率分别为2-3、2-2、2-1、…、212、213 Hz，共17个频率.测量范围为图 2a中深灰色区域，点距100 m，线距200 m，共861个物理测点，观测参数为Ex的实部和虚部.将理论模拟数据的Ex实部和虚部分别施以1×10-11的噪音后作为观测数据进行反演，网格区域略大于测量范围，在x、y和z三个方向被剖分为40×40×35共56000个网格，x、y方向的网格宽度均为200 m，z方向的网格，第一层厚度为20 m，其下各层厚度按系数1.10递增.反演初始模型采用两层背景模型.反演频率的选择分为两种方式：(1)所有观测频率；(2)隔一个抽取一个频率的组合成两组8个频率.所有观测频率为①号频率组合.两组8个频率组合分别为：②0.125+0.5+2+8+32+128+512+2048 Hz；③0.25+1+4+16+64+256+1024+4096 Hz. 2.2 反演结果分析

图 3~5为3种频率组合的反演结果，为了对比，将低阻异常体用方框画于图中.①号频率三维反演电阻率如图 3所示，从图可见，异常体的形状、规模与理论模型基本一致，对异常体的上界面反映明显，异常体的中心电阻率偏高，基本反映出异常体的特征.图 4为②号频率三维反演电阻率分布图，与①号对比异常较弱，异常体上界面有反映，下界面无反映，不能反映异常体的位置特征.图 5为③号频率三维反演电阻率分布图，与①号结果大致相同，明显反映出异常体的形状与位置，与理论模型接近.比较图 3~5的计算结果发现，反演结果的质量不是由频率的个数决定，而是由部分参加反演的频率决定，对反演结果起决定性作用.

图 3 ①号频率组合三维反演电阻率分布图. 图中黑线代表真实模型边界(下同)
(a)z=-950 m处平面图；(b)y=12000 m处x方向剖面图；(c)x=0 m处y方向剖面图.Fig. 3 3D resisitivity distribution at the end of inversion of No.1. Black line denotes the boundary of the model
(a)The plan view at z=-950 m;(b)The plan view at y=12000 m;(c)The plan view at x=0 m.

图 4 ②号频率组合三维反演电阻率分布图
(a)z=-950 m处平面图；(b)y=12000 m处x方向剖面图；(c)x=0 m处y方向剖面图.Fig. 4 3D resisitivity distribution at the end of inversion of No.2
(a)The plan view at z=-950 m;(b)The plan view at y=12000 m;(c)The plan view at x=0 m.

图 5 ③号频率组合三维反演电阻率分布图
(a)z=-950 m处平面图；(b)y=12000 m处x方向剖面图；(c)x=0 m处y方向剖面图. Fig. 5 3D resisitivity distribution at the end of inversion of No.3
(a)The plan view at z=-950 m;(b)The plan view at y=12000 m;(c)The plan view at x=0 m.

为验证关键频率在反演中的作用，设计④4+64+256 Hz；⑤4+16+256 Hz两组频率组合，反演结果见图 6和图 7.其中图 6为④号频率组合电阻率分布图，从图可见，几乎没有异常，只在异常体上界面有较弱的反映，且电阻率最低值约为750 Ω·m，远大于理论上的30 Ω·m.图 7为⑤号频率组合电阻率分布图，与④号频率对比，异常幅度增大，电阻率异常明显.为直观地比较异常变化，提取两组频率在y=12000 m，z=-958 m时的电阻率沿x剖面变化情况(如图 8所示)，从图可见，⑤号频率组合反演得到的异常电阻率比④号频率组合的异常电阻率小 200 Ω·m.由于两种频率组合只是16 Hz存在差别，说明16 Hz在反演中起到关键性作用.

图 6 ④号频率组合三维反演电阻率分布图
(a)z=-950 m处平面图；(b)y=12000 m处x方向剖面图；(c)x=0 m处y方向剖面图.Fig. 6 3D resisitivity distribution at the end of inversion of No.4
(a)The plan view at z=-950 m;(b)The plan view at y=12000 m;(c)The plan view at x=0 m.

图 7 ⑤号频率组合三维反演电阻率分布图
(a)z=-950 m处平面图；(b)y=12000 m处x方向剖面图；(c)x=0 m处y方向剖面图.Fig. 7 3D resisitivity distribution at the end of inversion of No.5
(a)The plan view at z=-950 m;(b)The plan view at y=12000 m;(c)The plan view at x=0 m.

图 8 y=12000 m，z=-958 m时x剖面电阻率曲线图Fig. 8 The profile map of resistivity at y=12000 m，z=-958 m

为了进一步说明关键频率在反演中的控制作用，用16 Hz+64 Hz的频率组合进行反演，三维反演电阻率分布如图 9所示，图中异常幅度最明显，异常电阻率最低值约400 Ω·m，对异常体的上、下界面均有明显反映.说明16 Hz和64 Hz为泥河模型可控源探测的可能最紧频率组合.

图 9 16 Hz+64 Hz三维反演电阻率分布图
(a)z=-950 m处平面图；(b)y=12000 m处x方向剖面图；(c)x=0 m处y方向剖面图.Fig. 9 3D resisitivity distribution at the end of inversion of 16 Hz+64 Hz
(a)The plan view at z=-950 m;(b)The plan view at y=12000 m;(c)The plan view at x=0 m.
2.3 异常响应分析

为了说明频率选择对反演的影响，了解频率对反演的控制，按如下公式计算异常体中心位置在地表投影(0比2小很多，0)点各频率异常响应相对变化，结果见表 1.

式中Ex0为背景场的电场值，Ex为异常体的电场值.

表 1 不同频率异常相对变化Table 1 Relative error of model having different frequency

由表可知，Ex实部的相对变化都小于20%，16 Hz异常相对变化最大，为-19.550%，其余频率实部的相对变化均小于10%，当频率大于128 Hz时小于1%.观察Ex虚部的相对变化，可发现较实部明显增大，0.5 Hz、16 Hz、32 Hz、64 Hz和128 Hz时Ex虚部的相对变化大于20%，64 Hz异常相对变化最大为618.43%.一般认为异常相对变化达到20%才能认为是可靠异常(蒋邦远，1998)，因此可以认为泥河模型能产生较为明显的异常，并且16 Hz附近的频率是模型的关键响应频率.这种变化说明，在野外工作中，结合地质情况，通过数值模拟，可对异常的可靠性进行分析，然后选择对野外工作有决定性影响的频率. 3 结 论

本文通过对泥河模型反演的研究，讨论了频率选择对三维可控源方法反演的影响.反演观测数据为Ex的实部和虚部，初始模型为两层介质背景模型，采用有限内存拟牛顿法反演，对比反演结果，得出这种条件下泥河模型反演最紧频率组合为16+64 Hz.根据泥河模型反演结果分析，得出以下结论：

(1)观测数据相同，初始模型相同的条件下，反演结果不是频率越多，效果越好.

(2)对于不同观测数据，在能获得合适的初始模型条件下，可能存在最紧频率组合使反演结果与实际情况接近.

(3)最紧频率组合可以通过数值模拟实现.在野外进行三维可控源电磁测深方法时，如果能获得合适的初始模型条件，可以先根据该地区的地质特点设计模型，进行数值模拟，确定最紧频率组合，提高三维可控源电磁测深方法反演的质量和效率.

(4)对于泥河模型这种低阻覆盖模型，可控源电磁测深方法的探测深度可达1000 m左右.这能大体上说，可控源探测深度可以达到1000 m，从而从另外一个侧面回答了可控源能用于大深度探测.

(5)虚部相对于实部，对地下异常电阻率变化更为敏感.因此，在可控源探测中，必须尽可能提高虚部或相位的观测质量.

参考文献
[1] Avdeev D B,Kuvshinov A V,Pankratov O V and Newman G A. 2002a.Three-dimensional induction logging problems. partⅠ. An integral equation solution and model comparisons[J].Geophysics, 67:413-426.
[2] Avdeev D B,Kuvshinov A V and Epova X A. 2002b.Three-dimensional modeling of electromagnetic logs from inclined-horizontal wells[J].Izvestiya.Phys.Solid Earth, 38:975-980.
[3] Avdeev D,Avdeeva A. 2009.3D magnetotelluric inversion using a limited-memory quasi-Newton optimization[J]. Geophysics. 74(3):F45-F57.
[4] Avdeeva A D, Avdeev D B. 2006.QN inversion of large-scale MT data[C].Progress in Electromagnetic Research Symposium.200-213.
[5] Commer M, Newman G. 2004.A Parallel Finite-Difference Approach for 3D Transient Electromagnetic Modeling with Galvanic Sources[J]. Geophysics, 69:1192-1202.
[6] Commer M,Newman G A. 2008.New advances in three-dimensional controlled-source electromagnetic inversion[J]. Geophys. J. Int., 172:513-535.
[7] Da U C.1995.A reformalism for computing frequency- and time- domain EM responses of a buried finite-loop[A].65th Ann .Internat. Mtg.,Soc.Expi.Geophys.Expanded Abstracts,1995:911-814.
[8] Di Q Y, Wang R. 2008. Controlled source audio-frequency magnetotellurics[M]. Science Precss, Beijing.
[9] He Z, Kurt S, Yu G, et al. 2008.On reservoir boundary detection with marine CSEM[J]. Applied Geophysics, 5(3): 181-188.
[10] Hu Z Z,Hu X Y. 2005.Review of three dimensional magnetotelluric inversion methods[J].Progress in Geophysics(in Chinese), 20(1): 214.
[11] Jiang B Y. 1998.Applied Near Zone Magnetic Source Transient Electromagnetic Exploration(in Chinese)[M].Beijing:Geological Publishing House.
[12] Liu Y H. 2011.Research on 3-D Controlled Source Electromagnetic Method Inversion Using Nonlinear Conjugate Gradients[Ph.D.thesis](in Chinese).Changchun:Jilin University.
[13] Liu Y H,Yin C C. 2013.3D inversion for frequency-domain HEM data[J]. Chinese J. Geophys. (in Chinese). 56(12): 4278-4287,doi:10.6038/cjg20131230.
[14] Liu Y H,Liu G X,Weng A H,et al. 2011.Calculation of low-frequency electromagnetic Green's functions using two-level approximate discrete complex image method[J]. Chinese J. Geophys. (in Chinese). 54(4):1114-1121.
[15] Mackie R L, Madden T R. 1993.Three-dimensional magnetotelluric inversion using conjugate gradients [J]. Geophys. J. Int., 115(1): 215-229.
[16] Nabighian M N. 1987.Electromagnetic Methods in Applied Geophysics[M].Vol.1,SEG,Tulsa,Oklahoma.
[17] Newman G A. Alumbaugh D L. 1997.Three-dimensional massively parallel electromagnetic inversion-I,Theory[J].Geophys.J.Int., 128:345-354.
[18] Newman G A, Alumbaugh D L. 2000.Three-dimensional magnetotelluric inversion using non-linear conjugate gradients[J]. Geophys. J. Int., 140: 410-424.
[19] Newman G A, Alumbaugh D L. 1995.Frequency-domain modelling of airborne electromagnetic responses using
[20] staggered finite differences[J]. Geophysical Prospecting, 43(8): 1021-1042.
[21] Newman G A, Commer M, and Carazzone J J. 2010.Imaging CSEM data in the presence of electrical anisotropy [J]. Geophysics, 75(2): F51-F61.
[22] Reddy I K, Phillips R J, Rankin D. 1977.Three-dimensional modelling in magnetotelluric and magnetic variational sounding[J]. Geophysical Journal International, 51: 313-325.
[23] Siripunvaraporn W, Sarakorn W. 2011.An efficient data space conjugate gradient Occam's method for three-dimensional magnetotelluric inversion[J].Geophysical Journal International, 186(2): 567-579.
[24] Smith J T, Booker J R. 1991.Rapid inversion of two-and three-dimensional magnetotelluric data[J]. Journal of Geophysical Research: Solid Earth (1978-2012), 96(B3): 3905-3922.
[25] Sun H F,Li X,Li S C. 2013.Three-dimensional FDTD modeling of TEM excited by a loop source considering ramp time[J]. Chinese J. Geophys. (in Chinese) . 56(3):1049-1064, doi:10.6038/cjg20130333.
[26] Tang J T,He J S. 2005.Controlled source audio-frequency Magnetotellurics[M].Central South University of Technology Press, Changsha.
[27] Tang J T,Ren Z Y,Hua X R. 2007.The forward modeling and inversion in geophysical electromagnetic field[J]. Progress in Geophysics(in Chinese), 22(4):1181-1194.
[28] Tan H D,Yu Q F,Wei W B. 2004.Magnetotelluric three-dimensional modeling using the staggered-grid finite difference method[J]. Chinese J. Geophys. (in Chinese),46(5): 705-711.
[29] Ting S C, Hohmann G W. 1981.Integral equation modeling of three-dimensional magnetotelluric response[J]. Geophysics, 46(2): 182-197.
[30] Wang X B,Liu Y,Gao Y C. 2009.Finite Element Modeling of 3D MT with Tetrahedron[A].The Chinses Geophysics.beijing.
[31] Wannamaker P E, Hohmann, G W, and SanFilipo, W A. 1984.Electromagnetic modeling of three-dimensional bodies in layered earths using integral equations[J]. Geophysics, 49, 60-74.
[32] Weng A H, Liu Y H, Jia D Y, Liao X D, Yin C C. 2012.Three-dimensional controlled source electroma- gnetic inversion using non-linear conjugate gradients[J]. Chinese J. Geophys. (In Chinese), 55(10):3506-3515,doi:10.6038/j.issn.0001-5733.2012.10.034.
[33] Weng A H, Wang X Q. 2003.Utilizingdirect integration to enhance calculation accuracy of 1d electromagnetic response for current dipole source[J]. Northwestern Seismological Journal, 25(3):193-197.
[34] Wu M A. 2011.Discovery of the Nihe Iron Deposit in Lujiang, Anhui and its Exploration Significance[J].Acta geologica sinica. 85(5):802-809.
[35] Xu Z F,Wu X P. 2010.Controlled source electromagnetic 3-D modeling in frequency domain by finite element method.[J] Chinese J. Geophys. (in Chinese) . 53(8): 1931-1939,doi:10.3969/j.issn.0001-5733.2010.08.019.
[36] Yang Y X,Sun W Y. 1997.Optimization theory and methods[M].Science Press, Beijing.
[37] Yee K S. 1966. Numerical solution of initial boundary value problems involving Maxwell's equations[J]. IEEE Trans. Antennas Propag, 14(3): 302-307.
[38] 底青云,王若. 2008.可控源音频大地电磁数据正反演及方法应用[M].北京:科学出版社.
[39] 胡祖志,胡祥云. 2005.大地电磁三维反演方法综述[J].地球物理学进展, 20(1): 214.
[40] 蒋邦远. 1998.实用近区磁源瞬变电磁法勘探[M].北京:地质出版社.
[41] 刘云鹤. 2011.三维可控源电磁法非线性共轭梯度反演研究[博士论文]长春:吉林大学.
[42] 刘云鹤, 殷长春.2013.三维频率域航空电磁反演研究[J]. 地球物理学报, 56(12): 4278-4287,doi:10.6038/cjg20131230.
[43] 刘云鹤,刘国兴,翁爱华,等.2011.基于二级近似离散复镜像法的低频电磁格林函数计算[J].地球物理学报,54(4):1114-1121.
[44] 孙怀凤, 李貅, 李术才,等. 2013.考虑关断时间的回线源激发 TEM 三维时域有限差分正演[J]. 地球物理学报, 56(3):1049-1064,doi:10.6038/cjg20130333.
[45] 谭捍东,余钦范,魏文博. 2004.大地电磁法三维交错采样有限差分数值模拟[J].地球物理学报,46(5): 705-711.
[46] 汤井田,何继善. 2005.可控源音频大地电磁法及其应用[M]. 长沙:中南大学出版社.
[47] 汤井田,任政勇,化希瑞. 2007.地球物理学中的电磁场正演与反演[J].地球物理学进展, 22(4):1181-1194.
[48] 王绪本,刘云,高永才.2009.大地电磁测深三维有限单元四面体剖分法数值模拟[A].中国地球物理学会第二十五届年会论文集.北京.
[49] 翁爱华,刘云鹤, 贾定宇, 廖祥东, 殷长春. 2012.地面可控源频率测深三维非线性共轭梯度反演[J].地球物理学报, 55(10):3506-3515,doi:10.6038/j.issn.0001-5733.2012.10.034.
[50] 翁爱华,王雪秋.2003.利用数值积分提高一维模型电偶源电磁测深响应计算精度[J]. 西北地震学报, 25(3):193-197.
[51] 吴明安. 2011.安徽庐江泥河铁矿的发现及意义[J].地质学报, 85(5):802-809.
[52] 徐志锋, 吴小平.2010.可控源电磁三维频率域有限元模拟[J]. 地球物理学报, 53(8): 1931-1939,doi:10.3969/j.issn.0001-5733.2010.08.019.
[53] 袁亚湘,孙文瑜. 1997.最优化理论与方法[M].北京:科学出版社.