地球物理学进展  2016, Vol. 31 Issue (4): 1513-1518   PDF    
基于多重网格的Helmholtz方程解算及简单地质模型的大地电磁正演模拟
杨振威1,2, 徐招峰1,2, 侯欣欣3, 赵秋芳1,2     
1. 河南理工大学资环学院, 焦作 454000
2. 中原经济区煤层(页岩)气河南省协同创新中心, 焦作 454000
3. 河南理工大学, 焦作 454000
摘要: 为了提高大地电磁正演模拟的计算效率,开展了基于多重网格算法的Helmholtz方程解算研究;对比Gauss-Seidel、Richardson迭代算法的收敛性,在达到同样的误差限值下,多重网格算法仅需十余次套迭代即可满足误差精度要求,而Gauss-Seidel和Richardson算法分别需要近500次和5000次迭代;多重网格算法的收敛速度极快,多重网格算法的收敛速度比一般数值算法快近3个数量级;由于其是一种基于套迭代技术,该方法完成一次循环的迭代效率明显比一般数值算法低;多重网格算法极快的收敛速度为将其用于提高大地电磁正演模拟效率提供了可能.
关键词多重网格     Helmholtz方程     收敛性     计算效率    
Solution of Helmholtz equation based on multigrid and magnetotelluric modeling based on simple model
YANG Zhen-wei1,2 , XU Zhao-feng1,2 , HOU Xin-xin3 , ZHAO Qiu-fang1,2     
1. School of Resources and Environment, Henan Polytechnic University, Jiaozuo 454000, China
2. Collaborative Innovation Center of Coalbed Methane and Shale Gas for Central Plains Economic Region, Henan Province, Jiaozuo 454000, China
3. Henan Polytechnic University, Jiaozuo 454000, China
Abstract: In order to improve computational efficiency for magnetotelluric forward simulation, solving helmholtz equation based on multigrid is studied. Comparing to iterative algorithm such as Gauss-Seidel、Richardson, reached the same error limit, multigrid algorithm need more than ten times, and Gauss-Seidel、Richardson algorithm need 500 times and 600 times respectively. The convergence rate of multigrid is very fast—an order of magnitude upper than the other algorithm. Because the multigrid is a nested iteration, the iteration efficiency is not high than the other numerical algorithm. The fast convergence rate of makes it possible to improve magnetotelluric forward simulation by multigrid.
Key words: multigrid     helmholtz equation     convergence     computational efficiency    
0 引 言

大地电磁法三维反演是计算地球物理研究的前沿和热点,而精度高且快速的三维正演计算是大地电磁三维反演走向实际应用的重要基础(谭捍东等,2003葛永斌等,2006李勇等,2012).目前,三维正演计算主要采用有限差分、有限元等数值模拟方法求解Helmholtz方程(Tan et al.,2006汤井田等,2007翁爱华等,2012),但最终都归结到求解大型线性方程组Au=f.如何快速、准确的解算大型线性方程组是实现正演计算的关键,一般的数值迭代算法(Gauss-Seidel、SOR等)解算线性方程组存在的一个重要问题是:计算量大、计算效率低(随着迭代次数的增加,收敛速度减慢),因此,开展可有效提高大型线性问题计算效率的数值算法的研究,具有重要的意义(刘卓等,1998田小波等,2004Pereira et al.,2006).

近年来,国内外学者开展了基于多重网格法求解椭圆形偏微分方程的的研究(Kwak,1999Aruliah and Ascher,2002葛永斌等,2006).欧阳洁(1994)开展了自适应控制层间转换的多重网格算法研究,结果表明该方法比高斯—赛德尔迭代(Gauss-Seidel)方法和逐次超松弛法(SOR)的计算量少的多;柳建新等(2011)为了准确预测和分析多重网格法用于大地电磁法正演计算的收敛效果,对多重网格法的收敛性进行了广义傅里叶谱分析;鲁晶津等(2009)开展了多重网格法求解三维泊松方程的数值模拟研究,结果表明:多重网格法的计算速度明显优于ICCG和Gauss方法;

本文在简要介绍多重网格算法的基础上,对比研究了多重网格、高斯—赛德尔、理查德森算法在求解Helmholtz方程过程中的迭代次数、迭代误差和计算时间,认为多重网格算法的收敛速度比一般迭代算法快近3个数量级,但单次循环所需时间较长.

1 多重网格算法

多重网格(MultiGrid method)算法是由苏联数学家Fedorenko于20世纪60年代首先提出的,其思想源于一种加速收敛技术——套迭代,即利用粗网格校正和迭代光滑高频误差分量,离散化为更精细的网格时,仍能保持较快收敛速度,可有效解决复杂数值计算中收敛速度慢的问题(Gupta et al.,1997唐学军等,2000;Chen et al.,2011).

采用数值算法求解偏微分方程的误差分量可以分为两类,一类是频率较高,振荡快的高频误差分量;另一类是频率变化较慢的低频误差分量.一般的迭代方法可以快速的消减振荡误差分量,但对低频误差分量,消除效率很低.高频误差分量和低频误差分量是相对而言的,在细网格上的低频误差分量,在粗网格上可被视为为高频分量.多重网格方法可以快速计算迭代求解由偏微分方程组离散后的代数方程组,其基本原理是一定尺度的网格容易消除波长与网格步长相对应的误差分量,采用不同尺度的网格消除不同波长的误差分量,即可实现快速迭代收敛,该算法首先在细网格上进行几次光滑迭代,然后将误差分量转移到较粗的网格上,消除与该层网格上相对应的较易消除的那些误差分量,这样逐层进行下去直到消除不同误差分量,求解粗网格方程后,再将粗网格上的解向量逐层投射到细网格上(Hemker et al.,1983柴玉珍等,1999Plaks et al.,2000),重复以上迭代过程,直至数值方程收敛(Zaslavsky and Schlick,1998;Briggs et al.,2011).一种简单的多重网格算法——二重网格V循环的迭代过程如图 1所示.

图 1 二重网格算法的V循环(Ωh代表细网格层,ΩH代表粗网格层) Figure 1 V-circle of double grid(Ωh represents fine grid,ΩH represents coarse grid)
2 求解Au=f的多重网格迭代步骤(以二重网格为例)

当多重网格算法的网格迭代层数为2时,多重网格算法即为二重网格法,二重网格的迭代步骤为:

① 前光滑化:即在细网格层上进行光滑化处理,进而消除高频误差分量;基于一般迭代算法(Richardson、Jacobi、Gauss-Seidel等)在细网格层上求解离散化的线性方程组Ahuh=fh,式中Ah为细网格层上的系数矩阵,uh为方程组在细网格层上的迭代解,f为方程组右端项;

② 计算残差:rh=fh-Ahuh,式中rh为细网格层上经过两次迭代计算后的残差;

③ 残量限制(图 2):r2h=Ih2h rh,将细网格层上计算得到的迭代误差通过限制算子I投射至粗网格层上,式中r2h为由细网格层投射到粗网格层上的残差;

图 2 误差向量由细网格限制至粗网格示意图 Figure 2 Schematic diagram on interpolation of error vector from fine grid to coarse grid

④ 求解粗网格方程:e2h=(A2h)-1r2h,在粗网格层上计算解的修正量e2h

⑤ 修正量插值回代:eh= P2hh e2h,通过插值算子P,将修正量从粗网格层转移至细网格层上,如图 3所示;

图 3 解的修正量由粗网格插值至细网格示意图 Figure 3 Schematic diagram on interpolation of error vector from coarse grid to fine grid

⑥ 计算细网格层上的修正解:uh=uh+eh,并多次重复上述过程,直至结果收敛.

其迭代过程如图 4所示.

图 4 二重网格法迭代流程示意图 Figure 4 Schematic diagram of double grid iteration
3 一维Helmholtz方程的多重网格算例

一维Helmholtz方程的Dirichlet边值问题为

(1)

首先,对控制方程在细网格层上进行离散化,在区间[-1,1]分成64等份,hi=1/32,在区间内节点xi,基于微分公式:

(2)

原控制方程可变形为

即:

(3)

在进行二重网格迭代之前,先基于Richardson算法进行2次光滑化,公式为

(4)

式中,ω=1/2是松弛因子,A是系数矩阵,h=1/64为细网格步长,F为细网格下方程右端项.二重网格法的解算步骤如下:

(1) 前光滑化,取迭代初始值u0=0,则

(5)
(6)

(2) 经过两次光滑化处理后,进入二重网格算法的第二步:残量(亏量)计算,公式为

(7)

(3) 残量限制:将残量从细网格限制到粗网格上,需要设定限制算子I2hh,考虑到 基于I2hh的映射:ΩhΩ2h,若定义vfvc分别是u在粗网格层和细网格层上的近似解,则,

I2hh×vh=v2h,且vhi2=$\frac{1}{4}$(v2i-1h+22ivh+v2i+1h),设I2hh为线性算子, RN-1⇒RN/2-1,对于N=4,则限制算子I=14[1 2 1 0 0;0 1 2 1 0],得到:

(8)

(4) 解算粗网格方程:

(9)

(5) 修正u的值:ufnew=uf-pec,式中,p为插值算子:p=$\frac{1}{2}$[1 2 1 0 0;0 1 2 1 0]T,p的取值方法与I的取值方法类似.

重复上述迭代过程,直至方程的解达到预设的误差限,此时,ufnew即是二重网格迭代的结果.

4 多种迭代算法对比 4.1 收敛速度对比

为了定量研究多重网格算法在求解Helmholtz方程的收敛特性,作者共采用了三种数值计算方法进行了对比计算:高斯-赛德尔(Gauss-Seidel)、理查德森(Richardson)和二重网格(Double Grid)算法,迭代次数分别设定为5次、50次、500次和5000次,迭代计算得到的Helmholtz方程曲线如图 5所示,图 5abcd分别显示了不同迭代次数的迭代计算曲线图,图中以+表示理查德森算法的Helmholtz方程迭代计算曲线,○表示高斯-赛德尔算法的Helmholtz方程解迭代计算曲线,☆表示二重网格算法的Helmholtz方程迭代计算曲线;从图中不难看出,二重网格算法的收敛速度最快,只需5次即可达到误差精度要求;达到同样的误差限值,Richardson算法的收敛速度最慢,需要5000次(选取的是最佳松弛因子),而Gauss-Seidel的收敛速度快于Richardson算法,但相比较于Double Grid算法,收敛速度还是慢近三个数量级.图 5直观的展示了多重网格算法在收敛速度上的优越性.

图 5 基于Gauss-Seidel、Richardson、Double-Grid算法不同迭代次数Helmholtz方程曲线图 Figure 5 The curve of numerical solution of helmholtz eqution by method of Gauss-Seidel、Richardson and Double Grid with different iteration number
4.2 迭代误差分析

迭代误差值设定为:Error=max(uiter-uex),其中uiter是迭代计算过程中的迭代解,uex为Helmholtz方程的精确解,当迭代误差小于误差限ep=1×10-6时,即Error<ep时,迭代终止,图 6a是基于Gauss-Seidel、Richardson算法求解Helmholtz方程的误差值随迭代次数的变化去曲线图,由图中可以看出经过不到500次迭代,误差值较快的从0.15左右减小至0.02,迭代次数为2000次时,误差值已基本达到极限,之后迭代次数继续增加至5000时,误差值基本不发生变化;

图 6 Gauss-Seidel、Richardson算法(a)和多重网格算法(b)迭代次数与迭代误差对比图 Figure 6 Contrast diagram of iterative number and iterative error with Gauss-Seidel、Richardson algorithm and Multigrid algorithm

图 6b为基于多重网格算法求解Helmholtz方程的误差值随迭代次数的变化曲线图,不难发现,初次迭代的误差值即减小至0.2×10-3至,经过不到500次迭代,迭代解与精确解的差值即达到误差限值,该图从另外一个角度呈现了多重网格的迭代收敛效率.

4.3 计算时间分析

表 1是不同数值算法Gauss-Seidel、Richardson和Multigrid求解Helmholtz方程在微型机(CPU2.20 GHz,内存2 GB)上的计算时间,相同的迭代次数(5次、50次、500次、5000次),三种数值算法迭代计算时间由小到大依次为:Richardson、Gauss-Seidel和Multigrid,由此可见,多重网格算法基于套迭代技术,虽迭代收敛速度快,但一次循环迭代的计算时间较长(迭代效率较低).

表 1 三种迭代算法不同迭代次数与所需时间对比表 Table 1 Contrast table of different iterative number and time needed with three iterative algorithm

表 1不难看出,随着迭代次数的增加,三种数值算法求解Helmholtz方程的计算时间存在一定差异,显然,在一定迭代次数的条件下,相较于Gauss-Seidel、Richardson数值算法,Multigrid算法的计算时间是最长的.

4.4 简单模型的大地电磁正演模拟

为了研究多重网格算法在大地电磁一维正演研究中的适用性,初步研究了大地电磁一维均匀半空间模型的大地电磁正演响应特征,设计模型电阻率为100 Ω·m,向下极限深度为1.0×105 m,模型下边界E=0,上边界H=1,考虑到Guass-Seidel迭代方法在迭代计算过程中,存在收敛性不稳定的缺陷,本次大地电磁正演计算中,采用了收敛性更加稳定的稳定双共轭梯度法(Bicgstab)算法,二重网格算法的单层网格迭代算法也采用Bicgstab算法.图 7图 8分别为基于Bicgstab和二重网格算法的大地电磁正演响应曲线,从图中不难看出,相同迭代次数条件下,二重网格算法得到的模型响应视电阻率值更接近真实电阻率.

图 7 (a)均匀半空间条件下基于Bicgstab算法的大地电磁一维正演模拟;(b)均匀半空间条件下基于Double grid算法的大地电磁一维正演模拟 Figure 7 (a) 1D magnetotelluric forward based on Bicgstab(uniformity half space model);(b)1D magnetotelluric forward based on Bicgstab(uniformity half space model)
图 8 均匀半空间模型的大地电磁一维正演响应 Figure 8 1D Forward modeling of magnetotelluric of uniformity half space model

图 7a存在的一个问题是在低频段,由于迭代次数偏少,计算误差较大,模型响应电阻率值几乎等于0,在高频段,模型响应视电阻率值在100 Ω·m附近快速上扬;图 7b中,在同样的迭代次数下(100次),模型计算误差相对较小,显示出二重网格算法收敛速度相较单层网格快.为了从另外一个侧面验证此说法的正确性,研究了直接解法求解模型的有限差分方程.图 8是网格单元个数为5000的直接解法求解大地电磁一维响应曲线,由图可知,直接解法得到的响应曲线在低频段f=10-3~10-2 Hz,视电阻率剧烈震荡,f=10--2~102 Hz处视电阻率曲线平直,视电阻率值稳定在100 Ω·m附近,当f>102 Hz时,视电阻率值逐渐增大,偏离了真实模型电阻率值.

5 结论和认识

5.1   本文以二重网格为例,基于一维Helmholtz方程的Dirichlet边值问题研究了多重网格算法相较于一般迭代算法的收敛速度以及计算时间.

(1) 多重网格算法的迭代收敛速度非常快,仅需3~5次迭代,即可收敛至精确解,而Gauss-Seidel和Richardson算法需要近5000次,相较于一般松弛迭代算法(Gauss-Seidel、Richardson),多重网格算法的收敛速度快近3个数量级(迭代次数);

(2) 多重网格算法的完成一次套迭代的计算时间较长,相同迭代次数条件下较一般松弛迭代算法(Gauss-Seidel、Richardson)所需时间更长;

(3) 初步研究了均匀半空间模型条件下基于多重网格法大地电磁正演模拟响应,认为多重网格法较单层网格算法在收敛速度方面有一定的优势.

5.2   上述研究证明:一方面,多重网格算法在收敛速度方面的优越性为将其应用于大地电磁正演模拟中提供了可能;另一方面,与一般大地电磁法正演模拟不同的是:上述迭代算法研究Helmholtz方程是在[0,1]区间网格均匀剖分,且不存在如何设定边界条件等问题,因此,将多重网格算法应用与大地电磁正演模拟研究,尚需进一步深入研究复杂模型的大地电磁正演模拟以及将自适应多重网格算法应用于大地电磁正演模拟的可能性.

参考文献
[1] Aruliah D A, Ascher U M.2002. Multigrid preconditioning for Krylov methods for time-harmonic maxwell's equations in three dimensions[J]. SIAM Journal of Scientific Computing, 24 (2) : 702–718. DOI:10.1137/S1064827501387358
[2] Briggs W L, Henson V E, McCormick S F.2011. A Multigrid Tutorial[M]. 2nd ed. Beijing: Tsinghua Publishing House .
[3] Chai Y Z, Qin X Y, Wang S L.1999. The application of the twogrid method on a simple example[J]. Journal of Taiyuan University of Technology (in Chinese), 30 (6) : 666–668.
[4] Chen X B, Zhao G Z, Tang J, et al.2005. An adaptive regularized inversion algorithm for magnetotelluric data[J]. Chinese Journal of Geophysics (in Chinese), 48 (4) : 937–946. DOI:10.3321/j.issn:0001-5733.2005.04.029
[5] Chen Z Y, Cheng D S, Feng W, et al.2011. A multigrid-based preconditioned Krylov subspace method for the Helmholtz equation with PML[J]. Journal of Mathematical Analysis and Applications, 383 (2) : 522–540. DOI:10.1016/j.jmaa.2011.05.054
[6] Ge Y B, Tian Z F, Ma H L.2006. A high accuracy multigrid method for the three-dimensional poisson equation[J]. Mathematica Applicata (in Chinese), 19 (2) : 313–318.
[7] Gupta M M, Kouatchou J, Zhang J.1997. Comparison of second and fourth-order discretizations for multigrid poisson solvers[J]. Journal of Computational Physics, 132 (2) : 226–232. DOI:10.1006/jcph.1996.5466
[8] Hemker P W, Kettler R, Wesseling P, et al.1983. Multigrid methods: Development of fast solvers[J]. Applied Mathematics and Computation, 13 (3-4) : 311–326. DOI:10.1016/0096-3003(83)90018-8
[9] Kwak D Y, Kwon H J, Lee S.1999. Multigrid algorithm for cell centered finite difference on triangular meshes[J]. Applied Mathematics and Computation, 105 (1) : 77–85. DOI:10.1016/S0096-3003(98)10094-2
[10] Li Y, Wu X P, Lin P R.2012. A study on three-dimension forward modeling and anomaly features of magnetotelluric sounding[J]. Progress in Geophysics (in Chinese), 27 (6) : 2452–2463. DOI:10.6038/j.issn.1004-2903.2012.06.020
[11] Liu J X, Guo R W, Tong X Z, et al.2011. General fourier analysis of multi-grid in magnetotelluric modeling[J]. Journal of Jilin University (Earth Science Edition)(in Chinese), 41 (5) : 1587–1595.
[12] Liu Z, Zeng Q C, Zhang L B.1998. The application of multigrid method in the adaptive mesh of acceleration generation[J]. Scientia Atmospherica Sinica (in Chinese), 22 (2) : 191–198.
[13] Lu J J, Wu X P, Spitzer K.2009. Multigrid method for 3D modeling of poisson equation[J]. Progress in Geophysics (in Chinese), 24 (1) : 154–158.
[14] Ouyang J.1994. The comparision between adaptive multiple grid method and successive overrelaxation method[J]. Communication on Applied Mathematics (in Chinese), 8 (1) : 27–33.
[15] Pereira F H, Verardi S L L, Nabeta S I.2006. A fast algebraic multigrid preconditioned conjugate gradient solver[J]. Applied Mathematics and Computation, 179 (1) : 344–351. DOI:10.1016/j.amc.2005.11.115
[16] Plaks A, Tsukerman I, Painchaud S, et al.2000. Multigrid methods for open boundary problems in geophysics[J]. IEEE Transactions on Magnetics, 36 (4) : 633–638. DOI:10.1109/20.877530
[17] Shahbazi K, Mavriplis D J, Burgess N K.2009. Multigrid algorithms for high-order discontinuous Galerkin discretizations of the compressible Navier-Stokes equations[J]. Journal of Computational Physics, 228 (21) : 7917–7940. DOI:10.1016/j.jcp.2009.07.013
[18] Tan H D, Tong T, Lin C H.2006. The parallel 3D magnetotelluric forward modeling algorithm[J]. Applied Geophysics, 3 (4) : 197–202. DOI:10.1007/s11770-006-4001-5
[19] Tan H D, Yu Q F, Booker J, et al.2003. Magnetotelluric three-dimensional modeling using the staggered-grid finite difference method[J]. Chinese Journal of Geophysics (in Chinese), 46 (5) : 705–711. DOI:10.3321/j.issn:0001-5733.2003.05.019
[20] 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. DOI:10.3969/j.issn.1004-2903.2007.04.025
[21] Tang J T, Wang F Y, Ren Z Y, et al.2010. 3-D direct current resistivity forward modeling by adaptive multigrid finite element method[J]. Journal of Central South University of Technology, 17 (3) : 587–592. DOI:10.1007/s11771-010-0527-z
[22] Tang X J, Wang J H, Kou X J, et al.2000. Multi-grid solution to adaptive nonlinear FEM and its application[J]. Rock and Soil Mechanics (in Chinese), 21 (2) : 102–107.
[23] Tian X B, Wu Q J, Zeng R S.2004. Multi-grid algorithm for numerical modeling of elastic wave field using finite-difference method[J]. Chinese Journal of Geophysics (in Chinese), 47 (1) : 81–87. DOI:10.3321/j.issn:0001-5733.2004.01.012
[24] Wang R, Wang M Y, Di Q Y.2006. Electromagnetic modeling due to line source in frequency domain using finite element method[J]. Chinese Journal of Geophysics (in Chinese), 49 (6) : 1858–1866. DOI:10.3321/j.issn:0001-5733.2006.06.035
[25] Weng A H, Liu Y H, Jia D Y, et al.2012. Three-dimensional controlled source electromagnetic inversion using non-linear conjugate gradients[J]. Chinese Journal of Geophysics (in Chinese), 55 (10) : 3506–3515. DOI:10.6038/j.issn.0001-5733.2012.10.034
[26] Zaslavsky L Y, Schlick T.1998. An adaptive multigrid technique for evaluating long-range forces in biomolecular simulations[J]. Applied Mathematics and Computation, 97 (2-3) : 237–250. DOI:10.1016/S0096-3003(97)10146-1
[27] 柴玉珍, 秦效英, 王淑丽.1999. 二重网格方法的应用简例[J]. 太原理工大学学报, 30 (6) : 666–668.
[28] 陈小斌, 赵国泽, 汤吉, 等.2005. 大地电磁自适应正则化反演算法[J]. 地球物理学报, 48 (4) : 937–946. DOI:10.3321/j.issn:0001-5733.2005.04.029
[29] 葛永斌, 田振夫, 马红磊.2006. 三维泊松方程的高精度多重网格解法[J]. 应用数学, 19 (2) : 313–318.
[30] 李勇, 吴小平, 林品荣.2012. 大地电磁测深三维正演模拟及异常特征研究[J]. 地球物理学进展, 27 (6) : 2452–2463. DOI:10.6038/j.issn.1004-2903.2012.06.020
[31] 刘卓, 曾庆存, 张林波.1998. 多重网格法在加速自适应网格生成中的应用[J]. 大气科学, 22 (2) : 191–198.
[32] 柳建新, 郭荣文, 童孝忠, 等.2011. 大地电磁法正演中多重网格法求解的广义傅里叶谱分析[J]. 吉林大学学报(地球科学版), 41 (5) : 1587–1595.
[33] 鲁晶津, 吴小平, SpitzerK.2009. 三维泊松方程数值模拟的多重网格方法[J]. 地球物理学进展, 24 (1) : 154–158.
[34] 欧阳洁.1994. 自适应多重网格法与超松弛法的比较[J]. 应用数学与计算数学学报, 8 (1) : 27–33.
[35] 谭捍东, 余钦范, BookerJ, 等.2003. 大地电磁法三维交错采样有限差分数值模拟[J]. 地球物理学报, 46 (5) : 705–711. DOI:10.3321/j.issn:0001-5733.2003.05.019
[36] 汤井田, 任政勇, 化希瑞.2007. 地球物理学中的电磁场正演与反演[J]. 地球物理学进展, 22 (4) : 1181–1194. DOI:10.3969/j.issn.1004-2903.2007.04.025
[37] 唐学军, 王建华, 寇新建, 等.2000. 自适应非线性有限元的多重网格法求解及应用[J]. 岩土力学, 21 (2) : 102–107.
[38] 田小波, 吴庆举, 曾融生.2004. 弹性波场数值模拟的隐式差分多重网格算法[J]. 地球物理学报, 47 (1) : 81–87. DOI:10.3321/j.issn:0001-5733.2004.01.012
[39] 王若, 王妙月, 底青云.2006. 频率域线源大地电磁法有限元正演模拟[J]. 地球物理学报, 49 (6) : 1858–1866. DOI:10.3321/j.issn:0001-5733.2006.06.035
[40] 翁爱华, 刘云鹤, 贾定宇, 等.2012. 地面可控源频率测深三维非线性共轭梯度反演[J]. 地球物理学报, 55 (10) : 3506–3515. DOI:10.6038/j.issn.0001-5733.2012.10.034