遙感反演2000—2020年青藏高原水儲量變化及其驅(qū)動機(jī)制

龍 笛，李雪瑩，李興東，韓鵬飛，趙凡玉，洪仲坤，王一鳴，田富強(qiáng)

(清華大學(xué)水沙科學(xué)與水利水電工程國家重點(diǎn)實(shí)驗(yàn)室，北京 100084)

青藏高原被譽(yù)為“亞洲水塔”、“地球第三極”，平均海拔約4 000 m，是亞洲主要大河的發(fā)源地，包括亞洲重要的跨境河流印度河、恒河、雅魯藏布江-布拉馬普特拉河、怒江-薩爾溫江、瀾滄江-湄公河以及中國長江和黃河等。同時，青藏高原在亞洲氣候系統(tǒng)的穩(wěn)定和高原生態(tài)系統(tǒng)(山地森林、高寒草甸、高寒草原、農(nóng)田等)的碳收支動態(tài)變化方面具有重要的生態(tài)安全屏障作用，是亞洲乃至北半球環(huán)境變化的調(diào)控器[1]。以南亞季風(fēng)及西風(fēng)帶為主導(dǎo)的大氣環(huán)流[2]和青藏高原獨(dú)特的高山地形，共同形成了該地區(qū)豐富的水資源儲量，區(qū)域內(nèi)湖泊、冰川眾多，湖泊群總面積約5萬km2，占中國湖泊面積的一半以上，其中大于1 km2的湖泊數(shù)量接近1 400個[3-4]，是世界上湖泊分布最密集的地區(qū)之一；共有冰川約7.7萬條，總面積達(dá)8.3萬km2，是除兩極和格陵蘭島之外冰川數(shù)目最多的地區(qū)[5]。其中，藏東南地區(qū)是中國最主要的海洋性冰川集中區(qū)，7 700余條冰川的面積超過7 200 km2。此外，青藏高原冬、春兩季的季節(jié)性積雪覆蓋范圍廣、深度較深、持續(xù)時間長，是中國三大穩(wěn)定積雪區(qū)之一[6-7]。青藏高原各類水儲量(包括湖泊、冰川、積雪、凍土等)變化不僅與區(qū)域能量、大氣循環(huán)密切相關(guān)，而且對亞洲主要大河的冰雪融水補(bǔ)給有重要影響，進(jìn)一步對泛第三極地區(qū)的20 多個國家和30 多億人口的生存發(fā)展具有重要意義[8]。

氣候變化對青藏高原具有顯著而深刻的影響。研究表明，青藏高原是全球范圍內(nèi)對氣候變化最敏感的地區(qū)之一[9-10]?？傮w而言，近年來青藏高原正經(jīng)歷顯著的變暖變濕趨勢[11-13]，持續(xù)性的氣候變化使“亞洲水塔”水儲量發(fā)生失衡[8]，主要表現(xiàn)為藏東南地區(qū)冰川快速消融[14-15]、喀喇昆侖地區(qū)冰川正平衡(即喀喇昆侖異常)[16-17]、內(nèi)流區(qū)湖泊顯著擴(kuò)張[18-19]等。這些現(xiàn)象表明，“亞洲水塔”水儲量及相關(guān)水文過程正發(fā)生深刻改變，冰湖潰決、冰崩、泥石流等自然災(zāi)害頻率和強(qiáng)度不斷增加，對“亞洲水塔”及泛第三極地區(qū)的水資源利用、水災(zāi)害防治和可持續(xù)發(fā)展帶來一系列重大挑戰(zhàn)[20]。

區(qū)域水儲量變化的地面觀測手段極其有限。近年來在美國對地觀測系統(tǒng)(Earth Observing System)、歐洲哥白尼計(jì)劃(Copernicus Programme)哨兵(Sentinel)衛(wèi)星系列、中國高分等重大對地觀測計(jì)劃開展下，遙感的時間、空間、光譜和輻射分辨率不斷提高，水循環(huán)要素的反演算法不斷發(fā)展，多衛(wèi)星組網(wǎng)觀測和數(shù)據(jù)融合技術(shù)不斷成熟，為水儲量變化的大范圍時空動態(tài)監(jiān)測提供了可能[21-22]。本研究在作者前期工作及數(shù)據(jù)集的基礎(chǔ)上[18，23-26]，集成衛(wèi)星測高、ICESat-2(Ice，Cloud and land Elevation Satellite)陸冰產(chǎn)品、青藏高原500 m分辨率雪深產(chǎn)品、重力恢復(fù)與氣候?qū)嶒?yàn)衛(wèi)星(Gravity Recovery and Climate Experiment，GRACE)及其后續(xù)衛(wèi)星GRACE-Follow On(GRACE-FO)等多源遙感信息，對部分水儲量變化數(shù)據(jù)進(jìn)行延長，并結(jié)合再分析數(shù)據(jù)，解析2000—2020年(以下簡稱近20 a)青藏高原湖泊水儲量、冰川水儲量、雪水當(dāng)量和總水儲量的變化及其氣候驅(qū)動機(jī)制。本研究將對認(rèn)識“亞洲水塔”各類水儲量對氣候變化的響應(yīng)和驅(qū)動機(jī)制提供關(guān)鍵科學(xué)支撐，并為“亞洲水塔”失穩(wěn)的適應(yīng)性對策制定提供科學(xué)依據(jù)。

1 研究區(qū)域

本研究以青藏高原及其典型冰川區(qū)藏東南地區(qū)(圖1)作為研究區(qū)域。青藏高原是世界上平均海拔最高的高原，是北半球冰凍圈的重要組成部分。青藏高原位于亞洲中部，地理范圍為68°E—104°E，26°N—40°N，總面積約300萬km2，其中中國境內(nèi)面積約250萬km2[27]。該地區(qū)是中國重要的水能、水資源戰(zhàn)略儲備區(qū)，主要發(fā)源于青藏高原的中國西南諸河(包括藏西和藏南諸河、雅魯藏布江、滇西諸河、怒江、瀾滄江和元江)年出境水量約5 800億m3，約占中國年用水總量6 100億方(2020年)的95%。高原內(nèi)降水、氣溫空間差異顯著，藏東南地區(qū)年平均氣溫達(dá)20 ℃，年降水量超過1 000 mm；而西北部地區(qū)年平均氣溫低于0 ℃，年降水量不足100 mm[28]。

圖1 青藏高原及藏東南地理位置Fig.1 Location information of the Tibetan Plateau and Southeastern Tibetan Plateau

藏東南地區(qū)主要包括帕隆藏布(雅魯藏布江一級支流)流域及其周邊地區(qū)，受豐沛的降水和高海拔帶來的低溫影響，該區(qū)域是中國海洋性冰川分布最為集中的地區(qū)，冰川水資源儲量十分豐富，同時也是青藏高原冰川消融最為劇烈的地區(qū)之一。藏東南地區(qū)的冰川多發(fā)育于山谷，冰舌較長，通常延伸至海拔非常低的區(qū)域，主要分布在念青唐古拉山中東段、橫斷山區(qū)的伯舒拉嶺、喜馬拉雅山東端的南迦巴瓦峰區(qū)域。該地區(qū)山高谷深，平均海拔超過4 000 m，最高峰為喜馬拉雅山東端的南迦巴瓦峰，海拔7 782 m，而最低的雅魯藏布江大峽谷海拔不到1 000 m。受地形影響，藏東南地區(qū)形成了罕見的水汽通道，降水豐沛，是整個青藏高原最為濕潤的地區(qū)，喜馬拉雅山脈及崗日嘎布山以南地區(qū)年降水量在1 000 mm以上；念青唐古拉山南麓的察隅、波密和易貢年降水量也在790 mm以上[29]。

2 基礎(chǔ)數(shù)據(jù)

2.1 湖泊水儲量

本研究在前期工作基礎(chǔ)上，使用了包括ICESat(2003—2009年)、Envisat(2002—2011年)、CryoSat-2(2010—2020年)和Jason-1/2/3(2002—2020年)系列等多源測高數(shù)據(jù)，補(bǔ)充延長了52個大、中型湖泊測高水位。利用Landsat 5/7/8影像(2002—2017年)提取上述湖泊岸線變化和面積變化，并分別用于生成光學(xué)水位和建立湖泊面積—水位關(guān)系[18]。不同測高衛(wèi)星覆蓋的湖泊數(shù)量不同，ICESat覆蓋42個湖泊，Envisat覆蓋35個湖泊，CryoSat-2覆蓋全部52個湖泊，Jason-1/2/3覆蓋12個湖泊，光學(xué)水位覆蓋全部湖泊。在各種測高衛(wèi)星中，Envisat和Jason-1/2/3是傳統(tǒng)的脈沖有限雷達(dá)測高衛(wèi)星，地面足跡通常在2～4 km，其中Jason-1/2/3系列衛(wèi)星擁有現(xiàn)有測高數(shù)據(jù)中最短的重訪周期(10 d)；CryoSat-2是采用了InSAR技術(shù)的新一代雷達(dá)測高衛(wèi)星，具有較小的地面足跡(沿軌道方向約300 m，垂直軌道方向約1.65 km)；ICESat是激光雷達(dá)測高衛(wèi)星，地面足跡約170 m，具有較高的精度。Envisat和CryoSat-2數(shù)據(jù)均可從歐洲航天局(ESA)官網(wǎng)公開下載(https:∥earth.esa.int/eogateway)，Jason-1/2/3數(shù)據(jù)可通過AVISO+(Archiving，Validation and Interpretation of Satellite Oceanographic data)官網(wǎng)下載(https:∥www.aviso.altimetry.fr/en/home.html)，ICESat數(shù)據(jù)可通過美國國家冰雪數(shù)據(jù)中心(NSIDC)官網(wǎng)獲取(https:∥nsidc.org/)。

2.2 冰川表面高程

本研究利用測高衛(wèi)星ICESat-2計(jì)算冰川表面高程變化[23]。ICESat-2(觀測時段為2018年10月至今)是ICESat的后續(xù)衛(wèi)星，于2018年9月15日成功發(fā)射。其有效載荷是先進(jìn)地形激光高度計(jì)系統(tǒng)(Advanced Topographic Laser Altimeter System，ATLAS)，是目前高程精度最高的星載激光雷達(dá)。ICESat-2的軌道高度為496 km，重訪周期為91 d。ICESat-2 發(fā)射3對(6束)激光束，每對光束包括一強(qiáng)一弱2束光。每個波束對之間的距離為3.3 km，同一波束對中強(qiáng)弱波束之間的距離為跨軌90 m，沿軌2.5 km。同一波束的相鄰足跡點(diǎn)之間的距離為 0.7 m。與ICESat相比，ICESat-2在數(shù)據(jù)覆蓋密度、數(shù)據(jù)精度和空間分辨率方面都有顯著提高。ICESat-2 共有 1 387 條軌道覆蓋整個地球表面，具有0.1 m的高程精度，可獲取精細(xì)的冰面地形。本研究使用L3A 陸冰高程數(shù)據(jù)(ATL06)，空間分辨率約為 20 m，數(shù)據(jù)可通過美國宇航局(NASA)數(shù)據(jù)平臺EARTHDATA(https:∥search.earthdata.nasa.gov/search)公開獲取。

2.3 積雪面積和雪水當(dāng)量

本研究利用Dai等[24]開發(fā)的青藏高原500 m分辨率雪深產(chǎn)品對青藏高原積雪時空變化規(guī)律進(jìn)行研究。青藏高原500 m雪深產(chǎn)品綜合利用了MODIS積雪面積、地表溫度和被動微波亮度溫度來反演雪深，利用地表溫度與被動微波K和Ka波段水平極化亮度溫度的關(guān)系，消除土壤對亮度溫度差(TBD)的貢獻(xiàn)，然后將被動微波格網(wǎng)上的TBD分配到有雪像元，最后利用TBD和雪深的關(guān)系提取雪深，得到的雪深產(chǎn)品分辨率和準(zhǔn)確度都較高。此外，本研究根據(jù)由歐洲中期天氣預(yù)報(bào)中心(European Centre for Medium-Range Weather Forecasts，ECMWF)官網(wǎng)(https:∥www.ecmwf.int/en/era5-land)公開獲取的ERA5-Land再分析資料，利用其月均積雪變量相關(guān)數(shù)據(jù)(雪深、雪密度、雪水當(dāng)量)和上述500 m分辨率雪深數(shù)據(jù)進(jìn)行對比，分析青藏高原積雪變化規(guī)律。ERA5-Land數(shù)據(jù)的空間分辨率為0.1°×0.1°，時間分辨率為1 h。

2.4 總水儲量

重力衛(wèi)星GRACE(觀測時段為2002年4月至2017年6月)及其后續(xù)衛(wèi)星GRACE-FO(觀測時段為2018年5月至今)的研制和發(fā)射，為陸地總水儲量變化(以下簡稱總水儲量變化，包括地表水儲量(如冰川、湖泊、積雪等)、土壤水儲量和地下水儲量的總變化)的反演開辟了新途徑。重力衛(wèi)星由2顆完全一樣的衛(wèi)星構(gòu)成，軌道距離地球約450 km，2顆衛(wèi)星的距離約220 km，將雙星之間的距離變化作為反演地球時變重力場的基本信息，可進(jìn)一步獲得由地球質(zhì)量重分布引起的總水儲量變化。重力衛(wèi)星可提供月尺度，空間分辨率約10萬km2(約300 km×300 km)的總水儲量變化信息[30]。本研究對比4種不同的GRACE三級產(chǎn)品在青藏高原的表現(xiàn)，包括2種基于球諧系數(shù)的網(wǎng)格產(chǎn)品(JPL-SH和CSR-SH)和2種基于質(zhì)量濃度斑塊(mascon)的網(wǎng)格產(chǎn)品(JPL-M和CSR-M)。其中，2種球諧系數(shù)產(chǎn)品為RL05版本，其名義分辨率為1°×1°；2種mascon產(chǎn)品為RL06版本，其名義分辨率分別為0.5°×0.5°(JPL-M)和0.25°×0.5°(CSR-M)。重建GRACE和GRACE-FO之間的缺失數(shù)據(jù)時，采用最新的GRACE JPL RL06版本。所有重力衛(wèi)星數(shù)據(jù)均為總水儲量相對于2004—2009年時段的距平值，可由NASA官網(wǎng)公開獲取(http:∥grace.jpl.nasa.gov)。

3 研究方法

3.1 湖泊水儲量變化

采用湖泊面積變化和湖泊水位變化觀測計(jì)算湖泊水儲量變化。湖泊面積的獲取主要通過遙感影像識別湖泊邊界，湖泊水位的獲取可通過測高衛(wèi)星直接觀測。此外，可以通過對湖泊岸坡局部地形的分析，建立特定區(qū)域岸線位置與水位之間的關(guān)系來推求湖泊水位。本研究采用了基于衛(wèi)星測高(測高水位)和基于光學(xué)影像的岸線變化(光學(xué)水位)2種方法來推求湖泊水位[18]。

基于云平臺GEE(Google Earth Engine)和水體指數(shù)法提取湖泊面積和岸線變化包括4個基本步驟：① 構(gòu)建包含湖泊邊界/岸線的緩沖區(qū)；② 構(gòu)建水體指數(shù)；③ 在緩沖區(qū)內(nèi)進(jìn)行水體指數(shù)計(jì)算，并選取合適的閾值將水體指數(shù)灰度圖二值化；④ 根據(jù)二值化圖像計(jì)算水體面積后導(dǎo)出。常用的水體指數(shù)有歸一化差異水體指數(shù)(NDWI)和改進(jìn)的歸一化差異水體指數(shù)(MNDWI)。其中，MNDWI可以更有效識別渾濁水體，更適用于夏季；NDWI可以一定程度區(qū)分湖岸的積雪和湖水，更適用于冬季。自動閾值的選取通?；谧畲箢愰g方差法，即大津法(Otsu′s Method)。

衛(wèi)星測高最初應(yīng)用于海平面和冰蓋監(jiān)測，近年來在內(nèi)陸水體以及河、湖、冰研究中得到廣泛應(yīng)用[18，31-33]。雷達(dá)測高衛(wèi)星數(shù)據(jù)通常還需要進(jìn)行波形重定，目的是消除復(fù)雜地形引起的雷達(dá)回波信號異常，從而提高水位的反演精度。Huang等[31]提出的50% Threshold and Ice-1 Combined(TIC)波形重定算法，在改進(jìn)閾值法的基礎(chǔ)上，優(yōu)化了測高足跡的選取，并針對不同季節(jié)使用不同的閾值，該方法可以顯著改善河湖水位的反演精度。在湖泊水位提取中，由于整體地形平坦、水位變動幅度較小，僅使用TIC算法的50%閾值模式即可滿足應(yīng)用需求。另一方面，不同的測高衛(wèi)星軌道密度、重訪周期以及精度各有不同，使用多源測高數(shù)據(jù)能夠延長和加密湖泊測高水位序列，但往往需要消除不同傳感器之間的系統(tǒng)誤差。在獲取多源雷達(dá)測高水位和光學(xué)水位后，可以通過比較數(shù)據(jù)重疊期的平均水位值，移除不同源測高水位之間的系統(tǒng)偏差[18]。經(jīng)與實(shí)測水位檢驗(yàn)以及理論推導(dǎo)，修正系統(tǒng)誤差后的光學(xué)和測高水位精度(均方根誤差)均可達(dá)到0.2 m以內(nèi)[18，34]。

通過建立湖泊岸線變化與同期測高水位之間的統(tǒng)計(jì)回歸關(guān)系，可將岸線變化轉(zhuǎn)化為水位信息，即光學(xué)水位。由于光學(xué)水位覆蓋的時段最長，以光學(xué)水位為參照，通過比較數(shù)據(jù)重疊期的平均水位值，可將多源測高衛(wèi)星數(shù)據(jù)之間的系統(tǒng)誤差全部移除，形成時空一致性良好的高時間分辨率水位序列。獲取湖泊水位和面積變化時間序列后，可以通過積分得到湖泊水量變化(常見的方法有臺體公式，面積—水位關(guān)系積分得到體積—水位關(guān)系等)。

3.2 冰川質(zhì)量平衡

本研究利用ICESat-2 L3A陸冰高程(ATL06 V3)產(chǎn)品，計(jì)算藏東南地區(qū)的冰川表面高程變化[23]。使用的參考DEM是最新發(fā)布的NASADEM(空間分辨率約30 m)，通過計(jì)算ICESat-2與NASADEM之間的高程差，即可得到2000—2019年間冰川表面高程變化值，進(jìn)而求得冰川厚度的變化速率和冰川物質(zhì)平衡。經(jīng)與野外考察獲取的GPS測點(diǎn)高程以及無人機(jī)獲取的帕隆四號冰川超高分辨率DSM(數(shù)字表面模型)相比，ICESat-2數(shù)據(jù)在相同位置的高程誤差在3 m以內(nèi)，具有極高的精度。

Δh=hICESat-2-hNASADEM+hp

(1)

式中：Δh為冰川表面高程變化值；hICESat-2和hNASADEM分別對應(yīng)ICESat-2足跡點(diǎn)的高程值和對應(yīng)NASADEM的高程值；hp為NASADEM對冰川的穿透深度， 在本研究中hp=1.5 m。

(2)

(3)

式中：R為冰川厚度變化速率；M為冰川物質(zhì)平衡；ρg和ρw分別為冰川和水的密度，ρg=850 kg/m3。

3.3 總水儲量變化

3.3.1 時間序列分解方法

本研究采用STL(Seasonal and Trend decomposition using Loess)時間序列分解方法，獲取總水儲量變化的長期趨勢(式(4) )，以消除季節(jié)性變化和殘差對評估總水儲量變化的影響。

Stot=Slt+Ssea+R

(4)

式(4) 將總水儲量距平(Stot)分解為長期趨勢項(xiàng)(Slt)、 季節(jié)性變化(Ssea)和殘差項(xiàng)(R)。 對分解后的長期趨勢項(xiàng)序列作線性擬合， 可獲取總水儲量變化的趨勢。 STL以魯棒局部加權(quán)回歸作為平滑方法， 被廣泛應(yīng)用于總水儲量時間序列分析。

3.3.2 重力衛(wèi)星缺失數(shù)據(jù)重建

GRACE和GRACE-FO之間的數(shù)據(jù)缺失(2017年7月至2018年5月)，造成總水儲量變化監(jiān)測連續(xù)性不足。近年來，機(jī)器學(xué)習(xí)模型在水文研究中的應(yīng)用日趨成熟，通過學(xué)習(xí)輸入層及輸出層變量之間的非線性相關(guān)關(guān)系，機(jī)器學(xué)習(xí)算法可以基于相關(guān)輸入變量對重力衛(wèi)星缺失數(shù)據(jù)進(jìn)行重建[35]。應(yīng)用于重建重力衛(wèi)星數(shù)據(jù)的機(jī)器學(xué)習(xí)算法經(jīng)歷了由相對簡單的自回歸移動平均模型(ARIMA)、多元線性及非線性回歸模型(multiple linear and nonlinear regressions)向較復(fù)雜的人工神經(jīng)網(wǎng)絡(luò)(ANN)、深度神經(jīng)網(wǎng)絡(luò)(DNN)及卷積神經(jīng)網(wǎng)絡(luò)(CNN)發(fā)展。相關(guān)研究中，Mo等[36]在全球尺度利用貝葉斯卷積神經(jīng)網(wǎng)絡(luò)對GRACE缺失數(shù)據(jù)進(jìn)行重建，該研究參照Long等[37]的方法將總水儲量變化信號分解為趨勢項(xiàng)和非趨勢項(xiàng)，可改進(jìn)前人研究中氣候因子較難捕捉總水儲量長期變化趨勢的這一缺陷，提升重建結(jié)果的可信度。本研究補(bǔ)充Mo等[36]的重建結(jié)果，用以分析2002—2020年完整的青藏高原水儲量變化，具體實(shí)施步驟如下：

(1) 確定輸入層變量。所有變量均由再分析資料ERA5-Land數(shù)據(jù)獲取，包括降水、氣溫、模擬的總水儲量距平(土壤水、雪水及冠層水儲量的加和)，以及累積總水儲量變化量(一段時間內(nèi)總降水減去總蒸散及總徑流)；

(2) 對輸入層所有變量去除線性趨勢，僅保留年際和季節(jié)性變化；

(3) 確定輸出層變量。針對不同的GRACE三級產(chǎn)品，將對應(yīng)的GRACE及GRACE-FO時間序列去除線性趨勢，保留年際和季節(jié)性變化信號作為神經(jīng)網(wǎng)絡(luò)的目標(biāo)變量；

(4) 基于貝葉斯卷積神經(jīng)網(wǎng)絡(luò)，依據(jù)輸入的氣候變量重建缺失時段總水儲量年際和季節(jié)性變化信號；

(5) 缺失時段最終總水儲量變化信號為第(3)步確定的線性趨勢與第(4)步確定的年際和季節(jié)性變化信號的和。

4 結(jié)果與討論

4.1 2000—2020年青藏高原各類水儲量變化

4.1.1 湖泊水儲量顯著擴(kuò)張

本研究研制了2000—2020年青藏高原52個大、中型湖泊的高時間分辨率水位、水量數(shù)據(jù)集。52個湖泊水儲量總體的不確定性約為±6%[18]，整體呈擴(kuò)張趨勢，2000—2020年間湖泊總水儲量增加約130 Gt，超過3個三峽水庫的總庫容。湖泊的水量變化可劃分為3個階段：① 2000—2012年水量增加速率較快，為6.35 Gt/a，其中內(nèi)流區(qū)湖泊水量增長速度約為5.95 Gt/a；② 2012—2017年整體增速放緩至1.42 Gt/a，其中內(nèi)流區(qū)湖泊為1.20 Gt/a；③ 2017年后湖泊再次出現(xiàn)快速擴(kuò)張，整體增速達(dá)到10.59 Gt/a，其中內(nèi)流區(qū)湖泊達(dá)到8.18 Gt/a(圖 2(a))。從內(nèi)流區(qū)的年降水量時間序列上可以看出，年降水量顯著高于多年平均降水量時，內(nèi)流區(qū)湖泊整體水量會快速增加，二者具有較高的一致性(圖 2(b))。

圖2 青藏高原52個大、中型湖泊累積水量與青藏高原內(nèi)流區(qū)年降水量變化時間序列Fig.2 Accumulated lake water storage change for 52 large and middle lakes and annual precipitation in the endorheic basin on the Tibetan Plateau

從湖泊水量變化趨勢的空間分布看(圖 3)，內(nèi)流區(qū)湖泊擴(kuò)張趨勢比外流區(qū)更為顯著，但在部分流域呈現(xiàn)出“大湖擴(kuò)張，小湖收縮”的情形[18]。以色林措流域?yàn)槔执胨靠焖贁U(kuò)張(0.93 Gt/a)，而其上游的吳如措(-0.007 Gt/a)、格仁措(-0.007 Gt/a)和錯愕湖(-0.004 Gt/a)水量則緩慢下降，這些湖泊與色林措之間有河道相連?？赡茉蚴请S著水流對上游湖泊下泄口的侵蝕，下泄口高程緩慢降低，湖泊蓄水能力下降，而時間分辨率較低的數(shù)據(jù)集難以反映這種現(xiàn)象。在外流區(qū)，黃河上游流域湖泊水量增加較為顯著，例如青海湖(0.54 Gt/a)、鄂陵湖(0.08 Gt/a)、扎陵湖(0.03 Gt/a)等。雅魯藏布江流域的主要湖泊出現(xiàn)水量下降情況，如羊卓雍措(-0.14 Gt/a)。

圖3 青藏高原52個大、中型湖泊2000—2020年水量變化趨勢Fig.3 Spatial pattern of water storage trends of 52 large and middle lakes for 2000—2020 in the Tibetan Platean

湖泊持續(xù)擴(kuò)張可能引發(fā)湖泊潰決或溢流洪水災(zāi)害，如青?？煽晌骼锏貐^(qū)的卓乃湖-庫賽湖-海丁諾爾湖-鹽湖系統(tǒng)。位于上游的卓乃湖于2011年發(fā)生潰決，湖水下泄到庫賽湖并引起庫賽湖溢流，進(jìn)而注滿了海丁諾爾湖，最后注入鹽湖。據(jù)估算，潰決后從卓乃湖下泄進(jìn)入庫賽湖的水量達(dá)到2.47 Gt，造成庫賽湖水位上升近8 m并引發(fā)庫賽湖溢流。鹽湖的面積近年來擴(kuò)大了數(shù)倍，也存在溢流的風(fēng)險(xiǎn)，其下泄點(diǎn)距離青藏鐵路沿線僅10 km左右，對鐵路運(yùn)輸安全可能構(gòu)成威脅[18]。

4.1.2 藏東南地區(qū)冰川后退

通過2020年的ICESat-2數(shù)據(jù)與2000年的NASADEM，可得到2000—2020年間20 a整個藏東南地區(qū)的平均冰川表面高程變化速率，該數(shù)據(jù)集的總體誤差約為±25%。20 a間該地區(qū)平均冰川表面高程變化速率為(-0.71±0.18)m/a，相當(dāng)于每年損失44.7億t的冰川水資源。分析2018年10月至今近3 a的ICESat-2數(shù)據(jù)，得到近3 a藏東南地區(qū)季節(jié)平均的冰川表面高程變化時間序列如圖4所示，其中逐月季平均指取該月前后3個月數(shù)據(jù)作平均；季平均指取每個季度的數(shù)據(jù)作平均；滑動平均是對季平均時間序列作周期為4的滑動平均；線性回歸是對季平均時間序列的線性擬合。可知近3 a藏東南地區(qū)冰川表面高程變化速率為(-0.97±0.59)m/a。Zhao等[23]得出藏東南地區(qū)10 a(2011—2020年)的冰川表面高程變化速率為(-0.87±0.13)m/a，故藏東南地區(qū)的冰川消融在不斷加速。通過高空間覆蓋范圍的ICESat-2數(shù)據(jù)，分析藏東南地區(qū)冰川物質(zhì)平衡的空間分布特征(圖5)。在0.5°×0.5°的網(wǎng)格尺度，藏東南地區(qū)冰川物質(zhì)損失最嚴(yán)重的區(qū)域位于橫斷山區(qū)域(藏東南地區(qū)的東北部)，2000—2020年間物質(zhì)平衡速率達(dá)到了(-1.43±0.36)m w.e./a(w.e.表示水當(dāng)量)。藏東南地區(qū)東部的冰川物質(zhì)損失比西部地區(qū)要快很多，東部區(qū)域的表面高程變化速率約為-1.0 m/a，而西部區(qū)域的表面高程變化速率在-0.5 m/a左右。此外，藏東南地區(qū)成片分布的大型冰川的表面高程變化速率，要低于周邊零星分布的小型冰川的速率[23]。

圖 4 近3 a藏東南地區(qū)冰川表面高程季節(jié)變化Fig.4 Seasonal time series of cumulative glacier elevation changes in the Southeastern Tibetan Plateau in recent three years

圖5 2000—2020年藏東南地區(qū)冰川表面高程變化空間分布(0.5°×0.5°網(wǎng)格)Fig.5 Spatial distribution of glacier elevation change during 2000—2020 in the Southeastern Tibetan Plateau

4.1.3 積雪深度和雪水當(dāng)量

2001—2015年青藏高原逐日500 m分辨率雪深產(chǎn)品季節(jié)性分布顯示(圖 6)，總體上積雪主要分布在山區(qū)，冬春季雪深較深，其次是秋季，夏季積雪零星分布在山區(qū)。雪深年內(nèi)變化一般從10月開始增加，2月達(dá)到峰值，而后逐漸消融減少，至6月完全融化。對比500 m分辨率雪深產(chǎn)品和ERA5-Land再分析產(chǎn)品(圖6)可知，2種產(chǎn)品的雪深空間分布一致性較高，然而ERA5-Land產(chǎn)品雪深數(shù)值明顯高于青藏高原500 m分辨率雪深產(chǎn)品。這與再分析數(shù)據(jù)存在對雪深或積雪面積的系統(tǒng)性高估有關(guān)[38]，而過大的降雪量可能是ERA5-Land再分析產(chǎn)品中雪深高估的主要因素。

圖6 青藏高原500 m分辨率雪深產(chǎn)品2001—2015年四季平均雪深空間分布Fig.6 Spatial distributions of multi-year (2001—2015) averaged snow depth derived from a 500 m snow depth product for four seasons

圖7 ERA5-Land產(chǎn)品2001—2015年四季平均雪深空間分布Fig.7 Spatial distributions of multi-year (2001—2015) averaged snow depth derived from the ERA5-Land product for four seasons

根據(jù)ERA5-Land產(chǎn)品的積雪密度數(shù)值，進(jìn)而通過雪深、雪密度與雪水當(dāng)量之間的關(guān)系，計(jì)算得到青藏高原500 m分辨率雪深產(chǎn)品的雪水當(dāng)量時間序列，進(jìn)一步對比其和ERA5-Land生成的雪水當(dāng)量時間序列可知(圖 8)，ERA5-Land產(chǎn)品數(shù)值明顯高于500 m雪深產(chǎn)品，但二者相關(guān)性很高。其中，雪深時間序列相關(guān)系數(shù)達(dá)0.83，雪水當(dāng)量時間序列相關(guān)系數(shù)達(dá)0.86。因此，ERA5-Land產(chǎn)品可一定程度反映青藏高原積雪時空變化規(guī)律，但其數(shù)值存在顯著高估。

圖8 2001—2015年青藏高原500 m分辨率雪深產(chǎn)品和ERA5-Land產(chǎn)品平均雪深和雪水當(dāng)量Fig.8 Time series of regionally averaged snow depth and snow water equivalent over the Tibetan Plateau during 2001—2015 for a 500 m snow depth product and the ERA5 product

對比2001—2015和2016—2020年間ERA5-Land產(chǎn)品反映的青藏高原多年月平均雪水當(dāng)量(圖9)可以發(fā)現(xiàn)，相比2001—2015年，2016—2020年間除7—10月積雪較少的月份，其余月份多年月平均雪水當(dāng)量呈增加趨勢。從年平均雪水當(dāng)量看，2001—2015年間青藏高原多年平均雪水當(dāng)量為13.6 cm，而2016—2020年間多年平均雪水當(dāng)量為14.1 cm，增加了3.7%。

圖9 2001—2015年和2016—2020年間ERA5-Land產(chǎn)品反映的青藏高原多年月平均雪水當(dāng)量Fig.9 Multi-year monthly averaged snow water equivalent derived from the ERA5-Land product for 2001—2015 and 2016—2020 over the Tibetan Plateau

4.1.4 總水儲量

本研究在獲取2002—2020年完整的總水儲量變化時間序列后，利用STL時間序列分解獲得水儲量多年變化趨勢(圖 10)。4種GRACE反演結(jié)果均顯示2002—2020年青藏高原外流區(qū)(印度河、恒河、雅魯藏布江、怒江、瀾滄江等)總水儲量呈顯著下降趨勢；內(nèi)流區(qū)(塔里木、柴達(dá)木和羌塘盆地)的總水儲量呈顯著增加趨勢。但GRACE不同產(chǎn)品在估算總水儲量變化趨勢的量級方面差異較大，2種球諧系數(shù)產(chǎn)品的結(jié)果比較一致，2種mascon產(chǎn)品的結(jié)果在外流區(qū)比較一致，但JPL-M比CSR-M能反映更豐富的水儲量空間變化信息。

圖10 2002—2020年GRACE不同三級產(chǎn)品反演的青藏高原水儲量變化趨勢Fig.10 Trends in total water storage over the Tibetan Plateau during 2002—2020 estimated from different GRACE solutions

基于GRACE JPL-M估算的總水儲量變化趨勢(圖 11)，2002—2020年間青藏高原總水儲量以-4.50 Gt/a的速率下降。其中，外流區(qū)(印度河-恒河-雅魯藏布江-怒江-瀾滄江-長江-黃河流域)總水儲量的下降速率達(dá)到-10.90 Gt/a，內(nèi)流區(qū)(阿姆河-塔里木-羌塘盆地-柴達(dá)木-河西區(qū)域)總水儲量則以6.40 Gt/a的速率上升。在冰川廣泛分布的印度河-恒河-雅魯藏布江流域(青藏高原海洋性冰川集中分布的地區(qū))，總水儲量以-9.50 Gt/a的速度下降，冰川快速劇烈消融是該區(qū)總水儲量下降的主要原因。例如，藏東南地區(qū)(念青唐古拉山中東段,橫斷山區(qū)的伯舒拉嶺,喜馬拉雅山東端的南迦巴瓦峰區(qū)域)的冰川消融速率達(dá)到-4.50 Gt/a，占該地區(qū)總水儲量下降速率(-5.90 Gt/a)的76%。印度河-恒河-雅魯藏布江流域的冰川是極為重要的固態(tài)水資源，冰川與積雪融水是下游印度河-恒河平原的主要灌溉水源，影響約1.3億人的生產(chǎn)生活用水[39]，冰川快速消融已成為該地區(qū)最嚴(yán)重的威脅之一。對于湖泊廣泛分布的羌塘盆地，總水儲量增加速率為5.20 Gt/a，主要由湖泊擴(kuò)張導(dǎo)致(湖泊水儲量增加速率為4.90 Gt/a，占總水儲量變化的94%)。盡管內(nèi)流區(qū)湖泊水儲量與下游地區(qū)的用水量聯(lián)系較弱，但顯著的湖泊擴(kuò)張可能導(dǎo)致湖泊潰決，進(jìn)而威脅周邊基礎(chǔ)設(shè)施和居民生命財(cái)產(chǎn)安全，如2011年9月青藏高原內(nèi)流區(qū)卓乃湖和庫塞湖的潰決事件[40-41]。

圖11 基于GRACE JPL-M獲取的2002—2020年總水儲量變化時間序列Fig.11 Time series of total water storage changes over the Tibetan Plateau，Indus-Ganges-Brahmaputra，and Inner Tibetan Plateau basins during 2002—2020 estimated from GRACE JPL-M

4.2 水儲量變化驅(qū)動機(jī)制

4.2.1 湖泊

降水是驅(qū)動湖泊水量變化最主要的因素。在定量解析氣候變量對湖泊水量平衡貢獻(xiàn)的研究中(包括色林措、納木措、當(dāng)若雍措、瑪旁雍措、佩古措、可可西里湖、勒斜武擔(dān)措等湖泊)[42-43]，降水對不同湖泊貢獻(xiàn)率為59%～95%。青藏高原內(nèi)流區(qū)是湖泊最集中的區(qū)域，近30 a降水量呈增加趨勢[44]。色林措、納木措和當(dāng)若雍措均位于青藏高原內(nèi)流區(qū)南部，這一區(qū)域的湖泊在2000—2010年前后經(jīng)歷了快速擴(kuò)張期，之后進(jìn)入緩慢擴(kuò)張/相對穩(wěn)定期?？煽晌骼锖屠招蔽鋼?dān)措位于內(nèi)流區(qū)的北部，呈現(xiàn)持續(xù)的水量上漲趨勢。整體上這些湖泊的水量變化與降水量的變化保持了較高的一致性。

冰川融水的貢獻(xiàn)可以通過模型模擬或遙感觀測獲得。冰川融水對湖泊水量的貢獻(xiàn)大多為10%～20%，對佩古措水量的貢獻(xiàn)較大，達(dá)到了30%[4]。從冰川分布看，內(nèi)流區(qū)北部的冰川面積大于南部地區(qū)，尤其喀喇昆侖地區(qū)冰川較為集中，但由于該區(qū)冰川質(zhì)量整體保持穩(wěn)定，冰川融水對其流域內(nèi)湖泊擴(kuò)張的貢獻(xiàn)相對有限。由于觀測資料匱乏，分析凍土對湖泊水量的貢獻(xiàn)通?；谀Ｐ湍M。針對色林措流域的模擬結(jié)果顯示，凍土對色林措水量變化的貢獻(xiàn)達(dá)到了28%[45]。

蒸散也是影響湖泊水量平衡的重要變量，是水離開內(nèi)流流域的主要途徑。在暖濕化背景下，青藏高原蒸散整體呈增加趨勢。溫度升高、凍結(jié)期縮短都將導(dǎo)致湖泊總蒸發(fā)量增大。同時，青藏高原湖泊面積擴(kuò)張也可能導(dǎo)致水面蒸發(fā)總量增加，但青藏高原的湖泊水量不可能無限增長，將達(dá)到某種平衡狀態(tài)[44]。這種平衡狀態(tài)所對應(yīng)的湖泊面積、湖泊水量、氣候條件以及對氣候變化的響應(yīng)規(guī)律有待進(jìn)一步探索。

4.2.2 冰川

氣溫和降水是影響冰川物質(zhì)平衡的重要因素。降水的多少決定了冰川物質(zhì)積累量的大小，氣溫的高低決定了冰川消融(積累)的快慢，降水和氣溫的變化是冰川變化的重要驅(qū)動因素[46]?；谟^測數(shù)據(jù)的統(tǒng)計(jì)結(jié)果表明，1993—2012年全球平均升溫速率約為0.14 ℃/(10 a)[47]；藏東南地區(qū)的氣象站點(diǎn)數(shù)據(jù)資料顯示，2003—2018年間，藏東南地區(qū)年均氣溫以0.23 ℃/(10 a)的速率上升，顯著高于全球平均水平。年降水量則基本保持穩(wěn)定，但降水的年際變化幅度明顯高于以往。氣溫升高使冰川消融速度加快，同時，降水相態(tài)也隨之發(fā)生改變[48-50]。當(dāng)越來越多的降雪轉(zhuǎn)變?yōu)榻涤陼r，冰川的物質(zhì)積累量減少，且降雨對冰川的沖刷作用增強(qiáng)。近些年來，氣溫高和降水少的年份多次出現(xiàn)(2009年、2014年和2018年)，這些年份的冰川消融尤為劇烈。氣溫升高和降水波動是近年來藏東南地區(qū)冰川后退的主要原因。

4.2.3 總水儲量

在人類活動影響相對較小的青藏高原，氣候變化是引起該地區(qū)近20 a來總水儲量顯著變化的主要原因。總體而言，降水和氣溫是影響總水儲量變化最重要的氣候因子。降水作為總水儲量的輸入，是引起冰川物質(zhì)正平衡、湖泊擴(kuò)張、土壤和地下水儲量增加的重要變量；氣溫與冰川融化及蒸散發(fā)密切相關(guān)，氣溫上升是引起總水儲量下降的主要原因。研究表明近年來南亞季風(fēng)逐漸減弱[51]，導(dǎo)致青藏高原南部(喜馬拉雅地區(qū))由南亞季風(fēng)帶來的降水減少，進(jìn)而引起該區(qū)總水儲量下降。此外，近年來西風(fēng)帶增強(qiáng)[52]導(dǎo)致帕米爾高原降水有所增加，是青藏高原西部阿姆河流域冰川物質(zhì)基本保持平衡的主要原因之一[2]。季風(fēng)變化深刻影響了青藏高原降水變化的時空格局，進(jìn)一步影響青藏高原冰川等固態(tài)水儲量和總水儲量的變化。

基于觀測數(shù)據(jù)統(tǒng)計(jì)結(jié)果，青藏高原1955—1996年升溫速率約為0.16 ℃/(10 a)，其中冬季升溫速率高達(dá)0.23 ℃/(10 a)[13]。顯著的升溫導(dǎo)致青藏高原冰川快速消融，進(jìn)一步引起總水儲量變化。例如，在海洋性冰川集中分布的藏東南地區(qū)，氣溫升高是冰川融化和總水儲量減少的主要?dú)夂蛟颉?/p>

5 結(jié) 論

本研究主要基于多源遙感反演(包括衛(wèi)星重力、衛(wèi)星測高和光學(xué)遙感影像等)對2000—2020年青藏高原各類水儲量變化規(guī)律及其氣候驅(qū)動機(jī)制進(jìn)行了研究，主要結(jié)論如下：

(1) 氣候暖濕化下降水量和冰川融水增加，引發(fā)了青藏高原湖泊水量整體上漲。52個大、中型湖泊水量變化過程分為3個階段：2000—2012年為平穩(wěn)增長階段(6.35 Gt/a)，2012—2017年為相對穩(wěn)定階段(1.42 Gt/a)，2017年后進(jìn)入快速增長階段(10.59 Gt/a)。

(2) 受氣溫上升及其帶來的降水相態(tài)改變等多種因素影響，藏東南地區(qū)的冰川在快速消融，達(dá)到-4.50 Gt/a，且有加快趨勢。

(3) 2001—2015年間青藏高原多年平均雪水當(dāng)量為13.6 cm，而2016—2020年間多年平均雪水當(dāng)量為14.1 cm，增加了3.7%，積雪變化主要受累積期平均氣溫和降水影響。

(4) 在氣溫上升、南亞季風(fēng)減弱和西風(fēng)帶增強(qiáng)綜合影響下，青藏高原外流區(qū)總水儲量下降速率達(dá)到-10.90 Gt/a，內(nèi)流區(qū)總水儲量則以6.40 Gt/a的速率上升。

在2000—2020年間，氣候變化顯著改變了青藏高原水儲量的狀態(tài)，由此帶來水資源量改變及其上下游分配等問題將產(chǎn)生安全、生態(tài)、經(jīng)濟(jì)、政治連鎖效應(yīng)，制定全球和區(qū)域氣候變化的減緩和適應(yīng)性對策并付諸實(shí)踐迫在眉睫。