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

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

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

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

1 加速度計(jì)校正基本理論

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

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

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

(1)

(2)

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

(3)

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

(4)

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

(5)

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

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

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

ACCnew=B+S·ACC1B

(6)

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

B=C0+C1×t

(7)

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

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

1.3 攝動力模型

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

表2 力模型

2 結(jié)果分析

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

2.1 軌道擬合結(jié)果

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖6 去除推進(jìn)器信號后的大地水準(zhǔn)面階方差

3 結(jié) 語

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