航空重力向下延拓病態問題的求解

蔣 濤，李建成，王正濤，張守建

1.武漢大學測繪學院，湖北武漢430079；2.地球空間環境與大地測量教育部重點實驗室，湖北武漢430079

航空重力向下延拓病態問題的求解

蔣 濤1，2，李建成1，2，王正濤1，2，張守建1，2

1.武漢大學測繪學院，湖北武漢430079；2.地球空間環境與大地測量教育部重點實驗室，湖北武漢430079

提出將廣義嶺估計用于求解航空重力向下延拓病態問題。研究求解逆Poisson積分問題的三種正則化方法：Tikhonov正則化、嶺估計和廣義嶺估計。利用EGM2008地球位模型設計模擬數值實驗，將飛行高度處含白噪聲的2.5′×2.5′重力擾動向下延拓至大地水準面上，與參考值作外部檢驗，全面檢驗、比較各向下延拓方法的可靠性、精度和穩定性。數值結果表明基于多個最優正則化參數的廣義嶺估計在延拓精度、穩定性和抗差性等方面表現良好。

航空重力；向下延拓；病態問題；正則化

1 引 言

隨著全球定位系統（GPS）定位精度的提高、慣性導航系統（INS）和重力測量系統的改進，航空重力測量在5～10km的空間分辨率尺度上能夠達到1～3mGal（1mGal＝10－5m／s2）的數據精度，所提供的重力場中、短波信息與地面重力測量和衛星重力測量形成有效互補，特別是在地面重力測量難以達到的地區（如兩極、森林和山區等），航空重力測量更是無可替代。航空重力測量觀測值經過數據預處理，可得到飛行高度處的重力異常或重力擾動，但大地測量領域需要的是大地水準面上的重力值，因此需將飛行高度處的重力異常或重力擾動向下調和延拓至大地水準面上。

航空重力向下延拓是一個數據噪聲放大的過程，屬于病態問題，若不采用合適的延拓方法，將無法得到穩定可靠的重力延拓解。在現有重力向下延拓方法中，以求解球外Dirichlet問題的逆Poisson積分法應用最為廣泛［1－6］，其他方法還包括最小二乘配置（LSC）［7］、快速傅里葉變換（FFT）［8－9］、直接代表法［10］和求解球內Dirichlet問題的直接法［11－12］。逆Poisson積分法涉及求解Poisson積分方程，當數據點間隔較小或延拓高度較大時法方程呈病態，解算時需引入正則化處理。以往研究大多采用基于單個正則化參數的Tikhonov正則化或嶺估計［4－6］，而理論上基于多個最優正則化參數的廣義嶺估計所得估值具有比普通嶺估計結果更小的平均均方誤差［13］，但目前尚未見到廣義嶺估計應用于重力向下延拓的研究結果發表，將基于廣義嶺估計的逆Poisson積分法用于向下延拓解算值得深入研究。本文以重力擾動為基本觀測量，介紹了基于逆Poisson積分法求解航空重力向下延拓病態問題的數學模型，給出逆Poisson積分問題的三種正則化解法——Tikhonov正則化、嶺估計和廣義嶺估計及確定正則化參數的詳細算法，利用EGM2008重力位模型設計模擬數值試驗全面比較分析了各正則化解法的可靠性、精度和穩定性。

2 數學模型

在球近似條件下，大地水準面可用半徑為R的地心參考球表面近似表示（R＝6 371km）。若已知大地水準面上的重力擾動δg（R，Ω′），求大地水準面外部空間某點的重力擾動δg（r，Ω），這一過程稱為重力擾動的向上延拓，可通過求解以大地水準面為邊界面、調和函數rδg所滿足的球外Dirichlet問題實現［14］。向上延拓值表達為Poisson積分如下

式中，σ表示全球積分區域；（r，Ω）表示球坐標（r，φ，λ）；K（r，ψ，R）表示Poisson積分核函數；L（r，ψ，R）表示點（r，Ω）與點（R，Ω′）之間的空間距離；ψ為兩點地心方向間的球面距離；dΩ′為積分面元。重力擾動的向下延拓是已知大地水準面外部空間的δg（r，Ω），求解Poisson積分的逆問題，得到大地水準面上的δg（R，Ω′）。此時式（1）成為第一類Fredholm積分方程，點（r，Ω）為觀測數據點，點（R，Ω′）為待求點，也可稱為流動點。

理論上Poisson積分式（1）應在全球范圍進行積分計算，但由于實際觀測數據有限，一般將球面積分區域劃分為近區和遠區。近區是一個以數據點（r，Ω）為中心，ψ0為半徑的球冠區域σ0；遠區是球面上的剩余部分σ－σ0。Poisson核K（r，ψ，R）隨著球面距離ψ增大而快速衰減。相對于近區流動點對數據點的貢獻，遠區流動點貢獻的量級很小，可由位系數和Poisson核截斷系數近似計算［1］

式中，Nf是所采用的位系數的最大階數；a為參考橢球長半軸；ˉCnm、ˉSnm為擾動位系數；ˉPnm（sinφ）為完全規格化勒讓德函數；δgn（Ω）為重力擾動面球諧展開式；Qn（r，ψ0，R）為Poisson核截斷系數。

式（1）的離散形式如下

式中，i＝1，2，…，M，M為觀測數據點個數；N為待求點個數，M≥N；Aij為M×N階設計矩陣A的元素。矩陣A的對角線元素為

A的非對角線元素為

式中，ΔΩj為第j個流動點所處格網的面積；Nc為處于球冠區σ0內的流動點個數。式（5）的解δg（R，Ωj）（j＝1，2，…，N）即為大地水準面上的重力擾動延拓值。

3 正則化方法

向下延拓的觀測方程（5）寫成標準Gauss－Markov模型形式如下

式中，y為M階觀測值向量；A為M×N階設計矩陣；x為N階未知數向量；e為觀測噪聲，其數學期望為0；D｛y｝為觀測值的方差－協方差陣；σ2是單位權方差；P是權陣，其最小二乘解為

向下延拓本質上是一個病態問題，延拓過程中會放大觀測噪聲和未模型化的信號。盡管將連續Poisson積分化為離散數值積分可在一定程度上抑制高頻噪聲，但隨著數據格網間隔的減小和延拓高度的增加，設計矩陣A的復共線性增強，有可能呈病態甚至奇異。若采用經典最小二乘方法求解，即使較小的觀測誤差也會導致解的嚴重失真甚至錯誤。解決該問題以得到穩定延拓解的一個途徑是對法方程作正則化處理。本文主要討論標準Tikhonov正則化、嶺估計和廣義嶺估計。不同正則化方法依據不同的正則化準則，正則化參數有所不同。Tikhonov正則化、嶺估計只有一個正則化參數α，廣義嶺估計則是基于多個參數αi（i＝1，2，…，N）的正則化方法。正則化參數起到平衡觀測誤差引起的估值誤差和正則化造成的估值偏差的作用。選擇最優的正則化參數可使得兩項誤差之和達到最小。α太小無法有效抑制數據觀測誤差，太大則使正則化解過于平滑，因此正則化方法的關鍵在于如何選擇合適的正則化參數α。

3.1 Tikhonov正則化

Tikhonov正則化采用如下最小化準則［15］

Tikhonov正則化解為

式中，I為單位陣；α＞0為正則化參數。Tikhonov正則化的關鍵是選擇合理的正則化參數α，本文采用廣義交叉檢驗（GCV）確定正則化參數。選取正則化參數的準則為滿足下述最小化條件［16］

式中，算子tr表示求矩陣的跡；Qα＝A（ATPA＋αI）－1ATP。正則化參數的確定采用搜索法，計算步驟如下：

（1）首先在α取較大的區間內（如［1×10－9，10］）按式（12）計算F（α）的值，為加快搜索速度α以10倍的量級遞增直至上限值，找出F（α）最小值對應的α0；

3.2 嶺估計

嶺估計的正則化解為［17－18］

式中，α＞0為嶺參數（正則化參數）。式（13）與Tikhonov正則化解（式11）是相同的，不同之處在于選擇正則化參數所依據的準則。嶺估計以使正則化解的平均均方誤差（MSE）最小作為確定正則化參數的準則。

對式（8）的法方程系數陣作特征值分解

式中，V為法方程系數陣的特征向量組成的正交矩陣，Λ為對角陣，其元素為特征值λ1≥λ2…≥λN＞0。式（13）可寫為

正則化解的平均均方誤差（MSE）定義為平均均方誤差矩陣（MSEM）的跡，寫成以特征值表示的譜分解式如下［17－18］

式中，zi為列向量Z＝VTx＾α的第i個元素值。將上式對參數α求偏導數得到

令上式為零，得到關于參數α的方程

滿足此方程的α即為能使平均均方誤差達到最小的最優正則化參數，可通過迭代計算得到α。由于式（18）是非線性方程，無法直接求解，應采用數值方法求解。具體迭代計算步驟如下：

（1）選取x＾α的迭代初始值。用最小二乘法求得的驗后單位權方差［KGσ＊5］⌒2代替σ2。

（3）若F（c）＞0，直接跳到步驟（4）；否則，令c＝c×s，s＞1，重復步驟（2）。

（4）利用二分法求得方程（18）的在區間［0，c］上的根α。

（5）利用步驟（4）中得到的參數α求得正則化解xα。

3.3 廣義嶺估計

與嶺估計相同，廣義嶺估計也以正則化解的平均均方誤差達到最小作為選擇最優正則化參數的標準，不同的是廣義嶺估計需要確定多個正則化參數αi（i＝1，2，…，N），而Tikhonov正則化和嶺估計只需確定單個正則化參數。由式（13）、式（14）可得廣義嶺估計的正則化解為［13，17］

式中，Μ是以αi（i＝1，2，…，N）為對角線元素的對角陣。正則化解的平均均方誤差為

式（20）對參數αi各自求偏導數

令各自的偏導數為零，解方程即求得使MSE（x＾α）達到最小的各正則化參數

最優正則化參數αi和對應的正則化解可迭代計算得到，具體迭代計算步驟如下：

（2）利用σ2和由式（22）計算得到參數αi；

（3）利用步驟（2）中的參數αi由式（19）得到新的正則化解

4 數值計算與分析

為評價上述向下延拓數學模型和正則化方法的數值精度、效率和穩定性，基于EGM2008地球位模型設計模擬試驗進行驗證。向下延拓采用移去—計算—恢復法（RCR），先由EGM2008位模型的2－2 190階位系數計算得到飛行高度H處（正高）的重力擾動δg（R＋H）觀測值，在觀測值中加入標準差為1.5mGal的白噪聲（模擬觀測值統計信息見表1），移去由GOCO01S衛星位模型［19］2－120階位系數計算的重力擾動得到殘余重力擾動δg′（R＋H），將δg′（R＋H）向下延拓至大地水準面得到延拓值δg′（R）；再恢復參考模型值得到大地水準面上的重力擾動延拓值δg（R），延拓過程中采用第二類Helmert凝集法計算地形質量影響。將由EGM2008的2－2190階位系數計算的大地水準面上的重力擾動作為真實值，比較重力擾動延拓值δg（R）與真實值，可以檢驗延拓模型和正則化解法的可靠性與穩定性并比較其優劣。為模擬真實航空重力測量條件，選定3個飛行高度，分別是H＝2km、3km、4km，數據點格網間隔取2.5′×2.5′，數據點范圍為緯線［30°，33°］和經線［107°，110°］圍成的區域，共包含5 184個數據點；Poisson積分球冠區半徑選取ψ0＝1°，為減弱邊界效應，計算點范圍取緯線［30.5°，32.5°］和經線［107.5°，109.5°］圍成的區域，共包含2 304個計算點。分別采用最小二乘法、Tikhonov正則化、嶺估計和廣義嶺估計求解Poisson積分方程（5），將飛行高度處的重力擾動觀測值向下調和延拓至大地水準面。采用EGM96位模型計算遠區流動點的貢獻，最大階數截斷至Nf＝360階，Tikhonov正則化由GCV法確定正則化參數，嶺估計和廣義嶺估計根據迭代算法確定，將觀測值作為待求參數的迭代初始值。

表1 重力擾動模擬觀測值統計表Tab.1 Statistics of simulated gravity disturbance mGal

表2給出了分別采用最小二乘法和三種正則化方法得到的向下延拓值與真實值的差值統計信息。可以看出，即使對于混入理想化噪聲的模擬重力數據，向下延拓解也是不穩定的，延拓高度越大，延拓值與真實值的偏差越大，且在相同的延拓高度上，不同求解方法的延拓結果相差很大。最小二乘解的均方根偏差（RMS）由3.963mGal（H＝2km）增大至14.880mGal（H＝4km），呈指數上升，表明白噪聲隨著延拓高度的增加被明顯放大，必須引入正則化方法才能獲取穩定延拓解。對比三種正則化算法發現：Tikhonov正則化的去噪效果不明顯，嶺估計能夠在一定程度上抑制白噪聲，但作用有限，在不同延拓高度處，廣義嶺估計所得延拓解的精度和穩定性都是最好的，RMS隨延拓高度的增大幅度也是最小的，當延拓高度為4km時，廣義嶺估計解的RMS僅為2.770mGal，僅占嶺估計所得延拓值RMS（6.325mGal）的44%，占Tikhonov正則化解RMS（10.627mGal）的26%。

以上數值結果都是在標準差為1.5mGal的白噪聲條件下得到的，為檢驗含較大量級噪聲時三種正則化方法的可靠性與精度，分別在重力擾動觀測值中加入標準差為2mGal和3mGal的白噪聲進行向下延拓計算。表3給出了不同噪聲量級時三種方法所得延拓值與真實值的差值統計信息，延拓高度取3km，總體上仍然是廣義嶺估計表現最優，具有最好的穩健性和抗差性。

表2 重力擾動延拓值與真實值的差值統計Tab.2 Statistics of the differences between continued values and true values mGal

表3 不同噪聲量級時重力擾動延拓值與真實值的差值統計Tab.3 Statistics of the differences between continued values and true values mGal

圖1是不同正則化解法所得重力擾動延拓值與真實值的差值的空間分布圖，延拓高度為4km，白噪聲標準差為1.5mGal，從中可以看出，基于廣義嶺估計的延拓值與真實值的偏差均在±10mGal以內，大部分計算點處的差值不超過±3mGal。相比之下，Tikhonov正則化和嶺估計所得重力延拓值與真實值存在較大偏差，特別是Tikhonov正則化的延拓結果嚴重失真，與真實重力場信號相差甚遠，說明延拓過程中噪聲放大嚴重，無法得到可靠的向下延拓解。數值比較結果表明基于多個正則化參數的廣義嶺估計能夠有效抑制測量噪聲，實現重力場信號的穩定向下延拓，在可靠性、精度和抗差性等方面都要顯著優于基于單個正則化參數的Tikhonov正則化和嶺估計。其原因是Tikhonov正則化和嶺估計通過在法方程系數陣主對角線元素上加上同一個常數α來改變法方程的病態特性，有可能出現待求參數被高估或被低估的情況，而廣義嶺估計在MSE（x＾α）最小化條件下迭代直至收斂得到每個對角線元素對應的正則化參數αi（i＝1，2，…，N），每個正則化解x＾αi都對應一個最優的正則化參數αi。

圖1 重力擾動延拓值與真實值的差值Fig.1 The difference between continued values and true values

5 結 論

基于EGM2008地球位模型設計模擬試驗，分別采用Tikhonov正則化、嶺估計和廣義嶺估計三種正則化方法求解逆Poisson積分問題，將飛行高度處含白噪聲的2.5′×2.5′重力擾動向下延拓至大地水準面上，與參考值作外部檢驗。數值比較結果表明，在延拓精度、穩定性和抗差性等方面，廣義嶺估計都要顯著優于Tikhonov法和嶺估計，其原因是廣義嶺估計需要確定多個最優正則化參數，所得估值具有比普通嶺估計估值更小的平均均方誤差，而Tikhonov法和嶺估計只需確定單個最優正則化參數，有可能出現部分待求參數被高估或被低估的情況。需要指出的是，多個正則化參數的迭代求解帶來了廣義嶺估計解法計算效率相對較低的數值問題，當數據量較大時將全部數據點分塊處理可解決這一問題，例如劃分成3°×3°的數據塊，各數據塊之間應保證有一定的重疊區域。

［1］ MARTINEC Z.Stability Investigations of a Discrete Downward Continuation Problem for Geoid Determination in the Canadian Rocky Mountains［J］.Journal of Geodesy，1996，70（11）：805－828.

［4］ LI Y C.Airborne Gravimetry for Geoid Determination［D］.Calgary：University of Calgary，2000.

［5］ KERN M.An Analysis of the Combination and Downward Continuation of Satellite，Airborne and Terrestrial Gravity Data［D］.Calgary：University of Calgary，2003.

［6］ 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.）

［7］ FORSBERG R，KENYON S.Evaluation and Downward Continuation of Airborne Gravity Data—the Greenland Example［C］∥Proc Int Symp Kinematic Systems in Geodesy，Geomatics and Navigation.Banff：［s.n.］，1994：531－538.

［8］ BáLHA T，HIRSCH M，KELLER W，et al.Application of a Spherical FFT Approach in Airborne Gravimetry［J］.Journal of Geodesy，1996，70（11）：663－672.

［9］ HWANG C，HSIAO Y，SHIH H，et al.Geodetic and Geophysical Results from a Taiwan Airborne Gravity Survey：Data Reduction and Accuracy Assessment［J］.Journal of Geophysical Research，2007，112（B4）：1－14.

［10］ SHI Pan，WANG Xingtao.The Frequence Domain Analysis for the Determination of Terrestrial Mean Gravity Anomaly Using Airborne Gravimetry［J］.Acta Geodaetica et Cartographica Sinica，1995，24（4）：301－308.（石磐，王興濤.空中測量地面平均重力異常的頻域分析［J］.測繪學報，1995，24（4）：301－308.）

［11］ SHI Pan，SUN Zhongmiao.The Solution to the Problem of the Spherical Interior Dirichlet and Its Application［J］.Acta Geodaetica et Cartographica Sinica，1999，28（3）：195－198.（石磐，孫中苗.球內Dirichlet問題解及其應用［J］.測繪學報，1999，28（3）：195－198.）

［12］ NOVáK P，HECK B.Downward Continuation and Geoid Determination Based on Band－limited Airborne Gravity Data［J］.Journal of Geodesy，2002，76（5）：269－278.

［13］ XU P，RUMMEL R.Generalized Ridge Regression with Applications in Determination of Potential Fields［J］.Manuscripta Geodaetica，1994，20：8－20.

［14］ HEISKANEN W A，MORITZ H.Physical Geodesy［M］.San Francisco：Freeman and Company，1967.

［15］ TIKHONOV A N.Solution for Incorrectly Formulated Problems and the Regularization Method［J］.Soviet Math Dokl，1963，4：1035－1038.

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

［17］ HOERL A，KENNARD R.Ridge Regression：Biased Estimation for Nonorthogonal Problems［J］.Technometrics，1970，12（1）：55－67.

［18］ XU P.Determination of Surface Gravity Anomalies Using Gradiometric Observables［J］.Geophys J Int，1992，110（2）：321－332.

［19］ PAIL R，GOIGINGER H，SCHUH W D，et al.Combined Satellite Gravity Field Model GOCO01S Derived from GOCE and GRACE［J］.Geophysical Research Letters，2010，37（20）：1－5.

Solution of Ill－posed Problem in Downward Continuation of Airborne Gravity

JIANG Tao1，2，LI Jiancheng1，2，WANG Zhengtao1，2，ZHANG Shoujian1，2
1.School of Geodesy and Geomatics，Wuhan University，Wuhan 430079，China；2.Key Laboratory of Geospace Environment and Geodesy，Ministry of Education，Wuhan University，Wuhan 430079，China

Generalized ridge estimate is suggested for solving the ill－posed problem in downward continuation（DWC）of airborne gravity.Three regularization methods are studied to solve the inverse Poisson problem，which are Tikhonov regularization，ridge estimate and generalized ridge estimate.By conducting a simulated numerical experiment designed based on EGM2008 Earth gravitational model，the2.5′×2.5′gravity disturbances with white noise at the flight levels are downward continued to the Geoid.The accuracy and stability of the DWC methods are tested and compared.The numerical results show that，in terms of continued accuracy，stability and robustness，the generalized ridge estimate based on multiple optimal regularization parameters performs significantly well.

airborne gravity；downward continuation；ill－posed problem；regularization

JIANG Tao（1984—），male，PhD candidate，majors in physical geodesy.

1001－1595（2011）06－0684－06

P223

A

國家自然科學基金（40637034；41074014；40804003；40904003）；地球空間環境與大地測量教育部重點實驗室測繪基礎研究基金（10－02－14）

叢樹平）

2010－11－01

2011－04－12

蔣濤（1984—），男，博士生，研究方向為物理大地測量。

E－mail：tjiang＠whu.edu.cn