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

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

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

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

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

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

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

1 引 言

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

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

2 基本原理

2.1 Tikhonov正則化方法

Tikhonov正則化方法是在經(jīng)典最小二乘準則的基礎(chǔ)上,加入?yún)?shù)向量范數(shù)最小的條件,同時引入正則化矩陣和正則化參數(shù),對于 Gauss-Markov線性模型有準則如下[12]:

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

式中,若正則化矩陣 K=I(為單位陣),則Tikhonov正則化就變?yōu)橛衅烙?如果正則化矩陣的對角線元素取重力場信號的Kaula階方差模型的倒數(shù),其他元素為零,就得到著名的Kaula正則化方法。

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

2.1.1 零次正則化矩陣

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

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

2.1.2 一次正則化矩陣

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

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

2.1.3 二次正則化矩陣

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

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

2.1.4 K aula正則化

正則化矩陣同樣可基于階方差模型,既可是經(jīng)驗性的重力信號階方差,也可基于現(xiàn)有的衛(wèi)星重力模型,根據(jù) Kaula規(guī)則[13],重力場信號階次方差滿足

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

2.2 最優(yōu)正則化參數(shù)的確定方法

2.2.1 L曲線法

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

2.2.2 GCV方法

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

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

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

3 數(shù)值結(jié)果與分析

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

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

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

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

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

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

圖3 正則化參數(shù)不同時求得模型的誤差譜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α

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

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

圖6 Kaula正則化矩陣對應(yīng)的GCV函數(shù)曲線Fig.6 The curve of GCV function according to Kaula matrix

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

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

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

4 結(jié) 論

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

[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].武漢:中國科學(xué)院測量與地球物理研究所,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衛(wèi)星重力梯度測量誤差分析及其模擬研究[J].大地測量與地球動力學(xué),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

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

2009-11-02

2010-03-15

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