吳 星,張傳定,劉曉剛

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

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

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

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

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

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

1 引 言

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

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

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

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

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

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

式中,l為觀測值向量;t是已測信號向量;s為待估信號向量;I為單位陣;v為隨機的觀測誤差向量。

隨機模型為

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

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

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

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

2.2 二階徑向擾動重力梯度的級數(shù)展開式

地球外部任意一點 P(r,θ,λ)處擾動位 T(r, θ,λ)的球諧函數(shù)級數(shù)表達式為[6]

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

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

并且可以通過遞推公式計算[6,17]。

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

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

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

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

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

其中,α1、α2、S1、S2為實常數(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),其他量是無量綱的數(shù)。將式(12)代入式(11)并代入式(10)后,完成級數(shù)求和,即得全球協(xié)方差的解析函數(shù)。

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

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

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

這說明了擾動引力位系數(shù)的協(xié)方差是一對角陣。

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

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

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

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

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

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

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

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

其中,

如此,對于 N2個子對角陣Dii′,若追蹤 m指標,則共有Nl個N×N維的滿秩矩陣

對這 Nl個Rm求逆,再按照對應項完成矩陣的乘法運算后,得到位系數(shù)的估值為

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

2.5 矩陣TSVD正則化求解

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

設矩陣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稱為截斷參數(shù),也就是正則化參數(shù)。

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

3 仿真試算與分析

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

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

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

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

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

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

為簡單起見,本文假設得到擾動引力梯度張量模擬觀測值是等精度觀測的,而且全球均勻分布,假定得到的40′×40′格網(wǎng)平均值的精度為1.0 mE,則根據(jù)式(37)可知,和50′×50′格網(wǎng)平均值的精度為0.84 mE,故數(shù)據(jù)B是在數(shù)據(jù)A的基礎上分別均加入了標準差為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)和仿真分析計算

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

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

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

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

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

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

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

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

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

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

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

4 結 論

通過本文的理論推導和數(shù)值試驗分析可得如下結論:

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

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

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

4.由于兩極地區(qū)子午線收斂問題,導致了全球協(xié)方差矩陣病態(tài)性問題,采用截斷奇異值分解正則化方法能夠有效地減弱法矩陣的病態(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ù)確定地球重力場的理論與方法[D].武漢:武漢測繪科技大學,1996.)

[6] ZHANG Chuanding.Satellite-borne Gravity Measurement: Foundation,Modeling Methods,and Data Processing Algorithms[D].Zhengzhou:Information Engineering University,2000.(張傳定.衛(wèi)星重力測量:基礎、模型化方法與數(shù)據(jù)處理算法[D].鄭州:信息工程大學,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].大地測量與地球動力學,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.(吳星.地球重力場調(diào)和分析方法研究[D].鄭州:信息工程大學,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è)大學出版社,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.(王振杰.測量中不適定問題的正則化解法[M].北京:科學出版社,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地球重力場模型在中國大陸適用性分析[J].測繪學報,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].測繪科學技術學報, 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].測繪學報,2009,38(2):101-107.)

(責任編輯:叢樹平)

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

國家自然科學基金(40774031);全國優(yōu)秀博士論文專項基金(200344);信息工程大學博士生創(chuàng)新基金(200707);中科院動力大地測量學重點實驗室開放基金(L06-06)。

2009-07-16

2010-03-02

吳星(1979—),男,博士,工程師,研究方向為大地測量學與測量工程。