基于GRACE-FO的微波與激光星間測距數據的重力場模型對比分析

朱紫彤 閆易浩 梁 磊 穆慶祿 王長青 馮 偉,4 鐘 敏,4

1 中國科學院精密測量科學與技術創新研究院大地測量與地球動力學國家重點實驗室，武漢市徐東大街340號，430077 2 中國科學院大學地球與行星科學學院，北京市玉泉路19號甲，100049 3 中山大學物理與天文學院，廣東省珠海市大學路2號，519082 4 中山大學測繪科學與技術學院，廣東省珠海市大學路2號，519082

隨著衛星重力測量技術的發展，特別是GRACE與GRACE Follow-On(GRACE-FO)衛星計劃的成功實施，為監測地球表面質量遷移提供了全新的視角，其產品廣泛應用于水文學、地震學、海洋學和地球物理學等相關學科研究[1]。GRACE-FO是GRACE的后續衛星，兩者原理基本一致，通過精確測量兩顆衛星之間的距離變化率來反演地球時變重力場模型。相較于GRACE，GRACE-FO除搭載K波段微波測距系統(KBR, K-band ranging)外，還搭載激光測距系統(LRI, laser ranging interferometer)，可實現nm級星間測距精度，相較于μm級K波段微波測距系統，其精度至少可提高2個量級[2]。

1 動力學法基本理論

1.1 反演流程

本文通過動力學法建立軌道及星間測距觀測方程，采用最小二乘方法解算重力場系數。數據處理流程主要分為以下步驟[3]：1)以GRACE-FO Level-1B 數據及背景場模型作為輸入數據，利用12階Gauss-Jackson積分器計算積分軌道，與GPS實測軌道作差得到軌道殘差；2)根據積分軌道計算星間距離變率參考值，與KBR或LRI觀測數據作差得到星間距離變率殘差；3)根據兩類觀測數據殘差的中誤差之比定權，再通過最小二乘方法求出時變重力場球諧系數。該過程使用的數據及模型見表1。

表1 重力場反演過程中使用的背景場模型

1.2 參數估計策略

與國內外多數重力衛星處理機構一樣，本文采取24 h弧長，利用動力學方法求解90階無約束地球時變重力場球諧系數。需要注意的是，GPS原始軌道數據為1 s采樣，為確保重力場高階部分的精度，本文將GPS觀測值降采樣至60 s。同時，對星間距離變率數據進行粗差探測，剔除RMS大于6倍中誤差的弧段，從而提高重力場解的可靠性。重力場反演過程中需要估計較多參數，包括衛星初始狀態參數、加速度計校正參數、星間測距經驗參數等。受限于GRACE-FO衛星加速度計的影響，其加速度計數據校正策略與GRACE并不一致。設置星間測距經驗參數與加速度計偏差參數每3 h估計1次，衛星初始狀態參數按弧段每24 h估計1次，加速度計尺度矩陣及重力場參數按月估計。本文采用的參數估計策略如表2所示，下文將簡單介紹引入的星間測距經驗參數及加速度計校正參數。

在星間測距觀測中，引入星間測距經驗參數來吸收未模型化的攝動力及觀測誤差，其數學模型為[10]：

C+Dt+Etcos(u)+Ftsin(u)

(1)

表2 重力場反演過程中參數估計策略

由于重力衛星加速度計數據產品中存在系統的儀器尺度因子和偏差漂移，因此在反演過程中必須對上述兩類參數進行確定[11]。加速度計數據的校正在科學參考框架(science reference frame, SRF)下進行，其表達式為：

acal=Saobs+b

(2)

式中，acal為科學參考框架下校準后的非保守力加速度，aobs為原始加速度計觀測值，S為3×3加速度計尺度矩陣，b為偏差參數。理想條件下，尺度矩陣S為單位陣。由于儀器缺陷導致加速度計軸之間相互影響，尺度矩陣還應包括非對角線元素。因此，本文給出尺度矩陣形式[12]：

(3)

該矩陣有3個對角線元素Sx、Sy、Sz，分別影響3個坐標軸加速度分量大小；α、β、γ為剪切參數，表示加速度計3軸之間在實際情況中由于非正交造成的影響；ζ、ε、δ為旋轉參數，表示加速度計坐標軸與科學參考框架不重合產生的影響。

加速度計偏差參數b會隨時間變化，本文采取如下校正形式：

b=C0+C1t

(4)

式中，C0、C1為待求參數，t為弧段起始處時間。在重力場反演過程中需要選擇合理的加速度計校正策略，從而獲得較優解。

2 結果分析與討論

為探討不同星間測距系統對重力場模型的影響，本文基于KBR和LRI數據分別計算2019-01～2021-08 GRACE-FO重力衛星時變重力場模型(APM_KBR與APM_LRI)時間序列，并與官方機構(center for space research, CSR)新發布的模型(基于KBR數據)進行對比分析，主要包括時變重力場模型、后驗殘差等。

2.1 時變重力場模型對比

當衛星機動或發生其他事件時均會引起LRI數據缺失，為評估相同條件(數據完整性及精度)下各模型解的精度，選取2019-05、2020-05、2021-05三個月(激光測距與微波測距均具有一致的數據連續性)結果的階方差(以大地水準面高表示)，并與CSR模型進行對比(圖1)。需要注意的是，時變重力場模型階方差前30～40階以時變信號為主，而高階部分以噪聲為主。由圖可見，2019-05各模型信號量級符合較好，高階部分的精度與CSR模型結果一致；2020-05各模型信號量級基本一致，但在40階之后APM_LRI的精度變差。這是由于該月利用LRI數據反演重力場時，在粗差探測過程中剔除了3 d星間距離變率精度較差的弧段，導致該月APM_LRI略差于APM_KBR與CSR模型；2021-05與2019-05類似，各模型在信號、噪聲等方面均基本處于同一水平。結果表明，在KBR和LRI數據都未缺失的條件下，反演得到的重力場模型精度水平一致。為方便比對，本文還選取2021-06重力場階方差進行對比，結果如圖2所示。如前文所述，LRI數據受衛星機動次數、載荷運行情況等影響較大，而該月GRACE-FO衛星機動及載荷異常事件較多，導致該月缺失13 d LRI數據[13]。從圖2可以看出，由于缺失數據較多，APM_LRI從30階之后噪聲變大。另外，APM_KBR與CSR模型在30～70階差別較大，說明在數據質量較差的情況下，重力場反演需要更精細的處理策略。

圖1 2019-05、2020-05和2021-05各模型階方差Fig.1 Degree variance of each model in 2019-05, 2020-05 and 2021-05, respectively

圖2 2021-06各模型階方差Fig.2 Degree variance of each model in 2021-06

此外，對2019-01～2021-08整個時間段重力場模型的質量異常趨勢進行統計，結果如圖3和4所示。在該過程中，由于空間分布存在較大的南北條帶誤差，因此本文采用如下后處理方式：300 km高斯平滑，扣除冰川均衡調整效應(glacial isostatic adjustment, GIA)，并替換重力場位系數中的質心項與C20項[14-16]。

從圖3(a)～3(c)可以看出，兩種模型在信號空間分布上與CSR模型基本一致，異常主要集中在南美洲地區、格陵蘭地區和非洲中南部區域。從圖3(d)和3(e)可以看出，APM_KBR和APM_LRI與CSR模型的差異主要集中在赤道附近，并且表現為條帶狀，該差異來源于后處理過程中未將條帶現象完全消除。圖4為兩模型與CSR模型之差及兩模型之差按地理緯度求得的標準差(standard deviation, STD)，可以看出，越靠近赤道附近，條帶誤差越大，模型的質量異常趨勢差異越大；而標準差量級較小，均在1 cm以內。同時統計分析顯示，APM_LRI的標準差略大于APM_KBR，推測其原因可能是APM_LRI中缺失較多星間距離觀測數據，導致模型噪聲較大。

圖3 CSR、APM_KBR與APM_LRI模型質量異常趨勢及CSR、APM_KBR、APM_LRI之間差異Fig.3 Mass abnormal trend of CSR, APM_KBR and APM_LRI, and differences among these models

圖4 按地理緯度統計的CSR、APM_KBR、APM_LRI之間質量變化趨勢之差的標準差Fig.4 Standard deviation of mass change trend of the difference between APM_KBR, APM_LRI and CSR by geographical latitude

為進一步分析重力場模型在區域內的水文信號，選取亞馬遜流域、格陵蘭地區、長江流域及撒哈拉沙漠4個典型區域作為研究對象，分別計算模型在各個流域內水儲量的時間序列，結果如圖5所示。在亞馬遜流域、格陵蘭島地區及長江流域各模型對應的水儲量變化曲線具有較高的一致性，反映出各個時變重力場模型的信號幅值基本一致。撒哈拉沙漠地區的質量變化較小，其信號可在一定程度上反映重力場解算結果的噪聲水平，對比該區域質量變化時間序列可以看到，APM_LRI在2021-06誤差較大，這是由于該月LRI缺失13 d數據，使得2021-06激光數據解算結果精度較差。同時，本文還給出亞馬遜流域、格陵蘭地區、撒哈拉沙漠3個典型區域不同水文模式的數據統計(以等效水柱高表示)，用于量化各模型的一致性水平，結果見表3。分析發現，各模型周年振幅幾乎一致，模型之間差異在0.1 cm以內；APM_KBR在格陵蘭地區的趨勢略小于另外2個模型，但差異也在0.1 cm左右；APM_LRI在撒哈拉沙漠的RMS略大于另外2個模型，結合圖5可知是由于LRI數據缺失所致。

2.2 星間距離變率后驗殘差分析

在實際重力場解算過程中，后驗殘差中仍存

圖5 各模型在4個典型區域內水儲量變化時間序列Fig.5 Time series of water storage changes of each model in four typical regions

表3 各模型在3個典型區域的水儲量變化統計

在背景模型中未充分描述的地球物理信號，其中包括比較典型的潮汐及非潮汐模型的混頻誤差等[3]。通過后驗殘差分析可以在一定程度上反映該月時變重力場解算質量及星間測距系統的測量精度，本文計算2021-05 KBR和LRI數據的后驗殘差，并分別在時域、頻域下對其進行對比分析。需要說明的是，由于不同天之間的后驗殘差結果基本一致，為清晰地展示兩類數據后驗殘差在時域下的特點，采用2021-05-01數據進行時域分析，結果如圖6所示。可以看出，兩類數據的后驗殘差量級均在±5×10-7m/s以內，但LRI精度明顯高于KBR。LRI觀測數據的后驗殘差整月RMS為 4.06×10-8m/s，KBR觀測數據的為8.63×10-8m/s(見表4，單位10-8m/s)。另外，KBR后驗殘差存在較多毛刺，即高頻噪聲較大，而LRI后驗殘差的時間序列更平滑。因此，相比于KBR測距系統，激光測距觀測數據具有更高的觀測精度。此外，為提高樣本選取的可信度，本文還隨機選取并統計了2個月(2019-10、2020-07)星間觀測數據后驗殘差的RMS(見表4，單位10-8m/s)。結果表明，APM_LRI后驗殘差的RMS比APM_KBR下降約50%。

圖6 2021-05-01 KBR與LRI的后驗殘差時間序列Fig.6 Post-fit residuals time series of KBR and LRI measurements on May 1, 2021

表4 KBR與LRI后驗殘差RMS統計

進一步給出2019-10、2020-07和2021-05后驗殘差在頻域下的對比情況(圖7)，值得說明的是，由于GRACE-FO重力衛星繞地周期約為90 min，因此在該周期對應的頻率下存在1 CPR(cycle-per-revolution)的非模型化攝動，而該影響會被星間測距經驗參數所吸收(見式(1))。圖7中黑色虛線為1 CPR頻率，約為0.18 mHz，可以看出，該頻段處2類觀測的后驗殘差振幅很低，即1 CPR處部分非模型化的攝動被星間測距經驗參數吸收。高頻處(大于等于60 CPR) KBR后驗殘差的噪聲大于LRI，這也與圖6的結論相符合。因此，從KBR及LRI星間距離變率后驗殘差的角度來看，LRI的測量精度顯著優于KBR。

最后對2021-05先驗殘差與后驗殘差的全球空間分布情況進行分析，結果如圖8所示。從圖8(a)和8(b)可以看出，APM_KBR與APM_LRI無論是信號分布還是量級差異均較小，這是因為先驗殘差中存在很強的非模型化攝動(1 CPR頻率處)，會將其他質量變化信號湮沒。而后驗殘差中(圖8(c)和8(d))經過星間距離變率經驗參數的擬合，非模型化的誤差大幅降低(圖7)，信號主要分布在亞馬遜流域、非洲中部及南太平洋等區域。另外，APM_KBR高頻噪聲較大，在空間分布中呈斑點狀，這也反映了LRI測量精度優于KBR。

圖7 2019-10、2020-07和2021-05 KBR與LRI的后驗殘差振幅譜Fig.7 The amplitude spectral of post-fit residuals of KBR and LRI measurements in October, 2019, July, 2020 and May, 2021

圖8 APM_KBR與APM_LRI的先驗殘差與后驗殘差空間分布Fig.8 The spatial distribution of pre-fit and post-fit residuals of APM_KBR and APM_LRI

3 結 語

本文通過動力學方法，基于KBR及LRI數據分別解算2019-01～2021-08時變重力場模型，通過對時變重力場模型的階方差、等效水柱高和后驗殘差進行對比分析，得出以下結論：

1)從單月解的精度來看，若KBR與LRI數據完整且質量較好，兩類數據反演得到的重力場精度水平一致，但由于LRI數據受衛星機動、載荷運行情況等的影響，研究時段內有較多激光測距數據缺失；質量異常趨勢項的空間分布及信號量級符合較好，其差異主要集中在赤道附近，應為重力場球諧系數中存在南北條帶誤差所致，統計顯示最大差異不超過1 cm。相較于APM_LRI，APM_KBR與CSR模型的質量異常趨勢一致性更好；各模型在亞馬遜流域、格陵蘭地區及長江流域的水儲量變化時間序列一致性較好，但由于LRI在2021-06缺失較多數據，導致APM_LRI模型在該月精度較差，在撒哈拉地區質量變化的RMS大于APM_KBR與CSR模型。

2)KBR與LRI兩類觀測數據在2021-05的后驗殘差量級在±5×10-7m/s以內，后驗殘差分析結果表明，LRI精度明顯優于KBR。同時，對時間序列進行分析發現，KBR后驗殘差中毛刺現象較多，而LRI的殘差更平滑，其原因為KBR的測量噪聲較大。從3個月的RMS統計結果來看，APM_LRI后驗殘差的RMS比APM_KBR下降約50%；振幅譜中兩者在低頻處差異較小，同時由于本文引入星間測距經驗參數，在1 CPR處兩者的后驗殘差被壓制，高頻處(大于等于60 CPR) KBR的后驗殘差噪聲顯著大于LRI；先驗殘差的空間分布中存在較大的非模型化攝動，導致兩模型差異不明顯，而后驗殘差中APM_KBR存在大量斑點，推測由KBR高頻噪聲所致。

通過對比分析可知，在數據完整且質量較好的條件下，基于LRI與KBR數據計算得到的重力場模型精度一致，激光干涉測距數據的高精度優勢在時變重力場模型反演過程中并無體現，這主要是由目前背景場模型不完備及儀器噪聲(主要源于加速度計)[17]共同導致，而參數策略的選取對重力場精度的影響同樣不可忽略。另外，LRI的測量精度顯著優于KBR，且高精度的星間測距系統對下一代重力衛星計劃意義重大。研究表明[18-19]，可直接利用LRI Level-1B數據探測一些瞬時的質量變化，如地震、洪水、海嘯等。此外，相比于KBR，LRI數據在高頻區域(大于等于60 CPR)具有更高的信噪比，其觀測精度優勢在反演高階地球靜態場模型過程中或可體現。