衛(wèi)星重力徑向梯度數(shù)據(jù)的最小二乘配置調(diào)和分析

吳 星,張傳定,劉曉剛

1.總裝備部工程設(shè)計(jì)研究總院,北京100028;2.海洋測(cè)繪研究所,天津300061;3.信息工程大學(xué)測(cè)繪學(xué)院,河南鄭州450052

衛(wèi)星重力徑向梯度數(shù)據(jù)的最小二乘配置調(diào)和分析

吳 星1,2,張傳定3,劉曉剛3

1.總裝備部工程設(shè)計(jì)研究總院,北京100028;2.海洋測(cè)繪研究所,天津300061;3.信息工程大學(xué)測(cè)繪學(xué)院,河南鄭州450052

研究利用衛(wèi)星重力梯度徑向分量確定地球引力場(chǎng)位系數(shù)的最小二乘配置(LSC)調(diào)和分析方法。首先論述最小二乘配置法的原理,推導(dǎo)擾動(dòng)引力梯度觀測(cè)量與球諧系數(shù)之間的協(xié)方差和自協(xié)方差矩陣。在擾動(dòng)引力梯度觀測(cè)數(shù)據(jù)為等經(jīng)差規(guī)則網(wǎng)格數(shù)據(jù)的情況下,引力位與擾動(dòng)引力梯度之間的協(xié)方差矩陣具有分塊Toeplitz循環(huán)陣的結(jié)構(gòu),有效地利用FFT變換技術(shù)將其降階。研究利用截?cái)嗥娈愔捣纸夥?TSVD)減弱協(xié)方差陣的病態(tài)性。最后得到引力梯度徑向分量的最小二乘配置調(diào)和分析的完整計(jì)算公式。模擬試算結(jié)果表明,基于TSVD的最小二乘配置調(diào)和分析方法,能夠以較高的精度還原全球重力場(chǎng),驗(yàn)證本文算法的有效性和實(shí)用性。

地球重力場(chǎng)模型;最小二乘配置;衛(wèi)星重力梯度;調(diào)和分析;協(xié)方差陣;Toeplitz矩陣;TSVD

1 引 言

國(guó)際上首顆重力梯度測(cè)量衛(wèi)星 GOCE已經(jīng)于2009年3月17日成功發(fā)射,其科學(xué)目標(biāo)是在100 km的空間分辨率下,實(shí)現(xiàn)大地水準(zhǔn)面精度達(dá)厘米級(jí)、重力異常精度達(dá)1～2 mGal的地球重力場(chǎng)(1 mGal=10-3cm/s2)。如何利用重力梯度觀測(cè)數(shù)據(jù)恢復(fù)地球重力場(chǎng)早已得到了廣泛研究,目前主要有兩大類方法:時(shí)域法和空域法,兩類方法各有優(yōu)缺點(diǎn)[1-2]?？沼蚍ㄖ械淖钚《伺渲梅ㄗ鳛榇_定全球或區(qū)域地球重力場(chǎng)的主要方法之一備受關(guān)注,其理論與方法也得到了完善和補(bǔ)充[3-4]。羅志才提出的頻域最小二乘配置法可以有效地提高計(jì)算速度[5],張傳定提出的最小二乘復(fù)數(shù)配置法解決了重力梯度張量水平分量的最小二乘配置問(wèn)題[6]。Tscherning和Sanso研究了能夠處理各類重力徑向分量的最小二乘配置法和快速配置法[7-8];汪海洪等針對(duì)不同分辨率的數(shù)據(jù)融合問(wèn)題,提出多分辨率最小二乘配置法[9];Arabelos等和Balmino分別研究了誤差協(xié)方差矩陣的計(jì)算問(wèn)題[10-11];Kotsakis提出具有協(xié)方差匹配約束的最小二乘配置法[12];Reguzzoni等提出最優(yōu)多步配置法[13]。這些研究工作都是利用了最小二乘配置法能綜合各類觀測(cè)數(shù)據(jù)確定地球重力場(chǎng),并且在數(shù)據(jù)處理中能顧及觀測(cè)量的誤差和計(jì)算待估量的誤差協(xié)方差矩陣的特點(diǎn)。但是,配置法要求有可靠的先驗(yàn)協(xié)方差模型,需要求解協(xié)方差矩陣的逆矩陣。在局部重力場(chǎng)逼近中,觀測(cè)量協(xié)方差矩陣維數(shù)不高,可直接求解,因此最小二乘配置在局部重力場(chǎng)逼近得到了廣泛應(yīng)用。在全球重力場(chǎng)逼近中由于協(xié)方差矩陣維數(shù)較大,必須對(duì)其進(jìn)行合理分解,才能給出穩(wěn)定解。

為此,本文主要研究衛(wèi)星重力徑向梯度規(guī)則網(wǎng)格觀測(cè)數(shù)據(jù)的最小二乘配置調(diào)和分析的算法,推導(dǎo)引力位系數(shù)與格網(wǎng)平均擾動(dòng)引力梯度間的協(xié)方差矩陣,研究協(xié)方差矩陣降階求逆方法及其數(shù)值穩(wěn)定性,給出仿真試驗(yàn)結(jié)果。

2 解算重力位系數(shù)的最小二乘配置法

衛(wèi)星重力梯度測(cè)量可以利用衛(wèi)星上的重力梯度儀直接測(cè)定重力位梯度張量的所有分量或部分分量。按照物理大地測(cè)量的傳統(tǒng)作法,首先選擇一合適的參考重力場(chǎng),將觀測(cè)量轉(zhuǎn)換為擾動(dòng)重力位梯度張量,然后將觀測(cè)數(shù)據(jù)歸算到以地球質(zhì)心為球心,衛(wèi)星標(biāo)稱球帶空間的平均球面上,并進(jìn)行網(wǎng)格化處理,將得到的球面規(guī)則網(wǎng)格平均擾動(dòng)引力梯度張量數(shù)據(jù)作為已測(cè)信號(hào),利用最小二乘配置法求解擾動(dòng)重力位系數(shù)。

2.1 最小二乘配置解法的原理

為了簡(jiǎn)便,假設(shè)不含非隨機(jī)參數(shù)向量,此時(shí)最小二乘配置的數(shù)學(xué)模型簡(jiǎn)化為[14]

式中,l為觀測(cè)值向量;t是已測(cè)信號(hào)向量;s為待估信號(hào)向量;I為單位陣;v為隨機(jī)的觀測(cè)誤差向量。

隨機(jī)模型為

其中,Ctt為已測(cè)信號(hào)協(xié)方差陣;Cvv為觀測(cè)誤差協(xié)方差陣。

利用最小二乘配置準(zhǔn)則,可得未測(cè)信號(hào)s的推估值和誤差協(xié)方差陣分別為

其中,Cst為已測(cè)信號(hào)和待估信號(hào)的互協(xié)方差陣。

顯然,對(duì)于由衛(wèi)星重力梯度數(shù)據(jù)解算地球重力場(chǎng)重力位系數(shù)的情況下,已測(cè)信號(hào)為衛(wèi)星重力梯度觀測(cè)信息,待估信號(hào)為所求的重力位系數(shù)。由此可見(jiàn),利用最小二乘配置法求解位系數(shù)的關(guān)鍵問(wèn)題是協(xié)方矩陣與誤差矩陣之和的求逆問(wèn)題。特別是已測(cè)信號(hào)的數(shù)據(jù)量較大時(shí),逆矩陣的求解方法及其穩(wěn)定性將直接影響配置結(jié)果的可靠性。

2.2 二階徑向擾動(dòng)重力梯度的級(jí)數(shù)展開(kāi)式

地球外部任意一點(diǎn) P(r,θ,λ)處擾動(dòng)位 T(r, θ,λ)的球諧函數(shù)級(jí)數(shù)表達(dá)式為[6]

當(dāng)且僅當(dāng)是規(guī)則球面點(diǎn)值或均值時(shí),協(xié)方差矩陣才具有可分解的結(jié)構(gòu)。為此,采用習(xí)用球面劃分[15-16],即把整個(gè)球面用經(jīng)、緯線劃分成經(jīng)度差為Δ θ、緯度差為Δ λ的 N×Nl個(gè)網(wǎng)格,其中, N、Nl分別為緯度方向和經(jīng)度方向的網(wǎng)格分化數(shù),Δ θ=π/N,Δ λ=2π/Nl。同時(shí)假設(shè)衛(wèi)星軌道面上的引力梯度數(shù)據(jù)經(jīng)過(guò)必要的歸算后得到全球覆蓋(非極軌衛(wèi)星,兩極空白地區(qū)用模型填充)的網(wǎng)格平均引力梯度觀測(cè)值ˉfi,k(i、k代表該網(wǎng)格所在的位置,i=0,1,…,N-1;k=0,1,…,Nl-1)。則有

式中,函數(shù)sinc(x)的表達(dá)式為

并且可以通過(guò)遞推公式計(jì)算[6,17]。

2.3 擾動(dòng)引力位系數(shù)與擾動(dòng)引力梯度的協(xié)方差

在研究地球重力場(chǎng)有關(guān)問(wèn)題時(shí),擾動(dòng)位 T之間的協(xié)方差是最基本的,其他擾動(dòng)場(chǎng)元之間的協(xié)方差函數(shù)都可以由它求得。設(shè) T(P)和 T(Q)是空間兩點(diǎn) P(r,θ,λ)和Q(r′,θ′,λ′)的擾動(dòng)位,ψ為P、Q兩點(diǎn)球面角距,則擾動(dòng)位空間協(xié)方差函數(shù)可表示為[14]

式中,Pn(cosψ)為勒讓德多項(xiàng)式;符號(hào)“＊”表示復(fù)數(shù)共軛;Kn為擾動(dòng)位階方差,一般可采用參考引力位模型系數(shù)計(jì)算,計(jì)算式為

另外,可用重力異常階方差模型求得,其關(guān)系式為

例如采用如下6參數(shù)重力異常階方差模型—T/R模型[18]

其中,α1、α2、S1、S2為實(shí)常數(shù);A、B為整常數(shù)。取α1=3.405,α2=140.03,S1=0.998 006,S2= 0.914 232,A=1,B=2。α1、α2的量綱是(mGal)2(1 mGal=10-3cm/s2),其他量是無(wú)量綱的數(shù)。將式(12)代入式(11)并代入式(10)后,完成級(jí)數(shù)求和,即得全球協(xié)方差的解析函數(shù)。

設(shè)Si為任一重力場(chǎng)元信號(hào)與擾動(dòng)位的線性泛函算子,則信號(hào)SiT(P)和 SjT(Q)間的協(xié)方差函數(shù)可由擾動(dòng)位的協(xié)方差函數(shù)表示[4],即有

根據(jù)球諧函數(shù)正交性,由式(5)可得

把式(9)代入式(16),并且令 dnm=Kn/((2n+1)· (2-δ0m)),則可得

這說(shuō)明了擾動(dòng)引力位系數(shù)的協(xié)方差是一對(duì)角陣。

由協(xié)方差函數(shù)傳播公式(13)和式(17),可得擾動(dòng)引力位系數(shù)與格網(wǎng)平均擾動(dòng)引力梯度徑向分量(Q)間互協(xié)方差為

格網(wǎng)平均擾動(dòng)引力梯度徑向分量間的自協(xié)方差為

從式(19)可以看出,對(duì)于等經(jīng)差規(guī)則網(wǎng)格數(shù)據(jù)而言,其第i個(gè)緯度圈與第i′個(gè)緯度圈間的協(xié)方差子矩陣是一個(gè)T oeplitz循環(huán)陣,利用傅里葉變換陣可將其化算為對(duì)角陣,降低矩陣求逆的維數(shù),提高計(jì)算結(jié)果的穩(wěn)定性和可靠性。正是利用這一特性,才使全球協(xié)方差矩陣的穩(wěn)定求逆成為可能。

2.4 最小二乘配置調(diào)和分析解

將規(guī)則擾動(dòng)引力梯度徑向分量觀測(cè)數(shù)據(jù)作以下排列

其中,ˉfi是第i(i=0,1,2,…,N-1)個(gè)緯度圈上所有觀測(cè)量所組成行向量

則協(xié)方差矩陣Ctt可分解成N2個(gè) Nl×Nl維的子矩陣,用表示(i,i′=0,1,…,N-1),其元素是所有第i行和第i′行平均引力梯度之間協(xié)方差函數(shù)值。另外,對(duì)于近圓極軌道衛(wèi)星,由于同一緯度圈上,每個(gè)子塊面積相同,衛(wèi)星重力梯度測(cè)量中一個(gè)完整覆蓋周期內(nèi)落入子塊內(nèi)的測(cè)量點(diǎn)數(shù)大致相同,子塊平均值的誤差方差大致相同[6]。第 i個(gè)緯度帶的誤差方差用表示,則誤差方差陣 Cvv同樣可分解成N2個(gè) Nl×Nl維的子矩陣,用表示(i,i′=0,1,…,N-1),并有

因此由(3)式可知,最小二乘配置調(diào)和分析解為

其中,

如此,對(duì)于 N2個(gè)子對(duì)角陣Dii′,若追蹤 m指標(biāo),則共有Nl個(gè)N×N維的滿秩矩陣

對(duì)這 Nl個(gè)Rm求逆,再按照對(duì)應(yīng)項(xiàng)完成矩陣的乘法運(yùn)算后,得到位系數(shù)的估值為

相應(yīng)地,由式(4)、式(22)與式(25),得位系數(shù)估值誤差的方差為

2.5 矩陣TSVD正則化求解

盡管利用分塊循環(huán)陣解決了協(xié)方差陣的降階求逆,但由于兩極地區(qū)子午線收斂問(wèn)題,即極區(qū)網(wǎng)格之間的球面角距ψ變化范圍較小,將導(dǎo)致矩陣Rm存在病態(tài)性問(wèn)題。實(shí)際數(shù)據(jù)表明,對(duì)于50′× 50′或40′×40′分辨率網(wǎng)格,矩陣 Rm的奇異值數(shù)值非常小,并接近于階梯型分布。為了解算該類病態(tài)方程,本文采用截?cái)嗥娈愔捣纸?TSVD)正則化方法[21]。該方法是把容易造成不穩(wěn)定的較小的奇異值直接截去,使原來(lái)的不適定問(wèn)題轉(zhuǎn)化為一個(gè)適定問(wèn)題來(lái)求解。

設(shè)矩陣Rm的奇異值分解式(SVD)為

其中,正交矩陣U=(u0,u1,…,uN-1),V=(v0,v1,…,vN-1),Σ=diag(σ0,…,σN-1),奇異值滿足σ0≥…≥σN-1>0。

則式(32)可表示為

由式(35),可得最小二乘配置的 TSVD正則化解為

其中,整數(shù)k稱為截?cái)鄥?shù),也就是正則化參數(shù)。

正則化參數(shù)的選取通常有容差原理法、廣義交叉驗(yàn)證法,還有L曲線法等[22]。容差原理法需要知道數(shù)據(jù)項(xiàng)中噪聲大小,而實(shí)際問(wèn)題中噪聲水平是未知的;廣義交叉驗(yàn)證法操作較為煩瑣,且結(jié)果不穩(wěn)定;L曲線法的關(guān)鍵是定位L曲線的最大曲率的那個(gè)點(diǎn),選擇所對(duì)應(yīng)的k值作為截?cái)鄥?shù)。L曲線的拐角點(diǎn)所對(duì)應(yīng)的解不但平衡了解范數(shù)和殘差范數(shù),而且趨于平衡正則化誤差和擾動(dòng)誤差。本文采用L曲線法確定正則化參數(shù),具體方法詳見(jiàn)[21,23]。

3 仿真試算與分析

3.1 協(xié)方差函數(shù)的計(jì)算

全球協(xié)方差函數(shù)模型采用兩種方法計(jì)算,第一種方法是采用當(dāng)前最新、精度最高的地球重力場(chǎng)模型EGM2008計(jì)算[23]。EGM2008模型是美國(guó)國(guó)家地理空間情報(bào)局在充分利用包括地面重力、衛(wèi)星測(cè)高、衛(wèi)星重力等最新數(shù)據(jù)基礎(chǔ)上研制出的新一代地球重力場(chǎng)模型。GPS水準(zhǔn)點(diǎn)外部檢核結(jié)果表明,EGM2008模型具有很高的精度,在全球范圍內(nèi)標(biāo)準(zhǔn)差達(dá)到13 cm[23],在我國(guó)高程異常的總體精度也高達(dá)20 cm[24]。第二種方法是采用T/R全球重力異常階方差模型轉(zhuǎn)換得到。擾動(dòng)位的階方差統(tǒng)計(jì)如圖1所示。

圖1 擾動(dòng)位階方差統(tǒng)計(jì)圖Fig.1 Chart of potential degree-variance

3.2 擾動(dòng)引力梯度徑向分量數(shù)據(jù)的仿真計(jì)算

利用截至300階次的 EGM2008模型,采用文獻(xiàn)[25]中球諧綜合方法分別模擬了全球分辨率分別為40′×40′和50′×50′的格網(wǎng)擾動(dòng)引力徑向梯度徑向分量數(shù)據(jù),衛(wèi)星軌道高度設(shè)為250 km。如表1所示,共產(chǎn)生了4套數(shù)據(jù),即數(shù)據(jù)1～4。數(shù)據(jù)1、2是由2～300階模型生成的,而數(shù)據(jù)3、4模擬生成時(shí),截掉了低階部分(2-50階)的影響。由文獻(xiàn)[25]可知,數(shù)據(jù) A1～A4的精度優(yōu)于10-7mE(1 mE=10-12s-2),其誤差可以忽略不計(jì)。為了使模擬觀測(cè)值具有實(shí)測(cè)性,需要向模擬梯度觀測(cè)值中加入觀測(cè)誤差。根據(jù)衛(wèi)星重力梯度的觀測(cè)精度、采樣間隔、飛行時(shí)間等不同,可得不同精度的重力梯度格網(wǎng)平均值,計(jì)算公式為[6]

其中,NT是整個(gè)飛行期的采樣點(diǎn)數(shù);Nik為落入第(i,k)網(wǎng)格的點(diǎn)數(shù)。

為簡(jiǎn)單起見(jiàn),本文假設(shè)得到擾動(dòng)引力梯度張量模擬觀測(cè)值是等精度觀測(cè)的,而且全球均勻分布,假定得到的40′×40′格網(wǎng)平均值的精度為1.0 mE,則根據(jù)式(37)可知,和50′×50′格網(wǎng)平均值的精度為0.84 mE,故數(shù)據(jù)B是在數(shù)據(jù)A的基礎(chǔ)上分別均加入了標(biāo)準(zhǔn)差為1.0 mE和0.84 mE的零均值白噪聲而產(chǎn)生的。

表1 EGM2008模擬的衛(wèi)星重力徑向梯度數(shù)據(jù)Tab.1 Radialsatellite gravity gradientssimulated using EGM2008

3.3 調(diào)和仿真分析計(jì)算

以計(jì)算位系數(shù)誤差的階均方根(error degree RMS)

為統(tǒng)計(jì)量對(duì)最小二乘配置調(diào)和分析方法進(jìn)行評(píng)估分析。式(38)中,(,)為位系數(shù)估值, (,)表示EGM2008模型位系數(shù)。

由數(shù)據(jù)1～4采用本文的最小二乘配置調(diào)和分析法(LSCHA),協(xié)方差函數(shù)利用截至360階次的EGM2008模型計(jì)算,所得統(tǒng)計(jì)結(jié)果見(jiàn)圖2。

圖2 最小二乘配置調(diào)和分析結(jié)果Fig.2 Least-squares collocation harmonic analysis of radial gravity gradients

在相同數(shù)據(jù)的情況下,采用截?cái)嘀?60階的T/R全球重力異常階方差模型轉(zhuǎn)換得全球協(xié)方差函數(shù),運(yùn)用LSCHA方法,所得結(jié)果與圖1幾乎完全一致。分析圖2可知:

1.對(duì)于截掉低階部分影響的數(shù)據(jù)3、4的計(jì)算結(jié)果要優(yōu)于數(shù)據(jù)1、2的計(jì)算結(jié)果。產(chǎn)生這一差異的原因是在數(shù)據(jù)1、2中考慮了重力場(chǎng)的低頻部分,由于低頻部分與高頻部分量級(jí)相差較大,這影響了高階位系數(shù)的計(jì)算結(jié)果。重力梯度數(shù)據(jù)反映了重力場(chǎng)的中頻特征,不能分辨重力場(chǎng)的低頻部分。

2.分辨率為40′×40′全球重力梯度數(shù)據(jù)比50′×50′的數(shù)據(jù)的調(diào)和分析結(jié)果有顯著的提高。主要原因是實(shí)際地球重力場(chǎng)是連續(xù)分布的,當(dāng)采用一定分辨率的離散格網(wǎng)數(shù)據(jù)表示時(shí),存在數(shù)據(jù)的離散化誤差。因此采用更高分辨率的觀測(cè)數(shù)據(jù)表示,可以減小離散化誤差,提高調(diào)和分析精度。

3.即使是采用精度較高的數(shù)據(jù) A1～A4,其最小二乘配置調(diào)和分析結(jié)果仍然均與EGM2008模型存在差異,即不能100%地一致恢復(fù)位系數(shù),這主要是由下列因素引起的:①由于協(xié)方差矩陣存在病態(tài)性問(wèn)題,其求逆勢(shì)必存在一定的計(jì)算誤差,即使采用正則化方法也不能完全避免;②需要采用更高分辨率的數(shù)據(jù),以進(jìn)一步減小離散化誤差;③配置方法是重力場(chǎng)數(shù)據(jù)處理的重要手段,但它畢竟是統(tǒng)計(jì)方法,能給出信號(hào)的估值,但未必能給出其真值,存在誤差是容許的。

在數(shù)據(jù)相同的情況下,采用了文獻(xiàn)[27]中的輪胎調(diào)和分析方法(THA)進(jìn)行了仿真試算,結(jié)果如圖3所示。

圖3 輪胎調(diào)和分析結(jié)果Fig.3 Torus harmonic analysis of radial gravity gradients

比較圖2和圖3可知,本文提出的基于 TS-VD的最小二乘配置調(diào)和分析結(jié)果略優(yōu)于輪胎調(diào)和分析結(jié)果,尤其是在觀測(cè)數(shù)據(jù)精度較高的情況下(數(shù)據(jù)A1～A4),更為明顯,而當(dāng)采用數(shù)據(jù)B1～B4時(shí),兩種方法所得結(jié)果幾乎一致,但輪胎調(diào)和分析的計(jì)算速度更快。

4 結(jié) 論

通過(guò)本文的理論推導(dǎo)和數(shù)值試驗(yàn)分析可得如下結(jié)論:

1.在重力梯度數(shù)據(jù)滿足一定條件(等經(jīng)差網(wǎng)格分布,以及同緯度帶具有相同誤差)的情況下,其協(xié)方差子矩陣是 Toeplitz循環(huán)矩陣,并可利用傅里葉變換陣將其化算為對(duì)角陣,降低矩陣求逆的維數(shù),從而提高計(jì)算結(jié)果的穩(wěn)定性和可靠性。正是這種特性才使得超大型矩陣的求逆成為可能,從而使最小二乘配置法真正用于恢復(fù)全球重力場(chǎng)。

2.由EGM2008模型和 T/R模型計(jì)算得到的擾動(dòng)位階方差在低階部分吻合較好。

3.本文基于 TSVD的最小二乘配置調(diào)和分析方法計(jì)算精度略優(yōu)于輪胎調(diào)和分析方法,但計(jì)算速度比輪胎調(diào)和分析方法慢。

4.由于兩極地區(qū)子午線收斂問(wèn)題,導(dǎo)致了全球協(xié)方差矩陣病態(tài)性問(wèn)題,采用截?cái)嗥娈愔捣纸庹齽t化方法能夠有效地減弱法矩陣的病態(tài)性。

[1] KOOP R.Global Gravity Field Modelling Using Satellite Gravity Gradiometry[R]. Amsterdam: Netherlands Geodetic Commission,1992.

[2] RUMMEL R,VAN GELDEREN M,KOOP R,et al. Spherical Harmonic Analysis of Satellite Gradiometry[R]. Amsterdam:Netherlands Geodetic Commission,1993.

[3] KRARUP T.A Contribution to the Foundation of Physical Geodesy[R].Copenhagen:Danish Geodetic Institute,1969.

[4] MORTIZ H.Advanced Physical Geodesy[M].England: Abacus Press,1980.

[5] LUO Zhicai.Theories and Method of the Determination of the Earth’s Gravity Field Using Satellite Gravity Gradient Data[D].Wuhan:Wuhan Technical University of Surveying and Mapping,1996.(羅志才.利用衛(wèi)星重力梯度數(shù)據(jù)確定地球重力場(chǎng)的理論與方法[D].武漢:武漢測(cè)繪科技大學(xué),1996.)

[6] ZHANG Chuanding.Satellite-borne Gravity Measurement: Foundation,Modeling Methods,and Data Processing Algorithms[D].Zhengzhou:Information Engineering University,2000.(張傳定.衛(wèi)星重力測(cè)量:基礎(chǔ)、模型化方法與數(shù)據(jù)處理算法[D].鄭州:信息工程大學(xué),2000.)

[7] TSCHERNING C C.Computation of Spherical Harmonic Coefficients and Their Error Estimate Using Least-squares Collocation[J].Journal of Geodesy,2001,75(1):12-18.

[8] SANSO F,TSCHERNING C C.Fast Spherical Collocation:Theory and Examples[J].Journal of Geodesy,2003, 77(1-2):101-112.

[9] WANG Haihong,LUO Zhicai,LUO Jia,et al.Preliminary Research on Multiresolution Least Square Collocation [J].Journal of Geodesy and Geodynamics,2006,26(1): 115-118.(汪海洪,羅志才,羅佳,等.多分辨率最小二乘配置法初探[J].大地測(cè)量與地球動(dòng)力學(xué),2006,26(1): 115-118.)

[10] ARABELOS D N,TSCHERNING C C.Error-covariances of the Estimates of Spherical Harmonic Coefficients Computed by LSC,Using Second-order Radial Derivative Functional Associated with Realistic GOCE Orbits[J]. Journal of Geodesy,2009,83(5):419-430.

[11] BALMINO G.Efficient Propagation of Error Covariance Matrices of Gravitational Models:Application to GRACE and GOCE[J].Journal ofGeodesy,2009,83(10): 989-995.

[12] KOTSAKIS C.Least-squares Collocation with Covariance-matching Constrains[J].Journal of Goedesy,2007, 81(10):661-677.

[13] REGUZZIONI M,TSELFES N.Optimal Multi-step Collocation:Application to the Space-wise Approach for GOCE Data Analysis[J].Journal of Goedesy,2009,83 (1):13-29.

[14] HEISKANEN WA,MORTIZ H.Physical Geodesy[M]. San Francisco:WH Freeman,1967.

[15] COLOMBO O L.Optimal Estimation from Data Regularly Sampled on a Sphere with Applications in Geodesy:Rep 291[R].Columbus:Ohio State University,1979.

[16] COLOMBO O L.Numerical Methods for Harmonic Analysison the Sphere[R]. Columbus:Ohio State University,1981.

[17] WU Xing.Research of Methods of Spherical Harmonic Analysis of the Earth’s Gravity Field[D].Zhengzhou: Information Engineering University,2005.(吳星.地球重力場(chǎng)調(diào)和分析方法研究[D].鄭州:信息工程大學(xué),2005.)

[18] TSCHERNING C C,RAPP R H.Closed Covariance Expressions for Gravity Anomalies,Geoid Undulations, and Deflections of the VerticalImplied by Anomaly Degree-variance Models:Rep 208[R].Columbus:Ohio State University,1974.

[19] XU Zhong,ZHANG Kaiyuan,LU Quan.Fast Arithmetic of Toeplitz Band Matrices[M].Xi’an:Northwestern Polytechnical University Press,2001:31-33.(徐仲,張凱院,陸全.TOEPLITZ矩陣類的快速算法[M].西安:西北工業(yè)大學(xué)出版社,2001:31-33.)

[20] HANSEN P C.The Truncated SVD as a Method for Regularization:NA-M-86-36[R].Stanford:Stanford University,1986,27:534-553.

[21] WANG Zhenjie.Regularization of Ill-posed Problems in Surveying[M].Beijing:Science Press,2006.(王振杰.測(cè)量中不適定問(wèn)題的正則化解法[M].北京:科學(xué)出版社,2006.)

[22] HANSEN P C.Analysis of Discrete Ill-posed Problems by Means of L-curve[J].SIAM Review,1992,34(4): 561-580.

[23] PAVILIS N K,HOLMES S A,KENYON S C,et al.An Earth Gravitational Model to Degree 2160:EGM2008 [R].Vienna:2008 General Assembly of the European Geosciences Union,2008.

[24] ZHANG Chuanyin,GUO Chunxi,CHEN Junyong. EGM2008 and ItsApplication Analysis in Chinese Mainland[J].Acta Geodaetica et Cartographica Sinica, 2009,38(4):283-289.(章傳銀,郭春喜,陳俊勇. EGM2008地球重力場(chǎng)模型在中國(guó)大陸適用性分析[J].測(cè)繪學(xué)報(bào),2009,38(4):283-289.)

[25] WU Xing,ZHANG Chuanding,YE Xiusong,et al. Simulation of Satellite Gravity Gradient Tensor Observations[J].Journal of Geomatics Science and Technology, 2008,25(6):391-395.(吳星,張傳定,葉修松,等.衛(wèi)星重力梯度數(shù)據(jù)的模擬研究[J].測(cè)繪科學(xué)技術(shù)學(xué)報(bào), 2008,25(6):391-395.)

[26] WUXing,ZHANG Chuanding,ZHAO Dongming. Generalized T orus Harmonic Analysis of Satellite Gravity Gradients Component[J].Acta Geodaetica et Cartographica Sinica,2009,38(2):101-107.(吳星,張傳定,趙東明.衛(wèi)星重力梯度分量的廣義輪胎調(diào)和分析方法[J].測(cè)繪學(xué)報(bào),2009,38(2):101-107.)

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

Least-squares Collocation Harmonic Analysis of the Radial Satellite Gravity Gradients

WU Xing1,2,ZHANG Chuanding3,LIU Xiaogang3
1.Center for Engineering Design and Research under the Headquarters of General Equipment,Beijing 100028,China;2.Institute of Hydrographic Surveying and Charting;Tianjin 300061,China;3.Institute of Surveying and Mapping,Information Engineering University,Zhengzhou 450052,China

The least-squares collocation harmonic analysis method,which is used to determine the earth geopotential coefficients from the radial satellite gravity gradient,is deeply studied.The principle of the least-squares collocation is firstly dissertated,and the covariance and self-covariance matrix between the disturbing gravity gradients and spherical harmonics are derived.When the disturbing gravity gradients are in the regularized equi-longitude grid,the covariance matrix between geopotential and disturbing gravity gradients has the configuration of blocked Toeplitz circulation matrix,and its degree can be lowered by using the fast Fourier transform technology effectively.Truncated singular value decomposition(TSVD)which is used to solve the ill-posed problem of covariance matrix is studied.The complete formula of least-squares collocation harmonic analysis of the radial gravity gradient is finally obtained.The simulation experiment results show that the least-squares collocation harmonic analysis based on TSVD can recover the global gravity field in a rather high precision and validity and practicability of algorithms in this paper are also testified.

Earth gravity field model;least-squares collocation;satellite gravity gradient;harmonic analysis; covariance matrix;Toeplitz matrix;TSVD

WU Xing(1979—),male,PhD,engineer, majors in Geodesy and Surveying Engineering.

E-mail:wuxing1979@163.com

1001-1595(2010)05-0471-07

P223

A

國(guó)家自然科學(xué)基金(40774031);全國(guó)優(yōu)秀博士論文專項(xiàng)基金(200344);信息工程大學(xué)博士生創(chuàng)新基金(200707);中科院動(dòng)力大地測(cè)量學(xué)重點(diǎn)實(shí)驗(yàn)室開(kāi)放基金(L06-06)。

2009-07-16

2010-03-02

吳星(1979—),男,博士,工程師,研究方向?yàn)榇蟮販y(cè)量學(xué)與測(cè)量工程。