衛星重力徑向梯度數據的最小二乘配置調和分析

吳 星,張傳定,劉曉剛

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

衛星重力徑向梯度數據的最小二乘配置調和分析

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

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

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

地球重力場模型;最小二乘配置;衛星重力梯度;調和分析;協方差陣;Toeplitz矩陣;TSVD

1 引 言

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

為此,本文主要研究衛星重力徑向梯度規則網格觀測數據的最小二乘配置調和分析的算法,推導引力位系數與格網平均擾動引力梯度間的協方差矩陣,研究協方差矩陣降階求逆方法及其數值穩定性,給出仿真試驗結果。

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

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

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

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

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

隨機模型為

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

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

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

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

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

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

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

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

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

2.3 擾動引力位系數與擾動引力梯度的協方差

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

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

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

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

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

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

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

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

這說明了擾動引力位系數的協方差是一對角陣。

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

格網平均擾動引力梯度徑向分量間的自協方差為

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

2.4 最小二乘配置調和分析解

將規則擾動引力梯度徑向分量觀測數據作以下排列

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

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

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

其中,

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

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

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

2.5 矩陣TSVD正則化求解

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

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

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

則式(32)可表示為

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

其中,整數k稱為截斷參數,也就是正則化參數。

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

3 仿真試算與分析

3.1 協方差函數的計算

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

圖1 擾動位階方差統計圖Fig.1 Chart of potential degree-variance

3.2 擾動引力梯度徑向分量數據的仿真計算

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

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

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

表1 EGM2008模擬的衛星重力徑向梯度數據Tab.1 Radialsatellite gravity gradientssimulated using EGM2008

3.3 調和仿真分析計算

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

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

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

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

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

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

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

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

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

圖3 輪胎調和分析結果Fig.3 Torus harmonic analysis of radial gravity gradients

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

4 結 論

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

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

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

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

4.由于兩極地區子午線收斂問題,導致了全球協方差矩陣病態性問題,采用截斷奇異值分解正則化方法能夠有效地減弱法矩陣的病態性。

[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.(羅志才.利用衛星重力梯度數據確定地球重力場的理論與方法[D].武漢:武漢測繪科技大學,1996.)

[6] ZHANG Chuanding.Satellite-borne Gravity Measurement: Foundation,Modeling Methods,and Data Processing Algorithms[D].Zhengzhou:Information Engineering University,2000.(張傳定.衛星重力測量:基礎、模型化方法與數據處理算法[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.(吳星.地球重力場調和分析方法研究[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].西安:西北工業大學出版社,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.(吳星,張傳定,葉修松,等.衛星重力梯度數據的模擬研究[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.(吳星,張傳定,趙東明.衛星重力梯度分量的廣義輪胎調和分析方法[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);全國優秀博士論文專項基金(200344);信息工程大學博士生創新基金(200707);中科院動力大地測量學重點實驗室開放基金(L06-06)。

2009-07-16

2010-03-02

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