位场向下延拓的CGNR法

张志厚¹，吴乐园¹，王瑞赛²，张金会³

(1. 浙江大学 地球科学系，浙江 杭州，310027；

2. 云南省有色地质 昆明勘测设计院，云南 昆明，650031；

3.安徽省勘察技术院，安徽 合肥，233005)

摘要：假定位场向下延拓的系数矩阵为正定的条件下，直接采用共轭梯度法求解位场向下延拓的第一类Fredholm型积分方程。理论模型试验表明：该方法收敛速度快，但抑噪能力较差。与位场向下延拓的积分迭代法相比，该方法收敛速度较快、发散也较快、适用性较差。为获得稳定近似解，将此不适定问题实施正则化过程，转化为最小二乘求极小值问题，再采用共轭梯度法迭代求解，实现抑噪能力较强的位场向下延拓的CGNR法。理论模型检验表明：位场向下延拓的CGNR法抑制噪声能力较强，并且相比同样具有抑制噪声能力较强的最小二乘最速下降法，CGNR法收敛速度很快，具有明显的计算优势。

关键词：向下延拓；CGNR；共轭梯度法；收敛速度；抑制噪声

中图分类号：P631          文献标志码：A          文章编号：1672-7207(2013)08-3273-09

CGNR method for potential field downward continuation

ZHANG Zhihou¹, WU Leyuan¹, WANG Ruisai², ZHANG Jinhui³

(1. Department of Earth Sciences, Zhejiang University, Hangzhou 310027, China;

2. Kunming Exploration and Design Institute, Yunnan Nonferrous Geology Bureau, Kunming 650031, China;

3. Geological Exploration Technologies institute of Anhui Province, Hefei 233005, China)

Abstract: The first kind of Fredholm integral equation in the downward continuation was solved by adopting conjugate gradient method with the assumption that the coefficient matrix is positive definite. The theoretical model test shows that downward continuation has a very fast convergence rate but suggests poor ability to compress the noise. Compared with integral iteration method, conjugate gradient method indicates a faster convergence rate, as well as divergence rate, and is not so applicable. In order to acquire stable solutions the regularization of the operator of the ill-posed problem was carried out and the problem was transformed to solve the minimum value of the least squares estimations. After using the conjugate gradient method to solve the least squares problem, a method was implemented for downward continuation based on CGNR method that has a better ability to suppress noise. Model test demonstrates that this method has strong ability in restraining noise. This method convergence fast and has obvious computation advantage by a comparison with least square steepest descent method.

Key words: downward continuation; CGNR; conjugate gradient method; convergence rate; anti-noise ability

位场(重磁)勘探作为非地震领域的重要分支，在研究地球内部结构、油气与固体矿产资源勘探、区域与深部构造研究、环境与工程勘探、考古等诸多领域发挥着重要的作用。位场的向下延拓运算是许多位场数据处理和反演计算中非常重要的步骤。此外，向下延拓后的场(重磁)可能分离出水平方向上“叠加”的弱异常，突出有用信息，有利于提高位场数据解释的可靠性；利用位场进行辅助导航^[1-6]，也使得位场向下延拓的方法研究备受我国诸多学者关注^[7-8]；位场的向下延拓归结为求解第一类算子方程，不适定性是其典型特征，因此，位场的向下延拓问题是难点，也是热点，研究精度高、效果好的大深度延拓方法有重要的意义。位场向下延拓的方法研究最早可追溯到20世纪40年代，Bateman^[9]和Peters^[10]采用泰勒级数法进行位场的向下延拓，延拓深度不超过2倍点距，为提高延拓深度，Pilkington等^[11]对原始数据进行低通滤波，并将航磁数据延拓至起伏地形线，Fedi等^[12]利用空间域有限差分并结合位场的拉普拉斯方程求解已知场的垂向各阶导数，实现位场向下延拓的泰勒级数展开法，该方法收敛速度快，但对噪声的放大也很快。20世纪70年代Dampancy^[13]和Syberg^[14]就已使用等效源法进行位场的二维和三维向下延拓，到目前为止等效源法仍是位场向下延拓的常用方法之一，如刘天佑等^[15]和黄翼坚等^[16]采用同样等效源法进行位场的向下延拓。实质上，等效源法最终也将求解一个第一类算子方程，最速下降法求解该方程计算时间长，使用价值不大^[17]，马奎特方法求得的等效源参量分布不稳定，下延效果差^[17]，非迭代直接算法的不稳定性和大计算量也是该方法的不利因素。在波数域利用FFT方法进行位场向下延拓不仅计算速度高，而且不会受到数据量大的限制，但是向下延拓极为不稳定。Clarke^[18]利用Weiner滤波设计向下延拓最佳滤波，消弱噪声的干扰，但延拓点距仍然有限。梁锦文^[19]讨论正则化方法位场向下延拓的波数响应及特点，简单的说，正则化波数响应曲线在低波数区与理论位场向下延拓的波数响应曲线是吻合的，在高波数区是抑制高频的。王邦华等^[20]利用引入位场频率，埋深与正则化因子等有关的校正函数，提出一种新的正则化方法，使得向下延拓过场源体不奇异，实质上，该校正函数有抑制高频成分的作用。姚长利等^[21]提出位场转换的抽样分组法，该方法是对离散数据进行等间距抽样重新分组，每1组做为独立的部分进行转换处理，最后将每组的结果合成，如果分组数选择合适，可以取得较好的效果。陈生昌等^[22]利用傅里叶变换矩阵的正交对称性，并结合矩阵的奇异值分解和广义逆原理，提出一种稳定的不需要进行求逆运算的位场向下延拓方法，一定程度上解决位场大深度向下延拓的不稳定性问题，该方法效率很高，但存在阻尼因子的选择问题。徐世浙提出的位场向下延拓的积分迭代法(IIM)，延拓深度大，可达到10倍甚至20倍点距以上^[23-26]，效果好^[25]，诸多学者讨论该方法的收敛性^[27-29]、噪声干扰^[30]、正则性^[31]以及地球物理含义^[32]等，使该方法得到极大发展。该方法不足之处是迭代过程中放大了噪声^{[30, 33]}。

共轭梯度法(CGM)是一种有效的反问题求解方法^[34]，已有学者将该方法应用在三维电阻率反演^[35]、地震偏移成像中^[36]。本文作者将该方法应用到位场的向下延拓中，模型试验的结果确表明该方法的适用性较差，为获稳定近似解，将向下延拓问题转化为最小二乘问题，实现位场向下延拓的CGNR法；在迭代计算中使用FFT法解决了数据量大的计算问题，提高计算效率；并且通过理论模型对本文所提出的方法进行检验，同时也与位场向下延拓的最小二乘最速下降法^[37](SDM)的收敛速度进行比对，最后对位场向下延拓的CGNR法在大噪声水平下的抑制噪声能力进行检验。

1 基本原理

1.1 空间域位场问题的表述

假设观测面为水平面(z=0)，该平面上的位场：

         (1)

计算观测面以下至场源以上空间(z=h＞0，z轴方向向下)平面的位场：

         (2)

根据向上延拓的公式，式(1)和(2)的关系为

 (3)

式(3)可以表示成褶积的形式：

        (4)

其中：

      (5)

数学上，位场的向上延拓可以表述为式(3)的1个第一类Fredholm积分方程。已知u₀，通过该积分方程求解u_h即为向下延拓。由于观测的位场数据与所求的向下延拓数据都是离散的，将式(3)离散化，并认为已知场的离散化数据与其正下方的待延拓点为一一对应的关系，离散化形式如下：

 (6)

其中：△x与△y分别为x与y方向的取样间距，M与N分别为x与y方向的取样总数，m=1，2，3，…，M，n=1，2，3，…，N。

式(6)用矩阵可以写成：

               (7)

其中：A为L×L(M×N)维矩阵，也称为位场向下延拓的系数矩阵，且为对称矩阵；U₀为已知位场组成的L×1维列向量；U_h为待求L×1维列向量。

1.2 CGM

对于式(7)，为表示方便，在不混淆的情况下，令X=U_h，B=U₀，则式(7)可以表示为

AX=B                  (8)

由于系数矩阵A对称矩阵，假设其正定，则式(8)对应的等价极值问题为

            (9)

采用CGM求解，其迭代过程描述如下：

            (10)

其中：

         (11)

式中：表示对向量求内积；利用梯度：

       (12)

构造P^(k)的共轭向量：

           (13)

其中：

       (14)

令X⁽⁰⁾=U₀开始迭代，且；式(10)~(14)就是求解式(7)的CGM，对于给定的精度ε，迭代到为止，其中ε是一个适当小的数，且存在，则认为。式(10)~(14)中的AP^(k)，AX^(k)与AR^(k) 均采用波数域的方法计算，可有效地提高计算效率。

1.3 CGNR法

考虑到不适定性是第一类算子的典型特性，这将导致观测数据中极小的误差都会使其解产生巨大的震荡现象^[38-39]，采用式(10)~(14)进行位场的向下延拓往往也可能会存在结果精度低或者过度迭代的缺点，且系数矩阵A是在假设正定的条件下才可以延拓计算，因此，有必要采用适用方法求解。为此，将式(8)实施正则化过程，转化为最小二乘求极小值问题：

      (15)

其中：||·||₂表示目标函数的二范数。式(15)实质上是求解方程组

             (16)

已知式(16)的系数矩阵A^TA为对称正定矩阵，此时采用CGM求解式(16)理论上是完备的。即用CGM方法求解与式(8)具有相同解的对称正定方程组(式(16))，这种方法称为CGNR方法，它属于一Krylov子空间

         (17)

的法方程组方法^[40]，其中 为初始残差。

由于反问题的不适定性，所以A的谱很“坏”，使得A^TA的谱更加分散，式(16)的病态性更为严重，从而使得迭代求解的收敛速度会变慢。然而大量的试验表明CGNR在许多场合下有较好的数值表现^[41-43]。长期以来，研究者过分看重了系数矩阵A^TA的“平方条件数”效应所带来的负面影响，而实际情况却并非如此。值得特别指出的是：CGNR不但能使残差范数单调递减^[43-44]，而且还具有使误差范数单调递减^[41,44]这一良好的理论特性；在稳健性方面，该方法有着长久的生命力^[44]。

鉴于位场向下延拓CGNR方法与CGM方法的区别是采用共轭梯度法分别求解式(16)与式(8)，故而将A^TA和A^TB分别替代式(10)~(14)中的A和B，则式(10)~(14)就转化为CGNR方法，下面给出CGNR方法的具体迭代过程。

步骤1：令，，，，对k=0，1，2，…进行以下迭代；

步骤2：；

步骤3：；

步骤4：；

步骤5：；

步骤6：；

步骤7：；

步骤8：。

2 理论模型试验

为检验本文提出的延拓方法的有效性，及验证直接采用CGM求解位场的向下延拓方程组的不稳定性，进行模型试验。采用球体在不同观察面的理论磁异常与向下延拓异常的对比来检验延拓效果，并且与位场向下延拓的积分迭代法、LS-SDM进行对比来检验CGNR法的收敛速度。球体模型参数如表1所示，其中：(x₀，y₀，h₀)为球体中心坐标；r₀为球体半径；M为磁化强度；I₀，A₀，I′和A′分别为地磁场倾角、剖面磁偏角、磁化强度倾角和剖面与磁北夹角。根据球体磁场的计算公式^[45](坐标系z轴方向向下)，分别计算z=1 km平面与z=0 km平面的理论磁异常，如图1和图2所示，其中计算网格数为401×401，数据资料点距为50 m。将图2所示的理论磁异常分别使用位场向下延拓的积分迭代法、CGM、LS-SDM以及CGNR法进行向下延拓20倍的点距，延拓结果如图3所示。由图3可见：位场向下延拓的积分迭代法与CGM的延拓结果的边部(见图3(a)和3(b))稍比LS-SDM与CGNR法的延拓结果的边部好(见图3(c)和3(d))。但总体而言，这4种方法的延拓结果除边部的零等值线略有偏差，其余部分都与理论场(图1)有很高的相似程度。

为定量的描述这4种方法延拓的精度、收敛性以及其与位场向下延拓的积分迭代法的收敛速度相比较，采用均方误差(MSE)来衡量随着迭代次数n的增加，其延拓结果与理论值的差异水平的变化，均方误差公式为：

   (18)

其中：u_con和u_theo分别为延拓值和理论值；l为用于计算的点数。

表1 球体模型参数表

Table 1 Parameters of globe model

图1 z=1 km平面理论磁异常等值线图(nT)

Fig. 1 Contour map of theoretical magnetic anomaly (in nT) at z=1 km

图2 z=0 km平面理论磁异常等值线图(nT)

Fig. 2 Contour map of theoretical magnetic anomaly (in nT) at z=0 km

图4所示为图2所示的理论磁异常分别采用位场向下延拓的积分迭代法、CGM、LS-SDM和CGNR方法随着迭代次数的增加，延拓结果的均方误差的变化。由图4可见：随着迭代次数的增加，均方误差都在急剧衰减，并且共轭梯度法的收敛速度比积分迭代法的收敛速度快，积分迭代法的收敛速度优于CGNR法，LS-SDM的收敛速度最慢。以上说明在无噪声模型下，CGM具有较好的收敛速度，及相对较好的计算精度。然而，在实际应用中，观测数据不可避免地带有各种噪声，在含噪声情况下，是否也存在这一计算优势？为此，有必要设计含噪声的模型试验。

假定该噪声具有可加性，则：

 (19)

其中：U₀为实测数据；△为(0，1)区间的噪声水平；rand(size(U_true))具有与U_true维数一致的随机噪声，U_true由无噪声模型计算生成，max(U_true)为U_true的最大值。

对图2所示的理论模型数据中加噪声水平为0.01的噪声，含噪声的磁异常如图5所示。将图5所示的含噪声的磁异常数据分别采用位场向下延拓的积分迭代法、CGM、LS-SDM以及CGNR方法向下延拓20倍的点距至z=1 km平面，各种延拓方法随着迭代次数的增加，延拓结果的均方误差的变化如图6所示，由图6可见：(1) CGM对含有噪声的数据进行延拓，其收敛的很快，发散的也很快，与积分迭代法相比较，积分迭代法的收敛的不是很快，发散的也不是很快。可认为在含噪声情况下直接采用CGM进行位场下延计算，其迭代的长期行为是一个不确定的，发散的结果。(2) LS-SDM与CGNR方法随着迭代次数的增加，延拓结果的均方误差在减小，可见最小二乘意义下的这2种梯度方法抗噪能力较好，从均方误差的衰减曲线上看，CGNR方法延拓结果的精度要高于LS-SDM延拓结果的精度。

图3 z=0 km平面数据向下延拓1 km得到的z=1 km平面磁异常等值线图(nT)

Fig. 3 Contour map of magnetic anomaly (in nT) at z=1 km from downward continuing magnetic (in nT) at z=0 km

图4 迭代次数与均方误差的关系

Fig. 4 Relationship between iteration number and MSE

图5 z=0 km平面含有噪声的磁异常等值线图(nT)

Fig. 5 Contour map of theoretical magnetic anomaly(in nT) with noise at z=0 km

如果通过有效的停机准则能够使CGM在均方误差达到极小时停机，不让其产生迭代不足或者过度迭代的现象，能否获得较好的延拓结果？本文作者给出CGM在含噪声情况下的理想迭代结果(迭代次数为4次)，如图7(b)所示。由图7(b)可见：CGM方法的理想迭代结果也并不理想，该方法迭代过程中累加噪声，影响延拓结果的精度。其余3种方法的延拓结果分别如图7(a)，7(c)和7(d)所示。由图7(c)和(d)可见：积分迭代法的抑制噪声能力差，LS-SDM与CGNR方法的延拓结果除去边部与零等值线外，其余部分效果较好。

图6 迭代次数与均方误差的关系

Fig. 6 Relationship between iteration number and MSE

图7 z=0 km平面数据向下延拓1 km得到的z=1 km平面磁异常等值线图(nT)

Fig. 7 Contour map of magnetic anomaly (in nT) at z=1 km from downward continuing magnetic (in nT) at z=0 km

综上所述，位场向下延拓的CGNR方法是抑制噪声能力较强的下延方法中的收敛速度快、延拓结果精度相对较高的位场大深度向下延拓方法。

为进一步检验位场向下延拓的CGNR方法的稳健性，对图2所示的理论磁异常分别叠加噪声水平△为10%，20%及30%的随机噪声，然后采用CGNR方法向下延拓20倍点距，延拓结果的均方误差与迭代次数的关系如图8所示。由图8可见：在高噪声水平下，位场向下延拓的CGNR方法仍具有良好的稳健性。

图8 迭代次数与均方误差的关系

Fig. 8 Relationship between iteration number and MSE

3 结论

(1) 位场的向下延拓问题是典型的褶积型线性反问题，由于其具有不适定性，为此将该问题转化为最小二乘问题，采用共轭梯度法进行求解，实现了位场向下延拓的CGNR方法。

(2) 在假设系数矩阵A为正定的条件下，直接采用共轭梯度法进行位场的向下延拓，发现该方法抑制噪声能力差。理论模型检验表明：位场向下延拓的共轭梯度法相比积分迭代法收敛速度较快，在含噪声的情况下，其发散的速度也很快，且位场向下延拓的共轭梯度法的迭代长期行为是一个发散的结果。

(3) LS-SDM与CGNR方法对高频噪声都有很好的抑制作用；不含噪声理论模型与含噪声理论模型试验表明，位场向下延拓的CGNR方法的收敛速度优于LS-SDM的收敛速度，延拓结果的精度相对较高，具有明显的计算优势。

(4) 高噪声情况下的模型试验进一步说明本文提出的CGNR下延方法具有良好的稳健性。

参考文献：

[1] 许大欣. 利用重力异常匹配技术实现潜艇导航[J]. 地球物理学报, 2005, 48(4): 812-816.

XU Daxin. Using gravity anomaly matching techniques to implement submarine navigation[J]. Chinese J Geophys, 2005, 48(4): 812-816.

[2] 彭富清. 地磁模型与地磁导航[J]. 海洋测绘, 2006, 26(2): 73-75.

PENG Fuqing. Geomagnetic model and geomagnetic navigation[J]. Hydrographic Surveying and Charting, 2006, 26(2): 73-75.

[3] 杨云涛, 石志勇, 关贞珍, 等. 地磁场在导航定位系统中的应用[J]. 中国惯性技术学报, 2007, 15(6): 686-692.

YANG Yuntao, SHI Zhiyong, GUAN Zhenzhen, et al. Application of geomagnetic field in navigation and localization system[J]. Journal of Chinese Inertial Technology, 2007, 15(6): 686-692.

[4] 郭才发, 胡正, 张士峰, 等. 地磁导航综述[J]. 宇航学报, 2009, 30(4): 1314-1319.

GUO Caifa, HU Zheng, ZHANG Shifeng, et al. A survey of geomagnetic navigation[J]. Journal of Astronautics, 2009, 30(4): 1314-1319.

[5] 杨昆, 康戈文, 李洪. 重力场和地磁场综合匹配在导航中的应用[J]. 航海工程, 2010, 39(1): 129-131.

YANG Kun, KANG Gewen, LI Hong, et al. Application of gravity and geomagnetism matching in navigation[J]. Ship & Ocean Engineering, 2010, 39(1): 129-131.

[6] 王向磊, 田艳峰. 基于地磁场的自主导航研究[J]. 地球物理学报, 2010, 53(11): 2724-2731.

WANG Xianglei, TIAN Yanfeng. Autonomous navigation based Geomagnetic research[J]. Chinese J Geophys, 2010, 53(11): 2724-2731.

[7] 徐世浙, 张秀达, 张昌达. 水下磁定位的若干问题[J]. 海军大连舰艇学院学报, 2007, 30(3): 4-6.

XU Shizhe, ZHANG Xiuda, ZHANG Changda. Several problems of underwater magnetic location[J]. Journal of Dalian Naval Academy, 2007, 30(3): 4-6.

[8] 陈龙伟, 张辉, 郑志强, 等. 水下地磁辅助导航中地磁场延拓方法[J]. 中国惯性技术学报, 2007, 15(6): 693-697.

CHEN Longwei, ZHANG Hui, ZHENG Zhiqiang, et al. Technique of geomagnetic field continuation in underwater geomagnetic aided navigation[J]. Journal of Chinese Inertial Technology, 2007, 15(6): 693-697.

[9] Bateman H. Some integral equations of potential theory[J]. J of Appl Physics, 1946, 17: 91-102.

[10] Peters L J. The direct approach to magnetic interpretation and its practical application[J]. Geophysics, 1949, 14: 290-320.

[11] Pilkington M, Roest W R. Draping acromagnetic data using the Taylor series[J]. Exploration Geophysics, 1992, 23: 373-380.

[12] Fedi M, Florio G. A stable downward continuation by using the ISVD method[J]. Geophys J Int, 2002, 151(1): 146-156.

[13] Dampancy C N G. The equivalent source technique[J]. Geophysics, 1969, 34(1): 39-53.

[14] Syberg F J R. Potential field continuation between general surfaces[J]. Geophysical Prospecting, 1972, 20(2): 267-282.

[15] 刘天佑, 杨宇山, 李媛媛, 等. 大型积分方程降阶解法与重力资料曲面延拓[J]. 地球物理学报, 2007, 50(1): 290-296.

LIU Tianyou, YANG Yushan, LI Yuanyuan, et al. The order-depression solution for large-scale integral equation and its application in the reduction of gravity data to a horizontal plane[J]. Chinese J Geophys, 2007, 50(1): 290-296.

[16] 黄翼坚, 王万银, 于长春. 等效源法三维随机点位场数据处理和转化[J]. 地球物理学进展, 2009, 24(1): 95-101.

HUANG Yijian, WANG Wanyin, YU Changchun. Processing and transform of 3-D arbitrarily distributed potential filed data using an equivalent source approach[J]. Progress in Geophysics, 2009, 24(6): 1095-1097.

[17] 钟本善, 周熙襄, 沙椿. 等效源法在重磁异常解释中应用的进一步探讨[J]. 物探与化探, 1986, 10(3): 201-206.

ZHONG Benshan, ZHOU Xixiang, SHA Chun. A further discussion on the application of equivalent source method to the interpretation of gravity and magnetic anomalies[J]. Geophysical & Geochemical Exploration, 1986, 10(3): 201-206.

[18] Clarke G K C. Optimum second-derivative and downward continuation filter[J]. Geophysics, 1969, 34: 424-437.

[19] 梁锦文. 位场向下延拓的正则化方法[J]. 地球物理学报, 1989, 32(5): 600-608.

LIANG Jinwen. Downward continuation of regularization methods for potential fields[J]. Chinese J Geophys, 1998, 32(5): 600-608.

[20] 王邦华, 王理. 重磁位场的正则化向下延拓[J]. 物探化探计算技术, 1998, 20(1): 30-35.

WANG Banghua, WANG Li. A new normalized method of downward extrapolation for potential field[J]. Computing Techniques for Geophysical and Geochemical Exploration, 1998, 20(1): 30-35.

[21] 姚长利, 管志宁, 黄卫宁. 位场转化的抽样分组法[J]. 石油地球物理勘探, 1997, 32(5): 696-702, 723.

YAO Changli, GUAN Zhining, HUANG Weining. Resampled grouping method for potential field transformation[J]. Oil Geophusical Prospecting, 1997, 32(5): 696-702, 723.

[22] 陈生昌, 肖鹏飞. 位场向下延拓的波数域广义逆算[J]. 地球物理学报, 2007, 50(6): 1816-1822.

CHEN Shengchang, XIAO Penfei. Wavenumber domain generalized inverse algorithm for potential field downward continuation[J]. Chinese J Geophys, 2007, 50(6): 1816-1822.

[23] 徐世浙. 位场大深度向下延拓[J]. 物探化探计算技术, 2006, 28(增刊): 29-31.

XU Shizhe. Downward continuation of long depth for potential fields[J]. Computing Techniques for Geophysical and Geochemical Exploration, 2006, 28(S): 29-31.

[24] XU Shizhe, YANG Jinyu, YANG Changfu, et al. The iteration method for downward continuation of a potential field from a horizontal plane[J]. Geophysical Prospecting, 2007, 55: 883-890.

[25] 徐世浙. 迭代法与FFT法位场向下延拓效果的比较[J]. 地球物理学报, 2007, 50(1): 285-289.

XU Shizhe. A comparison of effects between the iteration method and FFT for downward continuation of potential fields[J]. Chinese J Geophys, 2007, 50(1): 285-289.

[26] 徐世浙, 余海龙. 位场曲化平的插值—迭代法[J]. 地球物理学报, 2007, 50(6): 1811-1815.

XU Shizhe, YU Hailong. The interpolation—Iteration method for potential field continuation from undulating surface to plan[J]. Chinese J Geophys, 2007, 50(6): 1811-1815.

[27] 张辉, 陈龙伟, 任治新, 等. 位场向下延拓迭代法收敛性分析及稳健向下延拓方法研究[J]. 地球物理学报, 2009, 52(4): 1107-1113.

ZHANG Hui, CHEN Longwei, REN Zhixing, et al. Analysis on convergence of iteration method for potential fields downward continuation and research on robust downward continuation method[J]. Chinese J Geophys, 2009, 52(4): 1107-1113.

[28] 刘东甲, 洪天求, 贾志海等. 位场向下延拓的波数域迭代法及其收敛性[J]. 地球物理学报, 2009, 52(6): 1599-1605.

LIU Dongjia, HONG Tianqiu, JIA Zhihai, et al. Wave number domain iteration method for downward continuation of potential fields and its convergence[J]. Chinese J. Geophys, 2009, 52(6): 1599-1605.

[29] 王顺杰, 朱海, 栾禄雨. 水下地磁导航中位场积分迭代法收敛性分析[J]. 地球物理学进展, 2009, 24(3): 1095-1097.

WANG Shunjie, ZHU Hai, LUAN Luyu. Constringency analysis of the iteration method for continuation of potential fields in underwater geomagnetism navigation[J]. Progress in Geophysics, 2009, 24(6): 1095-1097.

[30] 于波, 翟国君, 刘燕春, 等. 噪声对磁场向下延拓迭代法的计算误差影响分析[J]. 地球物理学报, 2009, 52(8): 2182-2188.

YU Bo, ZHAI Guojun, LIU Yanchun, et al. Analysis of effect on the calculation error of downward continuation with iteration method[J]. Chinese J Geophys, 2009, 52(8): 2182-2188.

[31] 曾小牛, 李夕海, 刘代志, 等. 积分迭代法的正则性分析及其最优步长的选择[J]. 地球物理学报, 2011, 54(11): 2943-2950.

ZENG Xiaoniu, LI Xihai, LIU Daizhi, et al. Regularization analysis of integral iteration method and the choice of its optimal step-length[J]. Chinese J Geophys, 2011, 54(11): 2943-2950.

[32] 骆遥. 位场迭代法向下延拓的地球物理含义—以可下延异常逐次分离过程为例[J]. 地球物理学进展, 2011, 26(4): 1197-1200.

LUO Yao. Geophysical implications of the iteration method for downward continuation of potential field: a case by computiong regional anomaly can be continued downward[J]. Progress in Geophysics, 2011, 26(4): 1197-1200.

[33] 曾小牛, 李夕海, 韩绍卿, 等. 位场向下延拓三种迭代方法之比较[J]. 地球物理学进展, 2011, 26(3): 908-915.

ZENG Xiaoniu, LI Xihai, HAN Shaoqin, et al. A comparison of three iteration methods for downward continuation of potential fields[J]. Progress in Geophysics, 2011, 26(3): 908-915.

[34] 肖庭延, 于慎根, 王彦飞. 反问题的数值解法[M]. 北京: 科学出版社, 2003: 30-31.

XIAO Tingyan, YU Shengen, WANG Yanfei. Numerical methods for the solution of inverse problems[M]. Beijing: Science Press, 2003: 30-31.

[35] 吴小平, 徐果明. 利用共轭梯度法的电阻率三维反演研究[J]. 地球物理学报, 2000, 43(3): 420-427.

WU Xiaoping, XU Guoming. Study on 3-D Resistivity inversion using conjugate gradient method[J]. Chinese J Geophys, 2000, 43(3): 420-427.

[36] 王彦飞, 杨长春, 段秋梁. 地震偏移反演成像的迭代正则化方法研究[J]. 地球物理学报, 2009, 52(6): 1615-1624.

WANG Yanfei, YANG Changchun, DUAN Qiuliang. On iterative regularization methods for migration deconvolution and inversion in seismic imaging[J]. Chinese J Geophys, 2009, 52(6): 1615-1624.

[37] 陈龙伟, 徐世浙, 胡小平, 等. 位场向下延拓的迭代最小二乘法[J]. 地球物理学进展, 2011, 26(3): 894-901.

CHEN Longwei, XU Shizhe, HU Xiaoping, et al. The iterative Least Square method for downward continuation of potential fields[J]. Progress in Geophysics, 2011, 26(3): 894-901.

[38] Groetsch C W. The theory of tikhonov regularization for fredholm equations of the first kind[M]. Pibalan: Advanced Publising Program, 1984: 65-66.

[39] Engl H W, Hanke M, Neubauer A. Regularization of inverse problems[M]. Dordrecht: Kluwer, 1996.

[40] 李晓梅, 吴建平. Krylov子空间方法及其并行计算[J]. 计算机科学, 2005, 32(1): 19-20, 40.

LI Xiaomei, WU Jianping. Krylov subspace methods and parallel computation[J]. Computer Science, 2005, 32(1): 19-20, 40.

[41] Saad Y. Preconditioning techniques for nonsymmetric and indefinite linear systems[J]. J Comput Appl Matb, 1998, 24: 89-105.

[42] Saylor P E, Smolarski D C. Implementation of an adaptive algorithm for Richardsion’s method[M]. Lin Alg Appl: Springer, 1991: 154-156, 615-646.

[43] Nacbtigal N M, Reddy S C, Trefethen L N. How fast are nonsymmetric matrix derations[J]. J Matrix Anal Appl, 1992, 13: 778-795.

[44] 李春光, 徐成贤. 一个求解大型线性方程组的自适应CGNR算法[J]. 工程数学学报, 2001, 18(3): 71-77.

LI Chunguang, XU Chengxian. An adaptive CGNR algorithm for solving large linear systems[J]. Journal of Engineering Mathematics, 2001, 18(3): 71-77.

[45] 管志宁. 地磁场与磁力勘探[M]. 北京: 地质出版社, 2005: 183-184.

GUAN Zhining. Geomagnetic field and magnetic exploration[M]. Beijing: Geological Publishing House, 2005: 183-184.

(编辑 邓履翔)

收稿日期：2012-06-40；修回日期：2012-09-25

基金项目：国家自然科学基金资助项目(40774059，40874061)

通信作者：张志厚(1983-)，男，宁夏盐池人，博士，从事非震地球物理研究；电话：18683938811；E-mail：logicprimer@163.com