基于高斯牛顿法的二维直流电阻率法的 快速反演

赵东东，张钱江，戴世坤，陈龙伟，李 昆

(中南大学 地球科学与信息物理学院，三维电磁反演与成像实验室，长沙 410083)

直流电阻率法在水文地质、工程和资源勘探中得到了广泛的应用。随着仪器设备精度的提高和数据采集的自动化、智能化发展，直流电阻率法将获得更好的应用效果。由于地下构造是三维的，采用三维反演技术才能获得更可靠的反演结果，但考虑到实现三维数据采集和三维反演需要耗费大量的人工成本和计算资源，因此，在实际勘探中仍以二维数据采集和二维反演为主。

在二维直流电阻率法正反演中，国内外学者做了大量的研究。在正演数值模拟研究中，周熙襄等[1]将混合边界条件应用到点源二维直流电阻率法中，并采用交叉对称网格的剖分方法大大提升了有限元数值模拟的精度和数值结果的对称性。罗延钟等[2]提出一种求解二维线源起伏地形条件下中梯装置电阻率法正演问题的有限差分方法，得到比较满意的模拟精度和计算速度。MUNDRY等[3]采用有限差分法实现了点源二维直流电阻率法的求解。罗延钟等[4]讨论了二维构造上电阻率法有限单元模拟的边值条件、反傅氏变换算法及波数选取等问题，并提出了相应的算法。徐世 浙[5]详细介绍了有限元法在二维直流电阻率法中的求解过程。底青云等[6]从电性介质线源有限元模拟出发，深入探讨了改进复杂电性结构理论模型有限元正演资料仿真性和有效性的方法。阮百尧等[7]实现了电导率分块线性变化二维地电断面电阻率测深有限元数值模拟。柳建新等[8]提出了最优化离散波数的优化算法，使得点源二维直流电阻率法在计算量和计算精度方面均有很大的改善。汤井田等[9-10]介绍了基于非结构化网格的点源二维直流电阻率正演数值模拟算法，并提出一套具有更高精度和更大适用范围的傅里叶反变换离散波数。张东良等[11]采用有限差分方法模拟了二维起伏地表线源直流电阻率问题。PENZ等[12]提出了一种解决源奇异性问题的新方法，并在二维直流电阻率法中取得了良好的效果。麻昌英等[13]针对复杂地形下高密度激电法2.5维有限单元法数值模拟进行了研究，分析得出了山谷地形比山脊地形对极化率观测的影响大的结果。彭艳华等[14]利用边界校正法提高了源附近数值模拟精度，减小了边界条件对数值模拟结果的影响。肖晓等[15]和胡宏伶等[16]提出一种新的2.5D直流电阻率有限元-无限元耦合的方法和局部加密等级网格的方式，有效地提高了数值模拟精度和速度。

在二维直流电阻率法反演成像中，同样也有不少学者做了大量的研究工作。BARKER[17]利用ZOHDY[18]提出的用近似反演结果作为二维正演的初始模型，用实测测深曲线与预测测深曲线的比值来修改模型。LOKE等[19-20]基于最平滑约束进行了二维直流电阻率法反演，很大程度上减少了计算时间和存储空间。阮百尧等[21]采用一维直流电阻率测深的直接反演公式，结合二维正演和渐进迭代反演，推出一种精度相对较高的快速反演方法。徐海浪等[22]利用BP神经网络优化方法，实现了电阻率二维非线性反演，得到了比较理想的反演效果。刘海飞等[23]通过最小二乘法独立反演不同时刻采集的直流电阻率法实测数据，来重构间歇影像。程勃等24-25]基于遗传算法和统计学实现了电阻率测深二维反演，大大降低了计算量，并且取得了理想的效果。韩波等[26]结合全变差正则化方法和吉洪诺夫正则化方法的优点对二维分块常值电阻率反演进行了研究。戴前伟等[27]实现了基于混沌振荡PSO-BP算法的二维电阻率层析成像非线性反演。冯德山等[28-29]从超高密度电法的基本理论与数据采集方式入手，实现了全四极装置情形下的有限单元法正演与广义最小二乘反演，能够有效地模拟地形和复杂地质体。柳建新等[30]研究了直流电阻率测深二维自适应正则化反演，反演结果能比较准确地反映地下模型的真实电性结构。

虽然二维直流电法研究成果很多，大多集中在研究点源二维问题，对于线源二维问题的研究成果较少。在二维直流电阻率法中，假设选取波数为N个，则点源问题的计算量为线源问题计算量的N倍，因此在相同计算量前提下，线源问题能够采用更密集的网格剖分方式，从而达到精细勘探的目的。本文作者通过对二维线源问题的正反演算法进行研究，从而为直流电阻率法的数据处理和解释提供经验和理论支撑。首先，将直接解法引入到线性方程组的求解中，提升数值模拟的精度和速度；然后，将地震全波形反演[31-32]中常用的高斯-牛顿反演算法[33-36]的改进方法引入到二维线源直流电阻率法中，即在高斯-牛顿反演算法中，若近似海森矩阵的对角线元素占优，可用对角海森矩阵代替近似海森矩阵，大大简化了高斯-牛顿迭代方程组的求解；最后，对三极装置和偶极-偶极装置的反演效果和探测能力进行了对比分析。

1 正演问题

1.1 直流电阻率法二维线源变分问题

在直角坐标系中，设z方向垂直向上，x垂直于z水平向右，y沿着构造走向方向并垂直于zx平面。对于线源二维地电问题，电位应满足的微分方程为

式中： (,)xzσ 为电导率；fs为源项；U为电位。

单极线电流源和偶极线电流源的表达式如式(2)所示：

式中：A、B为线电流源的供电位置；IA和IB分别为A和B点的供电电流；xA和zA为A点坐标；xB和zB为B点坐标。

根据不同的收发装置，可以推导出相应的混合边界条件，从而得到与其边值问题等价的变分问题，混合边界条件如式(3)和(4)所示：

单极线电流源，

偶极线电流源，

式中：Ω为求解区域；Γ∞为无穷远边界；rA为测点至电源点A的矢量距离；rB为测点至电源点B的矢量距离； ),cos( nrA和 ),cos( nrB分别是矢径rA和rB与外法向n的夹角的余弦值。

1.2 数值模拟

采用有限单元法求解上述变分问题[5]，其中基单元为三角形，插值函数在笛卡尔坐标系上呈线性变化。令 0)( =UFδ ，线性代数方程组如式(5)所示：

式中：K为稀疏实系数矩阵；b为源向量。

求解上述线性方程组，可以得到每个节点的电位值U。

大型稀疏线性方程组的求解目前主要有两类方法：迭代解法和直接解法。迭代解法优点在于消耗计算资源低和在许多情况下收敛较快；缺点在于当地电模型复杂(系数矩阵条件数过大)时，迭代收敛很慢或者出现发散，并行性低，且在有限步内求解精度无法保证。直接解法通常采用针对大型稀疏线性方程组开发的直接解法求解器，其优点在于求解精度高，多源情况下能实现并行回代；缺点在于需要耗费时间对系数矩阵进行分解，且需要占用大量的内存资源。

对于线性方程组(5)的求解，采用Pardiso_64位直接解法求解器。测试计算机配置为测试计算机配置为CPU-Inter Core i7-2600，主频为3.40 GHz，内存为8.00 GB。图1所示为不同模型大小时Pardiso_64位求解器求解线性方程组时间。从图1中可以看出，随着线性方程组维数的增加，求解时间呈非线性增长特征。

图1 不同模型Pardiso_64位求解器求解线性方程组时间 Fig.1 Time for solving linear equations by Pardiso_64 Solver with different models

2 反演问题

直流电阻率法二维线源的反演问题的表达式如式(6)所示：

式中：m为模型参数矢量；d为数据矢量；e为误差矢量；F是正演模型函数；F(m)是预测数据。

取测点个数为N，反演模型节点个数为M，则野外观测数据表示为d=[d1,d2,…,dN]T，模型参数矢量表示为m=[m1,m2,…,mM]T。

由于反演问题本身存在多解性和不稳定性，因此采用Tikhonov正则化反演思想建立反演目标函数如式(7)所示：

数据目标函数表示为

式中：T为矩阵转置符号；Wd为数据权系数矩阵。

构建模型目标函数有多种方法，根据不同的约束条件可以分为：最小模型约束、最平缓模型约束、最光滑模型约束等。本文作者采用最平缓模型约束，即模型参数导数的平方和最小，约束条件如式(9)所示：

2.1 高斯-牛顿法

传统的牛顿法收敛速度很快，但却需要不断计算海森矩阵的逆矩阵，而当海森矩阵病态或者非正定时，牛顿法不能保证收敛。因此，有学者发展了近似构造海森矩阵的逆矩阵的高斯-牛顿法[33,37]。高斯-牛顿法继承了牛顿法在初始模型参数合理的条件下收敛速度较快的特点，同时，通过引入阻尼项保证了反演过程的稳定性。

对式(7)中反演目标函数 ()φm在m0进行二阶Taylor-Lagrange展开，省略二阶以上的高阶项，然后在式子两端同时对m求一阶导数，通过求解得到反演迭代方程式：

式中：J为 MN× 的偏导数矩阵，每个元素的表达式为

PRATT等[38]在二维声波方程频率域波形反演中对高斯-牛顿法进行了详细的讨论，并分别对比了高斯牛顿法和全牛顿法的反演效果。但是在实际计算过程中不管是计算全海森矩阵还是近似海森的逆都非常耗时。对于弱非均匀介质，海森矩阵具有主对角线严格占优的特点，非对角线元素对反演结果的贡献相对较弱。因此，在反演过程中采用近似海森矩阵的主对角线代替海森矩阵，能够提升25%的高斯-牛顿法的反演效率。则式(10)中迭代公式可以修改为

式中：diag表示取对角线元素。

2.2 反演算法流程

反演算法流程如下：

1) 给定初始模型m0；

3) 计算灵敏度矩阵J和梯度TΔJ d；

4) 计算近似海森矩阵的对角线元素diag(H) ；

5) 求解得到高斯-牛顿方向 kmΔ ；

6) 从mk出发，沿方向 kmΔ 做一维搜索，确定最优步长tk；

8) 返回第2)步再次循环，k=k+1，直到满足循环终止条件，得最优反演模型mk+1。

与常规高斯-牛顿反演算法相比，本文作者简化了近似海森矩阵，从而大大节省了高斯-牛顿迭代方程组的求解时间。由于模型参数矢量m表示电导率向量σ的对数，即 σm lg= 。在第k次迭代过程中模型参数修正公式如式(13)所示：

3 模型算例

在正演计算中，引入直接解法有效地保证了正演精度。模型测试中网格剖分大小均为nx=301，nz=151。线性方程组求解采用Pardiso_64位直接解法求解器，测试计算机配置为：CPU-Inter Core i7-2600，主频为3.40 GHz，内存为8.00 GB。

3.1 正演模型算例

为了验证该正演算法的正确性与可行性，设计了线源在地表时的均匀半空间模型，模型电导率为0.01 S/m。图2所示为均匀半空间模型不同装置数值模拟结果，其中蓝线表示解析解电位，红线表示有限元法数值模拟电位；□表示三极装置的视电阻率误差；*表示偶极-偶极装置的视电阻率误差。从图2可知，除了源附近极少数点之外，整个计算区域的数值模拟结果与解析解吻合较好，总体平均误差小于0.17%。假设线源可以放置在地下，同样在电导率为0.01 S/m的均匀半模型情形下测试。图3所示为均匀半空间模型不同装置数值模拟结果，从图3可以看出，整个计算区域的数值模拟结果与解析解吻合较好，精度很高，总体平均误差小于0.35%。

图2 均匀半空间模型不同装置数值模拟结果(源在地表) Fig.2 Numerical simulation results in homogeneous half space model (Source is on earth surface): (a) Pole-dipole array; (b) Dipole-dipole array; (c) Apparent resistivity error curves of both arrays

图3 均匀半空间模型不同装置数值模拟结果(源在地下) Fig.3 Numerical simulation results in homogeneous half space model (Source is under earth surface): (a) Pole-dipole array; (b) Dipole-dipole array; (c) Apparent resistivity error curves of both arrays

3.2 反演模型算例

3.2.1 算例1(水平地形下组合模型反演)

算例1是在背景电阻率是50 Ω·m的均匀介质中，有两个大小均为100 m×20 m、电阻率分别是200 Ω·m和10 Ω·m、顶部埋深为30 m的异常体。图4所示为水平地形电阻率模型与反演剖面图。从两种装置的反 演结果来看，基本都能使得目标体恢复到原来的形态和位置，可以清晰地反映地下真实电性结构。其中偶极-偶极装置的反演效果明显优于三极装置的反演效果，尤其在边界上与真实模型吻合程度更高。

图4 水平地形电阻率模型与反演剖面图Fig.4 Resistivity model and inversion profiles: (a) True model; (b) Dipole-dipole array; (c) Pole-dipole array

图5 所示为水平地形两种装置的反演拟合误差与迭代次数变化曲线关系。由图5可以看出，二者在反演过程中都能很快收敛且迭代稳定。初始模型选择均匀半空间模型，偶极-偶极装置的迭代误差均比三极装置的迭代误差大，偶极-偶极装置的视电阻率异常明显大于三极装置的视电阻率异常，即探测异常体的能力更强、分辨率更高。

3.2.2 算例2(起伏地形下组合模型)

算例2是在背景电阻率是50 Ω·m的均匀介质中，有两个电阻率分别是200 Ω·m和10 Ω·m的异常体组合模型。图6所示为起伏地形电阻率模型与反演剖面图。从两种装置的反演结果中都可以比较准确地反映出组合异常体的埋深和形态。对比图6(b)和(c)可以发现偶极-偶极装置反演得到的异常幅值大小，边界范 围更接近于真实模型。

图5 水平地形两种装置反演迭代收敛曲线 Fig.5 Iteration convergence curves of both arrays inversions

图7 所示为起伏地形下两种装置反演拟合误差与迭代次数变化曲线关系。由图7可以看出，同样在反演过程中能比较快地收敛且迭代稳定，但三极装置反演收敛速度比偶极-偶极装置收敛速度稍快。初始模 型选择均匀半空间模型，起伏地形情况下偶极-偶极装置的视电阻率异常也是明显大于三极装置的视电阻率异常，再次验证了偶极-偶极装置的分辨率比三极装置的分辨率高。

图6 起伏地形电阻率模型与反演剖面图Fig.6 Resistivity model and inversion profiles including rugged topography: (a) True model; (b) Dipole-dipole array; (c) Pole-dipole array

图7 起伏地形下两种装置反演迭代收敛曲线 Fig.7 Iteration convergence curves of both arrays inversions including rugged topography

在算例2(起伏地形下组合模型)反演过程中，研究了反演迭代过程中每个源对应的视电阻率迭代拟合误差，并结合两种装置进行了对比分析。图8所示为不同源对应的视电阻率迭代拟合误差曲线。初始模型均采用均匀半空间模型，由图8(a)可知每个源对应的初次迭代拟合误差分布，偶极-偶极装置的视电阻率拟合误差明显大于三极装置的视电阻率拟合误差，表明在相同的条件下，偶极-偶极装置相对三极装置而言，对异常体更敏感，异常响应明显较大。由图8(b)可知每个源对应的最终迭代拟合误差分布，偶极-偶极装置的整体异常响应依然比三极装置的异常响应大，尤其在异常体中心正上方位置表现最为明显。结合图8(a)和(b)可以看出，在整个反演迭代过程中，偶极-偶极装置反映异常的能力均强于三极装置，且随着源的位置靠近异常体，偶极-偶极装置的异常幅值明显增大，而三极装置的异常幅值相对较小，进一步表明偶极-偶极装置具有更好的分辨率。

图8 起伏地形不同源对应的视电阻率迭代拟合误差曲线 Fig.8 Apparent resistivity fitting error curves with different sources including rugged topography: (a) Initial fitting error of iteration; (b) Last fitting error of iteration

4 结论

1) 从二维线源问题出发，以有限单元法为基础，对二维地电断面模型进行了数值模拟，并采用快速、高精度的直接解法求解线性方程组，采用不同的装置验证了算法的稳定性和准确性。

2) 将地震全波形反演中常用的高斯-牛顿反演算法的改进方法引入到二维线源直流电阻率法中，简化了高斯-牛顿迭代方程组的求解；相比传统反演算法具有高精度、收敛稳定、速度快、且耗费计算资源少等优势。从模型算例可以看出，提出的算法具有非常高的数值精度和计算效率，合成数据的反演结果和实际模型吻合度非常高，并且能够适应复杂地形。

3) 通过对三极装置和偶极-偶极装置的探测能力和反演效果进行了测试，数值实验结果表明：偶极-偶极装置具有更好的探测能力、更高的分辨率。

[1] 周熙襄,钟本善,严忠琼,江玉乐. 电法勘探正演数值模拟的若干结果[J]. 地球物理学报,1983,26(5): 479-491.

ZHOU Xi-xiang,ZHONG Ben-shan,YAN Zhong-qiong,JIANG Yu-le. Some results of resistivity modelling[J]. Chinese J Geophys,1983,26(5): 479-491.

[2] 罗延钟,万 乐. 二维地形不平条件下均匀外电场的有限差分模拟[J]. 物化探计算技术,1984,6(4): 15-26.

LUO Yan-zhong,WAN Le. Finite difference simulation for 2D uniform external electric field with irregular topography[J]. Computing Techniques for Geophysical and Geochemical Exploration,1984,6(4): 15-26.

[3] MUNDRY E. Geoelectrical model calculations for two- dimensional resistivity distributions[J]. Geophysical Prospecting,1984,32: 124-131.

[4] 罗延钟,孟永良. 关于用有限单元法对二位构造作电阻率法模拟的几个问题[J]. 地球物理学报,1986,29(6): 613-621.

LUO Yan-zhong,MENG Yong-liang. Some problems on resistivity modeling for two-dimensional structures by the finite element method[J]. Chinese J Geophys,1986,29(6): 631-621.

[5] 徐世浙. 地球物理中的有限单元法[M]. 北京: 科学出版社,1994: 131-197.

XU Shi-zhe. The finite element method in geophysics[M]. Beijing: Science Press,1994: 131-197.

[6] 底青云,王妙月. 稳定电流场有限元法模拟研究[J]. 地球物理学报,1998(2): 252-260.

DI Qing-yun,WANG Miao-yue. The real-like 2D FEM modeling research on the field characteristics of direct electric current field[J]. Chinese J Geophys,1998(2): 252-260.

[7] 阮百尧,徐世浙. 电导率分块线性变化二维地电断面电阻率侧深有限元数值模拟[J]. 地球科学-中国地质学报,1998,23(3): 303-307.

RUAN Bai-yao,XU Shi-zhe. Fem for modeling resistivity sounding on 2D geoelectric model with line variation of conductivity within each block[J]. Earth Science Journal of China University of Geosciences,1998,23(3): 303-307.

[8] 柳建新,刘海飞. 计算最优化离散波数的优化算法[J]. 物探化探计算技术,2005,27(1): 34-38.

LIU Jian-xin,LIU Hai-fei. The optimum algorithm for the calculation of the optimized discrete wave-number[J]. Computing Techniques for Geophysical and Geochemical Exploration,2005,27(1): 34-38.

[9] 汤井田,王飞燕. 基于非结构化网格的2.5D直流电阻率模拟[J]. 物探化探计算技术,2008,5: 413-418.

TANG Jing-tian,WANG Fei-yan. 2.5D direct current resistivity simulation based on the unstructured mesh[J]. Computing Techniques for Geophysical and Geochemical Exploration,2008,5: 413-418.

[10] TANG Jing-tian,WANG Fei-yan,REN Zheng-yong. 2.5-D DC resistivity modeling by adaptive finite-element method with unstructured triangulation[J]. Chinese J Geophys,2010,53(3): 708-716.

[11] 张东良,孙建国,孙章庆. 2维和2.5维起伏地表直流电法有限差分数值模拟[J]. 地球物理学报,2011,54(1): 234-244.

ZHANG Dong-liang,SUN Jian-guo,SUN Zhang-qing. Finete- difference DC electrical field modeling on 2D and 2.5D undulate topography[J]. Chinese J Geophys,2011,54(1): 234-244.

[12] PENZ S,CHAURIS H,DONNO D,MEHL C. Resistivity modelling with topography[J]. Geophysical Journal International,2013,45(2): 679-688.

[13] 麻昌英,柳建新. 复杂地形下高密度激电法2.5维有限单元法数值模拟[J]. 物探化探计算技术,2014,36(5): 405-409.

MA Chang-ying,LIU Jian-xin. 2.5 dimensional numerical simulation of finite element method of IP high density in complex terrain[J]. Computing Techniques for Geophysical and Geochemical Exploration,2014,36(5): 405-409.

[14] 彭艳华,柳建新,刘海飞,孙丽影. 直流电阻率2.5维有限元数值模拟[J]. 地球物理学进展,2014,29(4): 1926-1931.

PENG Yan-hua,LIU Jian-xin,LIU Hai-fei,SUN Li-ying. 2.5D direct current resistivity modelling by finite-element method[J]. Progress in Geophysics,2014,29(4): 1926-1931.

[15] 肖 晓,原 源,汤井田. 基于有限元-无限元耦合的2.5D直流电阻率数值模拟[J]. 中南大学学报(自然科学版),2014,45(8): 2691-2700.

XIAO Xiao,YUAN Yuan,TANG Jing-tian. 2.5-D DC resistivity forward modeling by finite-infinite element coupling method[J]. Journal of Central South University (Science and Technology),2014,45(8): 2691-2700.

[16] 胡宏伶,肖 晓,潘克家,汤井田. 基于局部加密等级网格的2.5D直流电法有限元模拟[J]. 中南大学学报(自然科学版),2014,45(7): 2260-2267.

HU Hong-ling,XIAO Xiao,PAN Ke-jia,TANG Jing-tian. Finite element modeling of 2.5D DC resistivity based on locally refined graded mesh[J]. Journal of Central South University (Science and Technology),2014,45(7): 2260-2267.

[17] BARKER R. A simple algorithm for electrical imaging of the subsurface[J]. First Break,1992,10(2): 25-32.

[18] ZOHDY A A R. A new method for the automatic interpretation of Schlumberger and Wenner sounding curves[J]. Geophysics,1989,54(2): 245-253.

[19] LOKE M H,BARKER R D. Least-squares deconvolution of apparent resistivity pseudosetions[J]. Geophysics,1995,60(6): 1682-1690.

[20] LOKE M H,BARKER R D. Improvement to the Zohdy method for the inversion of resistivity sounding and pseudosection data[J]. Computer and Geosciences,1995,21(2): 321-332.

[21] 阮百尧,徐世浙. 二维直流电阻率测深曲线的快速反演[J]. 物探与化探,1996,20(6): 455-460.

RUAN Bai-yao,XU Shi-zhe. Rapid inversion of two-dimensional DC Resistivity sounding curve[J]. Geophysical and Geochemical Exploration,1996,20(6): 455-460.

[22] 徐海浪,吴小平. 电阻率二维神经网络反演[J]. 地球物理学报,2006,49(2): 584-589.

XU Hai-lang,WU Xiao-ping. 2D resistivity inversion using the neural network method[J]. Chinese J Geophys,2006,49(2): 584-589.

[23] 刘海飞,阮百尧. 直流电阻率最小二乘间歇反演[J]. 物探与化探计算技术,2007,29(4): 339-341.

LIU Hai-fei,RUAN Bai-yao. The least squares time-lapse inversion of DC resistivity[J]. Computing Techniques for Geophysical and Geochemical Exploration,2007,29(4): 339-341.

[24] 程 勃,底青云. 基于遗传算法和统计学的电阻率测深二维反演研究[J]. 地球物理学进展,2012,27(2): 0788-0795.

CHENG Bo,DI Qing-yun. 2D inversion of resistivity sounding base on GA and statistic method[J]. Progress in Geophysics,2012,27(2): 0788-0795.

[25] 程 勃,底青云. 基于偏导数和统计学方法的电阻率测深二维反演[J]. 石油地球物理勘探,2012,47(6): 1006-1013.

CHENG Bo,DI Qing-yun. 2D inversion of resistivity sounding based on partial derivative and statistics[J]. Oil Geophysical Prospecting,2012,47(6): 1006-1013.

[26] 韩 波,窦以鑫,丁 亮. 电阻率成像的混合正则化反演算法[J]. 地球物理学报,2012,55(3): 970-980.

HAN Bo,DOU Yi-xin,DING Liang. Electrical resistivity tomography by using a hybrid regularization[J]. Chinese J Geophys,2012,55(3): 970-980.

[27] 戴前伟,江沸菠. 基于混沌振荡PSO-BP算法的电阻率层析成像非线性反演[J]. 中国有色金属学报,2013,23(10): 2898-2904.

DAI Qian-wei,JIANG Fei-bo. Nonlinear inversion for electrical resistivity tomography based on chaotic oscillation PSO-BP algorithm[J]. The Chinese Journal of Nonferrous Metals,2013,23(10): 2898-2904.

[28] 冯德山,杨炳坤,戴前伟,王鹏飞. 超高密度电法全四极装置正反演[J]. 地球物理学进展,2014,29(2): 0718-0724.

FENG De-shan,YANG Bing-kun,DAI Qian-wei,WANG Peng-fei. The simulation and invertion for the full-four-pole array of ultra-high density resistivity method[J]. Progress in Geophysics,2014,29(2): 0718-0724.

[29] 冯德山,王鹏飞,杨炳坤. 超高密度电法有限单元法正演与广义最小二乘反演[J]. 中国有色金属学报,2014,24(3): 793-800.

FENG De-shan,WANG Peng-fei,YANG Bing-kun. Finite element forward simulation and least square generalized inversion of ultra-high density resistivity method[J]. The Chinese Journal of Nonferrous Metals,2014,24(3): 793-800.

[30] 柳建新,彭艳华. 直流电阻率测深二维自适应正则化反演[J]. 物探化探计算技术,2014,36(4): 385-388.

LIU Jian-xin,PENG Yan-hua. Two dimensional self-adaptive regularization inversion for DC resistivity data with active constraint balancing[J]. Computing Techniques for Geophysical and Geochemical Exploration,2014,36(4): 385-388.

[31] FICHTNER A,TRAMPERT J. Resolution analysis in full waveform inversion[J]. Geophysical Journal International,2011,187(3): 1604-1624.

[32] 胡光辉,贾春梅,夏洪瑞,贺剑波,宋 林,沈忠秋. 三维声波全波形反演的实现与验证[J]. 石油物探,2013,52(4): 417-425.

HU Guang-hui,JIA Chun-mei,XIA Hong-rui,HE Jian-bo,SONG Lin,SHEN Zhong-qiu. Implementation and validation of 3D acoustic full waveform inversion[J]. Geophysical Prospecting for Petroleum,2013,52(4): 417-425.

[33] BAE H S,PYUN S,CHUNG W. Frequency-domain acoustic- elastic coupled waveform inversion using the Gauss-Newton conjugate gradient method[J]. Geophysical Prospecting,2012,60(3): 413-432.

[34] JEONG W,LEE H Y,MIN D J. Full waveform inversion strategy for density in the frequency domain[J]. Geophysical Journal International,2012,188(3): 1221-1242.

[35] JOO Y,SEOL S J,BYUN J. Acoustic full- waveform inversion of surface seismic data using the Gauss-Newton method with active constraint balancing[J]. Geophysical Journal International,2013,61(1): 166-182.

[36] 高凤霞,刘 财,冯 晅. 几种优化方法在频率域全波形反演中的应用效果及对比分析研究[J]. 地球物理学进展,2013,28(4): 2060-2068.

GAO Feng-xia,LIU Cai,FENG Yuan. Comparisons and analyses of several optimization methods in the application of frequency-domain full waveform inversion[J]. Progress in Geophysics,2013,28(4): 2060-2068.

[37] RAVAUT C,OPERTO S,IMPROTA L,VIRIEUX J,HERRERO A,DELL’AVERSANA P. Multiscale imaging of complex structures from multifold wide-aperture Seismic data by frequency-domain full-waveform tomography: Application to a thrust belt[J]. Geophysical Journal International,2004,159(3): 1032-1056.

[38] PRATT R,SHIN C,HICKS G. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[J]. Geophysical Journal International,1998,133(2): 341-362.