石通， 劉新*， 穆大鵬， 李成名， 郭金運*， 邢云鵬

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

0 引言

GRACE重力衛(wèi)星計劃由美國宇航局(NASA)和德國空間飛行中心(DLR)聯(lián)合開發(fā)，于2002年3月成功發(fā)射，可以探測到地球重力場的時變特征，時間分辨率為10天至30天(Tapley et al.，2004).GRACE時變重力數(shù)據(jù)被廣泛應(yīng)用于地球表面質(zhì)量變化的研究，為陸地水儲量變化的反演提供了一種新的手段.目前，國內(nèi)外已經(jīng)有很多學(xué)者將GRACE重力衛(wèi)星應(yīng)用到如流域、區(qū)域陸地水儲量觀測(Chen et al.，2014, 2017；Guo et al.，2014，2016； Liu et al.，2020a)、冰川和兩極質(zhì)量的動態(tài)監(jiān)測(Luthcke et al.，2013；Velicogna and Wahr，2013；高春春等，2019)、洪澇干旱災(zāi)害的監(jiān)測(Long et al.，2014)、全球海平面變化(Jacob et al.，2012)等領(lǐng)域.這些研究表明，GRACE重力衛(wèi)星的觀測結(jié)果可以作為傳統(tǒng)地表水文觀測技術(shù)及重力觀測的有效補充，在近年來的水文學(xué)和地學(xué)研究中得到越來越廣泛的應(yīng)用.

經(jīng)過長達15年之久的運行，GRACE衛(wèi)星于2017年6月退役，GRACE Follow-On(GRACE-FO)重力衛(wèi)星于2018年5月22日發(fā)射，兩顆重力衛(wèi)星之間存在接近一年的間斷期.Velicogna等(2020)分析了格陵蘭島和南極洲的GRACE與GRACE-FO任務(wù)之間的數(shù)據(jù)連續(xù)性，發(fā)現(xiàn)在大陸和區(qū)域尺度上二者數(shù)據(jù)有較好的一致性.雖然現(xiàn)有的GRACE數(shù)據(jù)已經(jīng)能夠?qū)崿F(xiàn)監(jiān)測陸地水儲量變化，但是為了研究陸地水儲量長期變化，需要更長時間的數(shù)據(jù)積累.已有學(xué)者提出了對GRACE數(shù)據(jù)的重建方法以及間斷期的連接方法，如自回歸模型(Forootan et al.，2014)、線性回歸模型(Nie et al.，2016)、統(tǒng)計公式模型(Humphrey and Gudmundsson，2019)、奇異譜分析法(Li et al.，2019)、深度機器學(xué)習(xí)法(Sun et al.，2020；Li et al.，2021).

當前黃河流域生態(tài)環(huán)境脆弱，水資源保障形勢嚴峻，全球變暖背景下的氣候變化和人類活動使黃河流域的水循環(huán)發(fā)生顯著變化，對黃河流域水資源造成影響.已有學(xué)者針對黃河流域研究水儲量變化(尼勝楠等，2014；Jing et al.，2019；Liu et al.，2020b；張璐等，2020)，但其研究時段均較短，且未使用GRACE-FO衛(wèi)星觀測數(shù)據(jù).

本文基于多層感知機(Multi-Layer Perceptron, MLP)神經(jīng)網(wǎng)絡(luò)模型，結(jié)合GLDAS、CPC水文模型、降水、蒸散量以及土壤溫度數(shù)據(jù)，以黃河流域為研究區(qū)，重建GRACE與GRACE-FO之間缺失等效水高，重點分析黃河流域上中下游陸地水儲量變化與降水、蒸散發(fā)及土壤溫度的關(guān)系，并使用連續(xù)小波分析和交叉小波分析研究GRACE反演的陸地水儲量與GLDAS、CPC結(jié)果的周期特征.

1 研究區(qū)域與數(shù)據(jù)

1.1 研究區(qū)概況

黃河流域位于東經(jīng)96°—119°、北緯32°—42°之間，流域面積79.5萬km2.河口鎮(zhèn)以上為黃河上游，河道長3472 km，流域面積42.8萬km2；河口鎮(zhèn)至桃花峪為中游，河道長1206 km，流域面積34.4萬km2；桃花峪以下為下游，河道長786 km，流域面積只有2.3萬km2.黃河流域年降水量在100～1000 mm之間，中上游南部和下游地區(qū)降水大于650 mm，而在內(nèi)陸的西北寧夏、內(nèi)蒙古部分地區(qū)，年降水量不足150 mm，在秦嶺山脈北坡，由于受地形影響較大，降水多在700～1000 mm.流域整體春冬少雨干旱，秋夏雨量充足.圖1為黃河流域上中下游位置示意圖.

圖1 黃河流域位置Fig.1 Location of Yellow River basin

1.2 數(shù)據(jù)

1.2.1 時變重力場數(shù)據(jù)

GRACE及GRACE-FO月時變重力場模型的處理與發(fā)布主要由美國德克薩斯大學(xué)空間研究中心(CSR)、德國地學(xué)研究中心(GFZ)和美國噴氣動力實驗室(JPL)三個研究中心提供.本文采用CSR、JPL、GFZ提供的2002-04—2020-12共192個月(缺失2017-07—2018-05及部分月份)的GRACE level-2 RL06 GSM 60階月時變重力場模型，分別記為CSR-SH、JPL-SH、GFZ-SH.對于GRACE及GRACE-FO在其單獨任務(wù)期間的1個月或2個月缺失通過相鄰月份的平均值進行線性插值彌補(Long et al.，2015).

本文還采用了CSR、JPL發(fā)布的Mascon格網(wǎng)數(shù)據(jù)(Watkins et al.，2015；Save et al.，2016)，以下分別記為CSR-M與JPL-M.CSR-M與JPL-M格網(wǎng)大小分別為0.25°×0.25°與0.5°×0.5°.

1.2.2 水文模型

GLDAS水文模型與由GRACE反演的水儲量之間有很強的相關(guān)性(Long et al.，2013)，本文選取GLDAS-Noah水文模型(Rodell et al.，2004)，數(shù)據(jù)的時間跨度為2002-04—2020-12共225個月，時間分辨率為1月，格網(wǎng)大小為0.25°×0.25°.該水文模型是由美國宇航局戈達德空間飛行中心(GSFC)和美國國家環(huán)境預(yù)報中心(NCEP)聯(lián)合建立的全球水文模型，反映了陸地表面土壤水、冰和雪的變化.在GLDAS水文模型中提取0～2 m土壤濕度、積雪融水以及植被冠層水作為陸地水儲量.另外在GLDAS水文模型中提取土壤溫度和蒸散量，這兩個變量與陸地水儲量變化有一定的聯(lián)系(Li et al.，2021).

CPC水文模型(Fan and van den Dool，2004)全球土壤濕度也用于估計與氣候相關(guān)的水儲量變化(李圳等，2018； Zhong et al.，2019)，CPC水文模型由美國國家海洋和大氣管理局(NOAA)發(fā)布.選用該數(shù)據(jù)的時間跨度與GLDAS相同，格網(wǎng)大小為0.5°×0.5°，從CPC水文模型中提取0～1.6 m土壤濕度，并且包括了積雪融水等引起水平衡變化的因素.

1.2.3 降水數(shù)據(jù)

本文采用中國地面降水月值0.5°×0.5°格點數(shù)據(jù)集(V2.0)(趙煜飛等，2014)，該數(shù)據(jù)由中國氣象網(wǎng)提供的2472個地面雨量站點進行空間插值得到.臺站在中國東部地區(qū)密集且分布均勻，西北和東北地區(qū)相對稀疏.

數(shù)據(jù)來源信息如表1所示.

表1 數(shù)據(jù)來源信息Table 1 Data source information

2 方法

2.1 GRACE時變重力場數(shù)據(jù)處理

在GRACE數(shù)據(jù)后處理過程中，由于GRACE衛(wèi)星的參考框架原點為地球質(zhì)心，其時變重力場模型不包含一階項系數(shù)，為此需要使用SLR觀測的一階項來替換，從而顧及地心運動的影響(Swenson et al.，2008; Sun et al.，2016；孫玉等，2019)，同時GRACE對時變重力場的C20項也不敏感，此外從2016年10月開始，GRACE衛(wèi)星上一個加速度計出現(xiàn)故障，造成C30項精度較低(Loomis et al.，2020)，因此使用JPL及NASA GSFC發(fā)布的TN-13和TN-14數(shù)據(jù)分別替換一階項、C20與C30項.為了消除平均地球重力場的影響，將月時變重力場球諧系數(shù)減去其2004—2009年平均值，進行與Mascon相同的去平均方法，以研究時變重力場的特點.盡管RL06月重力場模型(Dahle et al.，2019)條帶誤差及高頻噪聲減小了很多，但剩余的誤差依然不可忽視.因此本文采用DDK3去相關(guān)濾波(Horvath et al.，2018)的方法對球諧系數(shù)進行處理.等效水高(Equivalent Water Height，EWH)(Wahr et al.，1998)為：

為了與GRACE數(shù)據(jù)進行對比，需要把GLDAS、CPC水文模型計算的陸地水儲量進行球諧展開到60階次，并進行300 km高斯濾波處理.

2.2 多層感知機重建等效水高

MLP是一種典型的前饋人工神經(jīng)網(wǎng)絡(luò)，適合對大量數(shù)據(jù)進行分類、建立復(fù)雜的非線性映射等問題.在設(shè)法獲得有限的輸入變量與相應(yīng)的輸出變量后，即可利用MLP網(wǎng)絡(luò)來模擬輸入與輸出之間的聯(lián)系，從而建立基于 MLP模型的預(yù)測方法(Pal and Mitra，1992；Kuremoto et al.，2014).由于本研究的重建方案是將若干輸入變量映射到輸出變量的回歸問題，所以使用MLP神經(jīng)網(wǎng)絡(luò)模型.

如圖2所示，MLP模型是由多層密集連接的信息處理節(jié)點所形成的網(wǎng)絡(luò)(Zhu et al.，2021)，這些信息處理節(jié)點稱為神經(jīng)元.MLP網(wǎng)絡(luò)由一個輸入層、一個或多個隱含層和一個輸出層組成.輸入層是用于輸入數(shù)據(jù)的節(jié)點，其數(shù)量與數(shù)據(jù)變量的個數(shù)一致，隱含層和輸出層分別由神經(jīng)元和輸出變量組成.相鄰層之間為全連接，而同一層的節(jié)點間未連接.本研究中使用三個隱含層，分別有512、256、128個神經(jīng)元.將時間、經(jīng)緯度、GLDAS水儲量、土壤溫度、蒸散量、降水量以及CPC水儲量當作神經(jīng)網(wǎng)絡(luò)輸入?yún)?shù)，輸出參數(shù)為GRACE等效水高.

MLP神經(jīng)網(wǎng)絡(luò)信息處理機制(Zhu et al.，2021)為：

y=f(Wx+b),(2)

其中，x是神經(jīng)元的輸入，W和b分別代表連接權(quán)值和偏置，f(·)為神經(jīng)元的激活函數(shù)，選擇tanh函數(shù)，以捕捉輸入信息的非線性信息.由于在本研究中MLP網(wǎng)絡(luò)被應(yīng)用到一個回歸問題并且輸出是任意的，因此不應(yīng)在最后一個隱含層中使用激活函數(shù)，所以輸出層的輸出機制為：

y=Wx+b.

(3)

通過神經(jīng)元對信息的處理及層層傳遞，MLP模型完成了從輸入到輸出的映射.MLP模型學(xué)習(xí)的本質(zhì)就是通過網(wǎng)絡(luò)的信息前向傳遞及誤差逆向傳播實現(xiàn)梯度下降算法，不斷調(diào)整網(wǎng)絡(luò)權(quán)值矩陣W和偏置矩陣b，使得網(wǎng)絡(luò)理想輸出與實際輸出的誤差達到最小.

在向MLP神經(jīng)網(wǎng)絡(luò)模型輸入?yún)?shù)之前，需要把數(shù)據(jù)進行標準化，即：

圖2 MLP網(wǎng)絡(luò)結(jié)構(gòu)，輸入層參數(shù)分別為時間、經(jīng)度、緯度、GLDAS與CPC水儲量、降水、蒸散發(fā)(ET)及土壤溫度(SoilTMP)，輸出參數(shù)為GRACE等效水高Fig.2 Structure of MLP network. The input layer parameters are time, longitude, latitude, GLDAS TWS (terrestrial water storage), CPC TWS, precipitation, evapotranspiration (ET) and soil temperature (SoilTMP),the output parameter is GRACE EWH

(4)

在月時間尺度上，由式(1)計算的等效水高與MLP網(wǎng)絡(luò)得到的等效水高重建精度評價選用皮爾遜相關(guān)系數(shù)r(式(5))、納什系數(shù)Ens(式(6))、平均絕對誤差MAE(式(7))與均方根誤差RMSE(式(8))：

(5)

(6)

(7)

(8)

相關(guān)系數(shù)r的取值范圍為|r|≤1，|r|越接近于1，則表明等效水高與其重建值相關(guān)程度越高.納什系數(shù)Ens反映擬合程度，Ens取值范圍為負無窮至1，Ens越接近1，表示重建質(zhì)量好，模型可信度高.MAE與RMSE越小，重建精度越高.

3 結(jié)果與分析

3.1 MLP神經(jīng)網(wǎng)絡(luò)重建結(jié)果驗證

對于已有GRACE及GRACE-FO的時變重力場數(shù)據(jù)，從2004—2015年逐次空缺一年進行實驗，以驗證等效水高重建精度.本文使用MLP神經(jīng)網(wǎng)絡(luò)的目的是重建GRACE與GRACE-FO兩代重力衛(wèi)星之間缺失的約一年數(shù)據(jù)，所以選擇缺失一年作為驗證，若扣除過多的年份則MLP神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)的樣本數(shù)據(jù)減少，且與真實缺失時間差異較大.基于MLP神經(jīng)網(wǎng)絡(luò)模型訓(xùn)練時間為GRACE運行時間(2002-04—2017-06)和GRACE-FO運行時間(2018-06—2020-12)，在訓(xùn)練時間中從2004—2015年逐年地扣除一年，共構(gòu)成12次實驗，模型具體訓(xùn)練時間如表2所示.

表2 MLP網(wǎng)絡(luò)訓(xùn)練時間Table 2 Training time of MLP network

模型訓(xùn)練后，再輸入空缺年份的土壤水含量、降水、蒸散發(fā)及土壤溫度等數(shù)據(jù)，用MLP模型重建出此缺失時段的GRACE EWH.利用訓(xùn)練樣本對MLP網(wǎng)絡(luò)進行訓(xùn)練，網(wǎng)絡(luò)的參數(shù)根據(jù)誤差梯度不斷修正，輸出結(jié)果逐步逼近訓(xùn)練集的輸出結(jié)果.圖3給出了黃河流域EWH與對應(yīng)重建值，圖3中數(shù)字對應(yīng)實驗序號.由圖3可以看出在逐年空缺的實驗中，等效水高重建值與真實值相符程度較高，相關(guān)系數(shù)r均大于0.93(p<0.01).但重建值在其拐點處與真實值有著一定的差距，而在拐點之間重建效果較好，其原因為極值點之間的數(shù)據(jù)多，MLP網(wǎng)絡(luò)訓(xùn)練程度高，而極值處的等效水高數(shù)據(jù)少，MLP網(wǎng)絡(luò)對極值不敏感，導(dǎo)致拐點處的重建效果差一些.

圖3 黃河流域等效水高重建Fig.3 Reconstructed EWH in the Yellow River basin

圖4給出了12次實驗重建等效水高與真實值之間各檢驗參數(shù)的平均值在空間上的分布.各格網(wǎng)點等效水高真實值與重建值的相關(guān)系數(shù)r均大于0.5(p<0.01)，且有97%以上的相關(guān)系數(shù)r大于0.8(p<0.01)，表明MLP網(wǎng)絡(luò)重建值與數(shù)據(jù)真實值之間的相關(guān)性較強.有82%以上Ens值大于0.5，表明MLP模型重建值的擬合程度高.至于(100°E，38°N)附近擬合程度較差，其原因為邊界效應(yīng)的影響，MLP模型訓(xùn)練該區(qū)域附近的格網(wǎng)點較少，導(dǎo)致其重建精度較差.MAE與RMSE最大值分別為1.52 cm和2.04 cm，誤差較大處均位于流域邊界處，原因也由于MLP網(wǎng)絡(luò)沒有訓(xùn)練黃河流域邊界之外的格網(wǎng)點數(shù)據(jù)，導(dǎo)致臨近邊界處誤差大，而流域內(nèi)部的MAE與RMSE均較小，滿足重建精度要求.以上實驗證明了MLP網(wǎng)絡(luò)對于重建GRACE與GRACE-FO兩代重力衛(wèi)星之間缺失等效水高數(shù)據(jù)的可行性.

圖4 黃河流域等效水高重建值與真實值的誤差檢驗Fig.4 Accuracy of the reconstructed and true EWH in the Yellow River basin

基于上述實驗，保持輸入?yún)?shù)不變，模型訓(xùn)練時間為2002-04—2017-06和2018-06—2020-12，將EWHMean-SH、EWHCSR-SH、EWHJPL-SH、EWHGFZ-SH、EWHCSR-M、EWHJPL-M逐一作為輸出參數(shù)進行多次實驗.模型訓(xùn)練后，輸入2017-07—2018-05的各輸入?yún)?shù)，MLP網(wǎng)絡(luò)將重建出該時間段的黃河流域各格網(wǎng)點等效水高值.圖5a給出了2002—2020年黃河流域等效水高，圖5b為MLP網(wǎng)絡(luò)重建的黃河流域在2017—2018年等效水高值.由于不同機構(gòu)對于GRACE觀測數(shù)據(jù)的處理方式不盡相同，導(dǎo)致計算的EWH值也有一定差異，造成MLP網(wǎng)絡(luò)建立的模型差異，但也可在預(yù)測期間體現(xiàn)出等效水高大致變化情況.EWHMean-SH與EWHCSR-M的RMSE為0.88 cm，GRACE 球諧系數(shù)反演結(jié)果與Mascon數(shù)據(jù)吻合程度較好，說明使用Mean-SH計算的EWH可信度較高，且避免了某些月份的異常值，所以在下文討論黃河流域水儲量變化時均使用EWHMean-SH.

3.2 GRACE與水文模型時間序列分析

GRACE反演的黃河流域EWHMean-SH與GLDAS、CPC水文模型得到的水儲量變化時間序列見圖6.圖6水儲量時間序列結(jié)果均表現(xiàn)出明顯的季節(jié)性，分別在秋季和春季達到最大值和最小值.對于圖6a黃河流域整體來說，GRACE數(shù)據(jù)反演的EWH在2015年之前與水文模型結(jié)果符合程度較高，GRACE與GLDAS相關(guān)系數(shù)r為0.69(p<0.01)，而在2015年以后二者相關(guān)系數(shù)r為0.37(p<0.01).在2015年之后，GLDAS得到的水儲量變化均高于GRACE數(shù)據(jù)反演的EWH，且GRACE反演的水儲量呈現(xiàn)明顯的下降趨勢，而水文模型得到的水儲量歷年結(jié)果大致持平甚至略微上升，具體原因可能為GRACE重力衛(wèi)星與水文模型監(jiān)測的具體內(nèi)容有所不同.GRACE反演結(jié)果包含多種因素，而GLDAS數(shù)據(jù)僅包含土壤水、積雪融水及植物冠層水含量，不包含深層水的變化，GLDAS、CPC水文模型僅能探測到淺層地表水，而黃河流域用水量大，對地下水超采嚴重，導(dǎo)致深層水儲量逐年減少，這部分能夠被GRACE重力衛(wèi)星探測到(馮偉等，2012).基于2004—2009年距平值，黃河流域EWH變化范圍為-10.8～7.9 cm.

由圖6b、c可知黃河流域上游與中游地區(qū)水儲量情況相似，在研究期間大致呈現(xiàn)持平趨勢，且上游與中游地區(qū)GRACE反演的EWH與水文模型相關(guān)性較好，上游與中游的GRACE與GLDAS相關(guān)系數(shù)r分別達到0.6與0.38(p<0.01)，進一步驗證了GRACE反演結(jié)果的可靠性.而黃河流域下游地區(qū)由GRACE反演的EWH呈現(xiàn)明顯的下降趨勢(見圖6d)，且年內(nèi)變化規(guī)律性較差，波動較多，GRACE與水文模型相關(guān)性均較小，而GLDAS與CPC相關(guān)系數(shù)r達到0.8(p<0.01)，說明下游地區(qū)GRACE探測到深層地下水虧損嚴重，人為因素影響大.

對于黃河流域整體及上中下游2002—2020年EWHMean-SH及GLDAS、CPC水文模型得到的水儲量時間序列按照(9)式(尼勝楠等，2014)進行最小二乘擬合，以研究水儲量在季節(jié)尺度上速率、振幅及相位的變化情況：

(9)

其中，Ai為振幅，φi為相位，Ti為周期，考慮周年及半年變化項，取a1作為EWH年變化速率.

由表3可見，由GRACE反演的黃河流域水儲量呈現(xiàn)以-0.51±0.03 cm·a-1的速率減少，而GLDAS及CPC水文模型變化速率大致相當，均為略微上升，GRACE反演結(jié)果與水文模型差距較大的原因為黃河流域地表土壤水保持情況較好，而深層地下水逐年減少，尤其以下游地區(qū)最為明顯.GRACE和GLDAS、CPC水文模型的周年相位較為一致，與GLDAS的周年振幅也較接近，而CPC的振幅明顯高于其他兩者，由于不同水文模型融合了不同的數(shù)據(jù)源，其包含的信息與精度存在一定的差異(李圳等，2018)，造成CPC水文模型的年周期振幅較大.

表3 黃河流域陸地水儲量變化速率、振幅及相位Table 3 Rate, amplitude and phase of EWH in the Yellow River basin

為研究黃河流域水儲量年際變化特征，采用13點滑動平均(13 points moving average，13PMA)(李武東等，2017)的方法來體現(xiàn)周年波動，結(jié)果見圖7.由圖7a可知，黃河流域水儲量整體上呈現(xiàn)下降趨勢，2003—2004年水儲量呈現(xiàn)明顯的上升趨勢后，在2004年春季達到最大值，此后直到2017年均大致呈現(xiàn)出不同程度的下降，于2017年后又呈現(xiàn)出上升趨勢.GLDAS與CPC水文模型結(jié)果顯示淺層水儲量大致呈現(xiàn)持平或上升趨勢.由圖7b、c可以看出黃河流域上游與中游水儲量在2002—2015年間變化不大，而2015—2016年間呈現(xiàn)出明顯下降趨勢，又在此后呈劇烈上升.而圖7d表明黃河流域下游地區(qū)在研究期間均呈現(xiàn)不同程度的下降趨勢，僅在2012年及2020年出現(xiàn)明顯上升趨勢.

圖7 采用13點滑動平均方法得到的等效水高年際變化趨勢(a) 黃河流域整體； (b) 黃河流域上游； (c) 黃河流域中游； (d) 黃河流域下游.Fig.7 EWH interannual trends using 13 points moving average method(a) Overall Yellow River basin； (b) The upper reach of Yellow River basin； (c) The middle reach of Yellow River basin；(d) The lower reach of Yellow River basin.

3.3 黃河流域陸地水儲量空間變化趨勢

為研究黃河流域在2002年至2020年陸地水儲量變化情況，利用GRACE球諧系數(shù)及Mascon計算了黃河流域水儲量年際變化趨勢.如圖8所示，黃河流域EWH變化趨勢上中下游差異明顯，其中上游水儲量保護情況較好，呈現(xiàn)略微上升趨勢，由圖8a—d可以看出球諧系數(shù)反演的結(jié)果對上游增長趨勢較為敏感；中游大致呈現(xiàn)持平或略微下降趨勢；而下游華北平原呈現(xiàn)嚴重虧損趨勢，以圖8e、f中Mascon對下游地區(qū)水儲量嚴重虧損的反映最為明顯.JPL-M格網(wǎng)大小雖然為0.5°×0.5°，但其許多相鄰格網(wǎng)的數(shù)值一致，導(dǎo)致圖8f中出現(xiàn)較大面積的相同數(shù)據(jù).

圖8 GRACE反演黃河流域等效水高變化率Fig.8 Rates of EWH from GRACE in the Yellow River basin

4 討論

4.1 黃河流域EWH影響因素分析

在流域水文循環(huán)過程中，陸地水儲量能夠綜合反映循環(huán)過程中降水、蒸散發(fā)等各種水文過程，其中降水能夠?qū)е铝饔蛩畠α吭黾?，而蒸散發(fā)等水文過程會使水儲量減少，土壤溫度的變化直接影響著蒸散發(fā)，進一步影響水文循環(huán).水儲量變化與降水、蒸散發(fā)及徑流的關(guān)系(尼勝楠等，2014)為：

(10)

其中S為水儲量，P為降水，E為蒸散發(fā)，R為徑流.

本節(jié)主要分析降水、蒸散及土壤溫度對黃河流域EWH的影響，并將EWH分別與降水量、蒸散量以及土壤溫度進行對比分析，如圖9所示，圖中降水量、蒸散量以及土壤溫度均已去除了2004—2009年的平均值，即距平值.

圖9 GRACE反演黃河流域EWH與降水、蒸散發(fā)及土壤溫度時間序列對比(a) 黃河流域整體； (b) 黃河流域上游； (c) 黃河流域中游； (d) 黃河流域下游.Fig.9 Time series of EWH by GRACE compared with precipitation, evapotranspiration and soil temperature in Yellow River basin(a) Overall Yellow River basin； (b) The upper reach of Yellow River basin； (c) The middle reach of Yellow River basin；(d) The lower reach of Yellow River basin.

對黃河流域整體及上中下游在研究期間由GRACE反演的水儲量變化與降水、蒸散量以及土壤溫度的相關(guān)系數(shù)以及滯后性進行分析，取最大滯后相關(guān)系數(shù)對應(yīng)的時間作為滯后時間，結(jié)果見表4.

表4 黃河流域EWH與降水、蒸散發(fā)及土壤溫度的滯后時間及相關(guān)系數(shù)Table 4 The lag time and coherence of the EWH with precipitation, evapotranspiration and soil temperature in Yellow River basin

4.1.1 降水對EWH的影響

降水是黃河流域水文循環(huán)過程中主要的水量來源，在水文循環(huán)過程中發(fā)揮著重要作用.黃河流域降水主要集中在汛期(6—10月)，汛期降水量占全年的70%以上，而冬季降水量僅占全年的2%～3.7%.由圖9可知黃河流域水儲量與月降水量呈現(xiàn)相同的季節(jié)性變化趨勢，在降水汛期水儲量顯著上升.水儲量在2003—2004年呈現(xiàn)上升趨勢，上升速率高達4.21±1.3 cm·a-1，這與2003年黃河流域降水量均較常年偏多，且出現(xiàn)了多年不遇的秋汛密切相關(guān)，導(dǎo)致黃河流域的水儲量距平值在2003年11月達到了研究期間的最大值+7.9 cm，并從此后2004年開始水儲量距平值呈現(xiàn)逐年減少趨勢，在2015—2016年間黃河流域EWH出現(xiàn)了明顯下降，速率達-2.52±0.65 cm·a-1，而2015年6—9月黃河流域平均降水量較常年同期偏少24%，水儲量變化總體趨勢與降水在一定程度上線性相關(guān).

由表4可以得出水儲量變化滯后于降水的時間為2～3個月，原因為降水量經(jīng)過一定時間的下滲等水文過程后轉(zhuǎn)化為水儲量從而被GRACE重力衛(wèi)星探測到，也包括了觀測月份之前的降雨.黃河流域上游地區(qū)水儲量與降水的變化趨勢基本一致，滯后2個月的EWH與降水相關(guān)性比中下游地區(qū)更強(r=0.31，p<0.01)，這主要是因為上游青藏高原、內(nèi)蒙古高原水儲量受自然因素影響較大，人為因素影響??；中游地區(qū)水儲量與降水的滯后相關(guān)系數(shù)較低(r=0.22，p<0.01)，原因為中游黃土高原生態(tài)環(huán)境脆弱，水儲量的影響因素較為復(fù)雜；而下游地區(qū)降水較為充沛，年際變化較大，而水儲量與降水基本不相關(guān)(r=0.08，p>0.1)，下游華北平原人口稠密，人類生活、工業(yè)用水量較大，且由《黃河水資源公報》可知部分地區(qū)由于長期過量開采地下水，黃河流域已形成多個淺層地下水降落漏斗，導(dǎo)致降水量對于下游水儲量的影響不占據(jù)主導(dǎo)因素.

4.1.2 蒸散發(fā)、土壤溫度對EWH的影響

由圖9a可以看出，黃河流域蒸散發(fā)量與土壤溫度在2002—2020年間均呈現(xiàn)輕微上升趨勢，而黃河流域年水儲量呈明顯下降趨勢，其主要原因為土壤溫度增高直接導(dǎo)致水分蒸散發(fā)的增多，而過高的蒸散發(fā)導(dǎo)致流域內(nèi)水儲量減少.流域蒸散發(fā)與土壤溫度呈明顯的季節(jié)性趨勢，年內(nèi)最大值出現(xiàn)在7—8月，此時的蒸散發(fā)與9—10月的水儲量相關(guān)性較強(r=-0.31，p<0.05)，此時土壤溫度與9—10月水儲量相關(guān)性也最為顯著(r=-0.4，p<0.05)，由此可見在夏季蒸散發(fā)與土壤溫度均達到年內(nèi)峰值，導(dǎo)致部分水儲量通過水循環(huán)回到大氣之中，而此水文過程大致需要2個月的時間.

由表4 還可知，對于黃河流域2002—2020年蒸散發(fā)、土壤溫度時間序列與滯后2～3個月的水儲量時間序列之間相關(guān)性最強，且蒸散發(fā)、土壤溫度對于上游地區(qū)的水儲量相關(guān)性影響較大，對中下游及整個黃河流域影響較小，這也和降水與水儲量的相關(guān)性較為相似，而對于黃河流域整體滯后相關(guān)系數(shù)為正的原因為其相關(guān)性均不顯著，且水儲量受多種因素影響.GLDAS水文模型未能夠充分考慮人類活動對蒸散發(fā)過程的影響，且由于其缺乏地下水模塊而未能考慮地下水等因素對流域蒸散發(fā)的影響作用，而GRACE重力衛(wèi)星則能夠充分探測到流域內(nèi)地下水和冰川等所有形式水資源的變化作用.總體而言，蒸散發(fā)與土壤溫度對于水儲量的影響是小于降水的，這與張璐等(2020)的結(jié)果一致.

4.2 EWH時間序列小波分析

相對于傅里葉變換，小波變換(Grinsted et al.，2004)在時域、頻域中都能表現(xiàn)出優(yōu)良的特性，利用該方法可以很好地將時間序列的周期性變化展現(xiàn)出來.本文采用連續(xù)小波分析(continuous wavelet transform，CWT)對2002—2020年共225個月GRACE反演的等效水高與GLDAS、CPC水儲量時間序列進行功率譜分析，對其周期性進行探測，圖10為連續(xù)小波分析的功率譜.從圖10可以看出，基于GRACE、GLDAS與CPC數(shù)據(jù)解算的水儲量都表現(xiàn)出了一定的周期變化特征，具體表現(xiàn)為均具有明顯的周年周期信號，CPC水文模型計算的水儲量周年信號最為明顯，而2014—2017年水儲量呈劇烈下降趨勢，降水偏少，導(dǎo)致此期間GRACE與GLDAS計算的水儲量周年信號不明顯.

圖10 GRACE、GLDAS及CPC計算的等效水高的小波分析Fig.10 The continuous wavelets transform of EWH by GRACE, GLDAS and CPC

圖11 GRACE分別與GLDAS、CPC的交叉小波分析與相干小波譜分析Fig.11 The cross wavelets transform and wavelet coherence of EWH by GRACE with GLDAS and CPC

圖12 GRACE與降水、蒸散發(fā)以及土壤溫度的交叉小波分析Fig.12 The cross wavelets transform of GRACE with precipitation, evapotranspiration and soil temperature

為研究GRACE反演的等效水高與GLDAS、CPC水文模型得到的水儲量之間的相關(guān)性和共振周期，分別對其時間序列進行了交叉小波(Cross Wavelet Transform，XWT)和小波相干譜分析(Wavelet Coherence，WTC)，結(jié)果如圖11所示，圖中的箭頭方向表示二者相位的相關(guān)性，向右表示同相位，向左表示反相位，向下表示相位超前90°變化，向上表示相位超前270°或者延后90°變化(郭金運等，2015).交叉小波譜的大能量區(qū)基本都出現(xiàn)在整個時間序列的周年位置上(見圖11a、b)，且相比于GRACE與GLDAS水儲量之間共振周期，GRACE與CPC反演的水儲量之間周年共振周期更為突出，說明基于GRACE反演的水儲量與CPC水文模型結(jié)果在周年周期上相關(guān)性程度較高；由圖11c、d相干小波譜分析可以看出，大部分的高值能量區(qū)都在影響錐線內(nèi)，說明GRACE反演的等效水高與GLDAS、CPC水文模型得到的水儲量較為一致，具有很強的相關(guān)性，且圖中箭頭大部分向右，說明兩種數(shù)據(jù)計算的等效水高具有大致相同的相位.

對GRACE反演的黃河流域水儲量分別與降水、蒸散量以及土壤溫度進行交叉小波分析(圖12)，以進一步分析其相關(guān)性.由圖12交叉小波功率譜可以看出，GRACE反演的水儲量與降水、蒸散量以及土壤溫度之間的周年周期最為突出；在小于1周年的周期中，GRACE與降水、蒸散量之間會間斷地出現(xiàn)一些高值能量帶，例如降水、蒸散發(fā)在半周年的周期位置上高能量區(qū)出現(xiàn)在2012—2014年，說明在半周年周期位置上GRACE與二者也具有一定的相關(guān)性，但不如周年周期明顯.周年周期的箭頭均朝上，表明GRACE 反演的水儲量相較于降水、蒸散量以及土壤溫度變化滯后大致3個月，此結(jié)果與表3的結(jié)果相吻合.

5 結(jié)論

本文基于多層感知機重建了GRACE與GRACE-FO之間缺失的等效水高，研究黃河流域在2002—2020年水儲量變化情況，并分析了降水、蒸散發(fā)、土壤溫度與水儲量的相關(guān)性，通過滑動平均方法與小波分析，討論了黃河流域上中下游等效水高的時變規(guī)律.分析結(jié)果表明：

(1)基于多層感知機對GRACE與GRACE-FO之間缺失等效水高的重建有一定的可行性，MLP神經(jīng)網(wǎng)絡(luò)通過訓(xùn)練相關(guān)的輸入?yún)?shù)后，可重建出精度較高的缺失等效水高，重建相關(guān)系數(shù)r均大于0.93(p<0.01).

(2)黃河流域水儲量變化呈現(xiàn)出明顯的季節(jié)性，水儲量分別在秋季與春季達到最大值與最小值，整體呈現(xiàn)以-0.51±0.03 cm·a-1速率下降的趨勢，2003—2004年水儲量呈現(xiàn)明顯的上升趨勢后，在2004年春季達到最大值，此后直到2017年均大致呈現(xiàn)出不同程度的下降， 2017年后又呈現(xiàn)出上升趨勢.黃河流域上游水儲量以0.01±0.02 cm·a-1的速率大致持平，中游以-0.11±0.02 cm·a-1的速率輕微下降，上游與中游地區(qū)水儲量與GLDAS水文模型相關(guān)系數(shù)r分別達到0.6與0.38(p<0.01)，而下游地區(qū)以-0.55±0.02 cm·a-1的速率顯著下降.

(3)流域水儲量變化受降水、蒸散發(fā)、溫度及人類活動等因素的共同影響，GRACE反演黃河流域水儲量變化相對于降水、蒸散發(fā)及土壤溫度滯后時間達2～3個月，降水與水儲量變化呈現(xiàn)出一致性，而水儲量變化與蒸散發(fā)、土壤溫度呈負相關(guān)，且降水對于水儲量的影響大于蒸散發(fā)和土壤溫度.上游與中游水儲量受降水、蒸散發(fā)及土壤溫度影響較大，而下游地區(qū)水儲量受氣象因素影響小，受人為因素影響較大.

(4)通過小波分析，驗證了GRACE與GLDAS及CPC得到的陸地水儲量具有很強的相關(guān)性，且均具有明顯的周年共振周期；GRACE反演結(jié)果與降水、蒸散量以及土壤溫度的周年共振周期同樣明顯.

致謝感謝CSR、JPL和GFZ提供的GRACE 衛(wèi)星數(shù)據(jù)，GSFC提供的GLDAS水文數(shù)據(jù)，NOAA提供的CPC水文數(shù)據(jù).感謝國家青藏高原科學(xué)數(shù)據(jù)中心(http:∥data.tpdc.ac.cn)提供中國區(qū)域基于降水重構(gòu)陸地水儲量變化數(shù)據(jù).感謝審稿專家對本文提出的寶貴修改意見.