GRACE-FO加速度計校正方法研究

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

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

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

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

1 加速度計校正基本理論

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

1.1 動力學基本方法

衛星在空間飛行滿足牛頓運動學定律，可寫成如下形式：

(1)

(2)

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

(3)

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

(4)

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

(5)

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

1.2 加速度計偏差參數校正策略

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

ACCnew=B+S·ACC1B

(6)

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

B=C0+C1×t

(7)

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

表1 加速度計偏差參數校正策略

1.3 攝動力模型

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

表2 力模型

2 結果分析

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

2.1 軌道擬合結果

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

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

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

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

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

2.2 時變重力場結果

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

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

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

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

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

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

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

進一步分析線性偏差參數形式(估計C0、C1項)對時變重力場解算結果的影響。表5為采用線性偏差參數形式(估計C0、C1項)對加速度計三軸分別采用不同校正頻率統計的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軸的校正頻率，以便從校正參數個數和不同方向的校正形式進行分析。從表5和圖4可以看出，校正參數個數為20、12、8時的計算結果基本相同，時變重力場解算結果未表現出明顯差異；而當校正參數個數為52時，在20、40和60階累計階方差處case22、case30的誤差較大，分別約為1.45 cm、1.47 cm，1.95 cm、1.94 cm和2.51 cm、2.58 cm，為所有策略中結果較差的2個。其中，case26在20、40、60階的大地水準面累計階方差為1.41 cm、1.84 cm、2.24 cm，為最優策略。不難看出，對于采用線性偏差參數形式進行加速度計校正時，參數個數的影響小于校正策略，如case22、case26和case30的結果表現出明顯差異，case26為所有策略中最優，即在X軸、Z軸上采用1 h校正頻率所計算的階方差會明顯變差，而對Y軸采用1 h校正頻率計算的階方差會明顯提高，這與采用常數偏差參數形式進行加速度計校正所得結論基本一致。

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

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

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

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

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

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

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

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

3 結 語

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