雅魯藏布江流域陸地水儲(chǔ)量變化特征及歸因

劉 瀏，王宣宣，牛乾坤，倫玉蕊，程湫雅，程 磊，徐宗學(xué)

（1.中國農(nóng)業(yè)大學(xué)水利與土木工程學(xué)院，北京 100083；2.中國農(nóng)業(yè)大學(xué)中國農(nóng)業(yè)水問題研究中心，北京 100083；3.武漢大學(xué)水利水電學(xué)院，武漢 430072；4.北京師范大學(xué)水科學(xué)研究院，北京 100875；5.城市水循環(huán)與海綿城市技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室，北京 100875）

0 引 言

陸地水儲(chǔ)量（Terrestrial Water Storage，TWS）是指儲(chǔ)存在地表或地表以下的各種形式的水，包括地表水、地下水、土壤水、冰川、湖泊、林冠截留水和生物質(zhì)水等[1]，在陸地和全球水循環(huán)中扮演著關(guān)鍵的角色，它的組成成分對(duì)水文、氣候系統(tǒng)有著重要影響[2]。作為水資源平衡的重要組成部分，陸地水儲(chǔ)量變化（Terrestrial Water Storage Change，TWSC）是指陸地水儲(chǔ)量（TWS）在時(shí)間和空間上的演變過程[3]。陸地水儲(chǔ)量及其變化在評(píng)估和預(yù)測(cè)氣候變化、農(nóng)業(yè)氣象、洪水和干旱等極端事件中起著不可替代作用[4]。

青藏高原位于中國西南部、亞洲中部，是世界上海拔最高、面積最大的高原，面積約 250 萬 km2，有地球的“第三極”之稱[5]。青藏高原地區(qū)是中國境內(nèi)眾多主要河流的發(fā)源地，湖泊面積約占中國湖泊總面積的一半左右，因此又被稱為“亞洲水塔”[6]。青藏高原作為全球氣候變化最為敏感的地區(qū)之一，直接關(guān)系到中國西南地區(qū)的水資源利用和流域的用水安全狀況[7]。近幾十年來，相對(duì)于全球平均溫度每10 a升高0.17 ℃，青藏高原的平均溫度每10 a升高0.35 ℃左右，其增溫幅度遠(yuǎn)高于全球平均水平[8]。青藏高原地區(qū)的溫度升高導(dǎo)致了冰川融化速度加快、徑流不斷增加等一系列的失衡現(xiàn)象，導(dǎo)致其水循環(huán)過程隨之發(fā)生巨大的轉(zhuǎn)變[9]。冰川體積不斷縮小導(dǎo)致青藏高原地區(qū)的固態(tài)水儲(chǔ)量減少[10]，與冰川的變化相反，青藏高原地區(qū)有超過三分之二的湖泊在擴(kuò)張，湖泊的水量和面積整體都呈現(xiàn)出明顯的增加趨勢(shì)[11]。同時(shí)，青藏高原地區(qū)的徑流也呈現(xiàn)出不同程度的增加，進(jìn)而導(dǎo)致了青藏高原及其周邊地區(qū)的水災(zāi)風(fēng)險(xiǎn)加大。開展青藏高原陸地水儲(chǔ)量變化規(guī)律的研究關(guān)乎到中國以及東南亞各國的用水安全，有助于進(jìn)一步加強(qiáng)對(duì)高寒地帶缺資料地區(qū)的水循環(huán)與水資源演變規(guī)律的了解，同時(shí)為科學(xué)合理開發(fā)利用青藏高原地區(qū)的水資源提供參考[12]。

青藏高原地區(qū)缺乏高精度的水文氣象實(shí)測(cè)資料[13]，限制了對(duì)其陸地水儲(chǔ)量的定量研究。2002年3月發(fā)射的GRACE（Gravity Recovery and Climate Experiment）重力衛(wèi)星能夠提供較高精度的水文信號(hào)，為水儲(chǔ)量的反演開辟了新的途徑[14]。然而，GRACE數(shù)據(jù)反演的水儲(chǔ)量的時(shí)間序列較短，不足20 a，難以對(duì)青藏高原地區(qū)過去幾十年來在全球氣候變暖的背景下水循環(huán)加強(qiáng)的進(jìn)行全面深刻的解析。因此，本研究選取青藏高原東南部的雅魯藏布江流域?yàn)檠芯繀^(qū)，基于2003—2017年的GRACE重力衛(wèi)星數(shù)據(jù)、ERA5（the fifth-generation reanalysis product of the European Centre for Medium Range Weather Forecasts）再分析數(shù)據(jù)和 GLDAS（Global Land Data Assimilation System）數(shù)據(jù)，采用水量平衡法和組分相加法兩種方法，將反演的4種陸地水儲(chǔ)量結(jié)果與GRACE重力衛(wèi)星數(shù)據(jù)反演的陸地水儲(chǔ)量在時(shí)間和空間尺度上進(jìn)行適應(yīng)性評(píng)估，獲取陸地水儲(chǔ)量反演的較優(yōu)方法，并基于該方法系統(tǒng)探究雅魯藏布江流域長序列（1948—2017年）陸地水儲(chǔ)量的時(shí)空演變規(guī)律及主導(dǎo)因素，以期為探究資料匱乏地區(qū)的生態(tài)水文過程提供可靠參考。

1 數(shù)據(jù)與方法

1.1 研究區(qū)概況

雅魯藏布江流域位于 27°N～32°N、82°E～97°E 之間，地處青藏高原東南部、喜馬拉雅山脈北麓，流域面積約24 萬km2，干流長度超過2 000 km[15]，流域地勢(shì)西高東低，海拔垂直差異大[16]。雅魯藏布江流域水資源儲(chǔ)藏豐富[17]，流域水資源狀況關(guān)乎到中國西藏以及下游地區(qū)的印度、孟加拉國等國家的用水的安全。上游地區(qū)平均海拔4 600 m以上，以高山谷地等地貌為主[18]；中游地區(qū)主要為屬于藏南谷地，主要為溫暖半干旱氣候[19]，是人類活動(dòng)集中的地區(qū)；下游地區(qū)包括高原溫帶濕潤、半濕潤氣候區(qū)和山地亞熱帶濕潤氣候區(qū)，氣候相對(duì)溫暖濕潤[20]。雅魯藏布江是世界上海拔最高的河流之一，受特殊的地形條件和海拔高度的影響，流域呈現(xiàn)出對(duì)氣候變化的響應(yīng)較為敏感和脆弱的特征[21]。尤其是近幾十年來，雅魯藏布江流域凍土不斷融化[22]，活動(dòng)層增厚[23]，冰川退縮幅度明顯[24]，這些都對(duì)流域生態(tài)環(huán)境和水資源保護(hù)產(chǎn)生了巨大影響。因此，準(zhǔn)確反演雅魯藏布江流域陸地水儲(chǔ)量變化，有助于加深在全球氣候變暖背景下對(duì)流域水文、生態(tài)系統(tǒng)變化對(duì)人類活動(dòng)以及氣候變化響應(yīng)機(jī)理的理解。研究區(qū)域概況如圖1所示。

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

本研究使用的兩種數(shù)據(jù)集分別是歐洲中期天氣預(yù)報(bào)中心第五代再分析資料ERA5和美國的全球陸面數(shù)據(jù)同化模型GLDAS數(shù)據(jù)集。利用這兩種數(shù)據(jù)集提供的降水、蒸散發(fā)、徑流、土壤含水量和雪水當(dāng)量來估算TWSC，再將得到的TWSC與GRACE數(shù)據(jù)反演的TWSC進(jìn)行比較。

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

GRACE重力衛(wèi)星是由美國宇航局（National Aeronautics and Space Administration，NASA）和德國航天局（German Aerospace Center，DLR）合作研發(fā)的雙星系統(tǒng)，于2002年3月成功發(fā)射[25]。目前有德克薩斯大學(xué)空間研究中心（Center for Space Research，CSR）、德國地學(xué)研究中心（German Research Centre for Geosciences，GFZ）以及美國噴氣動(dòng)力實(shí)驗(yàn)室（Jet Propulsion Laboratory，JPL）提供較為成熟的GRACE數(shù)據(jù)產(chǎn)品[26]。本研究使用的是JPL提供的 2003—2017年的月數(shù)據(jù)產(chǎn)品（https://grace.jpl.nasa.gov/data/get-data/ jpl_global_mascons/），具體版本為 JPL RL06M Version 2.0，空間分辨率為 0.5°×0.5°。

1.2.2 ERA5再分析數(shù)據(jù)

本研究使用的 ERA5數(shù)據(jù)是來自歐洲中期天氣預(yù)報(bào)中心（European Centre for Medium Range Weather Forecasts，ECMWF）最新發(fā)布的第五代全球大氣再分析資料 ERA5-Land的月平均數(shù)據(jù)（https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysi s-era5-land-monthly-means?tab=form），其空間分辨率較之前的版本提升至0.1°×0.1°，精度相對(duì)ERA-interim亦有所改善[27]。

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

本研究使用的 GLDAS數(shù)據(jù)是美國國家航空航天局（NASA）戈達(dá)德空間飛行中心（Goddard Space Flight Center，GSFC）和美國海洋和大氣局（National Oceanic and Atmospheric Administration，NOAA）國家環(huán)境預(yù)報(bào)中心（National Centers for Environmental Prediction，NCEP）聯(lián)合發(fā)布的基于衛(wèi)星、陸面模式和地面觀測(cè)數(shù)據(jù)的同化產(chǎn)品，能夠向用戶提供多種驅(qū)動(dòng)數(shù)據(jù)[28]。目前，GLDAS結(jié)合MOSAIC、NOAH（NCEP，OSU，Air Force and Office of Hydrology）、CLM（Community Land Model）、VIC（Variable Infiltration Capacity）4個(gè)陸面模式提供了大量的陸面數(shù)據(jù)產(chǎn)品。本研究采用的 GLDAS數(shù)據(jù)版本為 GLDAS_NOAH025_M（https://disc.gsfc.nasa.gov/datasets? keywords=GLDAS&page=1&temporalResolution=1%20mont）。與其他GLDAS產(chǎn)品相比，GLDAS-NOAH的空間分辨率更高，時(shí)間跨度更長，在青藏高原地區(qū)的應(yīng)用較為廣泛[29-31]。

1.3 研究方法

本研究采用水量平衡法和組分相加法兩種方法，基于ERA5和GLDAS兩種數(shù)據(jù)集，分別計(jì)算雅魯藏布江流域的TWSC，共生成4套TWSC數(shù)據(jù)集，并基于GRACE數(shù)據(jù)，定量評(píng)估不同數(shù)據(jù)源、不同方法對(duì)雅魯藏布江流域水儲(chǔ)量估算的適用性。為了便于比較，ERA5和GLDAS數(shù)據(jù)均采用雙線性插值，插值為 0.5°×0.5°，與 GRACE數(shù)據(jù)的空間分辨率保持一致。

1.3.1 水量平衡法

水量平衡法（PER）根據(jù)水量平衡方程[32]，用降水（Precipitation）、蒸散發(fā)（Evapotranspiration）、徑流（Runoff）來計(jì)算TWSC，公式如下：

式中S(t)為t時(shí)刻的TWS，mm；P(t)為總降水量，mm；E(t)為蒸散發(fā)量，mm；R(t)為該流域的徑流量，mm。

1.3.2 組分相加法

TWSC也可以用其主要的組成成分的變化來表示[33]，大量研究表明，在全球尺度和流域尺度上，用土壤含水量和雪水當(dāng)量的變化估算的TWSC與GRACE數(shù)據(jù)反演的TWSC有很好的一致性[34-36]。下式是組分相加法（SS）用土壤含水量（Soil moisture）和雪水當(dāng)量（Snow water equivalent）的變化近似表示TWSC：

式中So(t)和Sn(t)分別為某一區(qū)域和時(shí)間內(nèi)的土壤含水量和雪水當(dāng)量，mm。

1.3.3 評(píng)估指標(biāo)

本研究采用線性回歸方程來衡量水儲(chǔ)量的變化。

式中y代表水儲(chǔ)量的變化，x代表時(shí)間，a代表斜率，b代表截距[37]。

采用皮爾遜相關(guān)系數(shù)來衡量不同變量之間的相關(guān)性。

式中xi和yi分別為兩組序列的數(shù)據(jù)，n是變量的序列數(shù)目，和分別為兩組數(shù)列的平均值，rxy為兩個(gè)變量之間的相關(guān)系數(shù)[38]。

Mann-Kendall檢驗(yàn)方法是是世界氣象組織推薦的非參數(shù)檢驗(yàn)方法，能夠客觀反映時(shí)間序列趨勢(shì)，目前已經(jīng)被廣泛用于氣候參數(shù)和水文序列的分析中[39]。Mann-Kendall方法通過給定一個(gè)時(shí)間序列Kt=(k1,k2, …,km)，計(jì)算統(tǒng)計(jì)量Zc，計(jì)算公式如下：

式中kp和kq分別表示第p個(gè)和第q個(gè)樣本，1＜p＜q＜m，m為樣本容量；D(s)表示s的方差；β表示變化趨勢(shì)，β＞0表示增加趨勢(shì)，β＜0表示減少趨勢(shì)。|Zc|＞Z(1-α/2)表示在α顯著性水平下，變化趨勢(shì)顯著。

2 結(jié)果與分析

2.1 陸地水儲(chǔ)量反演結(jié)果的適應(yīng)性評(píng)估

2.1.1 陸地水儲(chǔ)量時(shí)間特征對(duì)比分析

圖2是基于ERA5和GLDAS數(shù)據(jù)，分別采用水量平衡法和組分相加法估算的2003—2017年雅魯藏布江流域TWSC逐月變化及其與GRACE數(shù)據(jù)反演結(jié)果的對(duì)比圖。如圖2a所示，GRACE數(shù)據(jù)反演的雅魯藏布江流域的TWSC在 2003—2017年間以 0.90 mm/月的速度減少（P＜0.01），相當(dāng)于 2.46 億 t/a。不同數(shù)據(jù)源以及不同方法反演的陸地水儲(chǔ)量存在較大差異性，PER（ERA5）、SS（ERA5）和 PER（GLDAS）的結(jié)果均呈現(xiàn)周期波動(dòng)的狀態(tài)，無顯著的增減趨勢(shì)。4種結(jié)果中，SS（GLDAS）以0.397 mm/月的速度減少（P＜0.01），在月尺度變化趨勢(shì)上與GRACE數(shù)據(jù)反演的結(jié)果最為一致。

如圖2b所示，PER（ERA5）和SS（ERA5）變化幅度分別為-44.77～77.15和-50.04～71.81 mm，PER（GLDAS）和SS（GLDAS）范圍分別為-62.22～82.97和-69.54～113.36 mm，均與GRACE數(shù)據(jù)反演的TWSC變幅（-159.23～208.19 mm）存在一定差異，但相較于其他3套結(jié)果，SS（GLDAS）反演的水儲(chǔ)量變幅與 GRACE結(jié)果更為接近。此外，本研究進(jìn)一步檢驗(yàn)了 4套反演結(jié)果2003—2017年間180個(gè)月的TWSC與GRACE數(shù)據(jù)反演結(jié)果的相關(guān)性，PER（ERA5）、SS（ERA5）、PER（GLDAS）和SS（GLDAS）與GRACE結(jié)果的相關(guān)系數(shù)分別為-0.03、0.36、-0.12和0.75。由此看出，組分相加法估算的 TWSC的逐月變化與 GRACE數(shù)據(jù)反演的TWSC逐月變化的相關(guān)性均優(yōu)于水量平衡法。此外，無論是水量平衡法還是組分相加法，基于GLDAS數(shù)據(jù)反演的月尺度TWSC均與GRACE的反演結(jié)果的相關(guān)性更高。綜上所述，SS（GLDAS）在月尺度變化趨勢(shì)和變化幅度上，與GRACE數(shù)據(jù)反演的結(jié)果更為一致。

2.1.2 陸地水儲(chǔ)量空間特征對(duì)比分析

圖3是基于 ERA5和 GLDAS數(shù)據(jù)的 TWSC與GRACE數(shù)據(jù)反演結(jié)果相關(guān)系數(shù)的空間分布?；贓RA5數(shù)據(jù)，采用水量平衡法（PER）和組分相加法（SS）估算的上游和中游地區(qū)的相關(guān)系數(shù)較大，最大值分別為0.37和 0.67，而相關(guān)系數(shù)較小的地區(qū)主要集中在下游的尼洋河流域。基于GLDAS數(shù)據(jù)，采用水量平衡法估算的大部分地區(qū)的相關(guān)系數(shù)較小，僅流域東南部的相關(guān)系數(shù)略高，最大值為僅 0.21；采用組分相加法估算的大部分地區(qū)的相關(guān)系數(shù)相對(duì)較大，大于0.50的面積超過61.3%，最大值為 0.82，而相關(guān)系數(shù)較小的地區(qū)主要集中在下游的尼洋河流域。由此看出，陸地水儲(chǔ)量反演結(jié)果較差的地區(qū)主要集中在下游的尼洋河流域，即帕隆藏布冰川分布地區(qū)[40]，該地區(qū)雪線高度的高程變化超過2 000 m[41]，受特殊的海拔和地形的影響，模型存在較大的測(cè)量誤差，導(dǎo)致 GLDAS數(shù)據(jù)估算的結(jié)果難以準(zhǔn)確反映該地區(qū)的TWSC。與時(shí)間變化特征對(duì)比分析結(jié)果類似，SS（GLDAS）與GRACE數(shù)據(jù)反演結(jié)果的空間相關(guān)性最好。

綜合兩套數(shù)據(jù)、兩種算法得出的四組TWSC的時(shí)空變化的結(jié)果表明：無論是基于ERA5數(shù)據(jù)還是GLDAS數(shù)據(jù)，采用組分相加法估算的TWSC與GRACE數(shù)據(jù)反演的TWSC在時(shí)間和空間上的一致性均優(yōu)于水量平衡法。其中，SS（GLDAS）與GRACE重力衛(wèi)星數(shù)據(jù)反演的結(jié)果最為一致。

2.1.3 組分相加法敏感性分析

為探究不同組分對(duì) TWSC反演結(jié)果的影響，本研究基于 GLDAS數(shù)據(jù)和組分相加法進(jìn)行了TWSC反演的敏感性分析。分別考慮以下組合：1）土壤含水量、雪水當(dāng)量、地表水、地下水和植被冠層截留（圖4a）；2）土壤含水量、雪水當(dāng)量和地下水（圖4b）；3）土壤含水量和雪水當(dāng)量（圖4c）。三種組合得到的相關(guān)系數(shù)空間分布特征相似，下游的尼洋河流域相關(guān)性較差，統(tǒng)計(jì)結(jié)果如表1所示?？紤]土壤含水量、雪水當(dāng)量、地表水、地下水和植被冠層截留時(shí)，SS（GLDAS）與GRACE數(shù)據(jù)反演的流域平均 TWSC的相關(guān)系數(shù)為0.50，48.4%的地區(qū)相關(guān)系數(shù)大于 0.50?？紤]土壤含水量、雪水當(dāng)量和地下水時(shí)，SS（GLDAS）與 GRACE數(shù)據(jù)反演的流域平均TWSC的相關(guān)系數(shù)為0.51，58.0%的地區(qū)相關(guān)系數(shù)大于 0.50。僅考慮土壤含水量和雪水當(dāng)量時(shí)，SS（GLDAS）與 GRACE數(shù)據(jù)反演的流域平均TWSC的相關(guān)系數(shù)最大，為0.53，且該條件下，61.3%的地區(qū)相關(guān)系數(shù)大于 0.50。表明綜合考慮地表水、地下水和植物冠層截留時(shí)反演的 TWSC結(jié)果與 GRACE反演結(jié)果的一致性并未有所提升。此外，已有研究表明，在雅魯藏布江流域，相對(duì)于其他組分，土壤含水量和雪水當(dāng)量的變化更加顯著[10,34,37]。因此，基于土壤含水量和雪水當(dāng)量，采用組分相加法反演雅魯藏布江流域陸地水儲(chǔ)量的變化是合理的。

表1 組分相加法敏感性分析結(jié)果對(duì)比Table 1 Comparison of sensitivity analysis results of the summation method

因此，基于GLDAS數(shù)據(jù)的組分相加法更適用于雅魯藏布江流域陸地水儲(chǔ)量的反演，可為準(zhǔn)確描述高寒地區(qū)水儲(chǔ)量變化提供重要參考。

2.2 1948—2017年雅魯藏布江流域陸地水儲(chǔ)量變化特征

鑒于 GLDAS數(shù)據(jù)的組分相加法反演的 TWSC與GRACE重力衛(wèi)星數(shù)據(jù)反演的結(jié)果最為一致，本研究基于1948—2017年的GLDAS數(shù)據(jù)，采用組分相加法進(jìn)一步探究 1948—2017年雅魯藏布江流域陸地水儲(chǔ)量的變化特征。

2.2.1 1948—2017年TWSC的時(shí)間變化特征

1948—2017 年間 TWSC呈現(xiàn)極顯著上升趨勢(shì)（P＜0.01），但TWSC在2002年前后發(fā)生突變。如圖5所示，2002年之前TWSC以0.024 mm/月的速度呈現(xiàn)極顯著上升的趨勢(shì)（P＜0.01），2002年之后TWSC以0.397 mm/月的速度呈現(xiàn)極顯著下降的趨勢(shì)（P＜0.01）。這與劉瀏等[42]在雅魯藏布江流域干濕轉(zhuǎn)換特征研究中的結(jié)果保持一致，即1982—2015年雅魯藏布江流域總體呈現(xiàn)濕潤化趨勢(shì)，但進(jìn)入21世紀(jì)以后流域表現(xiàn)出極顯著的干旱化趨勢(shì)。

本研究通過小波分析探究1948—2017年雅魯藏布江流域TWSC的長時(shí)間的變化特征，如圖6所示。從大尺度周期來看，對(duì)應(yīng)于時(shí)頻分布圖的上部，TWSC存在著40 a左右的大尺度震蕩周期，但鑒于時(shí)間序列總長度為70 a，因此對(duì)40 a左右的大尺度周期不予考慮。從中尺度周期來看，對(duì)應(yīng)于時(shí)頻分布圖的中部，TWSC存在24 a左右的中尺度振蕩周期，經(jīng)歷3次高低起伏。從小尺度周期來看，對(duì)應(yīng)于時(shí)頻分布圖的下部，TWSC在1948—1975年期間存在著5 a左右的小尺度振蕩周期，在1976—2017年期間存在著8 a左右的小尺度振蕩周期，TWSC隨著時(shí)間尺度的減小表現(xiàn)出更為復(fù)雜的周期振蕩。小波方差圖中顯示，在1948—2017年期間，陸地水儲(chǔ)量的年變化的主周期為35 a。從中尺度和小尺度來看，TWSC在2002年左右都發(fā)生了一次由大到小的突變。

2.2.2 1948—2017年TWSC的空間變化特征

1948—2017 年雅魯藏布江流域TWSC從上游到下游呈現(xiàn)遞增的規(guī)律，如圖7a所示，陸地水儲(chǔ)量最大減少值為-10.26 mm，位于雅魯藏布江上游西部，最大增加值為38.36 mm，位于下游東南部，中游和地區(qū)的陸地水儲(chǔ)量基本保持不變或增加。因此，1948—2017年間雅魯藏布江流域 TWSC呈現(xiàn)顯著增加趨勢(shì)主要與中下游地區(qū)的陸地水儲(chǔ)量增加有關(guān)。如圖7b，流域大部分地區(qū)的TWSC變化整體呈增加趨勢(shì)，但在中游的“一江兩河”（雅魯藏布江中游及其支流拉薩河和年楚河）地區(qū)和下游冰川分布較多的帕隆藏布地區(qū)的多年平均TWSC變化呈現(xiàn)減少的趨勢(shì)，進(jìn)一步表明人類活動(dòng)和冰川積雪集中地帶的冰雪消融對(duì)水儲(chǔ)量變化趨勢(shì)有著不可忽略的影響。

圖8是2002年前后雅魯藏布江流域多年平均TWSC和變化趨勢(shì)的空間分布圖。在1948—2002年，相對(duì)于上游地區(qū)，中下游地區(qū)的陸地水儲(chǔ)量較為豐富，但是中下游大部分地區(qū)TWSC呈現(xiàn)減少的趨勢(shì)，表明雖然中下游地區(qū)陸地水儲(chǔ)量在1948—2002年間有所增加，但增加的趨勢(shì)卻在變緩。在2003—2017年，流域陸地水儲(chǔ)量減少的地區(qū)主要分布在除下游東南部以外的地區(qū)，占流域總面積的83.7%。陸地水儲(chǔ)量增加的地區(qū)主要集中在下游東南部，但是增加的趨勢(shì)也明顯變緩。2002年前后，中游地區(qū)的陸地水儲(chǔ)量轉(zhuǎn)變最為明顯，陸地水儲(chǔ)量由增加變?yōu)闇p少。這可能是由于近些年來中游地區(qū)人類活動(dòng)強(qiáng)度較大，城市建設(shè)用地增加、過度放牧導(dǎo)致植被覆蓋面積減少，降低了該地區(qū)的持水能力，進(jìn)而導(dǎo)致了該地區(qū)的陸地水儲(chǔ)量急劇減少。

2.3 TWSC變化的歸因分析

本研究基于1948—2017年的GLDAS數(shù)據(jù)進(jìn)行70 a TWSC變化的歸因分析，圖9是1948—2017年土壤含水量和雪水當(dāng)量的逐月變化。在2002年前后，土壤水變化和雪水當(dāng)量的變化趨勢(shì)與TWSC的變化趨勢(shì)基本一致。1948—2002年土壤含水量的變化以 0.014 mm/月的速率呈現(xiàn)出極顯著上升的趨勢(shì)（P＜0.01），2003—2017年土壤含水量的變化以 0.394 mm/月的速率呈現(xiàn)出極顯著下降的趨勢(shì)（P＜0.01）。1948—2002年雪水當(dāng)量的變化以0.009 mm/月的速率呈現(xiàn)出極顯著上升的趨勢(shì)（P＜0.01），2003—2017年雪水當(dāng)量的變化以0.002 mm/月的速率呈現(xiàn)出極顯著下降的趨勢(shì)（P＜0.01）。根據(jù)Meng等[43]計(jì)算不同組分對(duì)青藏高原陸地水儲(chǔ)量變化貢獻(xiàn)率的方法，本研究計(jì)算2002年前后土壤含水量和雪水當(dāng)量的變化對(duì)TWSC的貢獻(xiàn)率。1948—2002年，土壤含水量和雪水當(dāng)量的變化對(duì)TWSC的貢獻(xiàn)率分別為61%和39%；2003—2017年，土壤含水量和雪水當(dāng)量的變化對(duì)TWSC的貢獻(xiàn)率分別為99%和1%。相較于雪水當(dāng)量，2002年之后土壤含水量的變化速率更為顯著，這可能是由于 2002年以后，青藏高原地區(qū)氣溫增加，帕隆藏布地區(qū)冰川和積雪消融轉(zhuǎn)化為土壤水[44]，多年凍土向季節(jié)性凍土轉(zhuǎn)化，導(dǎo)致土壤含水量的變化增加[45]，土壤含水量的變化對(duì)TWSC的貢獻(xiàn)率顯著增加。同時(shí)，2002年以后流域積雪消融，積雪的面積分布和深度都有所減少，導(dǎo)致雪水當(dāng)量的變化對(duì)TWSC的貢獻(xiàn)率急劇減少[46]。

圖10是1948—2017年TWSC、土壤含水量變化和雪水當(dāng)量變化量的空間分布?？梢钥闯?，土壤含水量和TWSC的空間分布規(guī)律更加一致，2002年以后，TWSC和土壤含水量減少的區(qū)域都在從上游向下游逐漸擴(kuò)張，其中中游地區(qū)變化最為明顯。中游地區(qū)的TWSC和土壤含水量在1948—2002年間基本保持不變，在2003—2017年間有所減少。1948—2002年和2003—2017年土壤含水量減少的地區(qū)主要分布在上游，土壤含水量增加的地區(qū)主要集中在下游，這可能是由于下游地區(qū)由于海拔相對(duì)較低且溫度較高，植被覆蓋密度相對(duì)較大[37]，促進(jìn)了下游地區(qū)土壤的持水能力[47]。1948—2002年雪水當(dāng)量減少的地區(qū)主要集中在上游地區(qū)，中下游大部分地區(qū)主要呈現(xiàn)增加的趨勢(shì)。2003—2017年雅魯藏布江流域大部分地區(qū)的雪水當(dāng)量與1948—2002年的雪水當(dāng)量空間分布整體相似，但下游的帕隆藏布地區(qū)的雪水當(dāng)量卻明顯減少，這再次表明了氣候變暖對(duì)冰雪消融的促進(jìn)作用[10]。

綜上所述，在1948—2002年和2003—2017年期間，土壤含水量變化和雪水當(dāng)量的變化的趨勢(shì)與TWSC時(shí)間和空間上的變化趨勢(shì)基本一致。在時(shí)間上，土壤含水量變化和雪水當(dāng)量的變化都在2002年前后呈現(xiàn)出相反的變化趨勢(shì)。在空間分布上，土壤含水量變化和雪水當(dāng)量的變化呈現(xiàn)出明顯的空間異質(zhì)性，主要體現(xiàn)在人類活動(dòng)強(qiáng)度較高的“一江兩河”地區(qū)和冰川分布集中的帕隆藏布地區(qū)。但是，相對(duì)于雪水當(dāng)量的變化，土壤含水量的變化對(duì)TWSC的變化起主導(dǎo)作用。

3 結(jié) 論

本研究利用 2003—2017年 180個(gè)月的 ERA5（the fifth-generation reanalysis product of the European Centre for Medium Range Weather Forecasts）和GLDAS（Global Land Data Assimilation System）的月數(shù)據(jù)，采用水量平衡法和組分相加法計(jì)算雅魯藏布江流域 TWSC（Terrestrial Water Storage Change）變化。基于 GRACE（Gravity Recovery and Climate Experiment）數(shù)據(jù)反演的陸地水儲(chǔ)量，進(jìn)行不同數(shù)據(jù)源、不同計(jì)算方法在雅魯藏布江流域的適應(yīng)性評(píng)估?；贕LDAS數(shù)據(jù)的組分相加法進(jìn)一步研究1948—2017年的雅魯藏布江流域TWSC時(shí)空演變規(guī)律及其主導(dǎo)因素。主要結(jié)論如下：

1）綜合兩套數(shù)據(jù)、兩種算法得出的四組TWSC的時(shí)空變化的結(jié)果表明：基于GLDAS數(shù)據(jù)的組分相加法反演的TWSC與GRACE數(shù)據(jù)反演的結(jié)果最為一致，該方法能夠有效識(shí)別雅魯藏布江流域陸地水儲(chǔ)量變化狀況，對(duì)進(jìn)一步探究高寒地區(qū)的水-植被-生態(tài)復(fù)雜互饋機(jī)理具有重要意義。

2）相對(duì)于水量平衡法，組分相加法估算的雅魯藏布江流域TWSC效果更好。但是組分相加法未考慮地下水的變化，因此基于 GLDAS數(shù)據(jù)的組分相加法估算的TWSC與 GRACE數(shù)據(jù)反演的結(jié)果存在一定的差異。例如，基于GLDAS數(shù)據(jù)的組分相加法估算的TWSC的變化幅度（-69.54～113.36 mm）比 GRACE數(shù)據(jù)反演的TWSC變化幅度（-159.23～208.19 mm）小，且二者之間存在一定的時(shí)空差異。

3）Mann-Kendall檢驗(yàn)和小波分析的結(jié)果表明，TWSC在2002年發(fā)生了突變。且1948—2017年間TWSC、土壤含水量變化和雪水當(dāng)量變化都在2002年前后呈現(xiàn)出相反的變化趨勢(shì)，即2002年之前呈現(xiàn)增加的趨勢(shì)，2002年之后呈現(xiàn)減少的趨勢(shì)。在空間分布上，TWSC、土壤含水量變化和雪水當(dāng)量變化都存在明顯的空間異質(zhì)性，主要體現(xiàn)在人類活動(dòng)強(qiáng)度較高的“一江兩河”地區(qū)和冰川分布集中的帕隆藏布地區(qū)。2002年前后，土壤含水量的變化對(duì)TWSC的貢獻(xiàn)率分別為61%和99%，表明相對(duì)于雪水當(dāng)量的變化，雅魯藏布江流域土壤含水量變化對(duì)TWSC起主導(dǎo)作用。