孫煒程,魏德宏,羅朱鍵,禤鍵豪,張興福

(1. 廣東工業(yè)大學(xué)土木與交通工程學(xué)院,廣東 廣州 510006; 2. 同濟大學(xué)測繪與地理信息學(xué)院,上海 200092)

GRACE和GRACE-FO重力衛(wèi)星為全球和區(qū)域陸地水儲量變化監(jiān)測提供了數(shù)據(jù)支撐,但GRACE與GRACE-FO之間存在約1年的銜接空白,且當(dāng)前GRACE-FO的D衛(wèi)星加速度計已出現(xiàn)故障。SWARM衛(wèi)星提供了新的對地觀測數(shù)據(jù),可作為GRACE和GRACE-FO數(shù)據(jù)填補或替換的有效手段。對于SWARM數(shù)據(jù)解算,ASU與IGG等科研機構(gòu)分別采用去相關(guān)加速度法和短弧積分法等方式反演SWARM時變模型[1],COST-G機構(gòu)則基于VCE方法對多機構(gòu)SWARM時變模型進行融合[2]?？梢妵H機構(gòu)足夠重視SWARM衛(wèi)星在時變重力場反演中的應(yīng)用。

近年來,SWARM時變模型逐漸得到廣泛應(yīng)用。文獻(xiàn)[2]組合多機構(gòu)SWARM時變模型,對亞馬孫、阿拉斯加等多個流域進行了陸地水儲量反演分析,發(fā)現(xiàn)SWARM與GRACE模型反演結(jié)果在多個流域的相關(guān)系數(shù)大于0.6。文獻(xiàn)[3]利用SWARM時變模型探測亞馬孫干旱事件,通過GRACE、水文模型與虛擬水文站數(shù)據(jù)進行比較,驗證了SWARM探測結(jié)果的可靠性。文獻(xiàn)[4]聯(lián)合SWARM衛(wèi)星觀測資料與ARIMA-MC預(yù)測方法,有效填補了GRACE與GRACE-FO在南極與格陵蘭島冰蓋質(zhì)量的間斷空白。由于SWARM在不同區(qū)域的反演精度存在差異性,因此需分析SWARM是否適用于區(qū)域陸地水儲量變化監(jiān)測,為后續(xù)開展SWARM的應(yīng)用提供參考。

松遼流域是我國主要的糧食產(chǎn)區(qū)之一,其中松花江流域曾發(fā)生嚴(yán)重的地下水匱乏[5],遼河流域則因農(nóng)業(yè)灌溉用水導(dǎo)致地下水被過度開發(fā)利用[6]。因此,開展松遼流域的陸地水儲量變化監(jiān)測研究對當(dāng)?shù)氐乃Y源管理與保護具有重要意義。本文選取該流域作為SWARM監(jiān)測陸地水儲量變化的研究對象,并綜合多源數(shù)據(jù)分析SWARM時變模型用于陸地水儲量變化監(jiān)測的適用性。

1 數(shù)據(jù)來源

2015年7月前,SWARM衛(wèi)星的星載設(shè)備尚未得到更新,受衛(wèi)星接收器通道較少與太陽活動等因素影響,其探測數(shù)據(jù)存在較大噪聲[2]。本文選取2015年7月及以后的SWARM衛(wèi)星數(shù)據(jù)進行反演分析。數(shù)據(jù)來源信息見表1。

表1 數(shù)據(jù)來源信息

2 原理與數(shù)據(jù)處理

2.1 時變模型階方差計算

由于采用的時變模型數(shù)量較多,為全面評估模型,計算模型各階次位系數(shù)相對于參考重力場模型位系數(shù)的大地水準(zhǔn)面階方差[7],公式為

(1)

2.2 時變模型反演陸地水儲量變化

利用時變重力場位系數(shù)計算陸地水儲量等效水高的公式為(參數(shù)詳見參考文獻(xiàn)[8])

(2)

計算區(qū)域陸地水儲量變化時間序列ΔTWS的公式為

(3)

式中,N為研究區(qū)域內(nèi)的網(wǎng)格數(shù)量;θ、λ分別為地心余緯與地心經(jīng)度。

本文利用TN-13一階項估計值補充時變模型的一階項位系數(shù)[9],同時結(jié)合SLR測定的C20項與C30項替換時變模型所對應(yīng)的位系數(shù)[10-11],并進行GIA改正[12];采用高斯濾波與Duan去相關(guān)濾波處理,以削弱時變模型的高頻噪聲與條帶誤差[13-14],并通過GLDAS模型計算尺度因子以恢復(fù)信號[15];利用奇異譜分析(SSA)填補GRACE數(shù)據(jù)空缺[16]。此外,由于GRACE與GRACE-FO模型的求解方法基本相同,本文將兩者統(tǒng)稱為GRACE。

2.3 確定SWARM時變模型與濾波方案

選取各機構(gòu)2015年7月—2016年12月的SWARM和GRACE時變模型計算階方差,如圖1所示。結(jié)果表明:ASU、COST-G、IGG的SWARM與GRACE階方差的前10階結(jié)果較接近,上述模型前10階的階方差均值分別為1.213、1.165、1.128、0.955 mm,可見時變模型的前10階位系數(shù)精度較為一致,但隨著階次的增加,SWARM時變模型位系數(shù)的噪聲逐漸增大。

圖1 不同機構(gòu)數(shù)據(jù)的大地水準(zhǔn)面階方差

為進一步分析各SWARM時變模型在松遼流域的精度差異,以GRACE反演結(jié)果為參考進行比較,并結(jié)合信噪比進行綜合評價,從而篩選模型并確認(rèn)濾波方案,對比時段為2015年7月—2020年12月。為了保證數(shù)據(jù)處理的一致性,對SWARM與GRACE時變模型均截斷至40階,并采用相同的濾波處理與泄露誤差改正,濾波半徑分別采用800、1000、1200 km。考慮松遼流域面積限制及信號削弱等因素的影響,本文未添加濾波半徑更大的方案參與比較。濾波結(jié)果統(tǒng)計見表2。

表2 時變模型濾波結(jié)果的數(shù)據(jù)統(tǒng)計

由表2可以看出,IGG-SWARM所得結(jié)果的均方根誤差(RMSE)與相關(guān)系數(shù)均優(yōu)于其他SWARM的濾波結(jié)果;且采用1200 km高斯濾波處理時,統(tǒng)計得到的RMSE與相關(guān)系數(shù)相對于其他濾波結(jié)果更佳。在信噪比方面,IGG-SWARM經(jīng)1200 km高斯濾波的信噪比數(shù)據(jù)相對于ASU和COST-G機構(gòu)的結(jié)果分別高62.47%、55.99%。因此,后續(xù)相關(guān)計算選擇IGG-SWARM時變模型(簡稱SWARM),濾波方案選用1200 km高斯濾波。

3 松遼流域陸地水儲量時空變化分析

分別選擇10、40階次的SWARM時變模型,以及10、40、60階次的GRACE時變模型,反演松遼流域陸地水儲量變化,并進行相關(guān)分析與對比。其中,10階模型未進行任何濾波處理;40階模型均采用1200 km高斯濾波與泄露誤差改正;而60階GRACE時變模型則采用300 km高斯濾波、Duan去相關(guān)濾波與泄露誤差改正。

圖2為松遼流域2015年7月—2020年12月的陸地水儲量及降水量的空間變化趨勢?？梢钥闯?①10階SWARM與GRACE反演松遼流域的陸地水儲量均在西南區(qū)域變化不大,在其他大部區(qū)域上升;②40階SWARM與GRACE均反映出流域陸地水儲量的上升中心位于中東部,且與降水量的空間變化趨勢符合度較好;③60階GRACE反演結(jié)果同樣顯示流域中東部呈顯著上升趨勢,相較于40階SWARM結(jié)果,60階GRACE結(jié)果的空間分辨率有所提高,能夠檢測到流域西南區(qū)域的陸地水儲量存在虧損。此前研究表明,2002—2018年,松花江流域的陸地水儲量自北向南逐漸減少,而遼河流域的陸地水儲量呈南少北多的形勢[17],與本文結(jié)果相吻合。綜合而言,SWARM和GRACE在同一尺度下所得松遼流域陸地水儲量的空間變化趨勢較一致,而模型在不同尺度下反演結(jié)果的空間分辨率存在一定的差異,但基本可以反映出流域陸地水儲量的空間變化情況。

圖2 松遼流域2015年7月—2020年12月陸地水儲量與降水量的空間變化趨勢

SWARM與GRACE反演松遼流域陸地水儲量變化的時間序列與降水量如圖3所示(虛線表示時變模型反演陸地水儲量的變化趨勢)。表3為變化趨勢與周年振幅?？梢钥闯?當(dāng)SWARM與GRACE同為10階時,兩者的趨勢項分別為0.59、0.48 cm/a,周年振幅差值為0.04 cm;當(dāng)SWARM與GRACE同為40階時,其趨勢項分別為0.96、0.71 cm/a,周年振幅差值為0.65 cm;60階GRACE的趨勢項為0.75 cm/a,綜合說明松遼流域的陸地水儲量整體呈上升趨勢。結(jié)合RMSE與相關(guān)系數(shù)計算得知,圖3中3組模型結(jié)果在數(shù)據(jù)重疊時段的RMSE分別達(dá)2.443、2.936、3.923 cm,相關(guān)系數(shù)依次為0.704、0.813、0.621,說明SWARM與GRACE反演結(jié)果具有較好的吻合度與相關(guān)性,其中同尺度模型反演結(jié)果的RMSE與相關(guān)系數(shù)較優(yōu)。此外,圖3反映出模型反演的陸地水儲量變化嚴(yán)重滯后于實際降水變化,通過相位差計算估計滯后約4～6個月。相關(guān)研究表明松花江流域在2002—2017年的GRACE反演結(jié)果與地下水實測數(shù)據(jù)具有長達(dá)8個月的滯后,這種滯后現(xiàn)象可能與流域氣候特征、流域面積和地下水超采等因素有關(guān)[5]。

圖3 時變模型反演松遼流域在2015年7月—2020年12月陸地水儲量變化的時間序列與降水量

表3 松遼流域陸地水儲量變化趨勢與周年振幅

由圖3可以看出,在SWARM與GRACE數(shù)據(jù)重疊時段,10、40階的SWARM與GRACE反演松遼流域的陸地水儲量于2016和2019年普遍達(dá)到較高水平。水文資料表明,2016年松遼流域雨季連續(xù)4月降水量大于9.0 cm,東遼河出現(xiàn)罕見洪水,吉林圖們江、黑龍江等干流全線超警[18],導(dǎo)致陸地水儲量顯著上升;2019年5—8月,松遼流域的降水量均值相對于2016—2020年同期高19.96%,引發(fā)了兩次編號洪水,共67條河流發(fā)生超警洪水,是范圍最廣的洪澇事件之一[18]。而在圖3(c)中,60階GRACE顯示2016年陸地水儲量達(dá)到2015年7月—2017年6月間的最高值,且2019年流域陸地水儲量水平高于2018年。整體而言,SWARM時變模型可有效識別松遼流域的顯著洪澇事件,其陸地水儲量的趨勢變化與GRACE反演的趨勢變化具有較好的一致性。

而在GRACE與GRACE-FO衛(wèi)星的銜接空白期,由于滯后現(xiàn)象,陸地水儲量變化與實際水文變化存在較大時間差。結(jié)果顯示,在10、40階尺度條件下,SWARM數(shù)據(jù)與GRACE填補數(shù)據(jù)的相關(guān)系數(shù)分別為0.565、0.733,具有較高相關(guān)性,并均于2017和2018年出現(xiàn)兩處明顯谷值;而在60階尺度條件下,GRACE填補數(shù)據(jù)的趨勢變化不明顯,這是由于當(dāng)GRACE大于40階時,其信號變異性的比例增大[16]。結(jié)合水文資料記載,松遼流域在2017年初連續(xù)4個月降水量低于1.0 cm,加上氣溫逐漸回暖使蒸散量增大,引起東北大部地區(qū)春夏連旱[18],該年地表水、地下水及水資源總量相對于2016—2020年的年均值分別低30.03%、15.75%和30.31%[19];而在2018年,松遼流域1—2月的降水量均值相對于2016—2020年同期降低10.48%,導(dǎo)致東北地區(qū)再次出現(xiàn)春旱[18]。綜上所述,在GRACE與GRACE-FO衛(wèi)星的間斷期,SWARM數(shù)據(jù)與GRACE空缺填補數(shù)據(jù)在10、40階尺度條件下的相關(guān)性較好,并與降水及水資源公報等水文資料記錄的水文變化一致。

4 結(jié) 論

本文選擇ASU、IGG和COST-G等機構(gòu)的SWARM時變模型反演松遼流域陸地水儲量變化,并與GRACE、GRACE-FO時變模型進行比較,分析了IGG-SWARM時變模型用于松遼流域陸地水儲量變化監(jiān)測的適用性,結(jié)論如下。

(1)ASU、IGG和COST-G等機構(gòu)的SWARM與GRACE時變模型的前10階位系數(shù)精度接近,其中IGG-SWARM時變模型更適用于松遼流域的反演研究,該模型采用1200 km高斯濾波后的信噪比相比ASU與COST-G模型結(jié)果分別高62.47%、55.99%。

(2)IGG-SWARM與GRACE/GRACE-FO時變模型在同尺度下均反映出在2015年7月—2020年12月,松遼流域西南區(qū)域的陸地水儲量變化較小,其他大部區(qū)域呈上升趨勢;而時變模型在不同尺度下的反演結(jié)果對流域陸地水儲量變化的探測能力存在一定的差異,但基本可以反映出流域陸地水儲量的空間變化情況。

(3)IGG-SWARM時變模型的時序結(jié)果可識別流域內(nèi)的顯著旱澇事件,與GRACE、GRACE-FO時變模型反演結(jié)果具有較好的一致性,兩者在數(shù)據(jù)重疊時段的相關(guān)系數(shù)可達(dá)0.6以上。綜合說明SWARM時變模型適用于監(jiān)測松遼流域的陸地水儲量變化研究。