GRACE/GRACE-FO空窗期的陸地水儲量變化數(shù)據(jù)間斷補償：以全球典型流域為例

徐鵬飛， 蔣濤， 章傳銀， 芮明勝， 劉宇

1 中國測繪科學(xué)研究院， 北京 100830 2 山東科技大學(xué)測繪與空間信息學(xué)院， 山東青島 266590

0 引言

陸地水是不可或缺的重要自然資源，與人類活動密切相關(guān)，由于水資源的分布不均，制約著國家和地區(qū)的高速發(fā)展，嚴重影響我國可持續(xù)發(fā)展戰(zhàn)略(Zhong et al.，2009).因此，陸地水儲量異常TWSA(Terrestrial Water Storage Anomaly)的長期可持續(xù)監(jiān)測，對于研究水循環(huán)過程，合理配置水資源，預(yù)防洪澇、干旱等自然災(zāi)害等有相當(dāng)重要的科學(xué)意義和參考價值(Chen et al.，2009).

2002年初，GRACE(Gravity Recovery and Climate Experiment)衛(wèi)星開始運行，基于GRACE獲取空間中長尺度的時變重力場信息，可反演TWSA，使TWSA的長期可持續(xù)監(jiān)測成為可能，國內(nèi)外學(xué)者基于GRACE數(shù)據(jù)對區(qū)域TWSA進行了大量研究(Tapley et al.，2004；Yeh et al.，2006；Tiwari et al.，2009；Chen et al.，2014；Long et al.，2016).直至2017年10月，GRACE任務(wù)結(jié)束，CSR(Center for Space Research)、GFZ(Helmholtz-Centre Potsdam-German Research Centre for Geosciences)、JPL(Jet Propulsion Laboratory)三所機構(gòu)最新發(fā)布的GRACE Level-2 RL06(Release 06)數(shù)據(jù)更新至2017年6月，但由于衛(wèi)星后期運行存在的質(zhì)量問題，一般2016年8月之后的數(shù)據(jù)不再使用.2018年5月，為延續(xù)GRACE任務(wù)，GRACE-FO(GRACE Follow-On)成功發(fā)射升空，該衛(wèi)星將繼續(xù)監(jiān)測全球重力場的變化.近期CSR、GFZ、JPL公布了最新的GRACE-FO Level-2 RL06球諧系數(shù)數(shù)據(jù)，下文簡稱GRACE-FO SH，截至2021年1月10日，數(shù)據(jù)的時間段為2018年6月至2020年11月，由此GRACE與GRACE-FO之間存在了20多個月的數(shù)據(jù)空窗期.為得到長期連續(xù)的TWSA監(jiān)測結(jié)果，更準確的研究TWSA的年際變化與長期趨勢特征，從而合理利用配置水資源，一種有效的針對GRACE/GRACE-FO數(shù)據(jù)空窗期的TWSA間斷補償方法是十分必要的.

數(shù)據(jù)間斷的補償可通過時間序列分析或機器學(xué)習(xí)等方法預(yù)測實現(xiàn)，國內(nèi)外學(xué)者針對水儲量變化的預(yù)測開展了大量研究工作.Mirzavand和Ghazavi(2015)利用ARMA(Autoregressive Moving Average)、ARIMA(Autoregressive Integrated Moving Average)、SARIMA(Seasonal Autoregressive Integrated Moving Average)對伊朗喀山平原的地下水位進行預(yù)測試驗，研究表明可提前60個月對地下水位取得較好的預(yù)測效果.Adamowski和Chan(2011)利用小波神經(jīng)網(wǎng)絡(luò)WNN(Wavelet Neural Network)預(yù)測加拿大魁北克的地下水位，并與ARIMA方法比較，驗證了該方法的有效性及潛力.Dos Santos和Pereira Filho(2014)利用ANN(Artificial Neural Network)，并引入天氣變量，預(yù)測圣保羅地區(qū)的用水需求，發(fā)現(xiàn)天氣對用水消耗存在反饋作用.另外還有WA-ANN(Wavelet Analysis-Artificial Neural Network)、SVR(Support Vector Regression)等其他方法(Adeli and Jiang，2006；Moosavi et al.，2013；王宇譜等，2013；Tiwari and Adamowski，2015；Al-Zahrani and Abo-Monasar，2015；Mukherjee and Ramachandran，2018).傳統(tǒng)的時間序列分析方法實現(xiàn)簡單，基于最小二乘原理構(gòu)造隨時間變化的函數(shù)，向后遞推，從而預(yù)測后續(xù)序列.然而該方法對于信號成分復(fù)雜的時間序列的預(yù)測效果較差，且易受到信號中噪聲成分的影響，從而預(yù)測失真.相較而言，ANN、WA-ANN、SVR算法的精度較高，但需額外的物理變量進行約束(如氣象數(shù)據(jù)等)，實現(xiàn)困難，計算效率較低.綜上所述，本文將引入奇異譜分析SSA(Singular Spectrum Analysis)+ARMA組合方法(Shen et al.，2018；牛余朋等，2020)，預(yù)測GRACE/GRACE-FO空窗期的TWSA，從而對TWSA的間斷進行補償.GRACE運行多年，積累了較為充足的樣本，且近期GRACE-FO數(shù)據(jù)的發(fā)布，為間斷補償?shù)木仍u定提供了可靠依據(jù).

SSA作為一種時間序列信號的處理方法，在GNSS(Global Navigation Satellite System)數(shù)據(jù)處理等領(lǐng)域應(yīng)用較多(王解先等，2013；Kondrashov and Berloff，2015；肖勝紅等，2016)，但在水儲量變化的預(yù)測方面應(yīng)用較少，從有限時間序列的動態(tài)重構(gòu)出發(fā)，結(jié)合特征向量分析(Eigenvector Analysis)，SSA方法可充分識別并提取信號中不同成分(如周期、趨勢、噪聲等)，且無需引入額外的物理變量，實現(xiàn)較為簡易，特別適合分析預(yù)測具有周期特征的TWSA時間序列.Li等(2019)利用SSA方法預(yù)測了中國多個地區(qū)的水儲量變化，取得了較好的效果，但關(guān)于如何選擇訓(xùn)練樣本、檢驗樣本提高預(yù)測模型精度的討論尚不充分.因此，本文將選取TWSA周期特征明顯的全球典型流域為例，驗證SSA方法在典型流域的適用性及間斷補償效果，并進一步討論如何選擇訓(xùn)練樣本、檢驗樣本構(gòu)建更適用于間斷補償?shù)念A(yù)測模型，從而提高GRACE/GRACE-FO數(shù)據(jù)空窗期TWSA間斷補償結(jié)果的精度.此外，雖然SSA方法可以有效對TWSA序列的周期項等主要成分進行分解與重構(gòu)，但剩余隨機序列中依然存在部分有效信息，而剩余隨機序列可以利用ARMA模型進行擬合遞推.為此，本文引入了SSA+ARMA組合方法進一步提高間斷補償結(jié)果的精度.從全球選取亞馬遜流域、多瑙河流域、恒河流域、密西西比河流域、尼日爾河流域、伏爾加河流域、葉尼塞河流域、長江流域這八個典型流域進行試驗，并分別利用SSA+ARMA、SSA、ARMA和WNN算法進行對比分析，驗證SSA+ARMA方法的適用性和有效性.后續(xù)采用GRACE-FO SH反演TWSA、Mascons、GLDAS(Global Land Data Assimilation System)數(shù)據(jù)，評定基于SSA+ARMA方法的TWSA間斷補償結(jié)果精度，驗證方法的可靠性.

1 數(shù)據(jù)

1.1 GRACE/GRACE-FO level-2 RL06

本文使用CSR的GRACE/GRACE-FO RL06 大地水準面球諧模型GSM(Geoid Spherical Harmonic Model)球諧系數(shù)文件，階數(shù)為60階，下文簡稱GRACE SH.選取GRACE SH的時間段為：2003年1月—2016年8月，共計149個月數(shù)據(jù)(部分月數(shù)據(jù)缺失).截至2021年1月10日，新發(fā)布的GRACE-FO SH時間段為：2018年6月—2020年11月，共計28個月數(shù)據(jù)(2018年8月、2018年9月數(shù)據(jù)缺失).與RL05數(shù)據(jù)一樣，這些數(shù)據(jù)扣除了固體潮汐、海洋潮汐、非潮汐的大氣和海洋信號(Bettadpur，2003).與RL05相比，RL06數(shù)據(jù)在參數(shù)選擇、算法設(shè)計、數(shù)據(jù)編輯、背景場模型構(gòu)建等方面進行了優(yōu)化，“條帶”誤差與截斷誤差有明顯的改善(Save，2018；G?ttl et al.，2018).數(shù)據(jù)預(yù)處理階段分別用SLR(Satellite Laser Ranging)提供的C20項和來自Swenson的一階項，替換GSM中的C20項、一階項(Swenson and Wahr，2002；Cheng and Tapley，2004).此外，冰川均衡調(diào)整(GIA)引起的重力場長期變化趨勢利用了三維壓縮地表負荷滯彈響應(yīng)模型進行扣除(Geruo et al.，2013).

1.2 Mascons結(jié)果

為驗證GRACE SH反演TWSA的可靠性及對間斷補償結(jié)果進行精度評定，本文選取JPL、CSR、GFZ的Mascons RL05以及JPL、CSR的Mascons RL06數(shù)據(jù)進行比較.Mascons RL05數(shù)據(jù)的分辨率為1°×1°，每月一值，數(shù)據(jù)時間段為2003年1月至2017年1月.JPL發(fā)布的Mascons RL06數(shù)據(jù)的分辨率為0.5°×0.5°，每月一值.CSR發(fā)布的Mascons RL06數(shù)據(jù)的分辨率為0.25°×0.25°，每月一值.截至2021年1月10日，基于GRACE-FO數(shù)據(jù)的Mascons RL06結(jié)果已經(jīng)更新至2020年10月.經(jīng)過泄漏誤差的修復(fù)、GIA改正等一系列處理后，具有更好的精度，因此廣泛應(yīng)用于GRACE SH反演TWSA結(jié)果精度的驗證分析(Scanlon et al.，2016).其中Mascons RL05采用的GIA改正模型為三維壓縮地表負荷滯彈響應(yīng)模型，與本文對GRACE SH的處理一致.而Mascons RL06采用的ICE6G-D模型(Purcell et al.，2016).二者模型存在差異，但差異主要體現(xiàn)在南極地區(qū)，在本文所選的研究流域中差異較小(馬超等，2016).

1.3 GLDAS水文模型

本文選用GLDAS V2.1的NOAH模型，時間段為2003年1月—2020年2月，其分辨率為0.25°×0.25°.該模型考慮了0～200 cm的土壤水、植物冠層地表水、雪水當(dāng)量，不包含地下水部分.本文將利用GLDAS數(shù)據(jù)求解單尺度因子，一定程度上修復(fù)球諧系數(shù)濾波引起的泄漏誤差，以及驗證補償結(jié)果的可靠性.

2 原理及方法

2.1 GRACE SH反演TWSA方法

GRACE SH可反演流域內(nèi)等效水高EWH(Equivalent Water Height)變化，然而GRACE重力場模型中存在高階噪聲及“條帶”誤差，因此本文引入了半徑300 km的FAN濾波和P3M15的去相關(guān)濾波的組合濾波方法，從而反演得到八個典型流域的TWSA，具體公式如下(Wahr et al.，1998； Han et al.，2005；Zhang et al.，2009；詹金剛等，2011)：

×[ΔCn mcos(mλ)+ΔSn msin(mλ)]

(1)

2.2 SSA原理

SSA本質(zhì)是針對一維時間序列的主成分分析方法，對于X=(x1,x2,…,xN)，SSA方法時間序列分析過程主要分為分解與重構(gòu)兩個方面，具體過程如下(Vautard et al.，1992；Chen et al.，2013)：

(1)軌跡矩陣

(2)

式中Xi=(xi,…,xi+M-1)T(1≤i≤K).X為M×K階軌跡矩陣.

(2)奇異值分解

X=X1+X2+…+Xd,

(3)

(3)分組

將d個初等矩陣分成m個互不相交的子集I1,I2,…,Im，設(shè)I={i1,i2,…,ip}，則對應(yīng)于子集I的合成矩陣可表示為：

XI=Xi1+Xi2+…+Xip,

(4)

將矩陣X以合成矩陣形式表示為：

X=XI1+XI2+…+XIm,

(5)

因此，分組的過程即是對I1,I2,…,Im的確定過程.

(4)時間序列重構(gòu)

(6)

重構(gòu)成分能夠反映時間序列的變化特征，各重構(gòu)成分功率譜密度的相對關(guān)系可用于有效信號提取與信噪分離.

2.3 SSA迭代預(yù)測方法與SSA+ARMA方法

本文采用的SSA迭代預(yù)測算法由SSA迭代插值法發(fā)展而來，二者的原理是一樣的(Schoellhamer，2001；Beckers and Rixen，2003；Kondrashov et al.，2010；郭金運等，2018).將GRACE SH反演的TWSA時間序列作為訓(xùn)練樣本，通過SSA分解對時間序列進行重構(gòu).迭代預(yù)測的具體步驟如下：

(1)設(shè)定預(yù)測時長n，將n個0值加到訓(xùn)練樣本(長度為N)的末尾處，構(gòu)建總長度為N+n的新序列.

(2)對長度為N+n的新序列進行SSA分解，將第一個重構(gòu)成分(RC1)末尾處的n個值替換新序列，作為預(yù)測數(shù)據(jù).之后循環(huán)該過程，最終前后兩次SSA分解的RC1之間的殘差小于閾值，根據(jù)經(jīng)驗閾值取0.01 mm(Shen et al.，2018；周茂盛等，2018).

(3)當(dāng)步驟(2)循環(huán)過程結(jié)束時，繼續(xù)通過SSA分解添加RC2，通過線性疊加RC1和RC2末尾處的n個預(yù)測值構(gòu)建新序列，循環(huán)該過程直至前后兩次分解的RC1+RC2之間的差值小于閾值.

(4)重復(fù)上述過程，直至用K個RC線性疊加(RC1+RC2+…+RCk)末尾處的n個預(yù)測值，且滿足閾值條件，則序列末尾n個值即為TWSA的預(yù)測值.

以預(yù)測時序與檢驗樣本符合精度最高為準則，結(jié)合W-correlation(Weighted Correlation)及FFT(Fast Fourier Transform)方法進行有效性和適用性檢驗(見3.2節(jié))，測試最佳窗口長度M和重構(gòu)階數(shù)K，并以此窗口長度M與重構(gòu)階數(shù)K，可以得到基于SSA方法TWSA預(yù)測結(jié)果.同時RC1+RC2+…+RCk的前N個值為與訓(xùn)練樣本對應(yīng)的SSA重構(gòu)序列，將訓(xùn)練樣本扣除SSA重構(gòu)序列后得到剩余隨機序列，建立ARMA(p，q)對剩余序列進行外推預(yù)測，SSA預(yù)測結(jié)果與ARMA外推結(jié)果之和為SSA+ARMA的間斷補償結(jié)果，此為SSA+ARMA方法(Shen et al.，2018；牛余朋等，2020).

2.4 精度評定指標(biāo)

RMSE(Root Mean Square Error)：

(7)

RMSE數(shù)值越小，表示模擬值序列與觀測值序列的誤差越小，精度越高.

相關(guān)系數(shù)R(Correlation Coefficient):

(8)

R取值在[-1,1]之間，數(shù)值越大，表明兩個序列的周期相位相關(guān)性越高.

納什效率系數(shù)NSE(Nash-Sutcliffe Efficiency Coefficient):

(9)

NSE取值在(-∞,1]之間，越接近1，表示預(yù)測模型可信度高.與R相比，NSE更能反映振幅相關(guān)性，因此被廣泛應(yīng)用于水文模型模擬結(jié)果的驗證(Merz and Bl?schl，2004).

3 結(jié)果分析

3.1 GRACE SH反演TWSA與Mascons比較

為方便表述，本文以相關(guān)代號表示各流域——亞馬遜流域(AZRB)、多瑙河流域(DNRB)、恒河流域(GNRB)、密西西比河流域(MSRB)、尼日爾河流域(NGRB)、伏爾加河流域(VGRB)、葉尼塞河流域(YNSRB)、長江流域(YZRB).八個流域邊界范圍如圖1所示，數(shù)據(jù)來自HydroSHEDS數(shù)據(jù)集(https:∥hydrosheds.org/page/overview).

首先利用CSR提供的2003年1月—2016年8月的GRACE RL06 SH反演各流域的TWSA序列.其中有15個月的球諧系數(shù)數(shù)據(jù)缺失，需對其進行插值補全.傳統(tǒng)方法為利用前后兩年或多年該月份球諧系數(shù)的均值代替，而本文將采用SSA迭代插值方法對TWSA結(jié)果進行插值補全(Shen et al.，2018；周茂盛等，2018).以MSRB為例，結(jié)果如圖2所示，圖中為傳統(tǒng)的前后兩年球諧系數(shù)均值替代法與SSA迭代插值方法的對比，可以看出SSA方法的插值結(jié)果更符合周期變化的特征，這有利于后續(xù)SSA方法的周期分解與重構(gòu).而傳統(tǒng)方法有明顯的跳變，引入了較大的誤差，特別2012年10月二者結(jié)果的差值達6.53 cm.為此，本文將采用SSA迭代插值方法對GRACE SH反演TWSA和Mascons結(jié)果的缺失數(shù)據(jù)補全.

圖2 基于SSA方法與SH平均方法的TWSA插值結(jié)果Fig.2 The interpolation results of TWSA using SSA and SH Mean

前文已知，本文采用組合濾波方法來抑制噪聲和“條帶”誤差，但也造成了信號幅值的減小，空間分布衰減，即泄漏誤差.為了恢復(fù)相對真實的TWSA，本文采用基于GLDAS數(shù)據(jù)的單尺度因子法，修復(fù)泄漏誤差，該方法將GLDAS格網(wǎng)數(shù)據(jù)球諧展開并截取至60階，采用與GRACE相同的處理流程得到區(qū)域TWSA序列，基于最小二乘原理求解TWSA時間序列與該區(qū)域原始GLDAS時間序列的比例系數(shù)(馮偉等，2012).由表1可知，不同流域的尺度因子存在差異，將研究流域GRACE SH反演的TWSA序列乘以對應(yīng)的尺度因子，作為最終的反演結(jié)果.為驗證GRACE SH反演TWSA結(jié)果的準確性和可靠性，本文將八個典型流域的TWSA分別與三所機構(gòu)的Mascons RL05和JPL、CSR的Mascons RL06結(jié)果進行比較.

表1 不同研究區(qū)域的單尺度因子Table 1 Single scale factor in different study areas

如圖3所示，GRACE SH反演的TWSA結(jié)果與Mascons結(jié)果相比，無論是振幅還是周期相位均有較強的一致性，這表明二者結(jié)果吻合性良好，驗證了本文GRACE SH反演TWSA結(jié)果的有效性和可靠性，且三所機構(gòu)不同版本間的Mascons結(jié)果差異也不明顯.八個流域TWSA序列均有一定的周期信號特征，其中AZRB、NGRB(圖3a，e)的TWSA序列周期特征最為明顯，規(guī)律性更強，而DNRB、MSRB、YZRB(圖3b，d，h)的信號成分復(fù)雜、周期特征相對較弱，序列中蘊含了更多的頻率信息.

圖3 GRACE SH反演研究區(qū)域的TWSA與Mascons結(jié)果比較Fig.3 Comparison of TWSA derived from GRACE SH with Mascons in study areas

本文通過計算八個流域的RMSE(cm)、NSE、R，定量評價GRACE SH反演TWSA結(jié)果與Mascons結(jié)果的符合精度.表2可知，GRACE SH反演TWSA結(jié)果與Mascons結(jié)果有較高的符合精度，大多數(shù)流域的RMSE在2 cm以內(nèi)，且NSE、R在0.95以上.其中在AZRB、NGRB的符合精度最高，NSE、R都達到了0.99，RMSE也在1.5 cm之內(nèi).由NSE可知，在MSRB中，GRACE SH反演的TWSA與Mascons結(jié)果有較明顯的差異，特別是RL06的Mascons結(jié)果.綜合考慮，CSR的Mascons RL05、Mascons RL06結(jié)果與GRACE SH反演結(jié)果吻合程度最高，這與本文選取CSR的GRACE SH數(shù)據(jù)反演TWSA有關(guān).由于Mascons RL05數(shù)據(jù)的時間范圍所限，后續(xù)將采用CSR的Mascons RL06數(shù)據(jù)與TWSA間斷補償結(jié)果對比驗證.

表2 GRACE SH反演研究區(qū)域的TWSA的精度評定Table 2 Accuracy assessment of regional TWSA derived from GRACE SH

3.2 SSA方法可行性分析

前文已知，在利用SSA方法預(yù)測前，需確定最佳窗口長度M和重構(gòu)階數(shù)K，本文以周期特征明顯的AZRB的TWSA為例，詳細說明該過程.為此將2003年1月至2012年12月的GRACE SH反演TWSA結(jié)果作為訓(xùn)練樣本，共計120個數(shù)據(jù).根據(jù)經(jīng)驗，GRACE SH反演的TWSA結(jié)果中可能存在3個月周期信號，因此在SSA分解時，將窗口長度M設(shè)置為3的倍數(shù)，連續(xù)測試窗口長度M，從而對分解的重構(gòu)成分進行周期對分離.以預(yù)測時序與檢驗樣本(2013年1月至2016年8月數(shù)據(jù))符合精度最高為準則，不斷測試，確定AZRB區(qū)域的最佳窗口長度M為48，重構(gòu)階數(shù)K為8.為定量描述各重構(gòu)成分之間的相關(guān)性，檢驗?zāi)Ｐ偷挠行院瓦m用性，本文采用W-correlation進行分析，結(jié)果如圖4所示，圖4a可知，前8階的兩兩相鄰的重構(gòu)成分相關(guān)系數(shù)均高于0.71，故前8階重構(gòu)成分中，相鄰的兩個階項可能為周期對(如RC1與RC2、RC3與RC4、RC5與RC6等均具有相同的周期).圖4b中，前8階特征值的貢獻率達98.29%.

圖4 TWSA時序SSA分解RCs的 W-correlation分析結(jié)果Fig.4 W-correlation analysis of the RCs of TWSA with SSA

本文采用FFT方法探測各重構(gòu)成分的周期特征，將具有相同周期信號特征的重構(gòu)成分進行分組(Kay and Marple，1981).圖5左側(cè)為相鄰兩個RC之和的重構(gòu)序列，右側(cè)為該序列的功率譜密度.由圖可知，RC1+RC2為具有年周期特征的信號，且其功率譜密度最大，對TWSA序列的貢獻率最高，其振幅為：-21.92至22.35 cm；RC3+RC4為具有36個月周期特征的信號；RC5+RC6為具有18個月周期特征信號；RC7+RC8為具有6個月周期特征信號.因此本文選取前8階重構(gòu)成分對AZRB區(qū)域2003年1月至2012年12月的TWSA序列進行重構(gòu)是適用的.

圖5 TWSA時序SSA分解RCs的FFT周期探測結(jié)果Fig.5 FFT period detection results of the RCs of TWSA with SSA

如圖6所示，將訓(xùn)練樣本TWSA序列與SSA重構(gòu)序列進行對比，兩種結(jié)果吻合很好，特別是在周期相位上有很強的一致性，殘差振幅也在2 cm左右，其符合精度分別為RMSE=1.71 cm、NSE=0.99、R=0.99.表明窗口長度M為48、重構(gòu)階數(shù)K為8的重構(gòu)模型在AZRB是有效的，利用該重構(gòu)模型繼續(xù)生成新的序列，即SSA的預(yù)測值，驗證了SSA迭代預(yù)測方法的可行性.同時，圖6中殘差序列即為SSA+ARMA方法中的剩余隨機序列，利用ARMA模型對其進行外推，并與SSA方法的預(yù)測值相加，即為SSA+ARMA的間斷補償結(jié)果.不同流域的SSA預(yù)測模型窗口長度M和重構(gòu)階數(shù)K存在差異，根據(jù)以上過程最終確定各流域的SSA預(yù)測模型參數(shù)詳見表4的預(yù)測模型A.

圖6 AZRB區(qū)域TWSA序列與其SSA重構(gòu)序列 結(jié)果對比Fig.6 Comparison of SSA reconstruction results with TWSA in AZRB

3.3 ARMA預(yù)測模型有效性檢驗

同樣以AZRB流域2003年1月至2012年12月的TWSA序列為樣本，描述預(yù)測模型ARMA(p，q)的構(gòu)建過程.首先對TWSA時間序列進行零均值處理，因為ARMA方法主要針對的是平穩(wěn)時間序列，所以本文采用了ADF(Augmented Dickey-Fuller)檢驗對TWSA時間序列的平穩(wěn)性進行檢測，若序列不滿足平穩(wěn)性要求則需要進行差分處理，直至序列平穩(wěn)(Luo et al.，2012；Mirzavand and Ghazavi，2015；牛余朋等，2020).利用AIC(Akaike Information Criterion)準則確定p、q，之后基于最小二乘原理求解模型的未知參數(shù).構(gòu)建模型后，還需有效性檢驗，有效性檢驗標(biāo)準為擬合樣本序列提取后的殘差序列應(yīng)與白噪聲相似.本文在AZRB流域所構(gòu)建ARMA(6，8)模型的有效性檢驗結(jié)果如圖7所示.圖7a為反映擬合序列與樣本序列偏差的殘差序列，振幅在-5.06 cm至7.02 cm之間.圖7b為殘差序列的正態(tài)QQ(Quantile-Quantile)圖，殘差的數(shù)值均分布在紅線附近，且該殘差序列滿足顯著性水平為0.05的Lilliefors假設(shè)檢驗，因而具有近似正態(tài)分布的特征.圖7c、圖7d中橫軸表示滯后階數(shù)，縱軸表示相關(guān)系數(shù)大小，當(dāng)滯后階數(shù)大于1時，殘差序列的自相關(guān)函數(shù)、偏自相關(guān)函數(shù)數(shù)值很小，均在隨機區(qū)間之內(nèi)，且不隨滯后階數(shù)增大而下降，具有獨立無關(guān)的隨機特性.綜上，殘差序列為服從近似正態(tài)分布的獨立無關(guān)隨機變量，滿足白噪聲特性，所以AZRB流域構(gòu)建的ARMA(6，8)是有效的，可以用于預(yù)測.

圖7 AZRB區(qū)域ARMA(6，8)模型檢驗結(jié)果Fig.7 Test results of the ARMA(6,8) model in AZRB

不同流域的ARMA(p，q)模型存在差別，本文對不同流域的參數(shù)p、q進行測試，并進行有效性檢驗，最終確定了八個流域的ARMA模型：AZRB采用ARMA(6，8)，DNRB采用ARMA(7，6)，GNRB采用ARMA(6，7)，MSRB采用ARMA(8，8)，NGRB采用ARMA(9，8)，VGRB采用ARMA(8，8)，YNSRB采用ARMA(10，10)，YZRB采用ARMA(9，10).在SSA+ARMA方法中，各流域剩余隨機序列的信號成分已發(fā)生變化，ARMA模型需要重新確定.

3.4 不同方法的TWSA預(yù)測結(jié)果比較

為比較SSA方法在典型流域的適用性，將前期120個月(2003年1月至2012年12月)的GRACE SH反演TWSA結(jié)果作為訓(xùn)練樣本，分別用SSA+ARMA、SSA、ARMA與WNN方法預(yù)測后續(xù)44個月(2013年1月至2016年8月)的TWSA(Chen et al.，2006).并以同期的GRACE SH反演TWSA作為真值檢驗，評估不同方法的預(yù)測效果，從而選擇適用于預(yù)測TWSA的方法.

如表3所示，四種方法TWSA預(yù)測序列均與真值序列有一定程度的吻合，但不同方法、不同流域之間還存在差異.SSA與SSA+ARMA方法的預(yù)測值與真值的吻合效果最好，且二者相差不大，ARMA方法效果最差.八個流域中，NGRB的TWSA預(yù)測結(jié)果精度是最高的，而DNRB、MSRB、YZRB預(yù)測結(jié)果精度相對較差.NGRB的TWSA時間序列的周期特征明顯，規(guī)律性強，因而四種方法的預(yù)測結(jié)果均與真值時間序列吻合得較好.相反，DNRB、MSRB、YZRB地區(qū)的TWSA序列更加復(fù)雜多變，蘊含更多的信息，因而預(yù)測結(jié)果與真值序列符合相對較差.綜合考慮三種精度指標(biāo)，SSA與SSA+ARMA方法的預(yù)測精度較高且相對穩(wěn)定，除AZRB外，其他多數(shù)流域的RMSE在4 cm以內(nèi)，所有流域的NSE值在0.7以上，R值在0.85以上.WNN方法預(yù)測結(jié)果也有較高的精度，在DNRB、VGRB流域的預(yù)測精度還略優(yōu)于SSA和SSA+ARMA方法，但在MSRB、YNSRB、YZRB預(yù)測精度較低，構(gòu)建的預(yù)測模型適用性較差.相較而言，除在NGRB外，ARMA模型的預(yù)測精度最差，特別是在YZRB中，ARMA模型預(yù)測結(jié)果的NSE只有0.46.

表3 基于不同方法的TWSA預(yù)測精度評定Table 3 Evaluation of TWSA′s prediction accuracy based on different methods

圖8為四種方法不同時段的預(yù)測精度，對多數(shù)流域來說，隨預(yù)測時間與訓(xùn)練樣本時間間隔越來越大，四種方法的預(yù)測精度逐漸下降，前兩年的預(yù)測精度優(yōu)于后兩年的預(yù)測精度.不同時段中，SSA與SSA+ARMA方法在各流域均有較高的預(yù)測精度，且受預(yù)測時間長度逐漸增加的影響較小，特別是在YNSRB和YZRB中，2016年SSA+ARMA方法的RMSE分別為2.72 cm和1.59 cm，明顯優(yōu)于ARMA與WNN方法.ARMA方法預(yù)測精度較差，且受預(yù)測時間長度的影響較大，2016年AZRB的RMSE達到了14.25 cm.在MSRB、YNSRB、YZRB中，WNN方法2016年的預(yù)測精度較差，而AZRB、MSRB中，2016年的預(yù)測精度較高.大部分流域中所有時段SSA+ARMA方法的預(yù)測精度均略優(yōu)于SSA方法，但二者差異較小.

圖8 不同方法的TWSA預(yù)測結(jié)果對比Fig.8 Comparison of TWSA prediction results by different methods

綜上所述，ARMA方法對于周期特征較弱的時間序列預(yù)測效果并不理想，易受到序列中噪聲成分的干擾，導(dǎo)致預(yù)測結(jié)果精度降低.WNN方法為機器學(xué)習(xí)方法，通過不斷調(diào)試學(xué)習(xí)速率、隱含層節(jié)點數(shù)等參數(shù)構(gòu)建預(yù)測模型，計算效率較低，且因為缺少相關(guān)約束，某些流域構(gòu)建的預(yù)測模型并不適用.在后續(xù)的研究中，若引入降水等氣象數(shù)據(jù)作為約束，WNN方法構(gòu)建模型的效率及預(yù)測精度可能會進一步提升，同樣具有很大的潛力.而結(jié)合特征向量分析，SSA方法可有效識別并提取時間序列中不同成分，受噪聲的影響小，能從復(fù)雜的信號成分中提取更多有用信息，特別適合分析預(yù)測具有復(fù)雜成分或信號周期特征明顯的TWSA時間序列.

由表3可知，SSA與SSA+ARMA方法相比，SSA+ARMA方法要略優(yōu)于SSA方法，且主要反映在RMSE上，其精度有0.01至0.1 cm量級的提升.這主要與各流域SSA重構(gòu)序列的重構(gòu)階數(shù)有關(guān)，如圖4b和圖6所示，AZRB的前8階重構(gòu)成分的貢獻率達98.29%，因此剩余時間序列包含的信息較少，精度提升并不明顯.但各流域的SSA預(yù)測模型存在差異，對于SSA重構(gòu)序列貢獻率較低的流域，SSA+ARMA方法會有較明顯的提升，如GNRB的重構(gòu)階數(shù)K為4(表4中預(yù)測模型A)，因而SSA重構(gòu)成分的貢獻率會較低，與SSA方法相比，其RMSE降低了0.35 cm，NSE、R均提升了0.2.

3.5 SSA+ARMA方法的預(yù)測模型構(gòu)建及補償結(jié)果精度評定

SSA+ARMA方法的預(yù)測精度主要還是由SSA方法的預(yù)測精度影響的.因而本文將討論如何選擇訓(xùn)練樣本、檢驗樣本構(gòu)建更適用的SSA預(yù)測模型，提高GRACE/GRACE-FO數(shù)據(jù)空窗期TWSA補償結(jié)果的精度.在之前的預(yù)測試驗中，SSA迭代預(yù)測方法構(gòu)建的各流域預(yù)測模型如表4中預(yù)測模型A所示，但該模型存在一定的局限性，訓(xùn)練樣本選取的2003年1月至2012年12月的120個樣本，并未充分考慮全部164個樣本，且訓(xùn)練樣本與GRACE/GRACE-FO間斷期的時間跨度較大.因此構(gòu)建的預(yù)測模型并不一定適用于數(shù)據(jù)間斷的補償.GRACE-FO數(shù)據(jù)作為GRACE的延續(xù)，本文對二者的數(shù)據(jù)處理過程一致，因此本文將GRACE-FO SH反演TWSA結(jié)果看作真值.為確定更準確、更適用于GRACE/GRACE-FO空窗期TWSA間斷補償?shù)腟SA模型，本文引入GRACE-FO SH反演TWSA及Mascons RL06結(jié)果作為檢驗樣本，提出了另外兩種構(gòu)建預(yù)測模型的方案，三種方案區(qū)別如圖9.預(yù)測模型A與預(yù)測模型B和C相比，分析了訓(xùn)練樣本對預(yù)測模型構(gòu)建的影響，其中預(yù)測模型B和C充分考慮了164個訓(xùn)練樣本.預(yù)測模型B和C的區(qū)別在于檢驗樣本的選擇，上文已知，隨預(yù)測時間與訓(xùn)練樣本時間間隔越來越大，預(yù)測精度逐漸下降.因而預(yù)測模型B引入GRACE-FO SH反演TWSA結(jié)果作為檢驗樣本，對間斷補償結(jié)果序列的末尾進行約束，從而得到更適用于間斷補償?shù)念A(yù)測模型.相較而言，預(yù)測模型C除了引入GRACE-FO SH反演TWSA結(jié)果，同時還考慮了Mascons RL06結(jié)果，可以對間斷補償結(jié)果序列的首尾兩端均進行約束.以預(yù)測序列與檢驗樣本符合精度最高為準則，確定了各流域的預(yù)測模型A、B、C的參數(shù)，即窗口長度M和重構(gòu)階數(shù)K.如表4所示，多數(shù)流域依據(jù)不同方案構(gòu)建的SSA預(yù)測模型參數(shù)存在差異，特別是重構(gòu)階數(shù)K，在預(yù)測模型B和C中，AZRB的K為4，表示該流域SSA的重構(gòu)序列中有4個重構(gòu)成分，因而剩余隨機序列中會包含較多的有效信息，相反NGRB的K為12，SSA方法對該流域TWSA序列的分解與重構(gòu)更加充分，剩余隨機序列中包含的有效信息較少.在VGRB、YZRB中，方案二和方案三構(gòu)建的預(yù)測模型B與C是一樣的.

圖9 三種方案的流程圖Fig.9 Flow chart of the three programmes

表4 基于不同方案的SSA預(yù)測模型Table 4 Prediction models of SSA based on different programmes

本文將GRACE SH反演的2003年1月至2016年8月TWSA序列作為樣本，基于SSA方法分別采用預(yù)測模型A、B、C，補償2016年9月至2020年2月的TWSA間斷結(jié)果.為檢驗補償結(jié)果的準確性及可靠性，受限于數(shù)據(jù)的包含范圍，本文選取2016年9月至2017年6月的CSR的Mascons RL06結(jié)果，2016年9月至2020年2月的GLDAS結(jié)果，2018年10月至2020年2月的GRACE-FO SH反演TWSA結(jié)果進行對比驗證.其中，GLDAS結(jié)果為對下載的格網(wǎng)序列球諧展開并截取至60階后，進行與GRACE數(shù)據(jù)一樣的處理流程后的結(jié)果.同樣的，各類數(shù)據(jù)均扣除了對應(yīng)數(shù)據(jù)2003年1月至2016年8月的均值，其中CSR提供的GRACE-FO RL06 SH扣除了該機構(gòu)GRACE RL06 SH的2003年1月至2016年8月均值.三種方案在間斷數(shù)據(jù)補償階段是一樣的，均采用2003年1月至2016年8月的164個TWSA結(jié)果為樣本，GRACE-FO SH結(jié)果與Mascons結(jié)果未參與實際預(yù)測過程，因此依據(jù)GRACE-FO SH結(jié)果與Mascons結(jié)果評定補償結(jié)果精度是有意義的.

圖10可知，除YZRB中預(yù)測模型A的補償結(jié)果外，其他間斷補償結(jié)果均取得了較好的效果.基于SSA方法的TWSA間斷補償結(jié)果與Mascons、GRACE-FO SH結(jié)果吻合得較好，這驗證了補償結(jié)果的可靠性.表5中RMSE可知，TWSA間斷補償結(jié)果與GLDAS結(jié)果在振幅上有明顯的差異，但由R可知，GLDAS結(jié)果與TWSA間斷補償結(jié)果在周期和相位變化上有較強的一致性，這亦能反映間斷補償結(jié)果的可靠性.振幅的偏差相對較大是因為GLDAS中只包含了0～200 cm土壤水、植物冠層地表水、積雪的變化，不包含地下水變化，而GRACE則反映了陸地總儲水量的變化和固體地球可能的物理遷移(Chao，2016).因而在后續(xù)主要依據(jù)GRACE-FO SH反演的TWSA及Mascons結(jié)果進行精度評定.

圖10 基于SSA方法TWSA間斷補償結(jié)果與GRACE-FO SH，Mascons和GLDAS結(jié)果的比較Fig.10 Comparison of TWSA gap compensation results by SSA with GRACE-FO SH，Mascons and GLDAS results

表5可知，各流域中預(yù)測模型A、B、C的間斷補償結(jié)果精度存在差異.預(yù)測模型A的符合精度相對較差，表明在構(gòu)建預(yù)測模型過程中，采用更充足的訓(xùn)練樣本有利于構(gòu)建更適用的預(yù)測模型.預(yù)測模型B、C均取得了較好的預(yù)測效果，但在細節(jié)上存在差異，預(yù)測模型B與GRACE-FO SH結(jié)果的符合精度最高，但可能會導(dǎo)致間斷補償結(jié)果序列的前部分失真，如DNRB的間斷補償結(jié)果與Mascons結(jié)果的NSE僅為-0.02；在GNRB中RMSE達10.60 cm.相較而言，預(yù)測模型C對間斷補償結(jié)果的首尾兩端均進行了約束，不會出現(xiàn)這種問題.但值得注意的是，在MSRB中，比較NSE可得，預(yù)測模型B的補償結(jié)果要優(yōu)于預(yù)測模型C.這與本文GRACE SH反演的TWSA在該流域與Mascons結(jié)果存在較大差異有關(guān)(見表2).因此當(dāng)引入Mascons RL06結(jié)果作為檢驗樣本約束間斷補償結(jié)果時，反而會降低補償精度.本文采用的GRACE/GRACE-FO SH均為CSR提供的RL06版本，且進行了相同的數(shù)據(jù)處理流程，可以看為同種數(shù)據(jù)，相反作為成熟產(chǎn)品的Mascons結(jié)果無法保證這一點.因而當(dāng)預(yù)測模型B、C的補償結(jié)果符合精度相當(dāng)時，本文優(yōu)先考慮預(yù)測模型B間斷補償結(jié)果.綜上所述，AZRB、DNRB、GNRB、NGRB采用預(yù)測模型C的SSA間斷補償結(jié)果，MSRB、YNSRB采用預(yù)測模型B的SSA間斷補償結(jié)果，在VGRB、YZRB構(gòu)建的預(yù)測模型B與C是一樣的.

表5 基于SSA方法的TWSA間斷補償結(jié)果精度評定Table 5 Accuracy assessment of TWSA gap compensation results by SSA

續(xù)表5

根據(jù)上文選取的基于SSA方法的間斷補償結(jié)果，本文進一步采用SSA+ARMA方法提升補償精度.表6可得，SSA+ARMA方法的間斷補償結(jié)果整體優(yōu)于SSA方法，且主要體現(xiàn)在RMSE上.其中AZRB的間斷補償結(jié)果的精度提升最為明顯，基于GRACE-FO SH結(jié)果的精度評定中，RMSE降低了1.35 cm，基于Mascons結(jié)果的精度評定中，RMSE降低了1.01 cm.除AZRB外，DNRB、NGRB、VGRB、YZRB的預(yù)測精度也有相對明顯的提升.

表6 基于SSA+ARMA方法的TWSA間斷補償結(jié)果 精度評定Table 6 Accuracy assessment of TWSA gap compensation results by SSA+ARMA

綜上所述，基于SSA+ARMA方法的TWSA間斷補償結(jié)果有較高的精度，與同期GRACE-FO SH結(jié)果相比，除AZRB的RMSE為4.96 cm外，其他流域的RMSE均在4 cm以內(nèi)；除YZRB的NSE、R分別為0.61、0.80外，其他流域的NSE值均在0.80以上，R值均在0.90以上.對TWSA序列相對復(fù)雜多變的DNRB、YZRB(圖3b，h)，間斷補償結(jié)果也取得了較好的效果，特別是YZRB(圖10h)的2018年10月—2019年7月周期特征不明顯，而SSA間斷補償結(jié)果與GRACE-FO SH結(jié)果卻有很強的一致性，這充分體現(xiàn)了SSA方法能識別、提取復(fù)雜成分信號中有用信息的特性.在八個典型流域中，NGRB的間斷補償結(jié)果精度最高，這可能與該流域受人類活動影響少，TWSA序列具有強周期特征有關(guān).

3.6 基于SSA+ARMA方法的TWSA間斷補償結(jié)果試驗

上文三種方案的預(yù)測模型構(gòu)建中，由于部分流域GRACE SH反演的TWSA與Mascon結(jié)果符合較差，因而預(yù)測模型B和C互有優(yōu)勢.但在實際應(yīng)用中，若采用同種數(shù)據(jù)作為訓(xùn)練樣本與檢驗樣本，方案三構(gòu)建的預(yù)測模型取得的間斷補償效果更好.截至2021年1月，基于GRACE-FO的CRS RL06 Mascon數(shù)據(jù)結(jié)果已經(jīng)更新至2020年10月，對基于SSA+ARMA方法的TWSA間斷補償?shù)膶嶋H應(yīng)用提供了充足的數(shù)據(jù)支持.

本文采用方案三的預(yù)測模型構(gòu)建思路，利用2003年1月至2016年8月的Mascon結(jié)果為訓(xùn)練樣本，2016年9月至2017年6月及2018年10月至2020年2月的Mascon結(jié)果作為檢驗樣本構(gòu)建SSA預(yù)測模型.構(gòu)建預(yù)測模型后，在實際應(yīng)用中為得到精度更高的間斷補償結(jié)果，無需考慮精度評定的問題，間斷補償階段應(yīng)盡可能選擇充足的樣本.因而采用2003年1月至2016年8月及2018年10月至2020年10月的Mascon結(jié)果作為樣本，利用SSA+ARMA方法對2016年9月至2018年9月的結(jié)果進行補償.如圖11所示，各流域的TWSA間斷補償結(jié)果均取得了較好的效果，補償結(jié)果與CSR RL06 Mascon結(jié)果連接較為緊密，且無論是趨勢項還是周期特征均與原始時序符合得較好.對于典型流域，利用SSA+ARMA方法可以對GRACE/GRACE-FO空窗期的TWSA間斷進行有效補償.表7為SSA+ARMA與SSA方法補償結(jié)果的統(tǒng)計信息對比.由標(biāo)準差和均方根值可知，SSA+ARMA方法的補償結(jié)果整體略優(yōu)于SSA方法，相較SSA方法，SSA+ARMA方法可以顧及到SSA方法未能充分重構(gòu)的剩余有效信息.

圖11 基于SSA+ARMA方法的TWSA間斷補償結(jié)果Fig.11 TWSA gap compensation results by SSA+ARMA

表7 基于SSA+ARMA和SSA的TWSA間斷 補償結(jié)果統(tǒng)計信息(單位：cm)Table 7 The statistics of TWSA gap compensation results by SSA+ARMA and SSA (Unit: cm)

4 討論

4.1 SSA預(yù)測模型參數(shù)確定

FFT及W-correlation方法均為提高構(gòu)建預(yù)測模型效率的輔助方法，本文最終的模型參數(shù)確定以預(yù)測序列與檢驗樣本的符合精度最高作為準則(包括WNN、ARMA方法).各典型流域均有明顯的周年信號，但半年信號在不同流域存在差異，在構(gòu)建各流域的SSA預(yù)測模型時，重構(gòu)階數(shù)K是否考慮半年信號還應(yīng)以預(yù)測序列與檢驗樣本的符合精度最高為準則.因而本文提出了三種方案來討論如何選取訓(xùn)練樣本、檢驗樣本，構(gòu)建更適用于間斷補償?shù)腟SA預(yù)測模型.引入GRACE-FO SH反演TWSA、Mascons結(jié)果作為檢驗樣本，不僅能使訓(xùn)練樣本更加充足，還能對間斷補償結(jié)果進行約束，防止間斷補償結(jié)果失真.方案一構(gòu)建的預(yù)測模型A與方案二和方案三構(gòu)建的預(yù)測模型相比表明，充足的訓(xùn)練樣本有助于構(gòu)建更適用于間斷補償?shù)腟SA預(yù)測模型.方案二的預(yù)測模型B的補償結(jié)果與GRACE-FO SH反演的TWSA符合得更好，但可能會導(dǎo)致補償結(jié)果前部分失真(如表5中DNRB、GNRB).方案三的預(yù)測模型C對補償結(jié)果的首尾兩端均進行約束，但在3.5節(jié)中，由于GRACE SH反演的TWSA與Mascons符合精度存在差異的問題，部分流域引入Mascons結(jié)果作為檢驗樣本反而可能會導(dǎo)致得不到最佳的預(yù)測模型(如表5中MSRB、YNSRB).而3.6節(jié)中，當(dāng)訓(xùn)練樣本及檢驗樣本均采用Mascon結(jié)果或同種數(shù)據(jù)時，采用方案三構(gòu)建預(yù)測模型的思路能取得更好的間斷補償效果.方案三的主要思想為引入檢驗樣本對間斷補償結(jié)果的首尾兩端進行約束，從而構(gòu)建更適用的預(yù)測模型.在實際應(yīng)用中，訓(xùn)練樣本與檢驗樣本的長度可以根據(jù)研究目的及需求的不同進行調(diào)節(jié).

4.2 SSA與SSA+ARMA方法

SSA+ARMA方法的間斷補償結(jié)果整體優(yōu)于SSA方法，且主要體現(xiàn)在RMSE上，其精度有0.01至1 cm量級的提升.表6可得，SSA+ARMA方法對AZRB、DNRB、VGRB、YZRB的精度提升相對明顯，特別是AZRB，基于GRACE-FO SH結(jié)果的精度評定中，該流域的RMSE降低了1.01 cm，基于Mascons結(jié)果的精度評定中，RMSE降低了1.35 cm.表4可知，AZRB的SSA預(yù)測模型C重構(gòu)階數(shù)K為4，扣除SSA重構(gòu)序列后的剩余隨機時間序列還保留了更多的有效信息，因而再利用ARMA進行外推可以有效提升最終的間斷補償結(jié)果的精度.在DNRB、VGRB、YZRB中，預(yù)測模型的重構(gòu)階數(shù)K同樣較小.相反，NGRB的預(yù)測模型的K為12，因為該流域的TWSA信號周期特征最為顯著，SSA方法能充分分解與重構(gòu)TWSA時間序列，剩余隨機時間序列包含的有效信號較少，所以SSA+ARMA方法的補償結(jié)果精度沒有明顯的提升.當(dāng)SSA方法無法充分的重構(gòu)TWSA信號，即預(yù)測模型的重構(gòu)階數(shù)K較小，剩余信號中保留更多的有效信息時，SSA+ARMA方法能更明顯的提升間斷補償結(jié)果的精度.如表7所示，由最大值、最小值及均值可知，SSA與SSA+ARMA補償結(jié)果的振幅差異較小，大部分差異在0.1 cm的量級.由標(biāo)準差和均方根值可知， SSA+ARMA方法的補償結(jié)果整體略優(yōu)于SSA方法，精度提升量級在0.01至0.1 cm，其中在AZRB的精度提升最為明顯，標(biāo)準差降低了0.78 cm，均方根值降低了0.79 cm.相較SSA方法，SSA+ARMA方法可以顧及到SSA方法未能充分重構(gòu)的剩余有效信息.在以后的實際應(yīng)用中，對于周期特征弱的非典型流域地區(qū)的TWSA序列進行間斷補償時，與SSA方法相比，SSA+ARMA方法可能會取得更好的補償精度提升效果.

5 結(jié)論

本文引入了SSA+ARMA方法對GRACE/GRACE-FO空窗期的TWSA間斷結(jié)果補償?shù)姆椒?聯(lián)合GRACE SH、GRACE-FO SH、Mascons、GLDAS數(shù)據(jù)，在AZRB、DNRB、GNRB、MSRB、NGRB、VGRB、YNSRB、YZRB八個流域進行了不同方法的預(yù)測試驗及SSA+ARMA方法的TWSA間斷補償試驗，分別驗證了SSA+ARMA方法的適用性及有效性.結(jié)論如下：

(1)傳統(tǒng)利用多年或前后兩年球諧系數(shù)均值替代缺失數(shù)據(jù)的方法會引入較大的誤差，而本文采用SSA迭代插值方法的插值結(jié)果更符合TWSA序列周期變化的特征，更利于后續(xù)SSA方法的周期分解與重構(gòu).本文GRACE SH反演的TWSA與CSR、JPL、GFZ的Mascons結(jié)果有較好的一致性，大多數(shù)流域的RMSE在2 cm以內(nèi)，NSE、R在0.95以上，驗證了反演結(jié)果的準確性及可靠性.CSR的Mascons RL05、Mascons RL06結(jié)果與GRACE SH反演TWSA結(jié)果吻合程度最高.

(2)通過四種方法的預(yù)測對比試驗發(fā)現(xiàn)，在多數(shù)流域中，隨預(yù)測時間長度的增加，四種方法的預(yù)測精度逐漸下降，但SSA與SSA+ARMA方法受其影響較小.SSA與SSA+ARMA方法的預(yù)測精度較高且相對穩(wěn)定，相較而言ARMA模型的精度最差.ARMA方法對于周期特征較弱的序列預(yù)測效果并不理想(如DNRB、MSRB、YZRB)，而SSA與SSA+ARMA方法取得了較好的效果，這充分體現(xiàn)SSA方法能識別、提取復(fù)雜成分信號中有用信息的特性.

(3)SSA+ARMA方法的補償精度主要受SSA方法補償精度的影響.引入GRACE-FO SH反演TWSA或Mascons結(jié)果作為檢驗樣本，采用更充足的訓(xùn)練樣本，可以構(gòu)建更適用于間斷補償?shù)腟SA預(yù)測模型，進而能有效提高最終TWSA間斷補償結(jié)果的精度.在實際應(yīng)用中，若引入同種數(shù)據(jù)作為檢驗樣本，對間斷補償結(jié)果的首尾兩端進行約束，構(gòu)建的預(yù)測模型能取得更好的間斷補償效果.

(4)基于SSA+ARMA方法的TWSA間斷補償結(jié)果有較高的精度，與GRACE-FO SH結(jié)果有很強的一致性.除AZRB的RMSE為4.96 cm外，其他流域的RMSE均在4 cm以內(nèi)；除YZRB的NSE、R分別為0.61、0.80外，其他流域的NSE值均在0.80以上，R值均在0.90以上.

(5)在實際應(yīng)用的試驗中，SSA+ARMA方法的間斷補償結(jié)果優(yōu)于SSA方法，且主要體現(xiàn)在標(biāo)準差和均方根值上，其精度提升可達0.1 cm量級.相較SSA方法，SSA+ARMA方法可以顧及到SSA方法未能充分重構(gòu)的剩余有效信息.

綜上，本文提出基于SSA+ARMA方法的TWSA間斷補償是可行且有效的，對實現(xiàn)TWSA的長期持續(xù)監(jiān)測有重要的參考價值和科學(xué)意義.

致謝感謝匿名審稿專家對本文提出的寶貴意見.