牛晗晗 王長青 鐘 敏 閆易浩,3 梁 磊

1 中國科學(xué)院精密測量科學(xué)與技術(shù)創(chuàng)新研究院大地測量與地球動力學(xué)國家重點實驗室，武漢市徐東大街340號，430077 2 中國科學(xué)院大學(xué)，北京市玉泉路19號甲，100049 3 華中科技大學(xué)物理學(xué)院，武漢市珞喻路1037號，430074

由于衛(wèi)星在空間受多種因素的影響，在重力場反演過程中需對加速度計觀測的非保守力進行校正[1-4]。不同機構(gòu)針對加速度計校正采用的方法不同。加速度計的校正主要是通過建立加速度計校正參數(shù)的設(shè)計矩陣與軌道或星間距變率的觀測方程來實現(xiàn)的，影響因素較多，是一個十分復(fù)雜的過程[5-10]，而時變重力場的解算精度和加速度校正策略具有一定的關(guān)聯(lián)。

相較于GRACE衛(wèi)星，GRACE-FO衛(wèi)星搭載了激光測距設(shè)備，對衛(wèi)星的姿態(tài)有較高的要求，會通過多次點火來調(diào)整衛(wèi)星姿態(tài)以維持激光測距，而加速度計能夠敏感捕捉到衛(wèi)星平臺頻繁機動引起的信號，并體現(xiàn)在加速度計觀測數(shù)據(jù)中。因此，在GRACE-FO時變重力場反演中，需要探索加速度計校正策略對時變重力場解算的影響，并綜合分析滿足重力場反演需要的加速度計校正方式。本文以GRACE-FO衛(wèi)星2019-01數(shù)據(jù)為例，探索不同加速度計偏差參數(shù)校正策略對時變重力場反演結(jié)果的影響。

1 加速度計校正基本理論

加速度計校正是通過建立加速度計參數(shù)設(shè)計矩陣與軌道觀測方程，并對觀測方程進行最小二乘求解，以實現(xiàn)對加速度計參數(shù)的校正。本文通過采用動力學(xué)方法建立觀測方程進行加速度計校正。

1.1 動力學(xué)基本方法

衛(wèi)星在空間飛行滿足牛頓運動學(xué)定律，可寫成如下形式：

(1)

(2)

式中，Φr、Sr為狀態(tài)轉(zhuǎn)移矩陣和參數(shù)敏感矩陣，是通過對r(t)在r0(t)處進行泰勒展開，取線性形式得到的，即

(3)

同理可得，星間距變率的觀測方程為：

(4)

此外，在重力場反演過程中需要聯(lián)合軌道和星間距變率觀測方程進行求解，其法方程可寫成如下形式：

(5)

式中，Nleo、Nkbrr為軌道和星間距變率觀測的法方程，bleo、bkbrr為觀測值，y為待求量。

1.2 加速度計偏差參數(shù)校正策略

加速度計的校正是在衛(wèi)星科學(xué)坐標系下進行的，其形式可寫成[11]：

ACCnew=B+S·ACC1B

(6)

式中，ACCnew為校正后的加速度計數(shù)據(jù)，B為加速度計偏差因子，S為尺度因子，ACC1B為L1B加速度計數(shù)據(jù)。由于本文只考慮加速度計的偏差校正，因此將加速度計的尺度參數(shù)設(shè)置為單位陣。對于加速度計的偏差因子B，設(shè)置的校正形式為：

B=C0+C1×t

(7)

式中，C0、C1為待求常數(shù)；t為時間，從每個弧段起始算起。設(shè)置加速度計校正的弧長為24 h，表1為偏差參數(shù)校正策略，包括常數(shù)偏差參數(shù)形式(估計C0項)和線性偏差參數(shù)形式(估計C0、C1項)，X、Y、Z表示的是SFR坐標系。

表1 加速度計偏差參數(shù)校正策略

1.3 攝動力模型

本文采用JPL發(fā)布的GNV1B軌道數(shù)據(jù)對GRACE-FO兩顆衛(wèi)星的加速度計偏差參數(shù)進行校正，在校正過程中需同時估計衛(wèi)星的初始狀態(tài)，所用數(shù)據(jù)主要為ACT1B、SCA1B、KBR1B、GNV1B，在構(gòu)建參考軌道時采用的力模型見表2。

表2 力模型

2 結(jié)果分析

本文通過分析表1中的33種偏差參數(shù)校正策略，對比軌道擬合結(jié)果及60階無約束重力場結(jié)果來探討所選取的加速度計校正策略對時變重力場解算的影響。

2.1 軌道擬合結(jié)果

表3為計算的每日軌道三軸殘差(RMS)結(jié)果和軌道3D殘差結(jié)果。結(jié)合表1和3，由case1～5(估計C0項)和case18～21(估計C0、C1項)可以看出，當只估計偏差參數(shù)常數(shù)項C0時，X、Y、Z方向校正的加速度計參數(shù)由72個減小到3個，計算的軌道三軸殘差逐漸增大，而軌道3D殘差也從1.51 cm逐漸增大到10 cm左右；同時估計加速度計常數(shù)項C0和線性項C1時，X、Y、Z方向校正的加速度計參數(shù)由48個減少到6個，軌道3D殘差由5 cm左右逐漸增大到10.97 cm。圖1(a)為計算的軌道3D殘差，圖中實線為采用常數(shù)偏差參數(shù)形式(估計C0項)，虛線為采用線性偏差參數(shù)形式(估計C0、C1項)，可以看出，計算的軌道3D殘差隨校正頻率的減小而逐漸增大；對比2種形式的計算結(jié)果發(fā)現(xiàn)，在校正頻率相同的情況下，加速度計校正參數(shù)越多軌道殘差越小，擬合的軌道越優(yōu)，24 h弧段每24 h校正1次加速度計參數(shù)時2種估計形式獲得的殘差沒有表現(xiàn)出明顯差異。

表3 采用不同策略擬合的軌道殘差

加速度計Y軸為非敏感軸，其測量精度相比于X軸和Z軸約低1個量級。進一步考察對加速度計三軸分別采用不同校正頻率進行計算獲取的軌道擬合結(jié)果，分析X、Y和Z軸采用每1 h、3 h、6 h、12 h、24 h校正1次計算所擬合的軌道殘差。結(jié)合表1和3中case5～17估計C0項和case21～33估計C0、C1項的結(jié)果可以看出：1)僅針對X軸，只估計C0項(case5～9)時，所校正的參數(shù)由1 h的26個減少到24 h的3個，計算的軌道3D殘差由6.75 cm增大到10.60 cm；當同時估計C0和C1項(case21～25)時，校正參數(shù)由1 h的52個減少到24 h的6個，計算的軌道3D殘差由5.27 cm增大到10.79 cm；2)僅針對Y軸，只估計C0項(case10～13)時，所校正的參數(shù)由1 h的26個減少到24 h的3個，計算的軌道3D殘差由8.71 cm增大到10.60 cm；當同時估C0和C1項(case26～29)時，校正參數(shù)由1 h的52個減少到24 h的6個，軌道3D殘差由8.70 cm增大到10.79 cm，并且對比發(fā)現(xiàn)，僅對Y軸進行不同校正頻率分析計算的殘差整體偏大；3)僅針對Z軸，只估計C0項(case14～17)時，所校正的參數(shù)由1 h的26個減少到24 h的3個，計算的軌道3D殘差由6.35 cm增大到10.60 cm；當同時估計C0和C1項(case30～33)時，校正參數(shù)由1 h的52個減少到24 h的6個，計算的軌道3D殘差由5.61 cm增大到10.79 cm。

另外，當校正的參數(shù)個數(shù)相同時(如case22、case26、case30)，所計算的軌道3D殘差在X軸和Z軸上基本相同(5.27 cm、5.61 cm)，但在Y軸上差異較大(8.70 cm)。將case21～33所計算的軌道3D殘差結(jié)果反映在圖1(b)中，由圖中case22～25可知，在對X、Y、Z軸采用不同的校正頻率時，加速度計的校正參數(shù)越多軌道殘差越小，軌道擬合結(jié)果越優(yōu)；校正參數(shù)個數(shù)相同時，擬合的軌道殘差在X軸和Z軸上基本接近，約為5.5 cm，在Y軸上結(jié)果較差，約為8.7 cm。

圖1 不同策略擬合的軌道3D位置殘差

2.2 時變重力場結(jié)果

2.2.1 三軸采用同樣的校正形式

圖2為對加速度計三軸采用相同的校正頻率所得的重力場結(jié)果，其中實線為采用常數(shù)偏差參數(shù)形式(估計C0項)計算的以大地水準面高表示的重力場階方差，虛線為采用線性偏差參數(shù)形式(估計C0、C1項)的計算結(jié)果，通常認為40階之前由時變重力場信號主導(dǎo)，大于40階誤差逐漸占主導(dǎo)。由圖可知，2～40階采用的幾種策略所計算的階方差基本相同；40階之后除了case1外，其他幾種策略計算的階方差結(jié)果差別不大；60階幾種策略的累計大地水準面階方差均為2.34 cm左右。對加速度計三軸采用相同校正形式不同策略計算出的重力場結(jié)果基本相同，需要指出的是，case1計算出的60階累計大地水準面階方差在2.49 cm左右，其原因可能是case1采用的是1 h校正頻率的偏差策略，所校正的加速度計參數(shù)過多或參數(shù)形式并不合適，出現(xiàn)過擬合現(xiàn)象，從而吸收了部分重力場信號，導(dǎo)致計算結(jié)果產(chǎn)生偏差。

圖2 不同策略反演的重力場

2.2.2X、Y、Z軸分別采用不同的校正參數(shù)個數(shù)

采用常數(shù)偏差參數(shù)形式(估計C0項)對加速度計X、Y、Z三軸分別采用不同校正參數(shù)個數(shù)計算的重力場大地水準面累計階方差結(jié)果見表4。由表可見，20階和40階的累計大地水準面階方差基本相同，約為1.42 cm和1.86 cm，反映在圖3則是紅線及藍線比較平緩。表4統(tǒng)計結(jié)果顯示，60階累計大地水準面階方差并未隨加速度計校正參數(shù)個數(shù)的減少變差，如case7～9、case11～13、case15～17等，加速度計校正參數(shù)由每天10個減少到4個，計算結(jié)果基本保持一致，且當三軸校正參數(shù)個數(shù)相同(參數(shù)個數(shù)為10)時，case7、case11、case15計算的 60階累計大地水準面階方差基本相同，約為2.34 cm(參數(shù)個數(shù)為6和4時，其結(jié)果與參數(shù)個數(shù)為10的結(jié)果基本一致)。而僅當參數(shù)個數(shù)為26時，分別對應(yīng)case6(X軸校正頻率為1 h，Y軸、Z軸為24 h)、case10(Y軸校正頻率為1 h，X軸、Z軸為24 h)和case14(Z軸校正頻率為1 h，X軸、Y軸為24 h)，計算結(jié)果差異明顯，約為2.38 cm、2.29 cm和2.40 cm，其中case10結(jié)果明顯優(yōu)于case6和case14。上述分析表明，采用常數(shù)偏差參數(shù)形式進行加速度校正時，參數(shù)個數(shù)并不能直接影響時變重力場解算的結(jié)果，但當不同方向的校正頻率不同時，其中Y軸1 h 校正1次、X軸和Z軸24 h校正1次，解算結(jié)果明顯優(yōu)于表4中其他策略。

圖3 20、40、60階重力場累計階方差(常數(shù)偏差形式)

表4 20階、40階、60階重力場累計階方差(常數(shù)偏差參數(shù)形式)

進一步分析線性偏差參數(shù)形式(估計C0、C1項)對時變重力場解算結(jié)果的影響。表5為采用線性偏差參數(shù)形式(估計C0、C1項)對加速度計三軸分別采用不同校正頻率統(tǒng)計的20階、40階和60階累計大地水準面階方差，其中case22～25為X軸分別采用1 h、3 h、6 h、12 h校正1次的加速度計線性偏差算例，Y軸、Z軸校正頻率固定為24 h；case26～29和case30～33采用同樣的處理方式，只分別改變Y軸和Z軸的校正頻率，以便從校正參數(shù)個數(shù)和不同方向的校正形式進行分析。從表5和圖4可以看出，校正參數(shù)個數(shù)為20、12、8時的計算結(jié)果基本相同，時變重力場解算結(jié)果未表現(xiàn)出明顯差異；而當校正參數(shù)個數(shù)為52時，在20、40和60階累計階方差處case22、case30的誤差較大，分別約為1.45 cm、1.47 cm，1.95 cm、1.94 cm和2.51 cm、2.58 cm，為所有策略中結(jié)果較差的2個。其中，case26在20、40、60階的大地水準面累計階方差為1.41 cm、1.84 cm、2.24 cm，為最優(yōu)策略。不難看出，對于采用線性偏差參數(shù)形式進行加速度計校正時，參數(shù)個數(shù)的影響小于校正策略，如case22、case26和case30的結(jié)果表現(xiàn)出明顯差異，case26為所有策略中最優(yōu)，即在X軸、Z軸上采用1 h校正頻率所計算的階方差會明顯變差，而對Y軸采用1 h校正頻率計算的階方差會明顯提高，這與采用常數(shù)偏差參數(shù)形式進行加速度計校正所得結(jié)論基本一致。

表5 20階、40階、60階累計階方差(線性偏差參數(shù)形式)

圖4 20、40、60階重力場累計階方差(線性偏差參數(shù)形式)

采用不同的加速度計校正策略，從計算的重力場階方差結(jié)果可以看出：1)對加速度計X、Y、Z三軸選用相同的校正頻率，采用常數(shù)偏差參數(shù)形式(估計C0項)所得結(jié)果與線性偏差參數(shù)形式(估計C0、C1項)基本相同，60階累計誤差均為2.34 cm左右；2)對加速度計X、Y、Z三軸分別采用不同的校正頻率，X軸、Z軸采用1 h校正頻率的結(jié)果最差，采用常數(shù)偏差參數(shù)形式和線性偏差參數(shù)形式計算的60階累計階方差分別約為2.38 cm、2.51 cm，2.40 cm、2.58 cm，對于Y軸采用1 h 校正頻率的結(jié)果最優(yōu)，采用常數(shù)偏差參數(shù)形式和線性偏差參數(shù)形式計算的60階累計誤差約為2.29 cm、2.24 cm，其余結(jié)果均為2.34 cm左右。綜合對比分析采用的不同校正策略發(fā)現(xiàn)，case26采用線性偏差參數(shù)形式(估計C0、C1項)，X軸、Z軸24 h校正1次、Y軸1 h校正1次時結(jié)果最優(yōu)。

圖5為計算的階方差結(jié)果和各模型的空間分布，由圖5(a)可見，重力場階方差與GFZ、JPL模型的結(jié)果較為接近，且從模型的空間分布可以看出，計算的時變重力場模型在信號強度、信號明顯區(qū)域(如亞馬遜流域、格陵蘭島、華北平原、中東地區(qū)及南極大陸等)都與官方結(jié)果一致。

圖5 階方差與模型空間分布

2.3 時變重力場關(guān)于加速度計校正的進一步分析與討論

采用常數(shù)和線性偏差參數(shù)形式，利用不同策略進行加速度計校正后發(fā)現(xiàn)，加速度計校正參數(shù)個數(shù)的多少對重力場解算結(jié)果的影響并不明顯，而當采用Y軸1 h校正1次、其他軸24 h校正1次的策略時，無論采用常數(shù)形式還是線性形式，結(jié)果均為最優(yōu)。由于Y軸為不敏感軸，且GRACE-FO的數(shù)據(jù)采集對姿態(tài)要求較高，而衛(wèi)星平臺頻繁機動，根據(jù)GRACE-FO推進器工作原理，其機動方向主要體現(xiàn)在Y軸和Z軸上，因此可通過處理加速度數(shù)據(jù)中推進器信號來解決此問題[12]。首先考慮在L1A-L1B數(shù)據(jù)中采用推進器模型(thruster)剔除影響加速度計數(shù)據(jù)的信號，并在L1B數(shù)據(jù)中采用剔除的推進器信號，然后對剔除部分進行插值來分析推進器因素對重力場反演結(jié)果的影響。圖6中紅線、藍線、黑線、綠線依次為直接采用原始數(shù)據(jù)、采用剔除插值方法(eliminate)、在L1A-L1B數(shù)據(jù)處理中未采用推進器模型(noTH)及加入了推進器模型(TH)反演的重力場結(jié)果，可以看出，對于在L1B數(shù)據(jù)中直接進行處理及在L1A-L1B數(shù)據(jù)處理中不加入推進器模型，20階之后所得的重力場結(jié)果較差；并且，在L1A-L1B數(shù)據(jù)處理中加入推進器模型后結(jié)果得到改善，但仍與原始結(jié)果有差距，說明對于加速度計數(shù)據(jù)，推進器信號的處理方法會直接影響時變重力場的解算精度。與該方法原理相同的case26為最優(yōu)策略，即Y軸的頻繁校正可能是推進器頻繁機動造成的。

圖6 去除推進器信號后的大地水準面階方差

3 結(jié) 語

本文通過動力學(xué)方法，利用2019-01 GRACE-FO衛(wèi)星實測數(shù)據(jù)，計算不同加速度計偏差參數(shù)校正策略，共計33個算例，并分析對比軌道擬合結(jié)果和相應(yīng)的時變重力場反演結(jié)果。軌道擬合作為時變重力場反演的中間步驟，是時變重力場解算的第1步，軌道擬合結(jié)果表明，計算的軌道3D殘差隨校正參數(shù)個數(shù)的減少而逐漸增大，并在加速度計校正頻率相同的情況下，校正參數(shù)個數(shù)越多軌道殘差越小。然而在對時變重力場反演策略進行分析時發(fā)現(xiàn)，反演的時變重力場精度與軌道擬合結(jié)果沒有直接關(guān)系，時變重力場解算結(jié)果的優(yōu)劣與加速度計校正策略有關(guān)，與加速度計校正參數(shù)個數(shù)無關(guān)。綜合分析不同的加速度計校正策略發(fā)現(xiàn)，對于加速度計三軸分別采用不同的校正策略會得到不同的結(jié)果，其差異主要體現(xiàn)在大于40階部分，其中Y軸校正參數(shù)個數(shù)較多時計算結(jié)果最優(yōu)，其原因可能是Y軸為相對不敏感軸，并且由于GRACE-FO衛(wèi)星裝載有激光測距儀，對姿態(tài)要求較高，在衛(wèi)星觀測中需頻繁對衛(wèi)星平臺進行姿態(tài)調(diào)整以保證測量精度和星間激光測距精度，而推進器信號耦合至加速度測量時主要作用于YZ平面，因此Y軸需要較多的校正參數(shù)。本文研究忽略了加速度計尺度因子，主要因為尺度因子和加速度計參數(shù)存在一定的相關(guān)性，在GRACE-FO的時變重力場反演中尺度因子通常為每月校正1次[13]，而對于尺度因子的選取需要在后續(xù)工作中進行進一步研究。