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

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

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

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

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

1 研究區(qū)域

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

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

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

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

2.1 湖泊水儲(chǔ)量

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

2.2 冰川表面高程

本研究利用測(cè)高衛(wèi)星ICESat-2計(jì)算冰川表面高程變化[23]。ICESat-2(觀測(cè)時(shí)段為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對(duì)(6束)激光束，每對(duì)光束包括一強(qiáng)一弱2束光。每個(gè)波束對(duì)之間的距離為3.3 km，同一波束對(duì)中強(qiáng)弱波束之間的距離為跨軌90 m，沿軌2.5 km。同一波束的相鄰足跡點(diǎn)之間的距離為 0.7 m。與ICESat相比，ICESat-2在數(shù)據(jù)覆蓋密度、數(shù)據(jù)精度和空間分辨率方面都有顯著提高。ICESat-2 共有 1 387 條軌道覆蓋整個(gè)地球表面，具有0.1 m的高程精度，可獲取精細(xì)的冰面地形。本研究使用L3A 陸冰高程數(shù)據(jù)(ATL06)，空間分辨率約為 20 m，數(shù)據(jù)可通過(guò)美國(guó)宇航局(NASA)數(shù)據(jù)平臺(tái)EARTHDATA(https:∥search.earthdata.nasa.gov/search)公開獲取。

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

本研究利用Dai等[24]開發(fā)的青藏高原500 m分辨率雪深產(chǎn)品對(duì)青藏高原積雪時(shí)空變化規(guī)律進(jìn)行研究。青藏高原500 m雪深產(chǎn)品綜合利用了MODIS積雪面積、地表溫度和被動(dòng)微波亮度溫度來(lái)反演雪深，利用地表溫度與被動(dòng)微波K和Ka波段水平極化亮度溫度的關(guān)系，消除土壤對(duì)亮度溫度差(TBD)的貢獻(xiàn)，然后將被動(dòng)微波格網(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)行對(duì)比，分析青藏高原積雪變化規(guī)律。ERA5-Land數(shù)據(jù)的空間分辨率為0.1°×0.1°，時(shí)間分辨率為1 h。

2.4 總水儲(chǔ)量

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

3 研究方法

3.1 湖泊水儲(chǔ)量變化

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

基于云平臺(tái)GEE(Google Earth Engine)和水體指數(shù)法提取湖泊面積和岸線變化包括4個(gè)基本步驟：① 構(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可以更有效識(shí)別渾濁水體，更適用于夏季；NDWI可以一定程度區(qū)分湖岸的積雪和湖水，更適用于冬季。自動(dòng)閾值的選取通?；谧畲箢愰g方差法，即大津法(Otsu′s Method)。

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

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

3.2 冰川質(zhì)量平衡

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

Δh=hICESat-2-hNASADEM+hp

(1)

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

(2)

(3)

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

3.3 總水儲(chǔ)量變化

3.3.1 時(shí)間序列分解方法

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

Stot=Slt+Ssea+R

(4)

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

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

GRACE和GRACE-FO之間的數(shù)據(jù)缺失(2017年7月至2018年5月)，造成總水儲(chǔ)量變化監(jiān)測(cè)連續(xù)性不足。近年來(lái)，機(jī)器學(xué)習(xí)模型在水文研究中的應(yīng)用日趨成熟，通過(guò)學(xué)習(xí)輸入層及輸出層變量之間的非線性相關(guān)關(guān)系，機(jī)器學(xué)習(xí)算法可以基于相關(guān)輸入變量對(duì)重力衛(wèi)星缺失數(shù)據(jù)進(jìn)行重建[35]。應(yīng)用于重建重力衛(wèi)星數(shù)據(jù)的機(jī)器學(xué)習(xí)算法經(jīng)歷了由相對(duì)簡(jiǎn)單的自回歸移動(dò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ò)對(duì)GRACE缺失數(shù)據(jù)進(jìn)行重建，該研究參照Long等[37]的方法將總水儲(chǔ)量變化信號(hào)分解為趨勢(shì)項(xiàng)和非趨勢(shì)項(xiàng)，可改進(jìn)前人研究中氣候因子較難捕捉總水儲(chǔ)量長(zhǎng)期變化趨勢(shì)的這一缺陷，提升重建結(jié)果的可信度。本研究補(bǔ)充Mo等[36]的重建結(jié)果，用以分析2002—2020年完整的青藏高原水儲(chǔ)量變化，具體實(shí)施步驟如下：

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

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

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

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

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

4 結(jié)果與討論

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

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

本研究研制了2000—2020年青藏高原52個(gè)大、中型湖泊的高時(shí)間分辨率水位、水量數(shù)據(jù)集。52個(gè)湖泊水儲(chǔ)量總體的不確定性約為±6%[18]，整體呈擴(kuò)張趨勢(shì)，2000—2020年間湖泊總水儲(chǔ)量增加約130 Gt，超過(guò)3個(gè)三峽水庫(kù)的總庫(kù)容。湖泊的水量變化可劃分為3個(gè)階段：① 2000—2012年水量增加速率較快，為6.35 Gt/a，其中內(nèi)流區(qū)湖泊水量增長(zhǎng)速度約為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ū)的年降水量時(shí)間序列上可以看出，年降水量顯著高于多年平均降水量時(shí)，內(nèi)流區(qū)湖泊整體水量會(huì)快速增加，二者具有較高的一致性(圖 2(b))。

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

從湖泊水量變化趨勢(shì)的空間分布看(圖 3)，內(nèi)流區(qū)湖泊擴(kuò)張趨勢(shì)比外流區(qū)更為顯著，但在部分流域呈現(xiàn)出“大湖擴(kuò)張，小湖收縮”的情形[18]。以色林措流域?yàn)槔?，色林措水量快速擴(kuò)張(0.93 Gt/a)，而其上游的吳如措(-0.007 Gt/a)、格仁措(-0.007 Gt/a)和錯(cuò)愕湖(-0.004 Gt/a)水量則緩慢下降，這些湖泊與色林措之間有河道相連。可能原因是隨著水流對(duì)上游湖泊下泄口的侵蝕，下泄口高程緩慢降低，湖泊蓄水能力下降，而時(shí)間分辨率較低的數(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個(gè)大、中型湖泊2000—2020年水量變化趨勢(shì)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ū)的卓乃湖-庫(kù)賽湖-海丁諾爾湖-鹽湖系統(tǒng)。位于上游的卓乃湖于2011年發(fā)生潰決，湖水下泄到庫(kù)賽湖并引起庫(kù)賽湖溢流，進(jìn)而注滿了海丁諾爾湖，最后注入鹽湖。據(jù)估算，潰決后從卓乃湖下泄進(jìn)入庫(kù)賽湖的水量達(dá)到2.47 Gt，造成庫(kù)賽湖水位上升近8 m并引發(fā)庫(kù)賽湖溢流。鹽湖的面積近年來(lái)擴(kuò)大了數(shù)倍，也存在溢流的風(fēng)險(xiǎn)，其下泄點(diǎn)距離青藏鐵路沿線僅10 km左右，對(duì)鐵路運(yùn)輸安全可能構(gòu)成威脅[18]。

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

通過(guò)2020年的ICESat-2數(shù)據(jù)與2000年的NASADEM，可得到2000—2020年間20 a整個(gè)藏東南地區(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é)平均的冰川表面高程變化時(shí)間序列如圖4所示，其中逐月季平均指取該月前后3個(gè)月數(shù)據(jù)作平均；季平均指取每個(gè)季度的數(shù)據(jù)作平均；滑動(dòng)平均是對(duì)季平均時(shí)間序列作周期為4的滑動(dòng)平均；線性回歸是對(duì)季平均時(shí)間序列的線性擬合。可知近3 a藏東南地區(qū)冰川表面高程變化速率為(-0.97±0.59)m/a。Zhao等[23]得出藏東南地區(qū)10 a(2011—2020年)的冰川表面高程變化速率為(-0.87±0.13)m/a，故藏東南地區(qū)的冰川消融在不斷加速。通過(guò)高空間覆蓋范圍的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月完全融化。對(duì)比500 m分辨率雪深產(chǎn)品和ERA5-Land再分析產(chǎn)品(圖6)可知，2種產(chǎn)品的雪深空間分布一致性較高，然而ERA5-Land產(chǎn)品雪深數(shù)值明顯高于青藏高原500 m分辨率雪深產(chǎn)品。這與再分析數(shù)據(jù)存在對(duì)雪深或積雪面積的系統(tǒng)性高估有關(guān)[38]，而過(guò)大的降雪量可能是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)而通過(guò)雪深、雪密度與雪水當(dāng)量之間的關(guān)系，計(jì)算得到青藏高原500 m分辨率雪深產(chǎn)品的雪水當(dāng)量時(shí)間序列，進(jìn)一步對(duì)比其和ERA5-Land生成的雪水當(dāng)量時(shí)間序列可知(圖 8)，ERA5-Land產(chǎn)品數(shù)值明顯高于500 m雪深產(chǎn)品，但二者相關(guān)性很高。其中，雪深時(shí)間序列相關(guān)系數(shù)達(dá)0.83，雪水當(dāng)量時(shí)間序列相關(guān)系數(shù)達(dá)0.86。因此，ERA5-Land產(chǎn)品可一定程度反映青藏高原積雪時(shí)空變化規(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

對(duì)比2001—2015和2016—2020年間ERA5-Land產(chǎn)品反映的青藏高原多年月平均雪水當(dāng)量(圖9)可以發(fā)現(xiàn)，相比2001—2015年，2016—2020年間除7—10月積雪較少的月份，其余月份多年月平均雪水當(dāng)量呈增加趨勢(shì)。從年平均雪水當(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 總水儲(chǔ)量

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

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

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

圖11 基于GRACE JPL-M獲取的2002—2020年總水儲(chǔ)量變化時(shí)間序列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 水儲(chǔ)量變化驅(qū)動(dòng)機(jī)制

4.2.1 湖泊

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

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

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

4.2.2 冰川

氣溫和降水是影響冰川物質(zhì)平衡的重要因素。降水的多少?zèng)Q定了冰川物質(zhì)積累量的大小，氣溫的高低決定了冰川消融(積累)的快慢，降水和氣溫的變化是冰川變化的重要驅(qū)動(dòng)因素[46]?；谟^測(cè)數(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)定，但降水的年際變化幅度明顯高于以往。氣溫升高使冰川消融速度加快，同時(shí)，降水相態(tài)也隨之發(fā)生改變[48-50]。當(dāng)越來(lái)越多的降雪轉(zhuǎn)變?yōu)榻涤陼r(shí)，冰川的物質(zhì)積累量減少，且降雨對(duì)冰川的沖刷作用增強(qiáng)。近些年來(lái)，氣溫高和降水少的年份多次出現(xiàn)(2009年、2014年和2018年)，這些年份的冰川消融尤為劇烈。氣溫升高和降水波動(dòng)是近年來(lái)藏東南地區(qū)冰川后退的主要原因。

4.2.3 總水儲(chǔ)量

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

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

5 結(jié) 論

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

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

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

(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ū)總水儲(chǔ)量下降速率達(dá)到-10.90 Gt/a，內(nèi)流區(qū)總水儲(chǔ)量則以6.40 Gt/a的速率上升。

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