李心如， 盧善龍， 王 勇， 李明陽， 方 純，杜 聰， 王 南

（1. 中國科學院 空天信息創(chuàng)新研究院 數(shù)字地球重點實驗室，北京100094； 2. 可持續(xù)發(fā)展大數(shù)據(jù)國際研究中心，北京 100094；3. 中國科學院大學，北京 100049； 4. 湖南科技大學，湖南 湘潭 411201； 5. 新疆大學，新疆 烏魯木齊 830046）

0 引言

全球氣候驅(qū)動影響下的冰川融化［1］、水文循環(huán)加速等問題［2-3］，使得全球很多湖泊呈擴張趨勢，這一趨勢的發(fā)展可能會導致內(nèi)流水系與外流水系重組等一系列重大變化［4-5］，進而引發(fā)湖泊潰堤、突發(fā)洪水等地質(zhì)危害。例如，氣候變化引起的冰川流失會減少漫灘和近海沉積［6-7］；河流流向的轉(zhuǎn)換會引發(fā)流域盆地的突然重組［8］；湖泊溢流會引起沖溝的發(fā)育、產(chǎn)生溯源侵蝕，使決口不斷擴大，最終導致潰壩［9］。當前全球范圍內(nèi)部分湖泊受氣候變化影響而產(chǎn)生的溢流情況，對河谷侵蝕、地貌改變都產(chǎn)生了一定的影響。比如，湖泊溢流洪水驅(qū)動引起Licus Vallis 山谷的河谷切割［10］、全球氣候變化使Owens湖在高穩(wěn)定水位時溢出［11］；Chicago 湖在湖水高位期間溢出［12］；Agassiz冰川湖的水位變化、溢流，向北大西洋注入淡水，淡水的流入削弱了經(jīng)向翻轉(zhuǎn)環(huán)流，阻礙了熱量向北緯的輸送，觸發(fā)了12.9—11.7千年前的較年輕的仙女木寒冷事件［13-15］；科羅拉多河橫截面發(fā)生溢流而導致峽谷切割［16］；秘魯Cordillera Blanca 山谷冰川湖暴發(fā)洪水導致谷底沉積物被搬運遷移［17］。

根據(jù)前人研究，全球范圍有101個大型內(nèi)-外流盆地具有潛在水力聯(lián)系，其中，青藏高原有26個［18］。青藏高原是世界上海拔最高的高原，也是全球氣候變化最為敏感的地帶之一。在過去幾十年，青藏高原經(jīng)歷了重大的環(huán)境變化［19］。而氣候變化對青藏高原的水儲量造成顯著影響［20-22］，研究表明，青藏高原水資源對氣候因子有明顯的響應特征［23］。最新的CMIP6 模擬試驗數(shù)據(jù)結(jié)果表明青藏高原年均地表氣溫在21世紀90年代上升2.5 ℃，年均降水在21世紀90 年代將增加12.8%［24］，青藏高原地表水資源量呈上升趨勢，相應的湖泊水位提高和湖泊面積擴大［25］。

由于青藏高原氣象和水文觀測站點少、數(shù)據(jù)獲取困難，衛(wèi)星遙感手段被廣泛應用于區(qū)域水文過程變化研究［26-28］。Li等［29］通過結(jié)合多個測高任務和光學遙感圖像，為2000 年至2017 年期間青藏高原上的52 個大型湖泊開發(fā)了高時間分辨率水位和存儲變化數(shù)據(jù)集。馬山木等［30］基于ICESat-2 衛(wèi)星陸地觀測產(chǎn)品數(shù)據(jù)覆蓋情況，對2018年10月—2021年4月期間青藏高原面積大于1 km2的473 個湖泊進行了高精度水位動態(tài)監(jiān)測。Yu等［31］研究古氣候表明，自末次間冰期以來持續(xù)干旱，可能導致青藏高原廣泛而漸進的湖面下降，因此青藏高原現(xiàn)代湖泊水位低于古湖泊水位。Zhang 等［32］研究發(fā)現(xiàn)，近幾十年來，青藏高原的水文循環(huán)顯著加劇，凈降水量的增加是湖泊水量增加的主要原因，其次是冰川質(zhì)量損失和多年凍土退化導致的地面冰融化。因此，青藏高原湖泊水位升高，會引起湖泊外溢、潰決等潛在風險。Cheng 等［33］通過衛(wèi)星和無人機遙感數(shù)據(jù)模擬了湖泊擴張、潰決時的風險威脅，并對實際湖泊附近的村莊搬遷提供了方案。

隨著全球氣候變暖，湖泊水位升高、水面范圍擴大，近幾十年來的氣候變化導致青藏高原內(nèi)流盆地大多數(shù)湖泊呈顯著擴張趨勢，擴大的湖泊淹沒區(qū)域可能對生態(tài)環(huán)境和當?shù)厝祟惿瞽h(huán)境造成不利影響和潛在威脅［33］。例如，青藏高原可可西里卓乃湖在2011 年發(fā)生潰決，連接下游庫賽湖、海丁諾爾湖和鹽湖［34-35］，導致鹽湖面積明顯擴大，直至2019年9月，鹽湖水位上升至應急排水工程引水口，重建了鹽湖與長江支流的水力聯(lián)系［3］。

為了探究青藏高原26 個有潛在溢出可能的湖泊的外溢潛在風險［18］，本文以衛(wèi)星遙感和氣象再分析數(shù)據(jù)為數(shù)據(jù)源，采用地統(tǒng)計分析方法，分析了這些湖泊的面積變化，以及與相應湖泊流域內(nèi)降水、氣溫和蒸發(fā)之間的相關(guān)性特征；在此基礎上，通過情景假設，分析了未來區(qū)域降水、氣溫等參量變化對湖泊水體面積的影響，并對這些湖泊未來水位上漲滿溢的可能性開展了推測與判斷，以期為區(qū)域未來氣候水文變化規(guī)律研究及致災風險應對策略的制定提供數(shù)據(jù)和理論參考。

1 研究區(qū)概況與數(shù)據(jù)獲取

1.1 研究區(qū)概況

青藏高原位于26°00′～39°47′ N，73°19′～104°47′ E，是我國面積最大、世界海拔最高的高原。青藏高原東西長約2 800 km，南北寬約300～1 500 km，總面積約250×104km2，地形上可分為羌塘高原、藏南谷地、柴達木盆地、祁連山地、青海高原和川藏高山峽谷區(qū)等6 個部分，包括中國西藏全部和青海、新疆、甘肅、四川、云南的部分以及不丹、尼泊爾、印度、巴基斯坦、阿富汗伊斯蘭共和國、塔吉克斯坦、吉爾吉斯斯坦的部分或全部［36-38］。青藏高原一般海拔在3 000～5 000 m之間，平均海拔4 000 m以上，為東亞、東南亞和南亞許多大河流發(fā)源地；青藏高原是我國湖泊分布最廣泛的地區(qū)，面積超過1 km2的湖泊數(shù)量超過1 400 個，總面積超過5×104km2［39］。青藏高原水資源分布呈東多西少，南多北少的分布態(tài)勢，水資源分布極其集中，主要集中在山南、林芝市，雅魯藏布江流域及藏南諸河流域［20］（圖1）。青藏高原正在變暖變濕，但是高原東側(cè)部分地區(qū)正在變暖變干，同時高原整體風速正在減小。升溫主要是夜間的最低溫度貢獻的，不同地區(qū)升溫速率有差異，中部地區(qū)高于東部地區(qū)。青藏高原平均溫度、最低和最高溫度的分布與高原海拔走勢一致。降水量空間分布上表現(xiàn)出從東南向西北逐級減少的分布特點。風速自東南向西北呈現(xiàn)出小—大—小的變化趨勢［40］。

圖1 青藏高原河流湖泊分布Fig. 1 Distribution of rivers and lakes in Qinghai-Tibet Plateau

1.2 數(shù)據(jù)獲取

本文的研究數(shù)據(jù)包含：①地表水遙感數(shù)據(jù)；②區(qū)域地面氣象要素驅(qū)動數(shù)據(jù)集；③青藏高原月平均地表蒸散發(fā)數(shù)據(jù)集；④數(shù)字高程數(shù)據(jù)；⑤其他矢量數(shù)據(jù)；⑥Landsat 8衛(wèi)星遙感數(shù)據(jù)。

（1）地表水遙感數(shù)據(jù)

本研究中的湖泊水面數(shù)據(jù)是來源于Google Earth Engine（GEE）平臺的全球水面分布數(shù)據(jù)集（JRC Yearly Water Classification History，v1.3），包含1984—2020年地表水的位置和時間分布圖、水面的范圍和變化的統(tǒng)計數(shù)據(jù)［41］。本研究中提取了青藏高原內(nèi)流盆地中主要湖泊1985—2020 年的永久水（一年四季都有水）和季節(jié)水（包括冰凍期，一年中在水下的時間不到12 個月），并對異常數(shù)據(jù)進行鄰近年份數(shù)據(jù)插值處理，從而獲取研究區(qū)近幾十年來的水體面積變化完整數(shù)據(jù)。

（2）區(qū)域地面氣象要素驅(qū)動數(shù)據(jù)集

氣溫和降水數(shù)據(jù)來源于國家青藏高原科學數(shù)據(jù)中心的“中國區(qū)域地面氣象要素驅(qū)動數(shù)據(jù)集”，時間段是1979—2018年，空間分辨率為0.1°。此數(shù)據(jù)集是以國際上現(xiàn)有的Princeton 再分析資料、全球陸面數(shù)據(jù)同化系統(tǒng)（Global Land Data Assimilation System，GLDAS）資料、Global Energy and Water Cycle Experiment-Surface Radiation Budget 輻射資料，以及熱帶降雨測量任務（Tropical Rainfall Measuring Mission，TRMM）降水資料為背景場，融合了中國氣象局常規(guī)氣象觀測數(shù)據(jù)制作而成［42］。原始資料來自于氣象局觀測數(shù)據(jù)、再分析資料和衛(wèi)星遙感數(shù)據(jù)。已去除非物理范圍的值，采用ANU-Spline（Australian National University 利用Fortran 開發(fā)的空間異相關(guān)插值模型，其利用薄盤平滑樣條函數(shù)對多變量數(shù)據(jù)進行分析和插值的工具）統(tǒng)計插值，精度好于國際上已有再分析數(shù)據(jù)的精度［43］。本研究根據(jù)此數(shù)據(jù)集提取并整理青藏高原內(nèi)流盆地的氣溫和降水數(shù)據(jù)。

（3）青藏高原月平均地表蒸散發(fā)數(shù)據(jù)集

經(jīng)過對湖泊面積進行突變點檢測，其突變時間在2000 年之后。使用來源于國家青藏高原科學數(shù)據(jù)中心的“青藏高原月平均地表蒸散發(fā)數(shù)據(jù)集”，可滿足構(gòu)建模型的數(shù)據(jù)需要，時間段是2001—2018年，空間分辨率為0.1°，數(shù)據(jù)存放格式為標準的NETCDF 格式，蒸散發(fā)量通過了青藏高原6 個湍流通量站的觀測數(shù)據(jù)的驗證［44］。此數(shù)據(jù)集主要以衛(wèi)星遙感數(shù)據(jù)（中分辨率成像光譜儀，Moderate-resolution Imaging Spectroradiometer）和再分析氣象數(shù)據(jù)（China Meteorological Forcing Dataset，CMFD）作為輸入，利用地表能量平衡系統(tǒng)模型（Surface Energy Balance System）計算得到［44］。

以內(nèi)流湖泊所處的湖泊盆地的矢量范圍裁剪氣象要素，對各盆地內(nèi)的降水、氣溫數(shù)據(jù)進行年度統(tǒng)計與分析，得到所對應的內(nèi)流盆地1979—2018 年各年份的年度降水、氣溫數(shù)據(jù)。采用同樣的方法得到各湖泊盆地2001—2018年各年份的年度蒸發(fā)數(shù)據(jù)。

（4）數(shù)字高程數(shù)據(jù)（DEM，Digital Elevation Model）

DEM 分別是來源于National Oceanic and Atmospheric Administration（NOAA）地形數(shù)據(jù)平臺的全球陸地海洋地形數(shù)據(jù)Earth Topography 1-minute（ETOPO1）DEM 和來源于Alaska Satellite Facility（ASF）EARTHDATA 平臺的Advanced Land Observing Satellite（ALOS）衛(wèi)星PALSAR 波段的高精度（12.5 m 分辨率）的DEM 數(shù)據(jù)。Earth Topography 1-minute ETOPO1 DEM 數(shù)據(jù)由National Geophysical Data Center（NGDC）美國地球物理中心發(fā)布，它是高程數(shù)據(jù)同時還包括海洋海底地形數(shù)據(jù)，ETOPO1 的分辨率為1′，模型校正時間為2011年［45］。ALOS（PALSAR）DEM 是日本的對地觀測衛(wèi)星獲取的數(shù)據(jù)，PALSAR 是ALOS 衛(wèi)星攜帶的一個L 波段的合成孔徑雷達傳感器，不受云層、天氣和晝夜影響，可全天候?qū)Φ赜^測，ALOS-PALSAR DEM兼顧高分辨率、直觀精細地形、受地表覆蓋影響小等多重特點，已被應用于地形分析等領域［46］。本研究中提取了青藏高原地區(qū)的兩種DEM數(shù)據(jù)。

（5）其他矢量數(shù)據(jù)

矢量數(shù)據(jù)包括全球內(nèi)陸盆地矢量范圍及內(nèi)外流盆地連接點，亞洲12 級流域盆地具體劃分，來源于美國地質(zhì)調(diào)查局官網(wǎng)https：//www. usgs. gov/。HydroBASINS 數(shù)據(jù)集對于青藏高原內(nèi)外流盆地劃分基本符合實際。

（6）Landsat 8衛(wèi)星遙感數(shù)據(jù)

Landsat 8 是美國陸地衛(wèi)星計劃（Landsat）的第八顆衛(wèi)星，于2013年2月11日在加利福尼亞范登堡空軍基地由Atlas-V 火箭搭載發(fā)射成功，最初稱為“陸地衛(wèi)星數(shù)據(jù)連續(xù)性任務”。Landsat 8上攜帶陸地成像儀（Operational Land Imager，OLI）和熱紅外傳感器（Thermal Infrared Sensor，TIRS），OLI 陸地成像儀包括9 個波段，空間分辨率為30 m，可在Google Earth Engine 平臺直接獲取應用?？芍苯永肎EE平臺“LANDSAT/LC08/C02/T1_RT_TOA”數(shù)據(jù)的夏季（6月）影像，提取所需各個時間段的水面數(shù)據(jù)。

2 研究方法

首先，利用ETOPO1 DEM 數(shù)據(jù)進行ArcGIS 水文徑流分析，獲取徑流線與內(nèi)外流盆地交界線的交點，并與當前全球內(nèi)流盆地溢出點進行對比驗證，即得到青藏高原地區(qū)內(nèi)流盆地湖泊潛在外溢點；其次，獲取內(nèi)外流盆地連接點所在的內(nèi)流盆地中的主要湖泊水體面積與氣象數(shù)據(jù)（年平均降水、年平均氣溫、年平均蒸發(fā)），并通過Pettitt 突變點檢驗、M-K趨勢分析及Pearson相關(guān)性分析等方法，分析近幾十年來青藏高原內(nèi)流盆地湖泊水體面積隨氣象數(shù)據(jù)的變化情況，建立并驗證湖泊水面面積-氣象參數(shù)線性擬合模型；最后，通過預測氣象數(shù)據(jù)的變化從而得到水體面積未來變化的可能，結(jié)合湖泊曾經(jīng)達到的最高水位，分析內(nèi)流湖泊一定時期內(nèi)發(fā)生外溢并與外流盆地重建水力聯(lián)系的可能性（圖2）。

圖3 青海湖2018年永久水和季節(jié)水水面分布Fig. 3 Permanent and seasonal surface water of Qinghai Lake in 2018

2.1 內(nèi)流盆地湖泊潛在外溢點檢測

利用ArcGIS 水文分析工具和DEM 數(shù)據(jù)，可以實現(xiàn)河流河網(wǎng)水系、出水口以及流域的提取。根據(jù)青藏高原DEM 數(shù)據(jù)徑流分析結(jié)果，設定閾值為10提取徑流河網(wǎng)，將得到的結(jié)果與HydroBASINS數(shù)據(jù)集的內(nèi)外流域分類結(jié)果進行空間疊加分析。Hydro-BASINS 數(shù)據(jù)集的圖層描述了全球范圍內(nèi)的流域邊界和子流域劃分。根據(jù)“Endo”字段劃分內(nèi)流區(qū)和外流區(qū)，這個字段共有3個值，其中，0表示不是內(nèi)流區(qū)的部分，1 表示內(nèi)流區(qū)的一部分，2 表示內(nèi)流區(qū)下游聚集部分。因此，當“Endo”字段值為0時，其區(qū)域為外流區(qū)，河流水系與外部海洋相連；當“Endo”字段值為1 或2 時，表示其區(qū)域為內(nèi)流區(qū)。從而獲取徑流河網(wǎng)與內(nèi)外流盆地交界線的交點［18］，并與當前全球內(nèi)流盆地溢出點進行對比驗證，即得到青藏高原地區(qū)內(nèi)流盆地湖泊潛在外溢點空間位置。

2.2 基于JRC數(shù)據(jù)的湖泊水面數(shù)據(jù)提取

湖泊面積通過JRC 全球水面分布數(shù)據(jù)集提取。使用JRC 數(shù)據(jù)集的“transitions”波段確定水體面積增大的湖泊。使用數(shù)據(jù)集中的“waterClass”波段提取各湖泊各年份永久水的年度面積。其中，紅色邊界線是青海湖歷史最大水痕，深藍色表示青海湖2018年永久水，淺藍色表示青海湖2018年季節(jié)水。

2.3 湖泊水面面積和湖泊流域內(nèi)氣象要素變化分析

2.3.1 Pettitt突變點檢驗

對時間序列湖泊年平均面積數(shù)據(jù)進行Pettitt 突變點檢驗［47］，獲取湖泊面積變化的突變年份信息。Pettitt檢驗對于序列中部的突變點更加敏感。Pettitt突變點檢驗法是一種非參數(shù)檢驗方法，直接利用秩序列來檢測突變點。對于包含N個樣本的一時間序列數(shù)據(jù)X，計算其統(tǒng)計量Uk：

式中：

式中：sgn是符號函數(shù)，定義如下：

取Ut中絕對值最大的值Kt，該點為最顯著的突變點，并計算Kt所對應的統(tǒng)計量P，如果P小于給定的顯著性水平（如P=0.05），表示存在統(tǒng)計顯著的突變點。

2.3.2 Mann-Kendall趨勢分析

在湖泊面積時間序列突變點檢測分析結(jié)果的基礎之上，對湖泊面積突變年份前后的時序湖泊面積數(shù)據(jù)進行M-K 趨勢分析［20］，獲得單個湖泊在不同時間段的面積變化趨勢特征。Mann-Kendall 法是一種非參數(shù)統(tǒng)計檢驗方法，但是不適用于檢測有多個突變點的序列。當Kendall系數(shù)為正，則表示湖泊面積為增加趨勢，系數(shù)越大，面積持續(xù)增加的可能性越大。

對于時間序列X，Mann-Kendall 趨勢檢驗的統(tǒng)計量為：

式中：xj為時間序列的第j個數(shù)據(jù)值；n為數(shù)據(jù)樣本的長度；sgn是符號函數(shù)，其定義如下：

Mann（1945 年）和Kendall（1975 年）證明，當n≥8 時，統(tǒng)計量S大致地服從正態(tài)分布，其均值為0，方差為：

式中：ti是第i組的數(shù)據(jù)點的數(shù)目。

標準化統(tǒng)計量，按照如下公式計算：

式中：Zc服從標準正態(tài)分布。

衡量趨勢大小的指標為：

式中：1

2.3.3 Pearson相關(guān)性分析

采用Pearson 相關(guān)性分析方法［48］對湖泊面積與同時期的氣象數(shù)據(jù)進行分析，判斷分析面積與氣象參量之間有強相關(guān)關(guān)系的湖泊。Pearson 相關(guān)系數(shù)衡量的是線性相關(guān)關(guān)系。相關(guān)系數(shù)越接近于1 或-1，相關(guān)度越強，相關(guān)系數(shù)越接近于0，相關(guān)度越弱。兩個變量之間的皮爾遜相關(guān)系數(shù)定義為兩個變量之間的協(xié)方差和標準差的商：

式（10）定義了總體相關(guān)系數(shù)，常用希臘小寫字母ρ作為代表符號。估算樣本的協(xié)方差與標準差，可得到Pearson相關(guān)系數(shù)，常用英文小寫字母r代表：

式中：r亦可由（Xi，Yi）樣本點的標準分數(shù)均值估計，得到與上式等價的表達式：

式中：（Xi-）/σX、及σX分別是對Xi樣本的標準分數(shù)、樣本平均值和樣本標準差。

計算出相關(guān)系數(shù)r之后，需要檢驗它是否具有統(tǒng)計學意義，即是否顯著，其計算公式為：

在t分布中找到對應的P值。若P值小于設定的顯著性水平，則拒絕零假設，認為X與Y之間存在顯著的線性相關(guān)；反之，則接受零假設，認為X與Y之間不存在顯著的線性相關(guān)。

2.4 湖泊水面面積-湖泊流域氣象參數(shù)關(guān)系模型構(gòu)建與精度驗證

2.4.1 模型構(gòu)建

氣候因素引起湖泊擴張是通過水量變化和流域水量平衡才能精確確定溢出可能性，但由于青藏高原地區(qū)的數(shù)據(jù)資料具有局限性，湖泊水量變化的計算需要以湖泊面積或水位為輸入進行計算，因此本研究直接使用了水面面積作為模型關(guān)系參量，構(gòu)建湖泊水面面積-氣象參數(shù)關(guān)系模型。

使用JRC 數(shù)據(jù)源的數(shù)據(jù)逐年提取湖泊面積，對所研究的湖泊水面面積與其所在湖泊流域內(nèi)的氣象數(shù)據(jù)（降水、氣溫、蒸發(fā)）進行相關(guān)性分析，根據(jù)相關(guān)性篩選影響湖泊水面面積變化的主要氣象因子，并將湖泊水面面積作為因變量，氣象因子作為自變量進行相關(guān)性關(guān)系模型擬合，得到湖泊水面面積-氣象參數(shù)關(guān)系模型。

根據(jù)歷史數(shù)據(jù)變化規(guī)律對該氣象因子進行時間序列預測，以氣象因子為因變量，時間為自變量進行關(guān)系擬合，得到氣象因子時間序列的預測模型。并根據(jù)此模型預測：隨著時間的變化，氣象因子相應數(shù)值的變化情況。再將其作為自變量帶入水面面積-氣象參數(shù)關(guān)系擬合模型中，預測因變量水面面積的變化情況。

將湖泊水面面積-氣象參數(shù)關(guān)系模型與氣象因子時間序列預測模型進行代入整合，使得氣象因子成為中間變量，從而得到面積-氣象因子-時間的擬合模型，可在此基礎上對湖泊未來水面面積變化推斷。

2.4.2 精度驗證

為了驗證湖泊水面面積-氣象參數(shù)關(guān)系擬合模型的可靠性，分別采用兩種方法進行精度驗證。方法一：根據(jù)擬合模型計算2019 年、2020 年的水面面積擬合數(shù)據(jù)，利用Landsat 8 影像數(shù)據(jù)提取2019 年、2020 年的湖泊面積數(shù)據(jù)，將計算得到的水面面積擬合數(shù)據(jù)與使用Landsat 8 衛(wèi)星提取的湖泊實際水面面積結(jié)果進行對比驗證，得到實際面積與模型擬合面積結(jié)果的誤差；方法二：利用突變年份2001—2015 年的歷史數(shù)據(jù)作為訓練集構(gòu)建擬合模型，將2016—2018 這三年的歷史數(shù)據(jù)作為驗證集進行對比驗證，將根據(jù)模型擬合得到的后一部分歷史數(shù)據(jù)的結(jié)果值與實際的歷史數(shù)據(jù)進行誤差對比分析，進行精度驗證。

2.4.3 基于Landsat影像的湖泊水面提取

選取Landsat 8 遙感影像（夏季影像）作為數(shù)據(jù)源，對遙感影像進行去云處理，結(jié)合歸一化水體指數(shù)（normalized difference water index，NDWI）［NDWI=（B（Green）-B（NIR））/（B（Green）+B（NIR））］設置水體掩膜。歸一化水體指數(shù)是根據(jù)水體光譜反射特性，基于近紅外波段與綠波段建立的歸一化比值指數(shù)。理論上，NDWI>0 表示地面有水或冰雪覆蓋，NDWI=0 表示地面有巖石或裸土等覆蓋，NDWI<0 表示地面有植被覆蓋。但在實際情況中，由于受到地物復雜性及噪聲等條件的干擾，區(qū)分水體與非水體的閾值往往不為0［49］。因此將計算出的NDWI 指數(shù)閾值設為0.005，即保留NDWI 大于0.005 的區(qū)域?qū)ρ芯繀^(qū)的水體進行提取。將得到的水體結(jié)果疊加在影像上，與事實基本一致。同時，將基于Landsat 影像的湖泊水面提取結(jié)果與基于JRC 數(shù)據(jù)的湖泊水面數(shù)據(jù)提取結(jié)果進行對比，二者基本一致。

2.5 湖泊外溢可能性量化分析

當湖泊面積持續(xù)增大，呈現(xiàn)擴張趨勢，可以使用湖泊歷史最高水位線來判斷湖泊是否會再次擴張到此位置，導致水位繼續(xù)升高，從而在湖泊流域邊界溢出、與外流區(qū)的河流相連的參考依據(jù)。結(jié)合各湖泊曾經(jīng)可能到達的歷史最大水位高度［50］和DEM數(shù)據(jù)，根據(jù)湖泊所在的湖泊流域等值線計算湖泊再次達到歷史最高水位時對應的理論面積。以理論最大面積為輸入值，利用湖泊水面面積-氣象參數(shù)擬合到模型中，計算得到理論溢出面積時對應的主要氣象因子的數(shù)值，并根據(jù)氣象因子時間序列的預測模型，計算得到湖泊達到理論溢出面積時對應的可能溢出時間，判斷短期內(nèi)是否有溢出的可能（圖4）。

圖4 湖泊外溢可能性分析流程Fig. 4 Possibility analysis process of lake overflow

3 結(jié)果分析

3.1 青藏高原內(nèi)外流盆地連接點分布

根據(jù)青藏高原DEM 數(shù)據(jù)徑流分析結(jié)果表明，本研究提取的徑流河網(wǎng)與內(nèi)外流矢量邊界線參考數(shù)據(jù)之間的連接點有26 個，其中，國內(nèi)青海有11個，西藏有12個；境外查謨-克什米爾有3個。這26個連接點位于國內(nèi)外五大流域中，如圖5所示。

圖5 青藏高原內(nèi)外流連接點Fig. 5 Connection point of internal and external flow in Qinghai-Tibet Plateau

3.2 典型湖泊面積變化趨勢研究

在26 個湖泊流域中，有23 個湖泊2018 年的水面面積相對于1984年的水面面積呈增加趨勢，其中有17 個在1984—2018 年間一直呈增加趨勢。一直增加是指1984—2018 年面積變化趨勢通過顯著性檢驗，這個時間段整體呈現(xiàn)增加趨勢，并非每一年都增加。17個一直增加的湖泊如表1所示。設定假設性檢驗閾值為0.05，對這17個湖泊的各年水面面積分別進行Pettitt檢驗以及M-K趨勢檢驗。結(jié)果表明，各湖泊在1987—2009年之間均有發(fā)生突變性變化，其中，有10 個湖泊的突變年份在2000—2003 年之間，且有6 個湖泊的面積在2001 年發(fā)生突變（表1）。根據(jù)M-K 趨勢檢驗結(jié)果，Kendall 系數(shù)值越大，增加的趨勢越明顯，設定Kendall 系數(shù)值為0.7，篩選出大于0.7的10個湖泊（表2）作進一步研究。

表1 17個典型湖泊的突變年份及增加趨勢Table 1 The 17 lakes with abrupt changes during the priod of 1987 to 2009

表2 10個湖泊水面面積與氣象參數(shù)相關(guān)性分析結(jié)果Table 2 Correlation analysis results between water body area and meteorological data of 10 lakes

對突變年份之后湖泊面積上升明顯的前10 個湖泊的水面面積與氣象參數(shù)數(shù)據(jù)進行相關(guān)性分析，設定顯著性水平閾值為0.05，當P值小于0.05，說明顯著相關(guān)，即湖泊水面面積與氣象參數(shù)之間有明顯相關(guān)性，反之，說明湖泊水面面積與氣象參數(shù)之間沒有明顯相關(guān)性（表2）。并對湖泊水面面積與不同氣象因子分別進行相關(guān)性分析，分析結(jié)果表明鹽湖、班公錯、澤錯、雀莫錯4 個湖泊的線性關(guān)系檢驗分析結(jié)果的相關(guān)系數(shù)R2大于0.5 且P值小于0.05，說明這4 個湖泊水面面積與氣象參數(shù)之間顯著相關(guān)，即這些湖泊水面面積的變化受到相應氣象因素的影響明顯。因此，進一步以鹽湖、班公錯、澤錯和雀莫錯為典型湖泊來研究氣象要素變化對湖泊水面面積變化的影響。

3.3 典型湖泊面積與關(guān)鍵氣象因子模型構(gòu)建

根據(jù)鹽湖、班公錯、澤錯和雀莫錯的水面面積和對應主要氣象因子數(shù)據(jù)，選擇突變年份之后至2018 年的數(shù)據(jù)構(gòu)建關(guān)系擬合湖泊水面面積-氣象參數(shù)關(guān)系模型，并對氣象因子進行時間推演，構(gòu)建氣象因子時間序列預測模型（表3）。將影響湖泊水面面積變化的主要氣象因子結(jié)合遙感圖像查看，結(jié)果表明，計算得到的各個湖泊受到的主要影響氣象因子與其所處地理位置的特點相符合。

表3 典型湖泊面積-氣象參數(shù)擬合模型Table 3 Fitting model of typical lake area meteorological parameters

位于青海省的鹽湖，由于2011 年卓乃湖潰決，將卓乃湖、庫賽湖、海丁諾爾湖、鹽湖進行連接［3］，之后鹽湖面積相較于之前驟增，因此主要研究時間段為潰決后的2012—2018 年。模型結(jié)果顯示湖泊與降水、氣溫的相關(guān)性較大，與鹽湖地勢較高，易受氣溫導致的冰川融水影響，以及連接其他湖泊受降水影響大的地理位置事實相符。位于西藏自治區(qū)較高海拔的班公錯和澤錯，位于青藏高原的西部，地勢高，周圍有明顯的雪山分布，湖泊以冰雪融水為主要補給來源，受氣溫影響比較明顯。位于青海省的雀莫錯，湖泊面積變化主要與氣溫和降水有關(guān)，該區(qū)域地勢較平緩，降水和氣溫數(shù)據(jù)起伏較大，總計年凈降水量為正值，湖泊的面積處于緩慢增加的趨勢。

3.4 模型精度驗證

為了驗證上述相關(guān)性模型的可靠性，采用兩種方法進行驗證。

（1）根據(jù)Landsat 8 影像提取的湖泊面積結(jié)果與模型擬合的結(jié)果進行對比，結(jié)果表明鹽湖、班公錯、澤錯、雀莫錯的實際面積與模型擬合面積結(jié)果的誤差范圍在6%以內(nèi)（表4）；并針對2019 年提取各個湖泊的JRC 水體面積，與同年Landsat 8 提取的水體面積進行對比，驗證Landsat 8 提取水體的精度（表5）。結(jié)果表明：Landsat 8 提取的面積誤差精度在5%以內(nèi)，此種方法提取水體的精度滿足要求。

表4 方法一：精度驗證結(jié)果Table 4 Method 1： accuracy verification results

表5 驗證Landsat 8 提取水體的精度Table 5 Verify the accuracy of water body extracted by Landsat 8

（2）將2000—2015 年數(shù)據(jù)作為訓練集擬合模型，并分別根據(jù)此模型計算班公錯、澤錯、雀莫錯的2016 年、2017 年、2018 年水體面積，與2016—2018年JRC 湖泊面積進行對比，結(jié)果表明，誤差百分比都在3%以內(nèi)（表6）。兩種精度驗證結(jié)果說明，構(gòu)建的4個湖泊的關(guān)系模型可靠性高。

表6 方法二：精度驗證結(jié)果Table 6 Method 2： accuracy verification results

3.5 湖泊滿溢及與外流水系重建水力聯(lián)系的可能性

以鹽湖、班公錯、澤錯、雀莫錯4 個湖泊理論溢出面積為輸入，根據(jù)各湖泊水面面積與主要氣象參數(shù)之間的相關(guān)性方程，可計算得到溢出位置對應的氣象參數(shù)值，結(jié)合各氣象因子隨時間變化的數(shù)值關(guān)系方程（表3），推算得到各湖泊達到歷史最高水位時的理論溢出面積的可能時間點。分析結(jié)果表明，除已于2019年外溢的鹽湖外，班公錯、澤錯、雀莫錯在近30年內(nèi)無溢出風險（表7）。

表7 典型湖泊溢出可能性分析Table 7 Possibility analysis of typical lake overflow

根據(jù)Google Earth 影像、DEM 數(shù)據(jù)以及河流水系分布可知，鹽湖在2019 年9 月溢出，沿著提前修建的引水河道流出。湖水從鹽湖流出，沿著東南方向的河道流入清水河，清水河整體走勢呈西北-東南方向，因此匯入清水河后繼續(xù)沿著東南方向的河道流入楚瑪爾河，楚瑪爾河作為長江源的北源，最終匯入東海。班公錯湖泊呈東南-西北走向，在地勢上東南高西北低，河流自東南向西北流，最終匯流到班公錯的西北角的河谷，一旦水流高過西北角出水河谷的最高高度，河流就會越過出水河谷，沿著歷史河道向西北方向流入什約克河，最終隨著什約克河注入印度河。當位于班公錯北部的澤錯湖泊水面升高至理論最大高度時，其從西南方向的出水口沿著河道最終向南匯入班公錯，最終跟隨班公錯的出水河流注入印度河。雀莫錯位于青海省，當湖泊水面面積增加至湖泊邊界溢出時，湖泊中的水可能會從湖泊邊界最低點沿著河道向西流入沱沱河，沱沱河作為長江的源頭，最終湖泊中的水會隨著沱沱河匯入長江、注入東海，如圖6所示。

圖6 典型湖泊的可能溢出河道Fig. 6 Possible overflow channels of typical lakes

4 結(jié)論與討論

本文使用GIS 方法對研究區(qū)的DEM 數(shù)據(jù)進行了水文徑流分析，獲取了徑流河網(wǎng)，以湖泊水面面積為因變量、氣象因子為自變量，構(gòu)建了時間序列上的線性擬合模型，分析了青藏高原26個具有潛在溢出可能性的內(nèi)流湖泊水面面積及其與所在湖泊流域主要氣象要素的變化趨勢及相關(guān)性特征，得到主要結(jié)論如下：

（1）青藏高原內(nèi)流盆地的湖泊水體整體擴張趨勢顯著。在26 個湖泊流域中，有23 個湖泊2018 年的水面面積相對1984年的水面面積呈增加趨勢，其中有17個在1984—2018年間一直呈增加趨勢。

（2）湖泊面積變化與氣候要素之間具有顯著的區(qū)域相關(guān)性。青藏高原東部青海省的鹽湖和雀莫錯的湖泊水面面積變化主要與氣溫和降水有關(guān)；而位于西藏自治區(qū)較高海拔的班公錯和澤錯的水面面積變化，主要與氣溫的變化有關(guān)。

（3）基于本文研究使用的數(shù)據(jù)，按照當前氣象數(shù)據(jù)的時間序列預測趨勢變化情況發(fā)展，根據(jù)短期內(nèi)四個典型湖泊水面面積隨氣象數(shù)據(jù)變化的線性擬合模型預測，在無突發(fā)事件的條件下，班公錯、澤錯、雀莫錯短期內(nèi)（30年內(nèi)）不會有外溢風險。

本文研究湖泊水面面積變化明顯、且與氣象因子具有顯著相關(guān)性的湖泊，構(gòu)建湖泊水面面積-氣象因子線性擬合模型，并對模型精度進行檢驗，結(jié)果顯示鹽湖的溢出時間與事實基本一致，而班公錯、澤錯和雀莫錯短期內(nèi)不會溢出，結(jié)果具有一定的可靠性，但同時也存在一些不足：

（1）主要選用了2000 年之后數(shù)據(jù)進行研究，基于趨勢分析檢測結(jié)果，相關(guān)氣象參數(shù)呈現(xiàn)的變化趨勢通過了假設檢驗，研究時間段較短，基于這個假設檢驗結(jié)果構(gòu)建模型并進行外推，推測研究湖泊未來30 年變化；并且本實驗所使用的數(shù)據(jù)集均來源于第三方數(shù)據(jù)，因此在數(shù)據(jù)的空間分辨率、不同模型得到的數(shù)據(jù)結(jié)果等數(shù)據(jù)精度方面存在一定的不確定性。

（2）僅考慮湖泊永久水體變化所致的溢出情況，未考慮極端降水、地質(zhì)災害（地震、滑坡和泥石流等）的影響，如果氣候系統(tǒng)再次出現(xiàn)極端變化，那么這種方法獲取的結(jié)果會不準確。

（3）研究結(jié)果是根據(jù)當前氣候變化趨勢通過簡單線性模型進行預測，未考慮冰期、間冰期、氣候周期變化、人為干預氣候、模型復雜性等方面的影響，因此對于長期準確預測有一定的局限性。

上述不足將在后續(xù)的研究中進一步完善，重點包括：

（1）選用長期（如近100 年）的時間序列進行分析，根據(jù)氣象站點數(shù)據(jù)以及多個氣象數(shù)據(jù)集得到更精準貼近實際情況的數(shù)據(jù)。

（2）增加地表徑流、地下水、冰川數(shù)據(jù)等指標，對不同補給來源的湖泊受到氣候因素的響應情況進行劃分。

（3）建立多要素的面積-氣象擬合模型，添加突發(fā)因素衡量指數(shù)，表明氣候周期變化，提高模型長期預測的準確性。