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

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

（1.中國農業大學水利與土木工程學院，北京 100083；2.中國農業大學中國農業水問題研究中心，北京 100083；3.武漢大學水利水電學院，武漢 430072；4.北京師范大學水科學研究院，北京 100875；5.城市水循環與海綿城市技術北京市重點實驗室，北京 100875）

0 引 言

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

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

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

1 數據與方法

1.1 研究區概況

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

1.2 數據來源

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

1.2.1 GRACE數據

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

1.2.2 ERA5再分析數據

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

1.2.3 GLDAS數據

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

1.3 研究方法

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

1.3.1 水量平衡法

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

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

1.3.2 組分相加法

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

式中So(t)和Sn(t)分別為某一區域和時間內的土壤含水量和雪水當量，mm。

1.3.3 評估指標

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

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

采用皮爾遜相關系數來衡量不同變量之間的相關性。

式中xi和yi分別為兩組序列的數據，n是變量的序列數目，和分別為兩組數列的平均值，rxy為兩個變量之間的相關系數[38]。

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

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

2 結果與分析

2.1 陸地水儲量反演結果的適應性評估

2.1.1 陸地水儲量時間特征對比分析

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

如圖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數據反演的TWSC變幅（-159.23～208.19 mm）存在一定差異，但相較于其他3套結果，SS（GLDAS）反演的水儲量變幅與 GRACE結果更為接近。此外，本研究進一步檢驗了 4套反演結果2003—2017年間180個月的TWSC與GRACE數據反演結果的相關性，PER（ERA5）、SS（ERA5）、PER（GLDAS）和SS（GLDAS）與GRACE結果的相關系數分別為-0.03、0.36、-0.12和0.75。由此看出，組分相加法估算的 TWSC的逐月變化與 GRACE數據反演的TWSC逐月變化的相關性均優于水量平衡法。此外，無論是水量平衡法還是組分相加法，基于GLDAS數據反演的月尺度TWSC均與GRACE的反演結果的相關性更高。綜上所述，SS（GLDAS）在月尺度變化趨勢和變化幅度上，與GRACE數據反演的結果更為一致。

2.1.2 陸地水儲量空間特征對比分析

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

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

2.1.3 組分相加法敏感性分析

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

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

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

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

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

2.2.1 1948—2017年TWSC的時間變化特征

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

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

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

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

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

2.3 TWSC變化的歸因分析

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

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

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

3 結 論

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

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

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

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