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

1 中國科學(xué)院精密測量科學(xué)與技術(shù)創(chuàng)新研究院大地測量與地球動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，武漢市徐東大街340號(hào)，430077 2 中國科學(xué)院大學(xué)地球與行星科學(xué)學(xué)院，北京市玉泉路19號(hào)甲，100049 3 中山大學(xué)物理與天文學(xué)院，廣東省珠海市大學(xué)路2號(hào)，519082 4 中山大學(xué)測繪科學(xué)與技術(shù)學(xué)院，廣東省珠海市大學(xué)路2號(hào)，519082

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

1 動(dòng)力學(xué)法基本理論

1.1 反演流程

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

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

1.2 參數(shù)估計(jì)策略

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

在星間測距觀測中，引入星間測距經(jīng)驗(yàn)參數(shù)來吸收未模型化的攝動(dòng)力及觀測誤差，其數(shù)學(xué)模型為[10]：

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

(1)

表2 重力場反演過程中參數(shù)估計(jì)策略

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

acal=Saobs+b

(2)

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

(3)

該矩陣有3個(gè)對(duì)角線元素Sx、Sy、Sz，分別影響3個(gè)坐標(biāo)軸加速度分量大小；α、β、γ為剪切參數(shù)，表示加速度計(jì)3軸之間在實(shí)際情況中由于非正交造成的影響；ζ、ε、δ為旋轉(zhuǎn)參數(shù)，表示加速度計(jì)坐標(biāo)軸與科學(xué)參考框架不重合產(chǎn)生的影響。

加速度計(jì)偏差參數(shù)b會(huì)隨時(shí)間變化，本文采取如下校正形式：

b=C0+C1t

(4)

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

2 結(jié)果分析與討論

為探討不同星間測距系統(tǒng)對(duì)重力場模型的影響，本文基于KBR和LRI數(shù)據(jù)分別計(jì)算2019-01～2021-08 GRACE-FO重力衛(wèi)星時(shí)變重力場模型(APM_KBR與APM_LRI)時(shí)間序列，并與官方機(jī)構(gòu)(center for space research, CSR)新發(fā)布的模型(基于KBR數(shù)據(jù))進(jìn)行對(duì)比分析，主要包括時(shí)變重力場模型、后驗(yàn)殘差等。

2.1 時(shí)變重力場模型對(duì)比

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

圖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

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

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

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

圖4 按地理緯度統(tǒng)計(jì)的CSR、APM_KBR、APM_LRI之間質(zhì)量變化趨勢之差的標(biāo)準(zhǔn)差Fig.4 Standard deviation of mass change trend of the difference between APM_KBR, APM_LRI and CSR by geographical latitude

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

2.2 星間距離變率后驗(yàn)殘差分析

在實(shí)際重力場解算過程中，后驗(yàn)殘差中仍存

圖5 各模型在4個(gè)典型區(qū)域內(nèi)水儲(chǔ)量變化時(shí)間序列Fig.5 Time series of water storage changes of each model in four typical regions

表3 各模型在3個(gè)典型區(qū)域的水儲(chǔ)量變化統(tǒng)計(jì)

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

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

表4 KBR與LRI后驗(yàn)殘差RMS統(tǒng)計(jì)

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

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

圖7 2019-10、2020-07和2021-05 KBR與LRI的后驗(yàn)殘差振幅譜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的先驗(yàn)殘差與后驗(yàn)殘差空間分布Fig.8 The spatial distribution of pre-fit and post-fit residuals of APM_KBR and APM_LRI

3 結(jié) 語

本文通過動(dòng)力學(xué)方法，基于KBR及LRI數(shù)據(jù)分別解算2019-01～2021-08時(shí)變重力場模型，通過對(duì)時(shí)變重力場模型的階方差、等效水柱高和后驗(yàn)殘差進(jìn)行對(duì)比分析，得出以下結(jié)論：

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

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

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