基于GRACE和GRACE-FO的黃河流域陸地水儲量及影響因素分析

任立良，王 宇，江善虎，衛林勇，王孟浩，張怡雅

(1.河海大學水文水資源與水利工程科學國家重點實驗室，江蘇 南京 210098；2.河海大學水文水資源學院，江蘇 南京 210098)

黃河流域是我國重要的生態屏障區和經濟發展區[1]。黃河流域上游水資源利用效率和經濟水平較低；中游水土流失嚴重，生態系統脆弱，水資源污染嚴重；下游用水需求大，水資源浪費嚴重[2-4]。受氣候變化和人類活動雙重影響，流域水循環過程發生了顯著變化，水資源保障形勢愈發嚴峻，研究變化環境下黃河流域水資源演變特征對黃河流域生態保護和高質量發展具有重大意義[5-6]。陸地水儲量(terrestrial water storage, TWS)變化綜合反映區域地表水、地下水以及河流湖泊水的變化，是衡量地區水資源變化的重要指標。由美國國家航空航天局和德國航空航天中心聯合開發的GRACE(gravity recovery and climate experiment)衛星計劃是監測大尺度TWS變化的有效方法。自2002年3月GRACE衛星發射以來，GRACE衛星數據被廣泛應用于研究陸地冰川消融、海平面與環流變化、干旱監測和TWS監測等[7-10]。張璐等[11]基于GRACE衛星遙感數據，發現黃河流域2003—2015年多年平均TWS西多東少，呈下降趨勢，其上游主要受降水和蒸散發的影響，中、下游受人類活動影響較大；李曉英等[12]發現黃河流域2003—2016年TWS變化呈下降趨勢，同一時段徑流對TWS變化的影響最直接，考慮時滯性，降水主要受森林、農田擴大，植被的葉面積指數和蒸散發增大的影響；Li等[13]發現2003—2016年黃河流域中、下游南部TWS下降幅度最大。GRACE衛星遙感數據截止到2017年6月，后續由GRACE Follow-On (GRACE-FO)衛星接替，GRACE-FO衛星遙感數據從2018年5月更新，兩個衛星遙感數據之間存在約1 a的數據缺失，不利于研究TWS長時間序列變化特征。針對這一問題，不同學者采用不同方法填補這一數據缺失，Sun等[14-15]采用深度神經網絡、BP神經網絡和徑向基函數神經網絡來填補GRACE與GRACE-FO間的數據缺失。

本文采用長短期記憶(long short-term memory, LSTM)神經網絡模型，重建GRACE和GRACE-FO間的TWS變化量，研究黃河流域2002年4月至2020年3月長時間序列TWS，分析氣候因子與TWS變化的相關性，為黃河流域水資源管理和利用提供科學依據，并為其他流域的TWS研究提供參考。

1 研究區概況和數據來源

1.1 研究區概況

黃河流域位于95°E～119°E、32°N～41°N之間，流域面積約79.5萬km2，地勢自西向東逐級下降(圖1(a))，流域橫跨青藏高原、內蒙古高原、黃土高原以及華北平原，流經青海、四川等9個省級行政區(下稱省區)。黃河流域2003—2016年平均TWS西部高、東部低，空間差異顯著(圖1(b)，圖中TWSA(terrestial water storage anomaly)為TWS的距平值，即為遙感TWS數據減去2004—2009年的月平均值)。2003—2019年流域降水量分布不均，從東南向西北遞減，降水集中在夏秋季節(圖1(c))；流域蒸發能力很強，東南部蒸散發量大于西北部(圖1(d))，最大蒸散發出現在7—8月，此時對TWS影響最大[16]。

(a) 高程

1.2 數據來源

1.2.1GRACE數據

GRACE數據主要由美國德克薩斯大學空間研究中心、美國太空總署噴氣動力實驗室和德國波茨坦地球科學中心3家科研機構對其進行解算和發布[17]。本文選取了2002年4月至2020年3月GRACE/GRACE-FO CSR RL06 Mascon V2產品TWSA，時間分辨率為1月，空間分辨率為0.25°。GRACE缺失數據是由于①GRACE衛星傳感器等問題導致某些月份數據缺失；②GRACE和GRACE-FO間存在數據缺失，時間跨度為2017年7月至2018年5月。本文采用六點樣條插值法填補GRACE衛星缺失月份數據，并利用LSTM填補GRACE和GRACE-FO間連續11月TWS變化量數據缺失，用于分析2002年4月至2020年3月黃河流域TWS的變化特征。

1.2.2GLDAS數據

全球陸面數據同化系統(global land data assimilation system, GLDAS)由美國航空航天局戈達德空間飛行中心和美國海洋與大氣局國家環境預報中心共同開發完成。GLDAS包括Noah、VIC(variable infiltration capacity)、CLM(community land model)、Mosaic等4個模型，通過汲取遙感衛星觀測數據和地面觀測數據，輸出全球地表狀態通量及變量數據[18]。本文采用最高空間分辨率為0.25°的Noah陸地表面模型，時間尺度為2002年4月至2020年3月，與GRACE/GRACE-FO CSR RL06 Mascon產品數據相匹配，并選取了Noah陸地表面模型中的雪水當量、土壤水分(土壤深度分別為0～10 cm、10～40 cm、40～100 cm、100～200 cm)和冠層水數據，累加得到GLDAS的TWSA。

1.2.3降水量和氣溫數據

本文采用中國氣象局發布的2002年4月至2020年3月黃河流域降水量和氣溫數據，分別來自中國地面降水月值(China gauge-based monthly precipitation analysis product, CPAP)0.5°×0.5°格點數據集V2.0和中國地面氣溫月值0.5°×0.5°格點數據集V2.0，時間分辨率為1月，降尺度為0.25°。

1.2.4蒸散發和徑流數據

本文蒸散發和徑流數據來源于Noah模型2002年4月至2020年3月的模擬結果，空間分辨率為0.25°，蒸散發數據時間分辨率為1 s，徑流數據時間分辨率為3 h，經累加計算將時間分辨率轉換為1月。

2 研究方法

2.1 LSTM神經網絡模型

LSTM模型通過改進隱藏層內部結構獲得長期記憶樣本數據中的歷史重要信息，解決了傳統循環神經網絡隨訓練時間加長或網絡層數增多引起的梯度爆炸或消失問題。LSTM單元中引入了輸入門(控制當前時刻神經單元的輸入)、遺忘門(以一定的概率控制上一時刻以及單元歷史信息的存儲)和輸出門(控制當前時刻記憶單元的信息輸出)，能較為準確地預測時間序列。Wei等[19]采用LSTM模型，結合GLDAS地表水、累積降水、氣溫數據重建了柴達木盆地GRACE和GRACE-FO間的缺失數據，研究了柴達木盆地TWS的時空變化特征及其原因。目前大多學者利用人工神經網絡重建的輸入層均采用降水、氣溫、蒸散發和GLDAS地表水等數據，直接模擬GRACE和GRACE-FO間的缺失數據。本文將GRACE和GRACE-FO衛星數據反演與GLDAS和水量平衡方程計算的TWS變化量作為輸入變量，采用LSTM模型重建黃河流域GRACE和GRACE-FO間的TWS變化量，反推GRACE和GRACE-FO間的缺失數據。

2.2 TWS變化量計算方法

基于GRACE和GRACE-FO衛星數據反演和LSTM模型重建的逐月TWSA序列，計算逐月TWS變化量時間序列[20]：

(1)

式中：CSt為流域t月的TWS變化量；ST,t+1、ST,t-1分別為t+1月和t-1月的TWSA。

基于GLDAS蒸散發及徑流數據，結合CPAP月降水數據，利用水量平衡方程計算黃河流域TWS變化量：

CSt=Pt-Et-Rt

(2)

式中Pt、Et、Rt分別為t月流域平均降水量、蒸散發量和徑流深。理論上，式(1)(2)得出的TWS變化量是相等的，因此，采用由水量平衡方法得到的TWS變化量來評估LSTM模型計算得到的TWS變化量是可行的。

2.3 地下水儲量計算方法

流域TWS包括地表水儲量、土壤水儲量和地下水儲量。為進一步分析流域TWS變化量的子成分，需對黃河流域地下水儲量進行計算，計算公式[21]為

SG=ST-SSM-SSW-SCI

(3)

式中：SG為流域地下水儲量；ST為流域TWS；SSM為流域土壤水儲量；SSW為流域雪水當量；SCI為流域冠層水儲量。

2.4 Spearman秩相關系數法

Spearman秩相關系數法常用于分析兩個變量之間關系的密切程度，計算公式為

(4)

3 結果與討論

3.1 GRACE和GRACE-FO間TWS變化量重建效果驗證

利用GRACE(2002年5月至2017年5月)和GRACE-FO(2018年7月至2020年2月)衛星數據反演的連續201月TWS變化量時間序列測試LSTM模型對TWS變化量的模擬效果，隨機抽取190月數據作為訓練樣本，剩余11月數據作為驗證樣本，結果如圖2所示。訓練期的Spearman秩相關系數為0.93，納什效率系數為0.87，標準化均方根誤差為0.06(圖2(a))；驗證期的Spearman秩相關系數為0.95，納什效率系數為0.83, 標準化均方根誤差為0.12(圖2(b)),表明LSTM模型對TWS變化量模擬性能良好，可用于重建GRACE和GRACE-FO間TWS變化量。

(a) 訓練期

GRACE和GRACE-FO衛星數據反演和LSTM模型重建的黃河流域2002年5月至2020年2月區域平均TWS變化量在月尺度上的變化過程如圖3(a)所示，與GLDAS和水量平衡方程計算結果相比可以發現三者的變化趨勢非常相似，峰值所處的時間段相近。如圖3(b)所示, 衛星反演和模型重建的TWS變化量與水量平衡方程計算結果的Spearman秩相關系數在GRACE時期為0.64、數據缺失時段為0.75、GRACE-FO時期為0.73、研究期為0.65；如圖3(c)所示，衛星反演和模型重建的TWS變化量與GLDAS計算結果的Spearman秩相關系數在GRACE時期為0.73、數據缺失時段為0.86、GRACE-FO時期為0.73、研究期為0.72。Spearman秩相關系數均通過了p<0.01的顯著性檢驗，比較結果驗證了LSTM模型重建GRACE和GRACE-FO間TWS變化量的可靠性。

圖3 TWS變化量在月尺度上的變化趨勢和相關性分析

水量平衡方程計算TWS變化量時只考慮了氣候因素的影響，而GRACE衛星遙感數據受人類活動的影響，兩者之間有一定的差異，且GRACE衛星遙感數據對于TWS虧損更加敏感。而GRACE衛星遙感數據和GLDAS得到的TWS變化量包含的水文變量不一致，GRACE衛星遙感數據包含了地下水和湖泊水等，兩者之間也存在一定差異性。

3.2 TWS趨勢分析

利用TWSA對TWS進行趨勢分析，圖4為黃河流域2002年4月至2020年3月的TWSA在月尺度上變化趨勢。從圖4可以看出，全流域TWSA呈下降趨勢，下降速率為-0.40 mm/月；上游TWSA呈平緩下降趨勢，下降速率為-0.16 mm/月；中游TWSA呈緩慢下降趨勢，下降速率為-0.67 mm/月；下游TWSA呈急劇下降趨勢，下降速率為-1.83 mm/月，表明黃河流域TWS空間差異顯著。由于黃河流域降水、氣溫、蒸散發等氣候因子具有季節性變化特征，TWS也表現出明顯的季節性變化特征，旱季TWS一般處于虧損狀態，濕季TWS一般處于盈余狀態。

圖4 TWSA在月尺度上的變化趨勢

如圖5所示，黃河流域2003—2019年多年月均TWSA、降水量、蒸散發量和氣溫均呈現明顯的季節性變化，黃河流域降雨集中在夏秋季節，多年最小月均TWSA和最大月均TWSA分別出現在6月和10月。6月之前雨季還未來臨，氣溫高，蒸散發量大，TWSA最低；7—8月降雨集中，TWSA逐漸增大，蒸散發量達到了最大值；TWSA峰值出現在10月，滯后降水量的峰值約2月。

圖5 2003—2019年各變量多年月均值

3.3 TWS各成分分析

圖6為黃河流域TWS各成分時間變化趨勢圖。從圖6可以看出，土壤水儲量以0.12 mm/月的速率平緩增長；雪水當量和冠層水儲量以0.001 8 mm/月的速率呈不明顯增長的趨勢；而地下水儲量以-0.52 mm/月的速率急劇下降，與TWSA變化速率-0.40 mm/月相近。

圖6 TWS成分在月尺度上的變化趨勢

對黃河流域年TWSA和年地下水儲量的相關性分析發現TWSA以-5.34 mm/a的速率下降，地下水儲量以-6.48 mm/a的速率快速下降，兩者變化趨勢相似。涂夢昭等[22]研究發現，黃土高原2005—2014年地下水儲量整體上以(-6.5±0.7) mm/a的速率急劇下降，與本文結果相近。年TWSA和年地下水儲量的Spearman秩相關系數為0.90(p<0.01)，兩者呈極顯著相關性，表明黃河流域TWS的急劇下降與地下水儲量急劇下降有關。

3.4 TWS影響因素分析

TWS受氣候因素和人類活動影響，本文基于2002年4月至2020年3月數據，重點分析氣候因素對TWS的影響。圖7為黃河流域TWS變化量與降水量和干燥度指數在年尺度上的變化關系。從圖7可以看出，年降水量上游和中游呈上升趨勢，下游呈下降趨勢；年干燥度指數上游呈下降趨勢，中游和下游呈上升趨勢。黃河流域上、中、下游年TWS變化量與年降水量和年干燥度指數分別呈極顯著的正相關關系和負相關關系，上、中、下游年TWS變化量與年降水量的Spearman秩相關系數分別為0.79、0.89和0.71；上、中、下游年TWS變化量與年干燥度指數的Spearman秩相關系數分別為-0.74、-0.93和-0.77，表明降水和蒸散發是黃河流域TWS變化量的重要影響因素。除氣候因素以外，人類活動(水庫調度、工農業用水等)對黃河流域TWS有較大影響[16]。

(a) 上游

4 結 論

a.基于LSTM模型重建GRACE和GRACE-FO間TWS變化量是可靠的，其與GLDAS和水量平衡方程計算的TWS變化量的相關系數分別為0.65和0.72。

b.研究期間，黃河流域TWSA分別以-0.40 mm/月和-5.34 mm/a的速率下降，上、中、下游TWSA分別以-0.16 mm/月、-0.67 mm/月和-1.83 mm/月的速率下降，表明黃河流域TWS變化空間差異顯著。TWSA最小值出現在6月，最大值出現在10月。黃河流域地下水儲量以-6.48 mm/a的速率快速下降，TWSA的急劇下降很可能是地下水急劇下降導致的，兩者的Spearman秩相關系數為0.90。

c.黃河流域上、中、下游年TWS變化量與年降水量和年干燥度指數分別呈極顯著的正相關關系和負相關關系。上、中、下游年TWS變化量與年降水量的Spearman秩相關系數分別為0.79、0.89和0.71；上、中、下游年TWS變化量與年干燥度指數的Spearman秩相關系數分別為-0.74、-0.93和-0.77，表明降水和蒸散發是黃河流域TWS變化量的重要影響因素。