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

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

0 引言

全球氣候驅(qū)動影響下的冰川融化［1］、水文循環(huán)加速等問題［2-3］，使得全球很多湖泊呈擴(kuò)張趨勢，這一趨勢的發(fā)展可能會導(dǎo)致內(nèi)流水系與外流水系重組等一系列重大變化［4-5］，進(jìn)而引發(fā)湖泊潰堤、突發(fā)洪水等地質(zhì)危害。例如，氣候變化引起的冰川流失會減少漫灘和近海沉積［6-7］；河流流向的轉(zhuǎn)換會引發(fā)流域盆地的突然重組［8］；湖泊溢流會引起沖溝的發(fā)育、產(chǎn)生溯源侵蝕，使決口不斷擴(kuò)大，最終導(dǎo)致潰壩［9］。當(dāng)前全球范圍內(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ā)生溢流而導(dǎo)致峽谷切割［16］；秘魯Cordillera Blanca 山谷冰川湖暴發(fā)洪水導(dǎo)致谷底沉積物被搬運(yùn)遷移［17］。

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

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

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

為了探究青藏高原26 個有潛在溢出可能的湖泊的外溢潛在風(fēng)險(xiǎn)［18］，本文以衛(wèi)星遙感和氣象再分析數(shù)據(jù)為數(shù)據(jù)源，采用地統(tǒng)計(jì)分析方法，分析了這些湖泊的面積變化，以及與相應(yīng)湖泊流域內(nèi)降水、氣溫和蒸發(fā)之間的相關(guān)性特征；在此基礎(chǔ)上，通過情景假設(shè)，分析了未來區(qū)域降水、氣溫等參量變化對湖泊水體面積的影響，并對這些湖泊未來水位上漲滿溢的可能性開展了推測與判斷，以期為區(qū)域未來氣候水文變化規(guī)律研究及致災(zāi)風(fēng)險(xiǎn)應(yīng)對策略的制定提供數(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，地形上可分為羌塘高原、藏南谷地、柴達(dá)木盆地、祁連山地、青海高原和川藏高山峽谷區(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ū)正在變暖變干，同時高原整體風(fēng)速正在減小。升溫主要是夜間的最低溫度貢獻(xiàn)的，不同地區(qū)升溫速率有差異，中部地區(qū)高于東部地區(qū)。青藏高原平均溫度、最低和最高溫度的分布與高原海拔走勢一致。降水量空間分布上表現(xiàn)出從東南向西北逐級減少的分布特點(diǎn)。風(fēng)速自東南向西北呈現(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)計(jì)數(shù)據(jù)［41］。本研究中提取了青藏高原內(nèi)流盆地中主要湖泊1985—2020 年的永久水（一年四季都有水）和季節(jié)水（包括冰凍期，一年中在水下的時間不到12 個月），并對異常數(shù)據(jù)進(jìn)行鄰近年份數(shù)據(jù)插值處理，從而獲取研究區(qū)近幾十年來的水體面積變化完整數(shù)據(jù)。

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

氣溫和降水?dāng)?shù)據(jù)來源于國家青藏高原科學(xué)數(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 輻射資料，以及熱帶降雨測量任務(wù)（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ù)進(jìn)行分析和插值的工具）統(tǒng)計(jì)插值，精度好于國際上已有再分析數(shù)據(jù)的精度［43］。本研究根據(jù)此數(shù)據(jù)集提取并整理青藏高原內(nèi)流盆地的氣溫和降水?dāng)?shù)據(jù)。

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

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

以內(nèi)流湖泊所處的湖泊盆地的矢量范圍裁剪氣象要素，對各盆地內(nèi)的降水、氣溫?cái)?shù)據(jù)進(jìn)行年度統(tǒng)計(jì)與分析，得到所對應(yīng)的內(nèi)流盆地1979—2018 年各年份的年度降水、氣溫?cái)?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 波段的合成孔徑雷達(dá)傳感器，不受云層、天氣和晝夜影響，可全天候?qū)Φ赜^測，ALOS-PALSAR DEM兼顧高分辨率、直觀精細(xì)地形、受地表覆蓋影響小等多重特點(diǎn)，已被應(yīng)用于地形分析等領(lǐng)域［46］。本研究中提取了青藏高原地區(qū)的兩種DEM數(shù)據(jù)。

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

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

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

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

2 研究方法

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

利用ArcGIS 水文分析工具和DEM 數(shù)據(jù)，可以實(shí)現(xiàn)河流河網(wǎng)水系、出水口以及流域的提取。根據(jù)青藏高原DEM 數(shù)據(jù)徑流分析結(jié)果，設(shè)定閾值為10提取徑流河網(wǎng)，將得到的結(jié)果與HydroBASINS數(shù)據(jù)集的內(nèi)外流域分類結(jié)果進(jìn)行空間疊加分析。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ū)下游聚集部分。因此，當(dāng)“Endo”字段值為0時，其區(qū)域?yàn)橥饬鲄^(qū)，河流水系與外部海洋相連；當(dāng)“Endo”字段值為1 或2 時，表示其區(qū)域?yàn)閮?nèi)流區(qū)。從而獲取徑流河網(wǎng)與內(nèi)外流盆地交界線的交點(diǎn)［18］，并與當(dāng)前全球內(nèi)流盆地溢出點(diǎn)進(jìn)行對比驗(yàn)證，即得到青藏高原地區(qū)內(nèi)流盆地湖泊潛在外溢點(diǎn)空間位置。

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

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

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

2.3.1 Pettitt突變點(diǎn)檢驗(yàn)

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

式中：

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

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

2.3.2 Mann-Kendall趨勢分析

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

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

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

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

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

標(biāo)準(zhǔn)化統(tǒng)計(jì)量，按照如下公式計(jì)算：

式中：Zc服從標(biāo)準(zhǔn)正態(tài)分布。

衡量趨勢大小的指標(biāo)為：

式中：1

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

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

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

式中：r亦可由（Xi，Yi）樣本點(diǎn)的標(biāo)準(zhǔn)分?jǐn)?shù)均值估計(jì)，得到與上式等價的表達(dá)式：

式中：（Xi-）/σX、及σX分別是對Xi樣本的標(biāo)準(zhǔn)分?jǐn)?shù)、樣本平均值和樣本標(biāo)準(zhǔn)差。

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

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

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

2.4.1 模型構(gòu)建

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

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

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

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

2.4.2 精度驗(yàn)證

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

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

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

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

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

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

3 結(jié)果分析

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

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

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

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

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

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

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

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

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

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

3.4 模型精度驗(yàn)證

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論與討論

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

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

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

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

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

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

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

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

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

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

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

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