地球物理学报  2014, Vol. 57 Issue (6): 1958-1967   PDF    
位场向下延拓的改进迭代维纳滤波法
曾小牛, 刘代志, 李夕海, 牛超, 杨晓君, 卢世坤    
第二炮兵工程大学, 西安 710025
摘要:根据维纳滤波理论导出的位场向下延拓滤波器为最佳下延滤波器,但因其实现需要已知待求位场和噪声的功率谱而在实际应用中受到限制.针对该问题,本文首先提出一种基于位场径向平均功率谱的位场噪声水平估计方法,进而利用偏差准则求取正则化参数,实现位场正则化向下延拓;然后将位场正则化下延结果的功率谱作为待求位场功率谱的估计初值,采用带修正项的迭代维纳滤波方法来更新对待求位场功率谱的估计,最后提出本文的位场向下延拓改进迭代维纳滤波方法.基于理论重力模型数据及航磁实测数据进行了向下延拓对比试验,结果表明,改进迭代法具有较好的收敛性,且下延精度优于Tikhonov正则化法和递增型维纳滤波法.
关键词位场     向下延拓     迭代维纳滤波     Tikhonov正则化     径向平均功率谱     偏差准则    
An improved iterative Wiener filter for downward continuation of potential fields
ZENG Xiao-Niu, LIU Dai-Zhi, LI Xi-Hai, NIU Chao, YANG Xiao-Jun, LU Shi-Kun    
Xi'an Research Institute of High Technology, Xi'an 710025, China
Abstract: The downward continuation filter deduced by Wiener filter theory is an optimum downward continuation. However, the realization of the Wiener filter is difficult because it needs the power spectrum of the noise and the unknown potential field. In this paper, we propose an improved iterative Wiener filter for downward continuation of potential fields. Firstly, a noise level estimation method is proposed based on the radially averaged power spectrum of the potential field. Then, we use the discrepancy principle to choose the regularization parameter. Next, the regularization downward continuation result is used as the initial value of the iterative method. At last, an iterative Wiener filter with a correction term is used to update the power spectrum estimation of the unknown potential field. The comparison analysis of a theoretical gravity model and real aeromagnetic data shows that the proposed iterative method yields a better downward continuation than the Tikhonov regularization method and the incremental Wiener filter.
Key words: Potential field     Downward continuation     Iterative Wiener filter     Tikhonov regularization     Radially averaged power spectrum     Discrepancy principle    

1 引言

地磁导航(胡小平和吴美平,2013Wu et al., 2013)和重力导航(Wang et al., 2012童余德等,2012)是近几年无源导航领域的研究热点.地磁/重力导航的实现,其中一个关键问题是不同高度高精度地磁/重力数据库的构建.经典场论中的重磁位场解析延拓理论,为解决该问题提高了途径(陈龙伟,2013).位场延拓方法中,又因向下延拓是一个典型的不适定线性反问题而一直是高精度重磁物探领域的研究难点和热点(张志厚,2013).对该不适定线性反问题的求解通常有两种策略:一是直接法;二是迭代法.直接法中最常用的是维纳滤 波法(Clarke,1969; Pawlowski,1995)和Tikhonov 正 则化法(梁锦文,1989; Pašteka et al., 2012; Abedi et al., 2013; Li et al., 2013)(或称广义逆法(陈生昌和肖鹏飞,2007)),虽然这两种方法的理论依据不同,但是它们却有很强的联系,即令维纳滤波中的噪声和信号功率谱之比为常数则可得Tikhonov正则化滤波算子.迭代法的种类很多,如积分迭代法(徐世浙,2006Xu et al., 2007)及其相关改进迭代法(张辉等,2009刘东甲等,2009曾小牛等, 2011b2013高玉文等,2012)、结合泰勒级数和位场导数的迭代法(Fedi and Florio, 2002; 王彦国等, 20112012Zhang et al., 2013; Ma et al., 2013)、等效源法(Dampancy,1969; Hansen and Miyazaki, 1984)、迭代正则化方法(Zhou et al., 2011; Zeng et al., 2013; 张志厚等,2013)、样条函数法(汪炳柱和王硕儒,1998Wang,2006)等.尽管迭代法并不能从根本上解决下延计算的不稳定性(姚长利等,2012李宏伟,2012),但根据文献(Zhou et al., 2011曾小牛等,2011bZeng et al., 2013Zhang et al., 2013)对一些迭代法和直接法进行的对比分析可知,一些改进的迭代法可以获得比直接法更好的向下延拓效果.

基于维纳滤波的向下延拓算子又被称为最佳下延算子(Clarke,1969),但它的具体实现受制于其需要已知噪声和待求位场的功率谱比.在实际应用中,如果将该功率谱比近似处理成一个常数就转化成了Tikhonov正则化法.Tikhonov正则化法下延结果稳定却过于平滑,究其原因是,正则化延拓算子在压制噪声的同时,同样压制了信号中的有用成份.张辉等(2009)采用递增型维纳滤波的方法更新对待求位场功率谱的估计,但对正则化参数和噪声水平并没有提出计算算法,而是令这两个参数都等于1.Li等(2013)采用L-曲线法求取正则化参数,并利用格林等效层的功率谱形式表示待求位场的功率谱,然后采用概率算法中的模拟退火方法来求解功率谱形式中的各个参数,最终实现了对待求位场功率谱的估计.

本文首先从位场径向平均功率谱(Spector and Grant, 1970)的物理特性出发,提出一种估计噪声水平的方法,并在此基础上基于后验策略的偏差准则(Morozov,1984)选取正则化参数,实现位场正则化向下延拓;然后,以正则化下延结果作为待求位场功率谱的估计初值,并采用带修正项的迭代维纳滤波方法来更新对待求位场功率谱的估计,进而提出本文的位场向下延拓改进迭代维纳滤波方法.最后,基于理论模型及实测数据的对比试验验证了改进迭代法的有效性.

2 位场维纳滤波向下延拓原理

若z=z1观测平面上的位场g(x,y)为已知,设z轴向下为正,则由g(x,y)向下延拓求z=z2(z2>z1,场源位于z=z2平面以下)平面上的位 场数据f(x,y)可表示成如下的第一类卷积型方程为

其中,为卷积符号;h(x,y)= (h=z2-z1>0 为延拓高度)为卷积核函数;n(x,y)为加性噪声,一般为高斯白噪声,设其方差为σ2.

对(1)式两边作Fourier变换得

其中,G(ωx,ωy)、H(ωx,ωy)、F(ωx,ωy)和N(ωx,ωy)分别为g(x,y)、h(x,y)、f(x,y)和n(x,y)的频域形式.维纳滤波亦称最小均方误差滤波,其对应的下延算子的频域形式如下(Clarke,1969)

其中,Pgx,ωy)、Pfx,ωy)和Pnx,ωy)分别为位场g(x,y)、f(x,y)和噪声n(x,y)对应的功率谱.采用经典的周期图法估计功率谱,则维纳滤波又可表示为

其中,*表示共轭.由于在实际应用中,σ2和F(ωx,ωy)一般是未知的,故在传统的维纳滤波法中通常将(4)式中的σ2/ |F(ωx,ωy)| 2用某个常量α来近似代替(亦被称为“参数化维纳滤波(Gonzalez et al., 2005)”),即

显然,这与Tikhonov正则化下延算子是一致的(梁锦文,1989陈生昌和肖鹏飞,2007).可见,Tikhonov正则化可视为维纳滤波的一种近视,只是在Tikhonov正则化中,α被称为正则化参数.

3 迭代维纳滤波原理及其改进 3.1 迭代维纳滤波原理

在参数化维纳滤波中,用常量α来代替σ2/ |F(ωx,ωy)| 2,虽然简单,但结果却使下延结果趋于平滑.如果能从已有的信息中提取或者估计出|F(ωx,ωy)|2,将会提高下延精度.由于实际应用中仅位场g(x,y)已知,基于此,我们可以通过计算|G(ωx,ωy)|2来大致估计|F(ωx,ωy)|2.同时,为了提高对|F(ωx,ωy)|2估计的准确度,可以考虑用迭代的方法来逼近其真实值,迭代维纳滤波即是基于此思想(Hillery and Chin, 1991).该迭代法的迭代过程可归纳如下:

设置初值为

迭代过程为

Hillery和Chin(1991)进一步指出(7)式不能收敛到理想的最小均方误差(除非误差n(x,y)为零),并提出添加修改项

其中,

式(9)为修改项,从上述迭代过程可知,虽然该迭代法解决了待求位场的功率谱估计问题,但依然存在需要估计位场噪声方差的问题.基于此,本文提出如下的噪声方差估计方法.

3.2 基于径向谱的噪声水平估计

Spector和Grant(1970)提出径向平均功率谱的概念,并被广泛应用于位场数据处理(Stefan and Vijay, 1995; Ravat et al., 2007; 郭良辉等,2012Guo et al., 2013).径向谱的定义如下:设位场数据的功率谱为P(ωx,ωy),其在径向方向φ上的取值为P′(ωr,φ),则

其中,ωr=为径向波数,φ为径向与ωx轴的夹角.当φ取定值时,P′(ωr,φ)仅是ωr的函数.求取径向平均功率谱的具体实施步骤为:(1)求位场功率谱;(2)以位场功率谱的中心作圆,圆的半径分别为基频的整数倍,从而形成一系列环带;(3)取环带内所有点的功率谱值作其平均值;(4)以径向波数ωr为横坐标,各频带内功率谱平均值作为纵坐标(取对数坐标),即可得径向平均功率谱.

一个假想的径向平均功率谱如图 1所示.从图 1可知,如果假设观测噪声为白噪声,则与位场信号不相关的白噪声的功率谱应该为常数,即对应于图 1中的水平部分.这样,对于位场的径向平均功率谱,存在一个截止频率ωc将位场信号谱和噪声谱大致分开.由此可以得到已知位场功率谱的等效噪声区域如图 2所示(仅显示四分之一区域).由以上讨 论可知,利用Parseval定理(Mitra and Kuo, 2006),我们可以通过如下式来近似估计噪声的方差为

图 1 径向平均功率谱示意图Fig. 1 Schematic diagram of the radially averaged power spectrum

图 2 1/4部分等效噪声区域示意图Fig. 2 Schematic diagram of one quadrant of the equivalent noise region
其中,gNoisex,ωy)表示等效噪声区域内相应位场g(x,y)的Fourier变换结果,Q为等效噪声区域内的数据总数值.

3.3 改进迭代维纳滤波原理

在实验中我们发现,当已知位场数据噪声干扰较大或向下延拓深度较大时,以已知位场作为迭代法迭代初值,3.1节中的迭代维纳滤波并不能获得理想的向下延拓结果,而且收敛缓慢,原因即出在迭代初值的设定上.因此,本文考虑用正则化向下延拓值作为迭代维纳滤波法的迭代初值.

正则化向下延拓实现的关键在于正则化参数的选取.正则化参数的选取,依据是否需要预先估计原始数据的噪声水平而分为先验和后验两类策略.L-曲线法(Abedi et al., 2013; Li et al., 2013; Zeng et al., 2013; 马涛等,2013)和拟最优法(或称C-范数法)(Pašteka et al., 20062012; Zhang et al., 2013)是最常用的先验策略;Morozov偏差原理(Morozov,1984)是典型的后验策略.以前的延拓文献中,通常都采用先验策略,未见有采用先验策略的例子.本文在上一节估计得到位场噪声方差的基础之上,采用偏差原理来实现位场向下延拓.

偏差原理的基本思想是,如果已知噪声的方差σ2,则最佳正则化参数α应该满足以下约束方程

其中,fα(x,y)为正则化参数取α时的正则化向下延拓结果,M为位场总的数据个数.式(12)的物理意义是显然的,即延拓误差应该同原始位场的噪声水平一致.在实际应用中,可构建如下的偏差函数为

式(13)可求解最佳正则化参数,值得指出的是,在噪声方差可以获取或近似得到的情况下,偏差准则是十分有效的正则化参数选择方法(王彦飞,2007).

综上所述,本文改进迭代维纳滤波方法的算法步骤如下:

(1)采用本文提出的基于径向平均功率谱的方法(式(11))估计出原始位场数据的噪声方差σ2

(2)采用偏差原理(式(13))求得正则化参数,进而求得位场的正则化向下延拓值;

(3)以第二步求得的正则化下延值作为迭代维纳滤波的迭代初值,即 F1(u,v)= H*x,ωy)H(ωx,ωy)2+α G(ωx,ωy);

(4)以迭代维纳滤波迭代式(7)式进行迭代,并在每次迭代过程中采用(8)式和(9)式对待求位场功率谱的估计进行修正.

4 数值试验及结果比较 4.1 理论模型试验

理论模型由两个不同深度层、不同大小的7个球体组合而成,图 3a显示了它们在直角坐标系(设z坐标向下为正)中的位置.各球体的具体参数见表 1.设点、线数同为512,点、线距都为50 m.h=0 m高度的重力异常如图 3b所示.向下延拓实验为将h=-1000 m高度的重力异常向下延拓到h=0 m高度.同时,为了模拟实际情况,在h=-1000 m高度重力异常中加入均值为零、方差为0.0011 mGal的白噪声,其等值线如图 3c所示.

图 3 理论模型重力数据向下延拓效果对比(a)球体位置坐标系;(b)h=0 m高度重力数据等值线图;(c)h=-1000 m高度含噪重力数据等值线图及其(d)径向对数功率谱;(e)求正则化参数的偏差原理函数曲线图;递增维纳滤波法和本文迭代法下延均方误差(f)和相对误差(g)随迭代次数变化的关系; 各方法的下延结果:(h)Tikhonov正则化法;(i)递增型维纳滤波法(迭代50次);(j)本文迭代法(迭代10次).Fig. 3 Effects comparison of downward continuation based on the model gravity data(a)Location of spheres coordinates system;(b)Contours of gravity data at h=0 m;(c)Contours of gravity data with noise at h=-1000 m and (d)its radially averaged power spectrum;(e)The function graph of the discrepancy principle for the choice of the regularization parameter;(f)Root mean square errors and (g)relative errors varying with the iteration number n for the incremental type Wiener filter method and the proposed iterative method in this paper; The downward continuation results for(h)Tikhonov regularization method,(i)the incremental type Wiener filter(50 iteration numbers), and (j)the proposed iterative method(10 iteration numbers).

表 1 仿真球体所用参数Table 1 Parameters of simulation spheres

首先,计算该模型h=-1000 m高度含噪重力异常数据的径向平均功率谱,如图 3d所示.通过功率谱形状分析,显然可将功率谱分为0~0.0008和0.0008~0.01两段,将前段视为信号部分,后一段近似视为噪声部分.如此,则可以确定功率谱的截止频率为ωc=0.0008.基于此,按公式(11)求得的噪声方差为0.001 14 mGal.求得噪声方差以后,再利用(13)式求取最优正则化参数,如图 3e所示,求得的正则化参数为0.0006.分别采用递增维纳滤波法和本文提出的方法将h=-1000 m平面处的含噪重力数据向下延拓1000 m(20倍点距)到h=0 m 高度,并用得到的延拓值fc与该高度处的真实重力 场值ft采用均方误差(Root Mean Square Error,RMSE)

和相对误差(Relative Error,RE)

来计算延拓误差.递增型维纳滤波法和本文方法的下延均方误差和相对误差随迭代次数n的变化关系分别如图 3f、3g所示(显然,依据本文改进迭代法的迭代步骤,本文迭代法在迭代次数n=1时对应Tikhonov法).Tikhonov方法、递增型维纳滤波方法、本文方法的延拓结果分别如图 3h图 3i图 3g所示.图 4a图 4b图 4c分别为三种方法下延结果与理论值的残差对比结果图.由图 3f—3j和图 4a—4c的对比可见:
图 4 各方法理论模型重力数据下延结果与理论值的残差对比(a)Tikhonov正则化法;(b)递增型维纳滤波法;(c)本文迭代法.Fig. 4 The difference between the downward continuation results and the theoretical value for the model gravity data(a)Tikhonov regularization method;(b)The incremental type Wiener filter ;(c)The proposed iterative method.

(1)本文提出的噪声水平估计方法是非常有效的,估计的噪声方差与真实值非常接近,其误差为3.64%,误差产生的原因应该是位场中的小部分高频分量和高频噪声混在一起,截断时被当作噪声,导致估计的噪声方差偏大;

(2)偏差函数φ(α)为α的单调函数(图 3e),且由于估计得到的噪声方差比真实值稍大,导致利用φ(α)得到的正则化参数较最优正则化参数也稍大,幸运的是,本文改进迭代法中,仅将正则化下延值作为迭代法的迭代初值,因此,该误差并不会影响迭代法的延拓精度;

(3)递增型维纳滤波法初始延拓误差很大,随着迭代的进行,虽然延拓误差迅速减少,但最终并不能得到相对较低的延拓误差,其本质原因在于不能有效估计噪声水平和正则化参数,而是简单地设噪声方差和正则化参数都为1;

(4)由于噪声水平估计较精确,且以Tikhonov正则化向下延拓结果作为迭代法迭代初值,本文改进迭代法的初始延拓误差较小,且延拓误差随着迭代次数的增加依然减少,能够得到较低的延拓误差.

本例中,Tikhonov法向下延拓的均方误差和相对误 差为0.2379 mGal、8.30%,递增型维纳滤波方法(迭代50次)向下延拓的均方误差和相对误差为 0.2181 mGal、6.37%,而本文所提出方法(迭代10次)向下延拓的均方误差和相对误差为0.1829 mGal、4.94%.

4.2 实测数据试验

航磁实测数据点距和线距均为100 m,其等值线图见图 5a.将原始数据向上延拓1000 m(10倍点距)以构建另一个平面的数据(图 5b);同时,考虑到向上延拓具有的平滑效应,为模拟实际情况,本文在向上延拓后得到的航磁数据中加入零均值、方差为23 nT的高斯白噪声,加噪后航磁数据等值线图如图 5c所示.计算该含噪航磁异常数据的径向平均功率谱,如图 5d所示.通过功率谱形状分析,可以确定功率谱的截止频率为ωc=0.0007.按公式(11)求得的噪声方差为23.1 nT.求得噪声方差以后,再利用(13)式求取最优正则化参数,如图 5e所示,求得的正则化参数为α=0.0011.

图 5 实测航磁数据向下延拓效果对比(a)原始航磁等值线图(nT);(b)向上延拓20倍点距后等值线图;(c)向上延拓航磁异常加噪数据等值线图及其(d)径向对数功率谱;(e)求正则化参数的偏差原理函数曲线图;递增维纳滤波法的下延均方误差(f)和相对误差(g)随迭代次数变化的关系;本文迭代法下延均方误差(h)和相对误差(i)随迭代次数变化的关系; 各方法的下延结果:(j)Tikhonov正则化法;(k)递增型维纳滤波法(迭代2次);(l)本文迭代法(迭代1次).Fig. 5 Effects comparison of downward continuation based on the real aeromagnetic data(a)Contours of the original aeromagnetic data(in nT);(b)Upward continuation 20 grid interval;(c)Contours of upward continuation aeromagnetic data with noise and (d)its radially averaged power spectrum;(e)The function graph of the discrepancy principle for the choice of the regularization parameter;(f)Root mean square errors and (g)relative errors varying with the iteration number n for the incremental type Wiener filter method;(h)Root mean square errors and (i)relative errors varying with the iteration number n for the proposed iterative method in this paper;The downward continuation results for(j)Tikhonov regularization method,(k)the incremental type Wiener filter(2 iteration numbers), and (l)the proposed iterative method(1 iteration numbers).

分别采用递增型维纳滤波法和本文提出的迭代法将向上延拓后加噪的航磁数据向下延拓1000 m到原来的高度.递增型维纳滤波法的延拓均方误差和相对误差随迭代次数n的变化关系分别如图 5f、5g所示.本文方法的延拓均方误差和相对误差随迭代次数n的变化关系分别如图 5h、5i所示(本文方 法在迭代次数n=1时对应Tikhonov法).Tikhonov方法、 递增型维纳滤波方法和本文方法的延拓结果分别如图 5j图 5k图 5l所示.图 6a图 6b图 6c分别为三种方法下延结果与理论值的残差对比结果图.由图 5f—5l和图 6a—6c的对比可见:

图 6 各方法实测航磁数据下延结果与理论值的残差对比(a)Tikhonov正则化法;(b)递增型维纳滤波法;(c)本文迭代法.Fig. 6 The difference between the downward continuation results and the theoretical value for the real aeromagnetic data(a)Tikhonov regularization method;(b)The incremental type Wiener filter;(c)The proposed iterative method.

(1)递增型维纳滤波法初始延拓误差很大,且只在第二步迭代得到了尚能接受的延拓误差,而后延拓误差急剧上升(图 5f、5g),究其原因是,实测数据的噪声方差远大于1,因此,简单的设噪声方差为1,会导致迭代法的严重半收敛性;

(2)由于噪声水平估计较精确,且以Tikhonov正则化向下延拓结果作为迭代法迭代初值,本文改进迭代法的初始延拓误差较小,且延拓误差随着迭代次数的增加依然减少,能够得到较低的延拓误差,迭代法的收敛性较好.本例中,Tikhonov法向下延拓的均方误差和相对误差为39.75 mGal、9.91%,递增型维纳滤波方法最好的(迭代2次)向下延拓的均方误差和相对误差为38.87 mGal、9.42%,本文方法(迭代10次)向下延拓的均方误差和相对误差为36.91 mGal、8.52%,可见本文方法有一定的优越性.

5 结论

(1)尽管迭代法不能从根本上消除位场向下延 拓的不稳定性,但是一些改进迭代法可以获得比直 接Tikhonov正则化法更好的延拓效果,因此,探讨一些改进的迭代法用于提高位场向下延拓的精度是有意义的.

(2)维纳滤波为最佳滤波,但其具体实现却存在困难;Tikhonov正则化法可视为维纳滤波的近似,但也正因为这种近似而导致其下延结果过于平滑.

(3)本文提出的位场噪声估计方法是基于位场功率谱的物理特性,虽然在估计噪声水平的过程中,因为位场高频信号和高频噪声不可能用一个截止频率完全地分开而导致估计到的噪声水平偏大,但可以应用到位场数据处理其它需要大致估计位场噪声水平的领域.

(4)递增型迭代维纳滤波的延拓误差和收敛性都与位场的噪声水平直接相关,由于设定正则化参数和噪声方差都为1,当噪声方差小于1时,该迭代法收敛性较好,但当噪声方差大于1时,该迭代法只能在个别迭代步上获得低于Tikhonov正则化法的延拓误差,此时,迭代法存在严重半收敛性.

(5)本文提出的改进迭代法由于噪声水平估计较精确,且以Tikhonov正则化向下延拓结果作为迭代法迭代初值,该迭代法的初始延拓误差较小,且延拓误差随着迭代次数的增加依然减少,能够得到较低的延拓误差,迭代法的收敛性较好,非常适合应用到实际位场数据向下延拓处理中.

参考文献
[1] Abedi M, Gholami A, Norouzi G H. 2013. A stable downward continuation of airborne magnetic data: A case study for mineral prospectivity mapping in Central Iran. Computers & Geosciences, 52: 269-280.
[2] Chen L W. 2013. Research on Spatial Domain Numerical Methods for Potential Field Continuation in Geomagmetic Navigation (in Chinese). Changsha: National University of Defense Technology.
[3] Chen S C, Xiao P F. 2007. Wavenumber domain generalize inverse algorithm for potential field downward continuation. Chinese J. Geophys. (in Chinese), 50(6): 1816-1822.
[4] Clarke G K C. 1969. Optimum second-derivative and downward-continuation filters. Geophysics, 34(3): 424-437.
[5] Dampancy C N G. 1969. The equivalent source technique. Geophysics, 34(1): 39-53.
[6] Fedi M, Florio G. 2002. A stable downward continuation by using the ISVD method. Geophysical Journal International, 151(1): 146-156.
[7] Gao Y W, Luo Y, Wen W. 2012. The compensation method for downward continuation of potential field from horizontal plane and its application. Chinese J. Geophys. (in Chinese), 55(8): 2747-2756,doi:10.6038/j.issn.0001-5733.2012.08.026.
[8] Gonzalez R C, Woods R E, Eddins S L. 2005. Digital Image Processing Using MATLAB. Beijing: Electronics Industry Press.
[9] Guo L H, Meng X H, Shi L, et al. 2012. Preferential filtering method and its application to Bouguer gravity anomaly of Chinese continent. Chinese J. Geophys. (in Chinese), 55(12): 4078-4088,doi:10.6038/j.issn.0001-5733.2012.12.020.
[10] Guo L H, Meng X H, Chen Z X, et al. 2013. Preferential filtering for gravity anomaly separation. Computers & Geosciences, 51: 247-254.
[11] Hansen R O, Miyazaki Y. 1984. Continuation of potential fields between arbitrary surfaces. Geophysics, 49(6): 787-795.
[12] Hillery A D, Chin R T. 1991. Iterative Wiener filters for image restoration. IEEE Transactions on Signal Processing, 39(8): 1892-1899.
[13] Hu X P, Wu M P. 2013. Technologies on Underwater Geomagnetic Field Navigation (in Chinese). Beijing: National Defense Industry Press.
[14] Li H W. 2012. Research on Iteration Method Using in Potential Field Transformations and Its Filter Characteristic (in Chinese). Beijing: China University of Geosciences (Beijing).
[15] Li Y G, Devriese S G R, Krahenbuhl R A, et al. 2013. Enhancement of magnetic data by stable downward continuation for UXO application. IEEE Transactions on Geoscience and Remote Sensing, 51(6): 3605-3614.
[16] Liang J W. 1989. The regularization method for downward continuation of Potential field. Chinese J. Geophys. (in Chinese), 32(5): 600-608.
[17] Liu D J, Hong T Q, Jia Z H, et al. 2009. Wave number domain iteration method for downward of potential fields and its convergence. Chinese J. Geophys. (in Chinese), 52(6): 1599-1605.
[18] Ma G Q, Liu C, Huang D N, et al. 2013. A stable iterative downward continuation of potential field data. Journal of Applied Geophysics, 98: 205-211.
[19] Ma T, Chen L W, Wu M P, et al. 2013. The selection of regularization parameter in downward continuation of potential field based on L-curve method. Progress in Geophys. (in Chinese), 28(5): 2485-2494,doi:10.6038/pg20130527.
[20] Mitra S K, Kuo Y. 2006. Digital Signal Processing: A Computer-based Approach. New York: McGraw-Hill.
[21] Morozov V A. 1984. Methods for Solving Incorrectly Posed Problems. New York: Springer-Verlag.
[22] Pašteka R, Richter F P, Karcol R, et al. 2006. Regularized derivatives of potential fields and their role in semi-automated interpretation methods. Geophysical Prospecting, 57(4): 507-516.
[23] Pašteka R, Karcol R, Kušnirák D, et al. 2012. REGCONT: A MATLAB based program for stable downward continuation of geophysical potential fields using Tikhonov regularization. Computers & Geosciences, 49: 278-289.
[24] Pawlowski R S. 1995. Preferential continuation for potential-field anomaly enhancement. Geophysics, 60(2): 390-398.
[25] Ravat D, Pignatelli A, Nicolos I, et al. 2007. A study of spectral methods of estimating the depth to the bottom of magnetic sources from near-surface magnetic anomaly data. Geophysical Journal International, 169(2): 421-434.
[26] Spector A, Grant F S. 1970. Statistical models for interpreting aeromagnetic data. Geophysics, 35(2): 293-302.
[27] Stefan M, Vijay D. 1995. Potential field power spectrum inversion for scaling geology. Journal of Geophysical Research, 100(B7): 12605-12616.
[28] Tong Y D, Bian S F, Jiang D F, et al. 2012. A new integrated gravity matching algorithm based on approximated local gravity map. Chinese J. Geophys. (in Chinese), 55(9): 2917-2924,doi:10.6038/j.issn.0001-5733.2012.09.011.
[29] Wang B Z, Wang S R. 1998. Spline function methods for upward continuation and downward continuation of 2D potential field. Computing Techniques for Geophysical and Geochemical Exploration (in Chinese), 20(2): 125-129
[30] Wang B Z. 2006. 2D and 3D potential-field upward continuation using splines. Geophysical Prospecting, 54(2): 199-209.
[31] Wang H B, Wang Y, Fang J, et al. 2012. Simulation research on a minimum root-mean-square error rotation-fitting algorithm for gravity matching navigation. Science China: Earth Sciences, 55(1): 90-97.
[32] Wang Y F. 2007. Computational Methods for Inverse Problems and Their Applications (in Chinese). Beijing: Higher Education Press.
[33] Wang Y G, Zhang F X, Wang Z W, et al. 2011. Taylor series iteration for downward continuation of potential field. OGP (in Chinese), 46(4): 657-662.
[34] Wang Y G, Wang Z W, Zhang F X, et al. 2012. Derivative-Iteration method for downward continuation of potential fields. Journal of Jilin University: Earth Science Edition (in Chinese),42(1): 240-245.
[35] Wu Z T, Hu X P, Wu M P, et al. 2013. An experimental evaluation of autonomous underwater vehicle localization on geomagnetic map. Applied Physics Letters, 103(10): 104102.
[36] Xu D X. 2005. Using gravity anomaly matching techniques to implement submarine navigation. Chinese J. Geophys. (in Chinese), 48(4): 812-816.
[37] Xu S Z. 2006. The integral-iteration method for continuation of potential field. Chinese J. Geophys. (in Chinese), 49(4): 1176-1182.
[38] Xu S Z, Yang J Y, Yang C F, et al. 2007. The iteration method for downward continuation of a potential field from a horizontal plane. Geophsical Prospecting, 55(6): 883-889.
[39] Yao C L, Li H W, Zheng Y M, et al. 2012. Research on iteration method using in potential field transformations. Chinese J. Geophys. (in Chinese), 55(6): 2062-2078,doi:10.6038/j.issn.0001-5733.2012.06.028.
[40] Zeng X N, Li X H, Han S Q, et al. 2011a. A comparison of three iteration methods for downward continuation of potential fields. Progress in Geophys. (in Chinese), 26(3): 908-915,doi:10.3969/j.issn.1004-2903.2011.03.016.
[41] Zeng X N, Li X H, Liu D Z, et al. 2011b. Regularization analysis of integral iteration method and the choice of its optimal step-length. Chinese J. Geophys. (in Chinese), 54(11): 2943-2950,doi:10.3969/j.issn.0001-5733.2011.11.024.
[42] Zeng X N, Li X H, Niu C, et al. 2013. The regularization-integral iteration method in wave number domain for downward continuation of potential fields. OGP (in Chinese), 48(4): 643-650.
[43] Zeng X N, Li X H, Su J, et al. 2013. An adaptive iterative method for downward continuation of potential-field data from a horizontal plane. Geophysics, 78(4): J43-J52.
[44] Zhang H, Chen L W, Ren Z X, et al. 2009. Analysis on convergence of iteration method for potential fields downward continuation and research on robust downward continuation method. Chinese J. Geophys. (in Chinese), 52(4): 1107-1113.
[45] Zhang H L, Ravat D, Hu X Y. 2013. An improved and stable downward continuation of potential field data: The truncated Taylor series iterative downward continuation method. Geophysics, 78(5): J75-J86.
[46] Zhang Z H, Wu L Y, Wang R S, et al. 2013. CGNR method for potential field downward continuation. Journal of Central South University (Science and Technology) (in Chinese), 44(8): 3273-3281.
[47] Zhang Z H. 2013. Numerical Computation Method of Downward Continuation of Potential Fields. Hangzhou: Zhejiang University.
[48] Zhou J, Shi G G, Ge Z L. 2011. Study of iterative regularization methods for potential field downward continuation in geophysical navigation. Journal of Astronautics, 32(4): 787-794.
[49] 陈龙伟. 2013. 面向地磁导航的位场延拓空间域算法研究[博士论文]. 长沙: 国防科学技术大学.
[50] 陈生昌, 肖鹏飞. 2007. 位场向下延拓的波数域广义逆算法. 地球物理学报, 50(6): 1816-1822.
[51] 高玉文, 骆遥, 文武. 2012. 补偿向下延拓方法研究及应用. 地球物理学报, 55(8): 2747-2756,doi:10.6038/j.issn.0001-5733.2012.08.026.
[52] 郭良辉, 孟小红, 石磊等. 2012. 优化滤波方法及其在中国大陆布格重力异常数据处理中的应用. 地球物理学报, 55(12): 4078-4088,doi:10.6038/j.issn.0001-5733.2012.12.020.
[53] 胡小平, 吴美平. 2013. 水下地磁导航技术. 北京: 国防工业出版社.
[54] 李宏伟. 2012. 位场转换计算中迭代法与其滤波特性的综合分析与研究[博士论文]. 北京: 中国地质大学(北京).
[55] 梁锦文. 1989. 位场向下延拓的正则化方法. 地球物理学报, 32(5): 600-608.
[56] 刘东甲, 洪天求, 贾志海等. 2009. 位场向下延拓的波数域迭代法及其收敛性. 地球物理学报, 52(6): 1599-1605.
[57] 马涛, 陈龙伟, 吴美平等. 2013. 基于L曲线法的位场向下延拓正则化参数选择. 地球物理学进展, 28(5): 2485-2494,doi:10.6038/pg20130527.
[58] 童余德, 边少锋, 蒋东方等. 2012. 一种新的基于局部重力图逼近的组合匹配算法. 地球物理学报, 55(9): 2917-2924,doi:10.6038/j.issn.0001-5733.2012.09.011.
[59] 汪炳柱, 王硕儒. 1998. 二维位场向上延拓和向下延拓的样条函数法. 物探化探计算技术, 20(2): 125-129.
[60] 王彦飞. 2007. 反演问题的计算方法及其应用. 北京: 高等教育出版社.
[61] 王彦国, 张凤旭, 王祝文等. 2011. 位场向下延拓的泰勒级数迭代法. 石油地球物理勘探, 46(4): 657-662.
[62] 王彦国, 王祝文, 张凤旭等. 2012. 位场向下延拓的导数迭代法. 吉林大学学报(地球科学版), 42(1): 240-245.
[63] 许大欣. 2005. 利用重力异常匹配技术实现潜艇导航. 地球物理学报, 48(4): 812-816.
[64] 徐世浙. 2006. 位场延拓的积分迭代法. 地球物理学报, 49(4): 1176-1182.
[65] 姚长利, 李宏伟, 郑元满等. 2012. 重磁位场转换计算中迭代法的综合分析与研究. 地球物理学报, 55(6): 2062-2078,doi:10.6038/j.issn.0001-5733.2012.06.028.
[66] 曾小牛, 李夕海, 韩绍卿等. 2011a. 位场向下延拓三种迭代方法之比较. 地球物理学进展, 26(3): 908-915,doi:10.3969/j.issn.1004-2903.2011.03.016.
[67] 曾小牛, 李夕海, 刘代志等. 2011b. 积分迭代法的正则性分析及其最优步长的选择. 地球物理学报, 54(11): 2943-2950,doi:10.3969/j.issn.0001-5733.2011.11.024.
[68] 曾小牛, 李夕海, 牛超等. 2013. 位场向下延拓的波数域正则-积分迭代法. 石油地球物理勘探, 48(4): 643-650.
[69] 张辉, 陈龙伟, 任治新等. 2009. 位场向下延拓迭代法收敛性分析及稳健向下延拓方法研究. 地球物理学报, 52(4): 1107-1113.
[70] 张志厚, 吴乐园, 王瑞赛等. 2013. 位场向下延拓的CGNR法. 中南大学学报(自然科学版), 44(8): 3273-3281.
[71] 张志厚. 2013. 位场向下延拓的数值计算方法[博士论文]. 杭州: 浙江大学.