改進高斯－牛頓法的位場向下延拓

曾小牛，劉代志，牛 超，齊 瑋

1.第二炮兵工程大學(xué)907教研室，陜西 西安 710025；2.清華大學(xué) 信息科學(xué)技術(shù)學(xué)院自動化系，北京 100084；3.第二炮兵裝備研究院，北京 100085

1 引 言

位場向下延拓在地質(zhì)勘探資料解析、磁性目標(biāo)探測定位及重力場和地磁場輔助慣性導(dǎo)航基準(zhǔn)圖制備等很多方面都有十分重要的作用。位場向下延拓是一個典型的線性不適定反問題，對于該類問題的求解一般有兩種策略［1］：一是正則化方法；二是迭代法。正則化方法方面，文獻［2］給出了航空重力向下延拓的正則化算法，并對重力數(shù)據(jù)處理中常用的幾種向下延拓方法進行了比較，得出正則向下延拓效果相對最好的結(jié)論；文獻［3］依據(jù)信噪比構(gòu)造正則化矩陣，以極小化均方誤差為目標(biāo)選取正則化參數(shù)構(gòu)造正則化矩陣，提出基于信噪比的正則化方法；文獻［4］提出將基于多個最優(yōu)正則化參數(shù)的廣義嶺估計用于航空重力向下延拓，并通過對比試驗說明，廣義嶺估計在可靠性、精度和穩(wěn)定性方面均優(yōu)于Tikhonov法和嶺估計；文獻［5］將系數(shù)矩陣相對較大和相對較小的奇異值區(qū)別對待，提出Tikhonov雙參數(shù)正則化法并驗證了該方法的有效性。迭代法方面，文獻［6—7］提出了位場向下延拓的積分迭代法，相較傳統(tǒng)的FFT法，該迭代法能夠穩(wěn)定向下延拓20倍點距，同時，有大量文獻進行了相關(guān)研究［8－15］。

本文首先將位場向下延拓視為一個線性優(yōu)化問題，然后采用高斯－牛頓法來實現(xiàn)位場向下延拓。其次，基于奇異值分解定理，探討了高斯－牛頓法的濾波函數(shù)及其濾波函數(shù)的濾波特性，并從濾波函數(shù)濾除高頻噪聲需要和目標(biāo)函數(shù)最小化的角度，提出了改進算法。改進算法采用遞增幾何級數(shù)計算迭代法正則參數(shù)，并引入最優(yōu)變步長策略來實現(xiàn)迭代法正則參數(shù)和迭代步長計算的完全自適應(yīng)。最后，通過理論重力模型和實測航磁數(shù)據(jù)驗證了改進迭代法的有效性。

2 位場正則向下延拓原理

位場向下延拓的數(shù)學(xué)模型可以簡化為

由此得位場向下延拓的最小二乘解

位場向下延拓本質(zhì)上是一個病態(tài)問題，延拓過程中會放大觀測噪聲。若采用最小二乘方法求解，即使較小的觀測誤差也會導(dǎo)致解的嚴(yán)重失真。一個解決該問題的有效途徑是在目標(biāo)函數(shù)中添加一個正則項。Tikhonov正則化即是一種廣泛應(yīng)用的方法，它是指最小化如下目標(biāo)函數(shù)

式中，λ為正則參數(shù)，用于平衡解的不穩(wěn)定性及光滑性。式（4）的解為

對系數(shù)矩陣A作奇異值分解

式（5）在頻域內(nèi)的形式為

3 位場改進高斯－牛頓法向下延拓原理

3.1 高斯－牛頓法基本原理

設(shè)xk是式（2）函數(shù)J（x）極小值的第k次近似，并將J（x）在xk處作二階Taylor展開，得

將J（x）的梯度ΔJ（x）＝AT（bδ－Ax）和Hessian矩陣代入式（10）得迭代法的迭代式為

由于反問題的病態(tài)性，直接用高斯－牛頓法求解，解出的值依然很不穩(wěn)定，效果很不理想。為解決這個問題，對于式（11），一般在ATA后面加上一個帶參數(shù)的單位陣，以使其具有良好的正則性質(zhì)，克服當(dāng)ATA奇異時，算法收斂到一個非駐點的情況。改進后的迭代格式如下

式中，rk＝bδ－Axk為殘差。

對比式（12）和式（5）可知，高斯 －牛頓法迭代得到的解xk＋1，都是上一步迭代得到的解xk加上用Tikhonov正則化方法處理殘差rk后得到的值的ak倍（即ak（ATA＋λI）－1ATrk）而得到的。每次迭代過程中，當(dāng)正則參數(shù)選擇合適，Tikhonov正則能夠有效利用殘差rk中有效信息進行迭代結(jié)果修正的同時，壓制殘差rk中的噪聲（來自殘差rk＝bδ－Axk中的bδ）。因此，定性地講，高斯－牛頓法能夠獲得比Tikhonov正則法更加精確的解，同時注意到每次迭代計算都僅包含一次Tikhonov正則計算，因此相較Tikhonov正則法，該迭代法并沒有增加過多的運算量。

值得注意的是，當(dāng)步長ak＝1，式（12）等價于迭 代 Tikhonov 法［18］、Levenberg－Marquardt法［19］或預(yù)處理 Landweber迭代法［20］。

3.2 濾波函數(shù)對比分析

當(dāng)ak＝1時，運用數(shù)學(xué)歸納法和系數(shù)矩陣A的奇異值分解可得

可見，ak＝1時的高斯－牛頓法的濾波函數(shù)為

當(dāng)k＝0時，濾波函數(shù)φ（）即轉(zhuǎn)化為Tikhonov正則濾波函數(shù)的對比如圖1所示（以4.1節(jié)理論重力模型的相關(guān)參數(shù)為例）。其中，設(shè)濾波函數(shù)所形成的低通濾波器的截止頻率為0.005（濾波函數(shù)的值降為0.5）。由圖1的比對可知，高斯－牛頓法濾波函數(shù)形成的濾波器的過渡帶明顯窄于（或稱下降陡度明顯大于）Tikhonov正則法，因此，從低通濾波器的濾波特性角度來講，高斯－牛頓法的濾波函數(shù)優(yōu)于Tikhonov正則法的濾波函數(shù)。

圖1 Tikhonov和高斯－牛頓法濾波函數(shù)特性的對比Fig.1 The characteristic comparison of filter functions for Tikhonov and Guass－Newton methods

3.3 正則參數(shù)的選取

針對步長ak＝1的情況，文獻［21］進一步考慮了式（12）的非穩(wěn)定隱式迭代形式

并相應(yīng)采用正則參數(shù)選擇的幾何級數(shù)法

式中，λ0為正則參數(shù)初始設(shè)定值，且λ0＞0，0＜q＜1。

文獻［22］建立了具有線性收斂速率的自適應(yīng)選擇方法

上述兩種正則參數(shù)計算方法是伴隨迭代的進行不斷地更新，具有自適應(yīng)的性質(zhì)，且計算簡單，因而具有相對較高的計算效率。但是，它們也存在嚴(yán)重的缺點，即存在嚴(yán)重的半收斂性：在迭代法迭代的過程中，伴隨迭代次數(shù)的增加，向下延拓的誤差先減后增。下面先定性分析其半收斂性產(chǎn)生的原因。

首先，將式（8）代入式（15）中可得迭代法非穩(wěn)定隱式迭代形式在頻域的形式

同樣以4.1節(jié)重力理論模型參數(shù)為例，Tik－h(huán)onov正則的濾波函數(shù)wλ）在不同正則參數(shù)λ條件下的變化如圖2所示?？梢?，wλ（）表現(xiàn)為一低通濾波，正則參數(shù)λ控制著該濾波函數(shù)的低通濾波特性：伴隨正則參數(shù)λ的減小，其低通濾波的效果越來越不明顯，即通過的高頻信息會越來越多。由此，當(dāng)用式（16）的幾何級數(shù)法或式（17）的自適應(yīng)法計算正則參數(shù)時，由于幾何級數(shù)法中公比0＜q＜1、自適應(yīng)法中，正則參數(shù)λk伴隨迭代法迭代的進行會不斷減小，則根據(jù)濾波函數(shù)wλ（）的濾波特性，式（18）右端第二項中的殘差rk＝bδ－Axk里面的高頻噪聲（來自bδ）伴隨迭代的進行會不斷地累積到最終延拓結(jié)果xk＋1中，使解的誤差不斷增大，并最終導(dǎo)致迭代法的半收斂性。由此可見，迭代法在計算正則參數(shù)λk時，應(yīng)使其在一定范圍內(nèi)不斷增加，以使迭代的最終結(jié)果穩(wěn)定。

圖2 奇異值和Tikhonov濾波函數(shù)wλ（）Fig.2 The singular valueand the filter function wλ（）for Tikhonov method

針對上述分析，本文相應(yīng)提出如下的幾何遞增級數(shù)法用于正則參數(shù)的計算

式中，λ為Tikhonov正則參數(shù)（通過廣義交叉校驗（GCV）準(zhǔn)則［23］或 L－曲線法［24］求得），而該方法的關(guān)鍵在于選取合適的大于1的公比q。同時，為了確保迭代法在迭代過程中，對上一步得到的結(jié)果進行充分修正，根據(jù)Tikhonov濾波函數(shù)的濾波特性分析，q不能過大，否則，λk的過快增長會使修正過程過快地結(jié)束而導(dǎo)致迭代法的最終結(jié)果改善有限，這在第四節(jié)的試驗中可以看出。

3.4 殘差最小步長準(zhǔn)則

首先，設(shè)迭代法第k次迭代時的殘差為rk＝bδ－Axk，則其第k＋1次迭代的殘差

3.5 算法具體實現(xiàn)步驟

綜上所述，本文改進高斯－牛頓法的具體實現(xiàn)步驟可描述如下：首先，令x0＝0，以廣義交叉校驗（GCV）準(zhǔn)則或L－曲線法（本文采用GCV準(zhǔn)則）求得Tikhonov正則法的正則參數(shù)λ；然后，用式（19）即λk＝λqk－1來計算迭代法的正則參數(shù)，且在每次迭代法計算過程中，采用由殘差最小準(zhǔn)則導(dǎo)出的步長計算公式（式（21））來更新迭代步長ak。從上面的描述可以看出，改進方法具有完全自適應(yīng)的性質(zhì)。下面的數(shù)值試驗表明，改進迭代法的收斂速度較快，一般迭代10～20次左右就可以獲得較好地收斂效果。

4 數(shù)值試驗及結(jié)果比較

4.1 理論重力模型

采用文獻［9］的雙球體重力模型。兩個球體的參數(shù)相同：剩余密度Δσ＝1.0t／m3、半徑r＝0.5km、中心埋深1.8km、球心x坐標(biāo)分別為10km和15km、球心y坐標(biāo)12.5km、z坐標(biāo)向下為正，線數(shù)M和每線點數(shù)N同為512、點距Δx和線距Δy都為50m。

試驗中，以h＝0m平面處的重力數(shù)據(jù)為觀測數(shù)據(jù)，并加入零均值、均方差為該高度理論重力異常絕對均值10%的高斯白噪聲以模擬實際情況，其結(jié)果如圖3所示。分別用Tikhonov法和本文的改進高斯－牛頓法將加噪數(shù)據(jù)向下延拓1000m（20倍點距），并用得到的延拓值xc與1000m深度處的真實重力異常值xt采用均方誤差（root mean square error，RMSE）

和相關(guān)性誤差（relative error，RE）

來計算、對比延拓誤差［25］，其結(jié)果如表1（迭代法迭代次數(shù)k＝10）和圖4（僅顯示主剖面）所示。

表1 Tikhonov法和本文方法的均方誤差及相關(guān)性誤差的比較Tab.1 The comparisons of root mean square errors and relative errors for Tikhonov and the method in this paper

圖3 加噪重力異常等值線圖Fig.3 Contours of the gravity anomaly data with noise

圖4 延拓結(jié)果在主剖面的對比Fig.4 The comparison of downward continuation results in main profile

本文迭代法的均方誤差和相關(guān)性誤差隨公比q和迭代次數(shù)的變化分別如圖5和圖6所示（迭代次數(shù)k＝0對應(yīng)Tikhonov法）。從圖5和圖6中的對比可知，延拓誤差隨q的變化可大致分為3個階段：當(dāng)公比q≤1時，迭代法存在非常顯著的半收斂性，即延拓誤差隨迭代次數(shù)的增加先減小，然后迅速增加；當(dāng)公比1＜q≤2，迭代法尚存在一定的半收斂性，但總體表現(xiàn)平穩(wěn)；當(dāng)q＞2時，迭代法完全穩(wěn)定，但是最終的延拓誤差比1＜q≤2時的延拓誤差要大。上述試驗結(jié)果與3.3節(jié)的理論分析是一致的。

4.2 航磁數(shù)據(jù)試驗

航磁實測數(shù)據(jù)點距和線距均為250m，其等值線圖見圖7。將原始數(shù)據(jù)向上延拓5000m（20倍點距）以構(gòu)建另一個平面的數(shù)據(jù)。同時，考慮到向上延拓具有的平滑效應(yīng)，為模擬實際情況，本文在向上延拓后得到的磁場數(shù)據(jù)中加入零均值、均方差為該延拓航磁數(shù)據(jù)絕對均值1%的高斯白噪聲，加噪后航磁數(shù)據(jù)等值線圖如圖8所示。與模型對比試驗一樣，分別采用Tikhonov法和本文方法分別將加噪向上延拓數(shù)據(jù)向下延拓20倍點距到原來的數(shù)據(jù)平面，并將得到的延拓值xc與原高度處的真實航磁異常值xt進行比較，兩種方法產(chǎn)生的均方誤差和相關(guān)性誤差對比如表2所示（迭代法迭代次數(shù)k＝20）。

圖5 延拓均方誤差隨公比q和迭代次數(shù)的變化Fig.5 Relationship between root mean square errors and iteration numbers for the method in this paper

圖6 延拓相關(guān)性誤差隨公比q迭代次數(shù)的變化Fig.6 Relationship between relative errors and iteration numbers for the method in this paper

表2 Tikhonov法和本文方法的均方誤差及相關(guān)性誤差的比較Tab.2 The comparisons of root mean square errors and relative errors for Tikhonov method and the method in this paper

圖7 原始航磁等值線Fig.7 Contours of the original aeromagnetic field

圖8 向上延拓航磁異常加噪等值線Fig.8 Contours of upward continuation of aeromagnetic field with noise

兩種方法的延拓結(jié)果分別如圖9（Tikhonov法）和圖10（本文方法）所示。本文迭代法延拓的均方誤差和相關(guān)性誤差隨迭代次數(shù)的變化分別如圖11和圖12所示（迭代次數(shù)k＝0對應(yīng)Tikhonov法）。綜合圖9～圖12分析可知：① 本文的改進迭代法向下延拓的過程穩(wěn)定，且相較Tikhonov法，改進迭代法的延拓結(jié)果具有更小的延拓誤差、帶有更多的細節(jié)；② 選擇公比1＜q≤2既可以確保迭代法延拓結(jié)果的收斂性，又可以保證迭代法的結(jié)果得到最好的修正。

圖9 Tikhonov法向下延拓結(jié)果等值線Fig.9 Contours of the downward continued data for Tikhonov regularization method

圖10 本文方法向下延拓結(jié)果等值線Fig.10 Contours of the downward continued data for the method in this paper

圖11 本文方法均方誤差隨迭代次數(shù)的變化Fig.11Relationship between root mean square errors and iteration numbers for the method in this paper

圖12 本文方法相關(guān)性誤差隨迭代次數(shù)的變化Fig.12 Relationship between relative errors and iteration numbers for the method in this paper

5 結(jié)束語

位場高精度穩(wěn)定向下延拓方法的研究一直是其諸多相關(guān)應(yīng)用領(lǐng)域的研究熱點。本文提出一種基于遞增幾何級數(shù)自適應(yīng)選擇正則參數(shù)和殘差最小步長準(zhǔn)則的改進高斯－牛頓法用于位場向下延拓，并給出了迭代法的具體實現(xiàn)步驟。將改進方法應(yīng)用到理論重力模型和實測航磁數(shù)據(jù)的向下延拓試驗中，結(jié)果表明，改進迭代法能夠有效地克服原迭代法的半收斂性，同時，相對傳統(tǒng)的Tikhonov正則法具有更高的延拓精度。

［1］ WANG Yanfei.Computational Methods for Inverse Problems and Their Applications［M］.Beijing：Higher Education Press，2007.（王彥飛.反演問題的計算方法及其應(yīng)用［M］.北京：高等教育出版社，2007.）

［2］ WANG Xintao，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.（王興濤，石磐，朱非洲.航空重力測量數(shù)據(jù)向下延拓的正則化算法及譜分析［J］.測繪學(xué)報，2004，33（1）：33－38.）

［3］ 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.（顧勇為，歸慶明.航空重力測量數(shù)據(jù)向下延拓基于信噪比的正則化方法的研究［J］.測繪學(xué)報，2010，39（5）：458－464.）

［4］ 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.（蔣濤，李建成，王正濤，等.航空重力向下延拓病態(tài)問題的求解［J］.測繪學(xué)報，2011，40（6）：684－689.）

［5］ DENG Kailiang，HUANG Motao，BAO Jingyang，et al.Tikhonov Two－parameter Regulation Algorithm in Downward Continuation of Airborne Gravity Data［J］.Acta Geodaetica et Cartographica Sinica，2011，40（6）：690－696.（鄧凱亮，黃謨濤，暴景陽，等.向下延拓航空重力數(shù)據(jù)的Tikhonov雙參數(shù)正則化法［J］.測繪學(xué)報，2011，40（6）：690－696.）

［6］ XU Shizhe.The Integral－iteration Method for Continuation of Potential Field［J］.Chinese Journal of Geophysics，2006，49（4）：1176－1182（徐世浙.位場延拓的積分－迭代法［J］.地球物理學(xué)報，2006，49（4）：1176－1182.）

［7］ XU S Z，YANG J Y，YANG C F，et al.The Iteration Method for Downward Continuation of a Potential Field from a Horizontal Plane［J］.Geophsical Prospecting，2007，55（6），883－889.

［8］ ZHANG Hui，CHENG Longwei，REN Zhixin，et al.Analysis on Convergence of Iteration Method for Potential Fields Downward Continuation and Research on Robust Downward Continuation Method［J］.Chinese Journal of Geophysics.2009，52（4）：1107－1113.（張輝，陳龍偉，任治新，等.位場向下延拓迭代法收斂性分析及穩(wěn)健向下延拓方 法 研 究 ［J］.地 球 物 理 學(xué) 報，2009，52（4）：1107－1113.）

［9］ LIU Dongjia，HONG Tianqiu，JIA Zhihai，et al.Wave Number Domain Iteration Method for Downward of Potential Fields and Its Convergence［J］.Chinese Journal of Geophysics.2009，52（6）：1599－1605.（劉東甲，洪天求，賈志海，等.位場向下延拓的波數(shù)域迭代法及其收斂性［J］.地球物理學(xué)報，2009，52（6）：1600－1605.）

［10］ 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.（于波，翟國君，劉雁春，等.利用航磁數(shù)據(jù)向下延拓得到海面磁場的方法［J］.測繪學(xué)報，2009，38（3）：202－209.）

［11］ YU Bo，ZHAI Guojun，LIU Yanchun，et al.Analysis of Noise Effect on the Calculation Error of Downward Continuation with Iteration Method［J］.Chinese Journal of Geophysics，2009，52（8）：2182－2188.（于波，翟國君，劉雁春，等.噪聲對向下延拓迭代法的計算誤差影響分析［J］.地球物理學(xué)報，2009，52（8）：2182－2188.）

［12］ ZENG Xiaoniu，LI Xihai，HAN Shaoqing，et al.A Comparison of Three Iteration Methods for Downward Continuation of Potential Fields［J］.Progress in Geophysics.2011，26（3）：908－915.（曾小牛，李夕海，韓紹卿，等.位場向下延拓三種迭代方法之比較［J］.地球物理學(xué)進展，2011，26（3）：908－915.）

［13］ ZENG Xiaoniu，LI Xihai，LIU Daizhi，et al.Regularization Analysis of Integral－iteration Method and the Choice of Its Optimal Step－length［J］.Chinese Journal of Geophysics，2011，54（11）：2943－2950.（曾小牛，李夕海，劉代志，等.積分迭代法的正則性分析及其最優(yōu)步長的選擇［J］.地球物理學(xué)報，2011，54（11）：2943－2950.）

［14］ 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.（姚長利，李宏偉，鄭元滿，等.重磁位場轉(zhuǎn)換計算中迭代法的綜合分析與研究［J］.地球物理學(xué)報，2012，55（6）：2062－2078.）

［15］ GAO Yuwen，LUO Yao.WEN Wu.The Compensation Method for Downward Continuation of Potential Field from Horizontal Plane and Its Application［J］.Chinese Journal of Geophysics，2012，55（8）：2717－2756.（高玉文，駱遙，文武.補償向下延拓方法研究及應(yīng)用［J］.地球物理學(xué)報，2012，55（8）：2747－2756.）

［16］ CHEN Longwei，HU Xiaoping，WU Meiping，et al.Research on Spatial Domain Continuation Method for Potential Field［J］.Progress in Geophysics，2012，27（4）：1509－1518.（陳龍偉，胡小平，吳美平，等.空間域位場延拓新方法研究［J］.地球物理學(xué)進展.2012，27（4）：1509－1518.）

［17］ YUAN Yaxiang，SUN Wenyu.Optimization Theory and Methods［M］.Beijing：Science Press，1997.（袁亞湘，孫文瑜.最優(yōu)化理論與方法［M］.北京：科學(xué)出版社，1997.）

［18］ KING J T，CHILLINGWORTH D.Approximation Generalized Inverse by Iterated Regularization［J］.Numerical Function Analysis and Optimization.1979，1（5）：499－513.

［19］ RANGANATHAN A.The Levenberg－Marquardt Algorithm［EB／OL］.2004［2012－03－08］.http：∥www.ananth.in／Notes＿files／lmtut.pdf，2004.

［20］ PIANA M，BERTERO M.Projected Landweber Method and Preconditioning［J］.Inverse Problems，1997，13（2），441－464.

［21］ HANKE M，GROETSCH C W.Nonstationary Iterated Tikhonov Regularization［J］.Journal of Optimization Theory and Applications，1998，98（1）：37－53.

［22］ PEREVERZEV S V，SCHOCK E.Brakhage’s Implicit Iteration Method and the Information Complexity of Equations with Operators Having Closed Range［J］.Journal of Complexity，1999，15（3）：385－401.

［23］ GOLUB G H，HEATH M，WAHBA G.Generalized Cross－validation as a Method for Choosing a Good Ridge Parameter［J］.Technometrics，1979，21（2）：215－223.

［24］ HANSEN P C，O＇LEARY D P.The Use of the L－curve in the Regularization of Discrete Ill－posed Problems［J］.SIAM Journal on Scientific Computing，1993，14（6）：1487－1503.

［25］ ZHOU J，SHI G G，GE Z L.Study of Iterative Regularization Methods for Potential Field Downward Continuation in Geophysical Navigation［J］.Journal of Astronautics，2011，32（4）：787－794.