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

任立良，王 宇，江善虎，衛(wèi)林勇，王孟浩，張怡雅

(1.河海大學(xué)水文水資源與水利工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210098；2.河海大學(xué)水文水資源學(xué)院，江蘇 南京 210098)

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

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

1 研究區(qū)概況和數(shù)據(jù)來(lái)源

1.1 研究區(qū)概況

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

(a) 高程

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

1.2.1GRACE數(shù)據(jù)

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

1.2.2GLDAS數(shù)據(jù)

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

1.2.3降水量和氣溫?cái)?shù)據(jù)

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

1.2.4蒸散發(fā)和徑流數(shù)據(jù)

本文蒸散發(fā)和徑流數(shù)據(jù)來(lái)源于Noah模型2002年4月至2020年3月的模擬結(jié)果，空間分辨率為0.25°，蒸散發(fā)數(shù)據(jù)時(shí)間分辨率為1 s，徑流數(shù)據(jù)時(shí)間分辨率為3 h，經(jīng)累加計(jì)算將時(shí)間分辨率轉(zhuǎn)換為1月。

2 研究方法

2.1 LSTM神經(jīng)網(wǎng)絡(luò)模型

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

2.2 TWS變化量計(jì)算方法

基于GRACE和GRACE-FO衛(wèi)星數(shù)據(jù)反演和LSTM模型重建的逐月TWSA序列，計(jì)算逐月TWS變化量時(shí)間序列[20]：

(1)

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

基于GLDAS蒸散發(fā)及徑流數(shù)據(jù)，結(jié)合CPAP月降水?dāng)?shù)據(jù)，利用水量平衡方程計(jì)算黃河流域TWS變化量：

CSt=Pt-Et-Rt

(2)

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

2.3 地下水儲(chǔ)量計(jì)算方法

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

SG=ST-SSM-SSW-SCI

(3)

式中：SG為流域地下水儲(chǔ)量；ST為流域TWS；SSM為流域土壤水儲(chǔ)量；SSW為流域雪水當(dāng)量；SCI為流域冠層水儲(chǔ)量。

2.4 Spearman秩相關(guān)系數(shù)法

Spearman秩相關(guān)系數(shù)法常用于分析兩個(gè)變量之間關(guān)系的密切程度，計(jì)算公式為

(4)

3 結(jié)果與討論

3.1 GRACE和GRACE-FO間TWS變化量重建效果驗(yàn)證

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

(a) 訓(xùn)練期

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

圖3 TWS變化量在月尺度上的變化趨勢(shì)和相關(guān)性分析

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

3.2 TWS趨勢(shì)分析

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

圖4 TWSA在月尺度上的變化趨勢(shì)

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

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

3.3 TWS各成分分析

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

圖6 TWS成分在月尺度上的變化趨勢(shì)

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

3.4 TWS影響因素分析

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

(a) 上游

4 結(jié) 論

a.基于LSTM模型重建GRACE和GRACE-FO間TWS變化量是可靠的，其與GLDAS和水量平衡方程計(jì)算的TWS變化量的相關(guān)系數(shù)分別為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最小值出現(xiàn)在6月，最大值出現(xiàn)在10月。黃河流域地下水儲(chǔ)量以-6.48 mm/a的速率快速下降，TWSA的急劇下降很可能是地下水急劇下降導(dǎo)致的，兩者的Spearman秩相關(guān)系數(shù)為0.90。

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