2. 北京航空航天大学 仪器科学与光电工程学院，北京 100191
2. School of Instrumentation Science and Opt-electronics Engineering,Beijing University of Aeronautics and Astronautics University,Beijing 100191,China
1 引 言
随着磁传感器技术的进步[1, 2]以及磁力测量系统的改进[3, 4, 5]，航空磁力测量已经日趋成熟，尤其是在海平面、地表和超低空等测量条件恶劣的地区，航空磁力测量更是不可或缺。采用航磁数据进行不同高度的向下延拓可以弥补条件恶劣地区磁场资料的不足，因此，深入研究稳定的向下延拓方法具有重要的理论意义和应用价值。
向下延拓问题是一个将噪声数据放大的过程，属于病态问题[6, 7, 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]
可见，式(1)右端为f(x,h)与h/π(x2+h2)的卷积。令fδx=f(x,0)、g(x)=f(x,h)以及φ(x)=h/π(x2+h2)，则式(1)可以表示为
式中，*代表卷积。对式(2)两边进行傅里叶变换，可得 式中，kx为角频率；F(kx)、G(kx)分别表示为F(kx)=∫－∞+∞fδ(x)e－ixkxdx,G(kx)=∫－∞+∞g(x)e－ixkxdx。
经计算可得φ(x)的傅里叶变换为Φ(kx)=e－kxh，即为波数域位场向上延拓算子；则Φ－1(kx)=ekxh为波数域位场向下延拓算子。由式(3)可得G(kx)=Φ－1(kx)F(kx)，然后再进行傅里叶反变换可得到向下延拓平面的位场
观测平面上的位场数据fδ(x)往往是带有高频噪声的数据，当位场向下延拓时，由于随着|kx|的增加，延拓算子(e|kx|h＞1)呈指数规律增加，所以求得的g(x)将非常不稳定，这就是位场向下延拓不稳定的本质。由此可见，当高频噪声存在时，向下延拓视为一个不适定问题，对其直接作傅里叶变换以及反变换所得到的结果，不能作为向下延拓平面位场的稳定近似解。
3 位场向下延拓的迭代法及相应的波数域算子同第一类Fredholm积分方程比较，式(2)可以表示为[19]
式中，K为Fredholm算子，对于形如式(5)的不适定问题求解一般采用正则化方法，其中Landweber迭代法是一种广泛应用的方法。 3.1 位场向下延拓迭代法 3.1.1 Landweber迭代法
Landweber迭代法是最速下降法的一种变形，其迭代式如下[20, 21]
式中，正则化参数α满足0＜α＜1/||K||2；K*为算子K的伴随算子，而且算子K显然为对称的线性紧算子，则有K*=K[22]。
对于给定的初值g0(x)，迭代式(6)可以展开为
继续按照式(7)中的步骤展开gn(x)，式(7)可以改写为下面的级数形式
当取迭代初值为g0(x)=αK*fδ(x)时，迭代格式可写为
式中，fδ(x)为观测面上受到干扰的位场数据，设观测面理想的位场数据为f(x)，且满足||f(x)－fδ(x)||≤δ，其中δ表示位场数据与理想数据的最大范数值。定义 则Landweber迭代法可以描述为下面的算法:
(1) 给定迭代初值T0=αK*。
(2) 计算Tn+1=Tn(I－αK*K)+αK*。
(3) 计算gn(x)=Tnfδ(x)，如果满足||Kgn(x)－fδ(x)||≤Cδ迭代终止；若不满足，n=n+1，转向(2)。其中，C一般为固定常数，它决定了Landweber迭代算法是否具有收敛性[19]。很显然，在δ确定后，C值越大使得迭代终止条件Cδ越大，则迭代次数变少，但是相应的结果精度会降低。
由上述算法可以看出，主要决定算法收敛速度的是第(2)步，在实际计算中，可以先对(2)迭代计算若干步之后，再进行第(3)步判断终止条件，可以在一定程度上加快计算速度。
3.1.2 加速Landweber迭代法通过对Landweber迭代法分析可知，如果足够大n(Tn的下标)能够快速计算出来，则迭代终止条件就能较易得到满足，从而使得计算效率提高。为此，对于Tn的计算，采用如下迭代形式
[23]
下面说明式(11)在计算速度方面的快速性：
由Tn+1=Tn(I－αK*K)+αK*可得
从而由式(12)可得 于是，设定
0=αK*时
1=0(2I－K0)=αK*+T0(I－αK*K)=T1
2=1(2I－K1)=T1+T1(I－KT1)=
从式(14)可以看出，用式(11)计算n次，相当于Landweber迭代法迭代2n-1次，这样可以使得计算效率明显提升。于是，加速Landweber迭代法步骤如下:
(1) 给定迭代初值0=αK*。
(2) 计算n+1=n(2I－Kn)。
(3) 计算gn(x)=nfδ(x)，如果满足||Kgn(x)－fδ(x)||≤Cδ迭代终止；若不满足，n=n+1，转(2)。
由式(9)可得
对式(15)两边同时作傅里叶变换可得
因此，位场向下延拓的Landweber迭代法对应的波数域算子为 3.2.2 加速Landweber迭代法波数域算子
由式(14)和式(10)可得
对式(18)右端作傅里叶变换可得加速Landweber迭代法波数域算子
比较式(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 |
设向下延拓深度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迭代法和原始波数域向下延拓算子滤波特性曲线 Fig. 3 Filtering characteristic curve of wave number domain downward continuation operator of Landweber iterative method,accelerated Landweber iterative method and original |
由图 3可以看出：① 原始波数域向下延拓算子具有高通滤波特性，而Landweber迭代法和加速Landweber迭代法波数域向下延拓算子为带通滤波，所以，在实际应用中，原始波数域向下延拓算子将放大原始位场中高频噪声，而后两种延拓算子则具有抑制高频噪声的作用；② 当迭代次数从5次增加到10次，Landweber迭代法波数域算子的中频放大作用变化不明显，加速Landweber迭代法波数域算子中频放大作用得到较大增强；③ Landweber迭代法和加速Landweber迭代法波数域算子对正则化参数α不敏感，即随α的改变其频率响应变化较小。
4.2 向下延拓结果分析根据球体磁场公式计算得到0 m平面的磁异常分布，如图 4所示，图中显示了0 m平面磁异常的等值线图，其中最大值是684.310 9 nT。为比较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(b)和6(d)对比可以看出，高斯白噪声对于地磁场的这种大尺度、缓慢变化的低频特性可以看作是小尺度、变化较为剧烈的高频噪声，所以可以通过这两种迭代算法对位场进行向下延拓。
|
| 图 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 |
通过上述分析，可以对两种迭代法的向下延拓能力进行仿真。将迭代终止条件中的C设定为2，用两种方法将其向下延拓500 m(5倍点距)，然后利用得到的延拓值同500 m深度的理论值进行误差计算，500 m深度理论值如图 7所示。计算延拓误差的公式为
|
| 图 7 500 m深度磁异常理论值 Fig. 7 The theoretical magnetic anomaly |
图 8为叠加了随机干扰噪声的主剖面数据，用两种方法将图 8中的含噪数据向下延拓500 m后，两种方法的延拓结果如图 9所示。图 9(a)和图 9(b)分别为两种迭代法延拓值与理论值的对比，红线表示理论值，黑线表示延拓值。从图 9(a)和图 9(b)可见，两种迭代方法都能稳定有效地进行向下延拓，显示出两种迭代法优秀的抗噪能力，而且在位场的极大值处下延误差最大。
|
| 图 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 nT和δ为0.684 3 nT，分别用两种迭代方法计算，得到的结果如表 1和表 2所示。
| 迭代法 | 迭代次数 | 误差/nT |
| Landweber迭代法 | 541 | 2.831 4 |
| 加速Landweber迭代法 | 11 | 2.672 6 |
| 迭代法 | 迭代次数 | 误差/nT |
| Landweber迭代法 | 5534 | 0.628 2 |
| 加速Landweber迭代法 | 16 | 0.602 7 |
从表 1和表 2计算结果可以看出，对于相同的扰动，加速landweber迭代法迭代次数比Landweber迭代法少得多，即收敛速度加快，计算效率明显提高，而且在误差上也满足要求，充分体现了加速Landweber迭代法的优越性。
在上述研究两种迭代法向下延拓的抗干扰性以及两种迭代法的计算效率中，将迭代终止条件中的C设定为2，显示了算法的收敛性。为了显示C值对加速效果的影响十分显著，选取不同的C值，分别采用两种迭代法进行向下延拓，噪声选取为0 m平面磁异常数据最大值的0.4%，即δ为2.737 2 nT得到的结果如表 3所示。
| 迭代法 | C值 | 迭代次数 | 误差/nT |
| Landweber迭代法 | 3 | 503 | 3.081 6 |
| 4 | 397 | 4.175 9 | |
| 加速Landweber迭代法 | 3 | 10 | 3.114 2 |
| 4 | 9 | 4.178 3 |
从表 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] | HUANG 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.) |