利用GOCE衛星觀測數據反演地球重力場模型

梁建青　沈云中　陳秋杰,2　張興福

1　同濟大學測繪與地理信息學院，上海市四平路1239號，200092 2　香港理工大學建設及環境學院，香港育才道11號，999077 3　廣東工業大學土木與交通工程學院，廣州市大學城外環西路100號，510006

?

利用GOCE衛星觀測數據反演地球重力場模型

梁建青1沈云中1陳秋杰1,2張興福3

1同濟大學測繪與地理信息學院，上海市四平路1239號，200092 2香港理工大學建設及環境學院，香港育才道11號，999077 3廣東工業大學土木與交通工程學院，廣州市大學城外環西路100號，510006

摘要：利用GOCE衛星約6個月的重力梯度數據和約1 a的幾何軌道數據，聯合解算250階次的地球重力場模型TJGOCE01。GOCE重力梯度數據的低頻誤差采用ⅡR數字濾波器處理，粗差采用閥值法和移動窗口閥值法組合探測與剔除。直接在梯度儀坐標系中建立GOCE衛星的重力梯度觀測方程，采用改進的短弧邊值法建立幾何軌道觀測方程。兩類觀測值的權根據其先驗精度確定，采用Kaula規則約束的正則化方法解算法方程。解算的TJGOCE01模型相對于EIGEN6C2模型在250階次的大地水準面誤差和大地水準面累積誤差分別為19.4 mm和177.9 mm。北美地區GPS水準觀測數據的檢驗結果表明，TJGOCE01模型的中誤差為0.544 m，略優于歐空局公布的同階次的第二代時域法和空域法解算的GOCE重力場模型。

關鍵詞：GOCE衛星；重力場模型；直接法；Kaula正則化；ⅡR 濾波器

GOCE計劃的預期目標是以優于100 km的空間分辨率，求得重力異常精度優于1～2 mGal、大地水準面精度優于1 cm的全球靜態重力場模型。然而，由于重力梯度儀本身的缺陷，無法給出全頻段的高精度梯度數據，且只有4個梯度分量在觀測頻段0.005～0.1 Hz范圍內能夠達到所需精度。在觀測頻段外，特別是低頻部分則包含很大的誤差，需要采用合適的濾波方法進行處理[1-2]。GOCE觀測數據存在粗差，需要分析這些粗差的來源和種類，并選擇合適的探測方法進行剔除[3]。GOCE衛星采用傾角約為96.7°的近圓形太陽同步傾斜軌道，因此出現極區數據空白，引起解算模型病態，需要進行正則化處理。此外，因為GOCE衛星的軌道數據對重力場的中長波信號敏感，重力梯度數據對重力場的中短波信號敏感，所以在聯合軌道數據和重力梯度數據解算高精度、高分辨率重力場模型的過程中需要合理地確定兩類觀測數據的權。目前，利用衛星重力梯度數據恢復重力場模型的代表性方法有直接法、時域法和空域法。空域法能夠快速解算出重力場模型，但是精度相對較低；時域法能夠解算得到獨立的GOCE重力場模型，但計算量較大；直接法解算重力場模型的精度最高，但計算量也最大。基于這些方法，ICGEM(International Centre for Global Earth Models)發布了利用GOCE衛星觀測數據解算的一系列地球重力場模型。

本文利用GOCE衛星2011-01-19～2011-06-30的重力梯度數據(約6個月)和2011-01-19～2011-12-31的軌道數據(約12個月)，采用無限脈沖帶通數字濾波器處理GOCE梯度數據的低頻誤差，采用閥值法和移動窗口閥值法相結合探測與剔除GOCE觀測數據的粗差，采用直接法解算250階次的重力場模型TJGOCE01，最后用EIGEN6C2模型和北美地區的GPS水準觀測數據對該重力場模型進行檢驗。

1 觀測模型和數據處理方法

1.1GOCE軌道數據恢復重力場的觀測方程

利用高低衛星跟蹤衛星數據恢復重力場模型的方法主要有Kaula線性攝動法、動力學積分法、短弧邊值法、加速度法和能量守恒法等。傳統的短弧邊值法的觀測方程為：

(1)

(2)

GOCE衛星在運行過程中受到地球中心引力、非球形引力、日月引力、固體潮和海潮等保守力與大氣阻力、太陽光壓等非保守力攝動。保守力可以根據高精度的模型計算得到，采用的保守力模型及其說明見表1[4]。GOCE衛星非保守力采用無阻力推進系統進行補償，但是由于衛星姿態調控等多種因素的影響，不能完全補償掉非保守力，因此通過估計3個方向的偏差參數和1周期的振幅參數來吸收殘余非保守力加速度[5-6]：

(3)

其中，(fx,fy,fz)為GOCE衛星在3個方向受到的非保守力，(ax,ay,az)和(bx,by,bz,cx,cy,cz)分別是非保守力加速度的偏差和振幅參數，v=2πt/T為真近點角，T=5 400 s為衛星的運轉周期。需要注意的是，式(3)只在局部軌道直角坐標系中成立。

根據改進短弧邊值法[4]，式(1)中所有衛星的位置向量，包括力模型中的位置向量均引入改正數。離散化式(1)，得：

(4)

式中，u0、p0分別是u、p的先驗值，δu、δp是要估計的未知參數，αk是積分系數，N是弧段的歷元總數。

消去殘余非保守力加速度的補償參數，并將所有弧段的法方程合并，得到用軌道數據求解位系數的總法方程：

(5)

式中，N1為法方程系數矩陣，w1為其常數項。

1.2GOCE梯度數據恢復重力場的觀測方程

GOCE衛星的重力梯度分量采用3對加速度計觀測得到，其9個重力梯度分量為：

(6)

由于Γ是對稱矩陣，且其3個主對角量滿足拉普拉斯方程，因此只有5個量是獨立的。由于重力梯度儀自身的缺陷，只有Γxx、 Γyy、Γzz、Γxz分量滿足測量精度要求。為避免精度較差的重力梯度分量對解算重力場模型的影響，只采用4個精度滿足要求的重力梯度分量，并在梯度儀坐標系(GRF)下直接建立觀測方程。

在地固坐標系(EFRF)下的引力梯度張量可以表示為：

(7)

將該引力梯度張量轉換至梯度儀坐標系，轉換矩陣為：

(8)

(9)

式中，E為觀測誤差組成的矩陣。因此，引入重力位系數的近似值u0代入式(9)右邊，將式(9)中滿足精度要求的4個觀測量表示成向量，并將各歷元的誤差方程合并在一起，表示為：

v=Aδu-d

(10)

式中，δu為位系數改正數向量，A為其設計矩陣，d為梯度觀測值與近似值u0求得的參考引力梯度值之差，v為殘差向量。

由于GOCE衛星的梯度測量頻帶為0.005～0.1 Hz,低頻部分包含較大的誤差，若采用有限脈沖帶通數字濾波器(FIR)處理低頻誤差，則計算量較大[7]。因此，本文采用無限脈沖帶通數字濾波器(ⅡR)處理梯度數據的低頻誤差。對誤差方程的系數矩陣和梯度觀測值進行同步濾波，從而避免濾波產生的相位漂移問題帶來的影響，具體濾波過程為：

(11)

式中，F{}為ⅡR濾波算子。經過ⅡR濾波之后的梯度儀坐標系下直接法求解重力場模型的法方程為：

N2δu=w2

(12)

式中，N2=(F{AGRF})TP(F{AGRF}),w2=(F{A})TPFg0gggggg。

1.3聯合GOCE軌道和梯度數據的解算模型

GOCE衛星的軌道和梯度數據分別反映了不同頻段的重力場信號，其中軌道數據對重力場的中長波信號比較敏感，梯度數據對重力場的中短波信號比較敏感。聯合這兩種類型的觀測數據解算高精度的重力場模型必須合理確定其權。本文根據軌道和梯度數據的先驗精度進行定權，其中GOCE軌道的精度約為0.02 m，梯度的測量精度約為2 mE。若計算式(5)和式(12)時，已經根據先驗精度進行了定權，則聯合軌道數和梯度數據解算重力場模型的法方程為：

(13)

由于利用GOCE觀測數據解算重力場模型是病態的，因此采用Kaula規則約束的正則化方法解算病態方程。GOCE病態方程的正則化解為：

(14)

式中，α為正則化參數，正則化矩陣K為對角陣，其對角線元素為Kaula階方差的倒數。根據Kaula 規則，n階系數對應的對角線元素為：

kn=2.0×1010n4

(15)

最優正則化參數α根據L曲線法確定[8-9]。

2 數據處理與重力場模型解算

2.1數據預處理

閥值法能夠準確識別時間較長區域的粗差，移動窗口的閥值法能夠較好地識別離散粗差或短歷元區域粗差，在探測粗差時兩者可以互補，所以本文采用閥值法和移動窗口的閥值法組合探測與剔除粗差[10]。由于重力梯度在觀測頻段外，特別是低頻部分包含大量的噪聲，因此需要對梯度數據進行濾波處理。本文采用9階巴特沃斯(Butter)的無限脈沖帶通數字濾波器(ⅡR)對誤差方程的系數矩陣和梯度擾動值進行同步濾波。

為驗證ⅡR數字濾波器處理梯度數據低頻誤差的效果，取2011-03-01～2011-03-10梯度數據的各個梯度分量進行濾波處理，濾波器的帶通頻率為0.005～0.1 Hz，根據式(11)得到濾波前后的誤差功率譜(圖1)。

由圖1可知，該濾波方法基本上能夠達到濾除低頻誤差的目的，同時也能夠有效保留觀測頻帶內的信號。

2.2重力場模型解算

采用GOCE衛星2011-01-19～2011-06-30的梯度觀測數據和2011-01-19～2011-12-31的幾何軌道觀測數據解算重力場模型，其中軌道數據只用于解算120階次的模型系數，梯度數據解算至250階次的模型系數。各梯度數據Γxx、Γyy、Γzz、Γxz的權相同，軌道與梯度數據的權根據其先驗精度確定，若軌道的權為1，則梯度的權為1.0×1020。參考引力梯度由參考模型EIGEN5C計算得到。

在重力場模型解算的過程中，按60 min的弧長單獨形成各個弧段的法方程，并按弧段消去其他參數，最后將所有弧段關于位系數的法方程疊加求解位系數。由于該解算模型是病態的，采用Kaula規則對二階以上的重力位系數約束后進行正則化解算，其最優正則化參數根據L曲線法確定(圖2)，為3.16×107，由式(14)計算重力場模型TJGOCE01。

3 重力場模型分析與檢驗

歐空局發布的第二代重力場模型TIM-R2、SPW-R2和DIR-R2只采用了GOCE衛星的觀測數據，其數據量與TJGOCE01模型接近。本文比較TJGOCE01、TIM-R2、SPW-R2和DIR-R2模型相對于EIGEN6C2模型的大地水準面誤差和大地水準面累積誤差，結果見圖3、圖4。

由圖3可知，TJGOCE01模型前120階次位系數精度均優于TIM-R2與SPW-R2模型。TJGOCE01模型前80階次位系數精度明顯低于EIGEN5C模型，是因為EIGEN5C模型采用了GRACE星間距離變率數據；但其在80～180階次部分精度優于EIGEN5C模型，說明梯度數據能有效改善中短波長重力位精度。TJGOCE01模型精度在160階次之前低于DIR-R2模型，因為其采用的梯度數據量相對較少，并且只采用了Kaula規則約束的正則化方法解算病態法方程；而DIR-R2模型采用ITG-GRACE2010S作為先驗信息的SCRA正則化方法和Kaula正則化方法共同解算法方程，因此其低階次的精度更接近于ITG-GRACE2010S模型的精度。TJGOCE01模型在180階次之后和TIM-R2、SPW-R2、DIR-R2模型的精度相當。

由圖4可知，TJGOCE01模型相對于EIGEN5C模型，在同階次的大地水準面累積誤差均小于TIM-R2和SPW-R2模型，主要原因是TJGOCE01模型采用了1 a的軌道數據。

為進一步驗證TJGOCE01模型的精度，選取2012年北美地區GPS水準網中的21 289個點對TJGOCE01模型進行檢核。由于各個模型的最大階次不一樣，將TJGOCE01、TIM-R2、SPW-R2和DIR-R2模型分別截斷至相同的階次并計算其大地水準面差距，最后將計算的大地水準面差距和GPS水準數據比較，以驗證TJGOCE01模型的外部符合精度(表2)。

由表2可知，各個模型計算的大地水準面差距與GPS水準數據間的標準差都超過0.5 m，主要原因是各模型最大階次僅為250，模型的截斷誤差較大。TJGOCE01模型的標準差與其他模型相當，驗證了TJGOCE01模型的可靠性。

4 結語

本文利用GOCE軌道數據和梯度數據，采用直接法解算了250階次的重力場模型TJGOCE01。采用無限脈沖帶通數字濾波器處理梯度數據的低頻誤差，并采用閥值法和移動窗口的閥值法組合探測與剔除粗差。在梯度儀坐標系中直接建立梯度觀測方程，采用改進短弧邊值法處理軌道數據，同時根據觀測數據的先驗精度確定軌道數據和梯度數據的最優權，并采用Kaula規則約束的正則化方法解算法方程，采用L曲線法求解最優正則化參數。

以EIGEN6C2模型為參考，TJGOCE01模型精度在120階次之前高于TIM-R2和SPW-R2模型，低于DIR-R2模型；在80～180階次高于EIGEN5C模型。相對于EIGEN6C2模型，TJGOCE01模型在250階次的大地水準面誤差和大地水準面累積誤差為19.4 mm和177.9 mm，略優于歐空局公布的第二代時域法和空域法解算的同階次的GOCE重力場模型。北美GPS水準觀測數據的檢驗結果表明，TJGOCE01模型的標準差與TIM-R2、SPW-R2和DIR-R2重力場模型的差異比較小，驗證了本文反演方法的正確性。

參考文獻

[1]Schuch W D. The Processing of Band-Limited Measurements:Filtering Techniques in the Least Squares Context and in the Presence of Data Gaps[J].Space Science Review，2003,108(1-2):67-68

[2]Reguzzoni M, Tselfes N. Optimal Multi-Step Collocation: Application to the Space-Wise Approach for GOCE Data Analysis[J].Journal of Geodesy, 2009, 83(1):13-29

[3]吳云龍，羅志才，李輝，等. 基于小波收縮閾值降噪的衛星重力梯度數據粗差探測方法[J].大地測量與地球動力學，2010，30(4):55-58(Wu Yunlong, Luo Zhicai, Li Hui, et al.Outlier Detection Algorithm for Satellite Gravity Gradient Data by Using Wavelet Shrinkage Denoising[J].Journal of Geodesy and Geodynamics, 2010,30(4):55-58)

[4]Chen Q J, Shen Y Z, Zhang X F, et al. Monthly Gravity Field Models Derived from GRACE Level 1B Data Using a Modified Short-Arc Approach[J].Journal of Geophysical Research: Solid Earth, 2015,120(3):1 804-1 819

[5]蘇勇, 范東明, 游為. 利用GOCE衛星數據確定全球重力場模型[J]. 物理學報, 2014,63(9):1-9(Su Yong, Fan Dongming, You Wei. Gravity Field Model Calculated by Using the GOCE Data[J]. Acta Phys Sin, 2014,63(9):1-9)[6]游為. 應用低軌衛星數據反演地球重力場模型的理論和方法[D].成都：西南交通大學，2008(You Wei. Theory and Methodology of Earth’s Gravitational Field Model Recovery by LEO Data[D]. Chengdu：Southwest Jiaotong University, 2008)

[7]Yi W Y. An Alternative Computation of a Gravity Field Model from GOCE[J]. Advances in Space Research, 2012, 50(3):371-384

[8]沈云中. 應用CHAMP衛星星歷精化地球重力場模型的研究[D]. 武漢：中國科學院測量與地球物理研究所,2000(Shen Yunzhong. Study of Recovering Gravitational Potential Model from the Ephemerides of CHAMP[D]. Wuhan：Institute of Geodesy and Geophysics, CAS, 2000)

[9]徐新禹，李建成, 王正濤, 等. Tikhonov正則化方法在GOCE重力場求解中的模擬研究[J].測繪學報,2010,39(5):466-470(Xu Xinyu, Li Jiancheng, Wang Zhengtao, et al. The Simulation Research on the Tikhonov Regularization Applied in Gravity Field Determination of GOCE Satellite Mission[J]. Acta Geodaetica et Cartographica Sinica,2010,39(5):466-470)

[10]周睿. GOCE衛星SGG數據處理及空域法恢復地球重力場模型研究[D].鄭州：信息工程大學，2013(Zhou Rui. Research on GOCE SGG Data Processing and the Space-Wise Approach to Recover Gravity Field Model[D]. Zhengzhou: Information Engineering University, 2013)

Foundation support:National Natural Science Foundation of China, No.41474017, 41274035.

About the first author:LIANG Jianqing, postgraduate, majors in satellite gravimetry, E-mail:liangjianqing9999@163.com.

Gravity Field Model Recovery Using the GOCE Data

LIANGJianqing1SHENYunzhong1CHENQiujie1,2ZHANGXingfu3

1College of Surveying and Geo-Informatics, Tongji University, 1239 Siping Road, Shanghai 200092, China 2Faculty of Construction and Environment, Hong Kong Polytechnic University,11 Yucai Street, Hong Kong 999077, China 3School of Civil and Transportation Engineering, Guangdong University of Technology, 100 West-Waihuan Road, Guangzhou 510006, China

Abstract:In this paper, a global gravity field model entitled TJGOCE01 up to the degree and order 250 is recovered from about 6 months of GOCE gravity gradient data and 12 months of GOCE orbit data,

using a direct approach. The gravity gradient data are filtered by an infinite impulse bandpass digital filter, and gross errors are detected and removed using the threshold and a moving window threshold methods. The gravity gradient observational equations are established in the gradiometer coordinate system directly and the kinematic orbit observation equations are established by using a modified short-arc approach. The weights of the two kinds of observations are determined based on their prior accuracies, and the regularization approach based on Kaula’s rule is applied to solve the normal equation. The geoid error and cumulative geoid error up to degree and order 250 of TJGOCE01 model with respect to EIGEN6C2 model are 19.4 mm and 177.9 mm, respectively. The test results with the GPS leveling data in North America show that the root mean square error of TJGOCE01 model is 0.544 m, which is superior to the second generation models developed by the European Space Agency using time-wise and space-wise approaches.

Key words:GOCE satellite; gravity model; direct approach; Kaula regularization; ⅡR filter

收稿日期：2015-06-19

第一作者簡介：梁建青，碩士生，主要研究方向為衛星重力測量，E-mail:liangjianqing9999@163.com。

DOI：10.14075/j.jgg.2016.06.001

文章編號：1671-5942(2016)06-0471-05

中圖分類號：P223

文獻標識碼：A

項目來源：國家自然科學基金(41474017, 41274035)。