Tikhonov正則化方法在GOCE重力場求解中的模擬研究

徐新禹,李建成,王正濤,鄒賢才

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

Tikhonov正則化方法在GOCE重力場求解中的模擬研究

徐新禹1,2,李建成1,2,王正濤1,2,鄒賢才1,2

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

給出四類可用于重力場解算的正則化矩陣(零次、一次、二次和Kaula),以及用于確定正則化參數的L曲線法和GCV方法的數學模型。基于SA方法,利用模擬數據分析討論了零次、一次以及Kaula正則化矩陣應用于GOCE全球重力場模型確定的有效性,并由Kaula正則化矩陣分析L曲線法和GCV方法確定正則化參數的可行性。數值結果表明三類正則化矩陣獲得的最優解(以大地水準面MSE最小為準則確定)的精度水平相近,關鍵在于相應正則化參數的確定,數值結果同時說明了GCV方法和L曲線法可用于確定正則化參數,且前者較后者具有更好的穩定性。

Kaula;正則化;GOCE;GCV;L曲線

1 引 言

衛星大地測量領域的線性問題往往是不適定的(病態的)[1-3],對衛星重力梯度測量任務來說,導致病態的主要原因可歸結為以下幾個方面:第一,重力信號向下延拓通常不穩定,該過程不僅放大觀測誤差,同時也會放大模型誤差(未能用模型描述的信號);第二,觀測數據的不規則分布破壞了原有球諧函數的正交性,例如 GOCE(gravity field and steady-state ocean circulation explorer)任務的PG(polar gap)問題和觀測值間斷等問題都使得觀測值不能全球均勻覆蓋;第三,重力觀測量本身不能完全包含所有頻譜的重力場信息,例如重力梯度張量觀測值的某個分量僅對重力場模型的某些頻段比較敏感,利用其恢復全頻譜的重力場模型往往導致不適定;第四,衛星重力梯度測量的有色噪聲特性也會使法方程系數陣的條件數增大,例如重力梯度儀并不是等精度獲得全頻段的重力信息。GOCE衛星已于2009年3月17日成功發射 (http:∥www.esa.int/SPECIALS/ GOCE),該任務的目標是實現厘米級精度的大地水準面,相應重力異常的精度為 1～2 mGal (1 mGal=10-3cm/s2),同時空間分辨率達到100 km。為實現該目標,研究利用其觀測值恢復重力場的相關理論與方法至關重要,存在許多問題有待于進一步研究和討論,其中病態問題就是一個需要深入研究的課題。

在物理大地測量領域,許多學者對病態問題都曾進行了深入細致的研究[4-9],提出很多計算正則化解和確定最優正則化參數的方法,在利用衛星重力觀測值(衛星跟蹤衛星和衛星重力梯度)求解重力場方面,主要有最小二乘配置法、Tikhonov正則化方法以及有偏估計、基于奇異值分解的方法等,正則化參數的選擇準則主要有均方誤差最小、基本跡估計器、Morozov不一致原理、方差分量估計等。本文主要研究物理大地測量領域廣泛采用的 Tikhonov正則化方法在GOCE全球重力場確定中的應用,利用模擬觀測值基于半解析法(semi-analytical,SA)討論不同正則化矩陣的優劣和最優正則化參數的確定。正則化矩陣主要討論零次、一次、二次和 Kaula正則化矩陣;正則化參數的選擇主要討論L曲線法[10]、廣義交叉檢驗(generalized cross-validation,GCV)[11]。

2 基本原理

2.1 Tikhonov正則化方法

Tikhonov正則化方法是在經典最小二乘準則的基礎上,加入參數向量范數最小的條件,同時引入正則化矩陣和正則化參數,對于 Gauss-Markov線性模型有準則如下[12]:

其中,α是正則化參數;K為正則化矩陣,其為正定對稱矩陣 K=UTU;x為待求參數,y為觀測值; A為觀測方程的設計矩陣;P=,其中為觀測值誤差協方差矩陣為單位權方差。由該準則可求得正則化解如下:

式中,若正則化矩陣 K=I(為單位陣),則Tikhonov正則化就變為有偏估計;如果正則化矩陣的對角線元素取重力場信號的Kaula階方差模型的倒數,其他元素為零,就得到著名的Kaula正則化方法。

式(1)中右邊第一項(殘差平方和最小)要求模型與觀測值有較好的一致性,第二項要求信號的范數最小,信號的范數最小可提高求解的穩定性,同時考慮兩項則需要找到這兩種標準的平衡點,不僅使模型與觀測值有較好的一致性,同時信號的范數較小,因而需要選擇合適的正則化矩陣和確定最優正則化參數,下面介紹幾種物理大地測量領域常用的正則化矩陣和正則化參數的確定方法[9]。

2.1.1 零次正則化矩陣

假設式(1)中正則化矩陣 K為單位陣,則右邊第二項球諧系數的范數可表示為引力位函數的平方在半徑為 R的球面上的積分,進一步推導可得

其中,I為單位矩陣。需要指出,式中的常數因子4π(GM)2在本文的研究中是將其包含在正則化參數中,后面兩個正則化矩陣也是將常數部分包括到正則化參數中。

2.1.2 一次正則化矩陣

根據零次正則化矩陣的推導過程,若將引力位的一階導數作為目標函數,取引力位面梯度的平方函數在半徑為R的球面上的積分,則有

K中的元素kij滿足kij=(δijn(i)(n(i)+1),δij為克羅內克符號;n(i)表示第i行對應的階數。

2.1.3 二次正則化矩陣

類似于一次正則化矩陣,將面Laplace算子作用于引力位函數,并在半徑為 R的球面上積分可得二次正則化矩陣表達式如下:

K中的元素kij滿足kij=δijn2(i)(n(i)+1)2。

2.1.4 K aula正則化

正則化矩陣同樣可基于階方差模型,既可是經驗性的重力信號階方差,也可基于現有的衛星重力模型,根據 Kaula規則[13],重力場信號階次方差滿足

將式(6)取倒數即得 Kaula正則化矩陣,將常數因子1010包括到正則化參數α中,則可得到正則化矩陣元素的形式為 kij=δijn4(i),式中符號的含義同上。從kij的形式可看出其與二次正則化矩陣比較接近,其高階趨近相同。另外,從數值上分析,正則化矩陣中相應低階次的元素乘以正則化參數的量級遠小于矩陣ATPA中的元素,對求解過程的影響基本相同,因而可認為 Kaula正則化矩陣與二次正則化矩陣是一致的。

2.2 最優正則化參數的確定方法

2.2.1 L曲線法

L曲線法的基本思路是將參數的信號范數η(α)=‖xα‖K與平差殘差的范數ρ(α)=‖Axα -y‖P分別取常用對數^η(α)=logη(α),^ρ(α)= logρ(α),并將其繪制成二維曲線圖,經證明對于不連續的病態問題,該曲線通常成L形,該方法定義最優正則化參數對應曲線的最大曲率點(即L曲線的頂點),曲線的曲率可表示為

2.2.2 GCV方法

GCV方法選擇正則化參數的基本思想是:首先將某一個觀測值yk從觀測值序列中除去,并利用剩下的觀測值求解某個正則化參數對應的解x[k]α,再用該解估計被剔除的觀測值,求其與真實觀測值之差Δyk=(Ax[k]α)k-yk,如果某個正則化參數能夠使所有的Δyk都較小,則認為該正則化參數是最優的,最優選擇準則如下[7]:

其中,Qα=A(ATPA+αK)-1ATP為影響矩陣,滿足Axα=Qαy;“tr”表示求矩陣的跡;n表示觀測值個數。考慮到Qα的計算中需要對法方程矩陣乘積和求逆,且Qα為n×n維矩陣,因此式(8)不適合大量觀測數據的情況,為了便于計算,將式(8)進一步變換,由等式tr(AB)=tr(BA)可得

其中,Tα=tr(U(ATPA+αK)-1UT),K=UTU,Tα為u×u的矩陣;u為未知參數個數。

3 數值結果與分析

為分析比較不同正則化矩陣和正則化參數確定方法在 Tikhonov正則化過程中的作用和效果,本文模擬了29 d的 GOCE觀測數據,數據采樣間隔為5 s,選擇EGM96為軌道模擬的參考模型。為使軌道模擬更具真實性,積分時模型的最高階取至250,軌道傾角約為96.7°,是一個近似重復軌道;重力梯度張量觀測值模擬時選擇EIGEN_CG03C模型,取至180階。本文下面的數值分析僅以重力梯度的徑向分量為例,在觀測值中加入根據GOCE任務解析形式誤差功率譜密度基于自回歸(auto regressive,AR)模型模擬的有色噪聲時間序列[14]。

首先分析比較幾類Tikhonov正則化矩陣對求解結果的影響。第二節中已經指明 Kaula正則化矩陣與二次正則化矩陣基本相同,因此這里僅比較零次、一次和 Kaula正則化矩陣。由模擬觀測數據,利用零次、一次和Kaula正則化矩陣,分別選擇不同的正則化參數(在某個區間內變化),基于半解析法(semi-analytical,SA)方法分別恢復180階重力場模型。為了評價解算模型的精度,由其計算緯度-80°×80°范圍內0.5°×0.5°大地水準面格網點值,并與“真”實重力場(這里指EIGEN_CG03C)對應值求差,計算大地水準面誤差RMS。圖1分別給出了零次、一次和Kaula正則化矩陣對應大地水準面誤差RMS隨正則化參數變化的曲線。

圖1 零次、一次和Kaula正則化矩陣對應的大地水準面誤差RMS隨正則化參數變化的曲線Fig.1 The curve of the geoid error RMS with the regularization parameters based on zero-order, first-order and Kaula regularization matrix

圖1中橫坐標表示Kaula正則化矩陣對應的正則化參數α,零次和一次正則化矩陣對應的參數α需要分別乘以106和103,其中三角形符號表示大地水準面誤差RMS最小值對應的點,Kaula正則化對應最小大地水準面誤差RMS的正則化參數為α=1010.7,圓點表示采用GCV方法確定的最優正則化參數所對應的點。從圖中可看出,基于這三類正則化矩陣求得的模型大地水準面誤差RMS隨正則化參數變化的曲線具有相同的趨勢,均隨著正則化參數的增大,先減小后增大;從不同曲線對應大地水準面誤差最小值的量級來看,三類正則化矩陣均能較好達到穩定求解目的,特別是一次和Kaula正則化矩陣,若能選擇最優的正則化參數,可使大地面誤差從幾十厘米降到10 cm左右。從圖中還可看出,采用GCV方法確定最優參數時,一次正則化矩陣獲得的參數值最接近大地水準面誤差最小的點(對應的正則化參數可看作是最優正則化參數),Kaula正則化方法確定的正則化參數對應大地水準面誤差RMS與一次正則化矩陣的最優值基本一致。此外,當正則化參數小于最優正則化參數時,Kaula正則化矩陣對應解的大地水準面誤差變化較慢,而正則化參數較大時,Kaula正則化矩陣對應解的大地水準面誤差變化較快,這說明在SA方法中,Kaula正則化矩陣對較大正則化參數的選擇較零次和一次正則化矩陣更為敏感。

圖1中誤差曲線隨正則化參數的變化趨勢說明,只有選擇合適的正則化參數才能使大地水準面誤差RMS最小。下面以 Kaula正則化矩陣分析不同正則化參數對求解結果的影響。從理論分析可知,正則化過程可看作是一個濾波過程,小正則化參數很難使求解結果趨于穩定,因而導致求解結果誤差較大,而大正則化參數會使求解過程過于穩定,即求解模型將會受到正則化帶來的過強“平滑效應”約束,從而使解過于平滑,尤其是高頻信號。圖2中給出了正則化參數分別選擇α=108、α=1010.7、α=1012時的大地水準面的誤差。從圖中可以看出,α=108時大地水準面誤差在近兩極地區明顯比α=1010.7時大,其他地區則差異不大,說明小的正則化參數對解的穩定性提高有限,仍然暴露出PG問題對求解結果的影響;當α=1012時,大地水準面在重力場高頻信號較明顯的地區(例如高山地區)的誤差較大,反映出高頻信號被過度平滑,同時大的正則化參數給接近兩極的地區也帶來較大的誤差;當α=1010.7時,求解的結果在比較的范圍內大地水準面誤差均較小。

圖2 不同正則化參數下求得模型的大地水準面誤差RMSFig.2 The geoid error RMS of the models estimated with differentα

除了從大地水準面誤差圖中可看出不同正則化參數對重力場模型不同頻段的影響外,從模型的誤差譜圖以及模型的階誤差RMS也可看出相同的特點。如圖3和圖4所示,這里誤差譜圖是將系數誤差的絕對值取常用對數得到的值,n、m分別表示階次。圖中α=108對應解的低次部分系數誤差明顯較大,α=1012時,圖3中(右圖)對應解的高頻部分的誤差普遍比α=1010.7偏大,圖4中模型的階誤差RMS在高于142階時比α=108對應模型的階誤差RMS大,這些都說明大的正則化參數對重力場模型的求解進行了強平滑濾波,從而損失了高頻信號,這與圖2(c)中的結果一致。

圖3 正則化參數不同時求得模型的誤差譜Fig.3 The error spectrum of the models estimated with differentα

圖4 α=108、α=1010.7、α=1012時求得模型的階誤差RMSFig.4 The degree error RMS of the models estimated with differentα

前面分析的結果均說明選擇合適的正則化參數對獲得較好求解結果的重要性,在數值模擬分析時可用大地水準面MSE最小準則采用測試法確定最優的正則化參數,例如前文討論不同正則化矩陣時的比較分析方法,但該準則無法直接應用于實測數據處理,因為無法獲知真實重力場模型,因此必須找到確定最優正則化參數的方法。目前大地測量領域討論較多正則化參數選擇方法的是L曲線法和 GCV方法,本文分別對這兩種方法在利用梯度觀測數據恢復重力場中應用的可行性和有效性進行數值分析討論。采用的模擬數據同前,選擇 Kaula正則化矩陣,分別采用L曲線法和GCV方法確定最優的正則化參數。L曲線如圖5所示,圖中五角星表示L曲線的最大曲率點,對應正則化參數為1011.3,從圖中可看出L曲線的頂點(最大曲率點)并不明顯。GCV函數曲線如圖6所示,圖中五角星表示函數最小值對應的點,正則化參數α為1010.3。

圖5 Kaula正則化矩陣對應的L曲線Fig.5 L-curve according to Kaula matrix

圖6 Kaula正則化矩陣對應的GCV函數曲線Fig.6 The curve of GCV function according to Kaula matrix

基于這兩種方法確定的正則化參數求解模型的大地水準面誤差的統計結果如表1所示。表中RMS的下標±90°和±80°分別表示在計算范圍為緯度 -90°～90°和-80°～80°,經度取0°～360°,表中還給出了由大地水準面誤差的RMS最小為準則確定的最優正則化參數對應模型的大地水準面誤差的統計結果。

表1 GCV方法和L曲線法求解模型的大地水準面誤差RMS統計Tab.1 The statistics of the geoid error RMS of model estimated with GCVand L-curve

從表1中可看出,L曲線方法確定的正則化參數α相對最優值1010.7偏大,GCV方法確定的偏小,兩者均不是理論上的最優。雖然 GCV方法確定模型的精度與最優值更為接近,但兩種方法獲得解的精度與最優解基本在同一量級,這說明了L曲線法和GCV方法在用Kaula正則化矩陣時確定最優正則化參數的可行性和有效性。從數值的穩定性來講,GCV方法相對更優,很容易找到GCV函數曲線的最小值,而L曲線法有些時候會由于曲線過于平滑而難于準確找到頂點,缺乏較好的穩定性。

4 結 論

本文在介紹 Tikhonov正則化方法、正則化矩陣的選擇以及正則化參數確定方法的基礎上,利用GOCE模擬數據基于SA方法分析比較了Tikhonov正則化方法中選擇零次、一次和 Kaula正則化矩陣對求解結果的影響,結果說明若能獲得最優正則化參數,三類矩陣均能獲得較好的結果;基于Kaula正則化矩陣分析比較了確定正則化參數的GCV方法和L曲線法,數值結果說明了兩種方法的可行性。需要指出,本文采用的SA方法是一種快速近似方法,本身具有模型誤差,這對求解結果會有一定影響,因此文中的分析更多的是定性分析不同正則化矩陣以及正則化參數選擇方法的優劣,但總的來說,模擬結果已經說明在利用GOCE觀測數據求解高階全球重力場時采用 Kaula正則化矩陣的有效性,以及采用GCV方法確定最優正則化參數的可行性。

[1] RUMMEL R.Determination of Short-wavelength Components of the Gravity Field from Satellite-to-Satellite Tracking or Satellite Gradiometry:An Attempt to an Identification of Problem Areas[J].Manuscripta Geodaetica,1979,4(2): 107-148.

[2] BOUMAN J,KOOP R.Quality Differences between Tikhonov Regularization and Generalized Biased Estimation in Gradiometric Analysis[J].DEOS Progress Letters,1997,97(1): 42-48.

[3] WANG Zhenjie.Research on the Regularization Solutions of Ill-posed Problems in Geodesy[D].Wuhan:Institute of Geodesy and Geophysics,Chinese Academy of Sciences, 2003.(王振杰.大地測量中不適定問題的正則化解法研究[D].武漢:中國科學院測量與地球物理研究所,2003.)

[4] XU P L.Determination of Surface Gravity Anomalies Using Gradiometric Observables[J].Geophys J Int,1992,110(2): 321-332.

[5] IL K K H.Regularization for High Resolution Gravity Field Recovery by Future Satellite Techniques[C]∥ANGER G,et al.Inverse Problems:Principles and Applications in Geophysics, Technology and Medicine:74.Berlin:Akademie Verlag, 1993:189-214.

[6] KUSCHE J,KLEES R.On the Regularization Problem in Gravity Field Determination from SatelliteGradiometric Data[C]∥áDáM J,SCHWARZ K.Vistas for Geodesy in the New Millennium:IGA 2001 Scientific Assembly.Budapest:Springer,2001:175-180.

[7] KUSCHEJ,KLEES R.Regularization of the Gravity Field Estimation from Satellite Gravity Gradients[J].Journal of Geodesy,2002,76(6-7):359-368.

[8] YANG Y X,SONG L,XU T H.Robust Estimator for Correlated Observations Based on Bifactor Equivalent Weights[J].Journal of Geodesy,2002,76 (6-7): 353-358.

[9] DITMAR P,KUSCHE J,KLEES R.Computation of Spherical Harmonic Coefficients from Gravity Gradiometry Data to be Acquired by the GOCE Satellite:Regularization Issues[J].Journal of Geodesy,2003,77(7-8):465-477.

[10] HANSEN P C.The L-curve and Its Use in the Numerical Treatment of Inverse Problems[C]∥JOHNSTON P. ComputationalInverse Problems in Electrocardiology. Southampton:WIT Press,2001:119-142.

[11] 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.

[12] TIKHONOV A N.Regularization of Ill-posed Problem [J].Dokl Akad Nauk SSSR,1963,151(1):49-52.

[13] KAULA W M.Theory ofSatelliteGeodesy[M]. Waltham:Blaisdell Publishing Company,1966:98.

[14] XU Xinyu,WANG Zhengtao,ZOU Xiancai,et al.The Research on theGOCE SatelliteGravity Gradiometry Error Analysis and Simulation[J].Journal of Geodesy and Geodynamics,2010,30(1):1-5.(徐新禹,王正濤,鄒賢才,等.GOCE衛星重力梯度測量誤差分析及其模擬研究[J].大地測量與地球動力學,2010,30(1):1-5.)

(責任編輯:叢樹平)

The Simulation Research on the Tikhonov Regularization Applied in Gravity Field Determination of GOCE Satellite Mission

XU Xinyu1,2,LI Jiancheng1,2,WANG Zhengtao1,2,ZOU Xiancai1,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

The Tikhonov regularization is widely applied in the geodesy,the principle of which is discussed in this paper,including the mathematical models of four types of regularization matrices(zero-order,first-order,secondorder and Kaula)and the regularization parameter selection methods:L-curve and GCV.The validation of zeroorder,first-order and Kaula regularization matrices applied in the gravity field determination with GOCE simulated data is analyzed based on the SA method.And the applicability of L-curve and GCV is also discussed using the simulated data.The results show that the accuracies of the optimized solutions(selected by minimizing geoid MSE) with the three types of regularization matrices are at the same level.The key point is the selection of the corresponding regularization parameter.The results also show that GCV and L-curve can be applied in the regularization parameter estimation,and the former method is more stable than the latter one.

Kaula;regularization;GOCE;GCV;L-curve

XU Xinyu(1978—),male,PhD,lecturer, majors in satellite geodesy.

E-mail:xyxu@sgg.whu.edu.cn

1001-1595(2010)05-0465-06

P223

A

國家自然科學基金(40904003,40637034,40704004);國家863計劃(2006AA12Z309)

2009-11-02

2010-03-15

徐新禹(1978—),男,博士,講師,研究方向為衛星大地測量。