朱占龙,杨功流,,杨淑洁

1.东南大学仪器科学与工程学院,江苏南京 210096;2.北京航空航天大学仪器科学与光电工程学院,北京 100191

位场向下延拓的加速Landweber迭代法

朱占龙1,杨功流1,2,杨淑洁2

1.东南大学仪器科学与工程学院,江苏南京 210096;2.北京航空航天大学仪器科学与光电工程学院,北京 100191

利用航空测量数据向下延拓得到不同高度的位场数据可以提高测量成果的综合利用率。Landweber迭代法是一种有效解决位场向下延拓的实用方法。鉴于Landweber迭代法的收敛速度比较慢,提出采用加速Landweber迭代法,推导得到两种迭代法对应的波数域算子并通过仿真分析算子的滤波特性,最后结合模型实例验证了所提出的加速Landweber迭代法不仅能有效进行位场向下延拓,而且比Landweber迭代法更高效。

位场;向下延拓;正则化方法;加速Landweber迭代法;收敛速度

1 引 言

随着磁传感器技术的进步[1-2]以及磁力测量系统的改进[3-5],航空磁力测量已经日趋成熟,尤其是在海平面、地表和超低空等测量条件恶劣的地区,航空磁力测量更是不可或缺。采用航磁数据进行不同高度的向下延拓可以弥补条件恶劣地区磁场资料的不足,因此,深入研究稳定的向下延拓方法具有重要的理论意义和应用价值。

向下延拓问题是一个将噪声数据放大的过程,属于病态问题[6-8]。若采用不合适的延拓方法,将无法得到稳定的解。传统的向下延拓方法是快速傅里叶变换法,但是利用该方法向下延拓得到的解极其不稳定[9-10],可靠的延拓深度一般不会超过2～3倍点距[11]。随着数学的研究发展,许多新理论和方法被成功地应用到此类病态问题的计算中来。文献[11]采用的迭代法能够稳定向下延拓20倍点距,在一定程度上解决了此类问题。近些年来,采用正则化方法[12]进行向下延拓成为该工程技术领域的热点,文献[13]考虑了正则化参数和地形对延拓结果的影响;文献[14]利用该方法研究了在航空重力测量数据向下延拓问题并分析了两种正则化方法的差异[15],正则化方法中的Landweber迭代法也在一定程度上解决了此类问题[16]。

以上方法均为向下延拓问题的解决提供了一定思路,但是工程化的应用需要计算更快速、效率更高效,因此本文着重研究向下延拓方法的运算效率问题,提出采用加速Landweber迭代法[17],结合模型数据分析Landweber迭代法和加速Landweber迭代法波数域算子的滤波特性,展示了加速Landweber迭代法在位场向下延拓中的计算效率与延拓精度。

2 位场向下延拓原理

如图1所示,水平面z＝0和z＝h之间是无源空间,设z轴向下为正,场源位于z＝h以下(h＞0)。观测面为水平面(z＝0),该平面上的位场f(x,0)为已知,计算观测面以下至场源以上平面的位场f(x,h)称为位场的向下延拓。

图1 向下延拓示意图Fig.1 Schematic diagram of downward continuation

根据位场向上延拓公式,得观测平面位场f(x,0)与向下延拓位场f(x,h)之间的关系[18]

观测平面上的位场数据fδx()往往是带有高频噪声的数据,当位场向下延拓时,由于随着的增加,延拓算子呈指数规律增加,所以求得的g x()将非常不稳定,这就是位场向下延拓不稳定的本质。由此可见,当高频噪声存在时,向下延拓视为一个不适定问题,对其直接作傅里叶变换以及反变换所得到的结果,不能作为向下延拓平面位场的稳定近似解。

3 位场向下延拓的迭代法及相应的波数域算子

同第一类Fredholm积分方程比较,式(2)可以表示为[19]

式中,K为Fredholm算子,对于形如式(5)的不适定问题求解一般采用正则化方法,其中Landweber迭代法是一种广泛应用的方法。

3.1 位场向下延拓迭代法

3.1.1 Landweber迭代法

Landweber迭代法是最速下降法的一种变形,其迭代式如下[20-21]

对于给定的初值g0x(),迭代式(6)可以展开为

继续按照式(7)中的步骤展开gnx(),式(7)可以改写为下面的级数形式

(3)计算gn(x)＝Tnfδ(x),如果满足迭代终止;若不满足,n＝n＋1,转向(2)。其中,C一般为固定常数,它决定了Landweber迭代算法是否具有收敛性[19]。很显然,在δ确定后,C值越大使得迭代终止条件Cδ越大,则迭代次数变少,但是相应的结果精度会降低。

由上述算法可以看出,主要决定算法收敛速度的是第(2)步,在实际计算中,可以先对(2)迭代计算若干步之后,再进行第(3)步判断终止条件,可以在一定程度上加快计算速度。

3.1.2 加速Landweber迭代法

通过对Landweber迭代法分析可知,如果足够大n(Tn的下标)能够快速计算出来,则迭代终止条件就能较易得到满足,从而使得计算效率提高。为此,对于Tn的计算,采用如下迭代形式[23]

从式(14)可以看出,用式(11)计算n次,相当于Landweber迭代法迭代2n－1次,这样可以使得计算效率明显提升。于是,加速Landweber迭代法步骤如下:

3.2 位场向下延拓波数域算子

3.2.1 Landweber迭代法波数域算子

由式(9)可得

因此,位场向下延拓的Landweber迭代法对应的波数域算子为

3.2.2 加速Landweber迭代法波数域算子

由式(14)和式(10)可得

比较式(17)和式(19)右端部分可以看出,这两种迭代法的向下延拓算子具有形式上的相似性。

4 数值试验

假设在一定深度处分布4个球体,建立如图2的坐标系:以球1的球心在平面的投影点作为坐标原点O,Y轴指向北方向,X轴垂直于Y轴指向东,Z轴垂直向下,XOY面为水平面。图2中1—4个球体球心坐标分别为(0,0,－1000)、(2000,2000,－900)、(2000,－2000,－1200)、(－2000,－2000,－1200),选定仿真区域范围为X(－4000,4000),Y(－4000,4000),网格间距100 m×100 m,点数81×81。4个球体的磁化倾角分别是55°、50°、60°、52°,4个磁化偏角分别是8°、10°、5°、7°。球体半径都为250 m,磁化强度200 A/m。

图2 球体位置坐标系Fig.2 Location of spheres coordinates system

4.1 波数域向下延拓算子滤波特性

设向下延拓深度h＝500 m(即5倍点距),以向下延拓原始波数域算子ekxh的频率响应为对比参考,Landweber迭代法波数域算子和加速Landweber迭代法波数域算子的频率响应如图3所示。图3中正则化参数α分别为0.5、0.75、1;图3(a)、3(b)显示了Landweber迭代法在n＝5、n＝10时的频率响应;图3(c)、3(d)显示了加速Landweber迭代法在n＝5、n＝10时的频率响应。

由图3可以看出:①原始波数域向下延拓算子具有高通滤波特性,而Landweber迭代法和加速Landweber迭代法波数域向下延拓算子为带通滤波,所以,在实际应用中,原始波数域向下延拓算子将放大原始位场中高频噪声,而后两种延拓算子则具有抑制高频噪声的作用;②当迭代次数从5次增加到10次,Landweber迭代法波数域算子的中频放大作用变化不明显,加速Landweber迭代法波数域算子中频放大作用得到较大增强;③Landweber迭代法和加速Landweber迭代法波数域算子对正则化参数α不敏感,即随α的改变其频率响应变化较小。

4.2 向下延拓结果分析

根据球体磁场公式计算得到0 m平面的磁异常分布,如图4所示,图中显示了0 m平面磁异常的等值线图,其中最大值是684.310 9 n T。为比较Landweber迭代法和加速Landweber迭代法的抗噪能力,在0 m高度磁异常数据中引入最大值的0.5%的正态分布随机干扰,并以该高度的含噪数据为观测平面位场数据进行向下延拓,0 m平面含噪磁异常的等值线图如图5所示。为方便起见,只图示y＝0 m剖面(这里称主剖面)的数据,通过4.1节的分析可知,两种算法对高频噪声具有抑制作用。所以首先判断噪声相对地磁异常数据为高频噪声才能进行向下延拓,为此,对随机干扰以及y＝0剖面的数据进行谱分析,如图6所示。图6(a)为对主剖面数据所加的随机干扰曲线,6(b)为理想的主剖面(图4中y＝0)数据,图6(c)为图6(a)所对应的谱分析,图6(d)为6(b)所对应的谱分析。在图6(c)中,不同的频率对应着一定大小的振幅,而在图6(d)中只有0～20 Hz范围内有较大振幅对应,其余频率段振幅很小或近乎为零,由6(c)和6(d)对比可以看出,高斯白噪声对于地磁场的这种大尺度、缓慢变化的低频特性可以看作是小尺度、变化较为剧烈的高频噪声,所以可以通过这两种迭代算法对位场进行向下延拓。

通过上述分析,可以对两种迭代法的向下延拓能力进行仿真。将迭代终止条件中的C设定为2,用两种方法将其向下延拓500 m(5倍点距),然后利用得到的延拓值同500 m深度的理论值进行误差计算,500 m深度理论值如图7所示。计算延拓误差的公式为

图8为叠加了随机干扰噪声的主剖面数据,用两种方法将图8中的含噪数据向下延拓500 m后,两种方法的延拓结果如图9所示。图9(a)和图9(b)分别为两种迭代法延拓值与理论值的对比,红线表示理论值,黑线表示延拓值。从图9(a)和图9(b)可见,两种迭代方法都能稳定有效地进行向下延拓,显示出两种迭代法优秀的抗噪能力,而且在位场的极大值处下延误差最大。

图3 Landweber迭代法、加速Landweber迭代法和原始波数域向下延拓算子滤波特性曲线Fig.3 Filtering characteristic curve of wave number domain downward continuation operator of Landweber iterative method,accelerated Landweber iterative method and original

图4 观测平面磁异常Fig.4 Magnetic anomaly of observation plane

图5 观测平面含噪磁异常Fig.5 Magnetic anomaly of observation plane with noise

图6 随机噪声、主剖面数据与相应的谱分析Fig.6 Random noise,potential field data and their spectral analysis respectively

图7 500 m深度磁异常理论值Fig.7 The theoretical magnetic anomaly

图8 主剖面含噪声位场值Fig.8 Potential field value with noise

图9 延拓值与理论值的比较Fig.9 The comparison of continuation value and theoretical value

下面同样采取仿真的方法验证加速Landweber迭代法具有更高的计算效率。选取迭代终止条件中的δ为0 m平面磁异常数据最大值的0.5%和0.1%,即δ为3.421 6 n T和δ为0.684 3 n T,分别用两种迭代方法计算,得到的结果如表1和表2所示。

表1 δ＝3.421 6 nT时两种迭代法的迭代次数与误差Tab.1 Iterative times and error of two methods whenδ＝3.421 6 nT

表2 δ＝0.684 3 nT时两种迭代法的迭代次数与误差Tab.2 Iterative times and error of two methods whenδ＝0.684 3 nT

从表1和表2计算结果可以看出,对于相同的扰动,加速landweber迭代法迭代次数比Landweber迭代法少得多,即收敛速度加快,计算效率明显提高,而且在误差上也满足要求,充分体现了加速Landweber迭代法的优越性。

在上述研究两种迭代法向下延拓的抗干扰性以及两种迭代法的计算效率中,将迭代终止条件中的C设定为2,显示了算法的收敛性。为了显示C值对加速效果的影响十分显著,选取不同的C值,分别采用两种迭代法进行向下延拓,噪声选取为0 m平面磁异常数据最大值的0.4%,即δ为2.737 2 n T得到的结果如表3所示。

表3 不同C值情形下两种迭代法的迭代次数与误差Tab.3 Iterative times and error of two methods in different value of C

从表3结果可知,随着C值的增加会使得两种迭代法迭代次数相应变少,但是相应的误差会随之增加,使得延拓精度降低。如果在延拓的精度要求不是很高的情形下,可以适当考虑增加C值以增快迭代进程,反之则相应减小C值提高延拓精度。

5 结 论

Landweber迭代法能够处理位场向下延拓问题,但是缓慢的收敛速度使其广泛应用产生了困难。本文在分析Landweber迭代法的基础上,提出采用加速Landweber迭代法高效解决位场向下延拓问题。推导得到两种迭代法的波数域向下延拓算子,详细分析了两种算子的滤波特性,发现两种算子都能够抑制随机噪声。最后给出一则实例表明加速Landweber迭代法能够提高计算速度,而且满足一定的延拓精度要求,能够对其在位场向下延拓问题的广泛应用提供一定参考。

[1] PANINA L V.Asymmetrical Giant Magneto-impedance (AGMI)in Amorphous Wires[J].Journal of Magnetism and Magnetic Materials,2002,249(1):278-287.

[2] HRISTOFOROU E,CHIRIAC H,NEAGU M,et al.New Load Cells and Torque Meters Based on Soft Magnetic Amorphous Alloy Wires[J].Sensors and Actuators,1998,68:307-315.

[3] H UANG Xuegong,WANG Jiong.Error Analysis and Compensation Methods for Geomagnetic Signal Detection System[J].Acta Armamentarii,2011,32(1):33-36.(黄学功,王炅.地磁信号检测系统误差分析与补偿方法研究[J].兵工学报,2011,32(1):33-36.)

[4] LIU Jianjing,ZHANG He,DING Libo.Static Calibration of Geomagnetic Sensors for Attitude Measurement[J].Journal of Nanjing University of Science and Technology, 2012,36(1):127-131.(刘建敬,张合,丁立波.姿态检测地磁传感器静态校正技术[J].南京理工大学学报,2012, 36(1):127-131.)

[5] ALONSO R,SHUSTER M.Complete Linear Attitude Independent Magnetometer Calibration[J].The Journal of the Astronautical Sciences,2002,50(4):477-490.

[6] COOPER G.The Stable Downward Continuation of Potential Field Data[J].Exploration Geophysics,2004,35: 260-265.

[7] YAO Changli,LI Hongwei,ZHENG Yuanman,et al.Research on Iteration Method Using in Potential Field Transformations[J].Chinese Journal of Geophysics, 2012,55(6):2062-2078.(姚长利,李宏伟,郑元满,等.重磁位场转换计算中迭代法的综合分析与研究[J].地球物理学报,2012,55(6):2062-2078.)

[8] YU Bo,ZHAI Guojun,LIU Yanchun,et al.The Downward Continuation Method of Aeromagnetic Data to the Sea Level[J].Acta Geodaetica et Cartographica Sinica,2009,38(3):202-209.(于波,翟国君,刘雁春,等.利用航磁数据向下延拓得到海面磁场的方法[J].测绘学报,2009,38 (3):202-209.)

[9] NING Jinsheng,WANG Haihong,LUO Zhicai.Downward Continuation of Gravity Signals Based on the Multiscale Edge Constraint[J].Chinese Journal of Geophysics,2005, 48(1):63-68.(宁津生,汪海洪,罗志才.基于多尺度边缘约束的重力场信号的向下延拓[J].地球物理学报, 2005,48(1):63-68.)

[10] SCHWARZ K P,SIDERIS M G,FORSBERG R.The Use of FFT Techniques in Physical Geodesy[J].Geophysical Journal International,1990,100:485-514.

[11] XU Shizhe.A Comparison of Effects between the Iteration Method and FFT for Downward Continuation of Potential Fields[J].Chinese Journal of Geophysics,2007,50(1): 285-289.(徐世浙.迭代法与FFT法位场向下延拓效果的比较[J].地球物理学报,2007,50(1):285-289.)

[12] JIANG Tao,LI Jiancheng,WANG Zhengtao,et al.Solution of Ill-posed Problem in Downward Continuation of Airborne Gravity[J].Acta Geodaetica et Cartographica Sinica,2011,40(6):684-689.(蒋涛,李建成,王正涛,等.航空重力向下延拓病态问题的求解[J].测绘学报, 2011,40(6):684-689.)

[13] WANG Xingtao,SHI Pan,ZHU Feizhou.Regularization Methods and Spectral Decomposition for the Downward Continuation of Airborne Gravity Data[J].Acta Geodaetica et Cartographica Sinica,2004,33(1):33-38.(王兴涛,石磐,朱非洲.航空重力测量数据向下延拓的正则化算法及其谱分解[J].测绘学报,2004,33(1): 33-38.)

[14] GU Yongwei,GUI Qingming.Study of Regularization Based on Signal-to-noise Index in Airborne Gravity Downward to the Earth Surface[J].Acta Geodaetica et Cartographica Sinica,2010,39(5):458-464.(顾勇为,归庆明.航空重力测量数据向下延拓基于信噪比的正则化方法的研究[J].测绘学报,2010,39(5):458-464.)

[15] GU Yongwei,GUI Qingming,BIAN Shaofeng,et al.Comparison between Tikhonov Regularization and Truncated SVD in Geophysics[J].Geomatics and Information Science of Wuhan University,2005,30(3): 238-241.(顾勇为,归庆明,边少锋,等.地球物理反问题中两种正则化方法的比较[J].武汉大学学报:信息科学版,2005,30(3):238-241.)

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

[17] MARTIN H.Accelerated Landweber Iterations for the Solution of Ill-posed Equations[J].Numerische Mathematik,1991,60:341-373.

[18] GUAN Zhining.Geomagnetic Field and Magnetic Exploration[M].Beijing:Geological Publishing House,2005.(管志宁.地磁场与磁力勘探[M].北京:地质出版社,2005.)

[19] WANG Yanfei.Computational Methods for Inverse Problems and Their Applications[M].Beijing:Higher Education Press,2007.(王彦飞.反演问题的计算方法及其应用[M].北京:高等教育出版社,2007.)

[20] SCHERZER O.A Modified Landweber Iteration for Solving Parameter Estimation Problems[J].Applied Mathematics Problems,1998,35:48-65.

[21] JIN Q N,AMATO U.A Discrete Scheme of Landweber Iteration for Solving Nonlinear Ill-posed Problems[J].Journal of Mathematical Analysis and Applications,2001, 253:187-203.

[22] XU Tianzhou.Appliance Functional Analysis[M].Beijing: Science Press,2002.(许天周.应用泛函分析[M].北京:科学出版社,2002.)

[23] WANG Guorong.Matrix and Operator of Generalized Inverse [M].Beijing:Science Press,1998.(王国荣.矩阵与算子广义逆[M].北京:科学出版社,1998.)

(责任编辑:陈品馨)

Accelerated Landweber Iteration Method for the Downward Continuation of Potential Field

ZHU Zhanlong1,YANG Gongliu1,2,YANG Shujie2
1.School of Instrumentation Science and Engineering,Southeast University,Nanjing 210096,China;2.School of Instrumentation Science and Opt-electronics Engineering,Beijing University of Aeronautics and Astronautics University,Beijing 100191,China

Aeromagnetic data downward continuation to different altitudinal data of field can improve the utilization percent of survey data.Landweber iteration method is an impactful way to ravel out the downward continuation of potential field.In view of the slow convergence speed of Landweber iteration method,an accelerated Landweber iteration method is used.Then,two wave number domain operators of two iteration methods are deduced and the filtering characteristics of operators are analyzed.Finally,the results of the numerical simulations show that the accelerated Landweber iteration method can not only deal with downward continuation of potential field,but also more efficient than Landweber iteration method.

potential field;the downward continuation;regularization method;accelerated Landweber iteration method;convergence speed

ZHU Zhanlong(1984—),male,PhD candidate,majors in geomagnetic data processing.E-mail:zhuzhanlong1984＠163.com

P223

A

1001-1595(2014)05-0458-08

中央高校基本科研业务费专项(YWF-10-01-B30)

2013-03-14

朱占龙(1984—),男,博士生,研究方向为地磁数据处理。

ZHU Zhanlong,YANG Gongliu,YANG Shujie.Accelerated Landweber Iteration Method for the Downward Continuation of Potential Field[J].Acta Geodaetica et Cartographica Sinica,2014,43(5):458-465.(朱占龙,杨功流,杨淑洁.位场向下延拓的加速Landweber迭代法[J].测绘学报,2014,43(5):458-465.)

10.13485/j.cnki.11-2089.2014.0088

修回日期:2013-09-02