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

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

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

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

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

1 數(shù)據(jù)來(lái)源

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

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

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

2.1 時(shí)變模型階方差計(jì)算

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

(1)

2.2 時(shí)變模型反演陸地水儲(chǔ)量變化

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

(2)

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

(3)

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

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

2.3 確定SWARM時(shí)變模型與濾波方案

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

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

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

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

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

3 松遼流域陸地水儲(chǔ)量時(shí)空變化分析

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

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

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

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

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

表3 松遼流域陸地水儲(chǔ)量變化趨勢(shì)與周年振幅

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

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

4 結(jié) 論

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

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

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

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