Tikhonov正則化方法在GOCE重力場(chǎng)求解中的模擬研究

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

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

Tikhonov正則化方法在GOCE重力場(chǎng)求解中的模擬研究

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

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

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

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

1 引 言

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

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

2 基本原理

2.1 Tikhonov正則化方法

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

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

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

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

2.1.1 零次正則化矩陣

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

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

2.1.2 一次正則化矩陣

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

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

2.1.3 二次正則化矩陣

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

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

2.1.4 K aula正則化

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

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

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

2.2.1 L曲線法

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

2.2.2 GCV方法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié) 論

本文在介紹 Tikhonov正則化方法、正則化矩陣的選擇以及正則化參數(shù)確定方法的基礎(chǔ)上,利用GOCE模擬數(shù)據(jù)基于SA方法分析比較了Tikhonov正則化方法中選擇零次、一次和 Kaula正則化矩陣對(duì)求解結(jié)果的影響,結(jié)果說明若能獲得最優(yōu)正則化參數(shù),三類矩陣均能獲得較好的結(jié)果;基于Kaula正則化矩陣分析比較了確定正則化參數(shù)的GCV方法和L曲線法,數(shù)值結(jié)果說明了兩種方法的可行性。需要指出,本文采用的SA方法是一種快速近似方法,本身具有模型誤差,這對(duì)求解結(jié)果會(huì)有一定影響,因此文中的分析更多的是定性分析不同正則化矩陣以及正則化參數(shù)選擇方法的優(yōu)劣,但總的來說,模擬結(jié)果已經(jīng)說明在利用GOCE觀測(cè)數(shù)據(jù)求解高階全球重力場(chǎng)時(shí)采用 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.(王振杰.大地測(cè)量中不適定問題的正則化解法研究[D].武漢:中國(guó)科學(xué)院測(cè)量與地球物理研究所,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)星重力梯度測(cè)量誤差分析及其模擬研究[J].大地測(cè)量與地球動(dòng)力學(xué),2010,30(1):1-5.)

(責(zé)任編輯:叢樹平)

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

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

2009-11-02

2010-03-15

徐新禹(1978—),男,博士,講師,研究方向?yàn)樾l(wèi)星大地測(cè)量。