氣候變化對可可西里鹽湖流域湖泊水量變化的影響分析

李 琳，譚德寶，文雄飛，王 瑩，劉希勝，王 崗

(1.長江科學(xué)院 空間信息技術(shù)應(yīng)用研究所，武漢 430010； 2.武漢市智慧流域工程技術(shù)研究中心，武漢 430010；3.長江水利委員會長江科學(xué)院，武漢 430010； 4.青海省水文水資源測報中心，西寧 810001)

1 研究背景

湖泊在地球水體中僅占0.4%，卻在自然界的水循環(huán)中發(fā)揮著重要的作用[1]。湖泊是氣候環(huán)境變化的指示器，對氣候變化反應(yīng)迅速，是連接地球表層系統(tǒng)冰凍圈、大氣圈、水圈和生物圈的紐帶，湖泊面積、水位及水量變化是其所在流域水量平衡綜合作用的結(jié)果[2]。青藏高原分布著大量的湖泊，這些湖泊較少受到人類活動的影響，對氣候暖濕化感應(yīng)強烈，可以反映自然狀態(tài)下的區(qū)域氣候變化規(guī)律，是氣候變化的前哨[3]。近20多年來，湖泊水量的變化尤為顯著，其原因大多為氣候變化的影響[4]。

鹽湖流域湖泊位于青藏高原腹地高寒地帶，受氣候變化影響較大。杜玉娥等[5]認(rèn)為2011年前氣候暖濕化是鹽湖面積緩慢擴大的主要原因，冰川和凍土融水可能是次要原因。姚曉軍等[6]認(rèn)為可可西里地區(qū)湖泊動態(tài)變化主要與氣候變化有關(guān)，降雨增多、蒸發(fā)減少是湖泊擴大的主要原因，而氣候變暖引起的冰川融水增加、凍土水分釋放是次要原因。Qiao等[7]利用1976—2019年多時相Landsat遙感影像與地形數(shù)據(jù)估算了可可西里湖泊蓄水量的變化，利用五道梁氣象站數(shù)據(jù)分析多年平均、最大最小氣溫和多年降雨量的變化趨勢，將氣溫和降雨的變化相比較，認(rèn)為降雨增加可能主導(dǎo)湖泊的擴張。陳強等[8]統(tǒng)計分析了1961—2018年可可西里主要氣候要素變化趨勢，并從水量平衡角度探討了水文情勢變化后的鹽湖水體補排、泄、蓄關(guān)系。

以往的研究主要是基于遙感衛(wèi)星等技術(shù)手段研究不同氣候條件下湖泊面積的變化，但是由于湖泊周邊地形的差異，面積變化不能充分反映湖泊對氣候變化的響應(yīng)[9]。研究湖泊水量變化需要考慮流域及周邊相連的水系，常常借助數(shù)學(xué)模型解析湖泊水量變化的影響因素[4]。VIC(Variable Infiltration Capacity)模型是基于SVATS(Soil Vegetation Atmospheric Transfer Schemes)機制的大尺度分布式水文模型。VIC模型以混合產(chǎn)流作為產(chǎn)流機制，在模擬產(chǎn)匯流過程中考慮了融雪和凍土因素的影響，并在徑流模擬方面有很多研究與進展，適用于鹽湖流域湖泊水量模擬與分析[10]。

由于實測數(shù)據(jù)難得，對于鹽湖流域水文模型模擬研究并不多，且之前未將湖泊水量變化量如冰川融水、凍土和地下水補給量化。因此本文結(jié)合多源遙感數(shù)據(jù)和氣象站數(shù)據(jù)，在利用統(tǒng)計方法計算長時間氣象要素與湖泊面積關(guān)系的基礎(chǔ)上，通過VIC模型模擬的水量和利用冰川數(shù)據(jù)集[11-12]估算的冰川融水，得到2015—2018年鹽湖流域湖面降雨、陸面凈補給、凍土釋水和地下水、冰川融水徑流組成與占比。探討氣候變化對湖泊水量變化的影響機制，辨識影響湖泊水量的主導(dǎo)因素，為青藏高原的水資源調(diào)查及合理利用和開發(fā)提供科學(xué)依據(jù)和理論依據(jù)[13]。

2 研究區(qū)域概況

鹽湖流域位于青海省玉樹州治多縣西部、昆侖山脈南側(cè)，屬可可西里腹地，范圍在35°19′N—35°54′N，91°21′E—93°39′E，海拔在4 400 m以上。研究區(qū)屬于半干旱氣候區(qū)，為蒸發(fā)大于補給的水量收支負(fù)平衡區(qū)域，降雨量由東南向西北逐漸減少[10]。土溫低、結(jié)凍期長，下伏深厚的多年凍土，地表寒凍冰緣地貌形態(tài)多種多樣。自西向東依次分布有卓乃湖、庫賽湖、海丁諾爾湖和鹽湖。2011年以前，4個湖泊各自水量基本平衡。卓乃湖外溢后，4個湖泊自西向東建立了水力聯(lián)系，鹽湖來水量不僅包括原有的湖面降雨、四周河流、溝道匯水補給等，還增加了從上游卓乃湖、庫賽湖、海丁諾爾湖向其輸送的水量。研究區(qū)地理位置見圖1。

3 數(shù)據(jù)與方法

3.1 數(shù)據(jù)準(zhǔn)備

3.1.1 氣象要素年際分析數(shù)據(jù)

五道梁氣象站位于93°03′E、35°08′N，距鹽湖子流域中心直線距離約30 km，區(qū)域大氣熱力和動力條件相似，采用五道梁氣象站數(shù)據(jù)具有一定的空間代表性[8]。因此選用五道梁氣象站年降雨、年平均氣溫、年蒸發(fā)數(shù)據(jù)用于1989—2018年鹽湖流域氣象要素年際間變化分析。

3.1.2 遙感數(shù)據(jù)

研究區(qū)12—5月份湖泊受冰封影響，選取1989—2018年10—11月無云遮擋的國產(chǎn)高分一號衛(wèi)星和國外Landsat衛(wèi)星數(shù)據(jù)。高分一號空間分辨率為16 m，數(shù)據(jù)來源于中國資源衛(wèi)星應(yīng)用中心。Landsat衛(wèi)星包括Landsat5TM、Landsat7ETM+、Landsat8OLI，空間分辨率為30 m，數(shù)據(jù)來源于美國地質(zhì)調(diào)查局[14]。具體遙感數(shù)據(jù)時間見表1。

3.1.3 VIC模型模擬所用數(shù)據(jù)

基于筆者之前的研究[10]可認(rèn)為VIC模型適用于鹽湖流域水文模擬，因此采用30 m分辨率的SRTM(Shuttle Radar Topography Mission)數(shù)字高程數(shù)據(jù)，并將可可西里鹽湖流域劃分為分辨率0.1°×0.1°的子網(wǎng)格120個、FAO世界糧農(nóng)組織的5′分辨率的土壤質(zhì)地分布圖、2009年歐空局基于MERIS數(shù)據(jù)發(fā)布的全球陸地覆蓋數(shù)據(jù)(ESA GlobCover)和陽坤研發(fā)的中國區(qū)域高時空分辨率地表氣象驅(qū)動數(shù)據(jù)集[15-17](CMFD)用于模型模擬。

由于研究區(qū)附近只有五道梁氣象站，將單一氣象站數(shù)據(jù)插值到整個流域可能難以滿足模型模擬精度，因此選用CMFD代替。將1989—2018年五道梁氣象站的年降雨量、年平均氣溫的觀測值與CMFD

表1 衛(wèi)星影像數(shù)據(jù)

數(shù)據(jù)集五道梁氣象站所在格點數(shù)據(jù)進行比較，評估結(jié)果見圖2。

圖2 1989—2018年五道梁氣象站觀測值與CMFD數(shù)據(jù)對比Fig.2 Comparison between observation and CMFD data of Wudaoliang meteorological station from 1989 to 2018

年降雨量的觀測值與CMFD數(shù)據(jù)的決定系數(shù)R2為0.952 6，RMSE為8.313 9；年平均氣溫的觀測值與CMFD數(shù)據(jù)的決定系數(shù)R2為0.986 8，均方根誤差RMSE為6.201 2。決定系數(shù)R2>0.95，且RMSE在誤差范圍內(nèi)，證明CMFD數(shù)據(jù)集在鹽湖流域具有較高的合理性和適用性，可用于流域的氣候變化與水文模擬研究。

3.1.4 冰川水量分析數(shù)據(jù)

估算冰川融水?dāng)?shù)據(jù)源于第二次冰川編目數(shù)據(jù)集[11](http：//westdc.westgis.ac.cn)與2017—2018年中國西部冰川編目數(shù)據(jù)集[12](http：//www.dx.doi.org/10.11922/sciencedb.j00001.00227)。

3.2 方法介紹

本文總體研究路線見圖3。

圖3 總體技術(shù)路線Fig.3 General technical route

3.2.1 氣象要素年際變化分析方法

利用線性傾向估計法計算分析氣候要素的變化速率。使用累積距平法來分析氣候要素的階段性變化特征。Mann-Kendall(M-K)檢驗用于氣象要素突變檢驗。

3.2.2 湖泊面積計算方法

對高分一號和Landsat衛(wèi)星數(shù)據(jù)進行輻射定標(biāo)、大氣校正、幾何校正、影像裁剪等預(yù)處理。針對陸地水體，基于衛(wèi)星影像多波段反射率數(shù)據(jù)歸一化比值運算得到的水體指數(shù)，能夠自動消除地形起伏的影響，并通過突出水體、抑制其他地物信息，是區(qū)分水體和非水體地表覆蓋的一種常用和有效辦法[18]。針對高分一號和Landsat衛(wèi)星的光譜特征，本文選擇歸一化水體指數(shù)(Normalized Difference Water Index，NDWI)作為特征參數(shù)用于提取湖泊面積。

具體步驟如下：

(1)使用大尺度均值漂移算法[19](Large-scale Mean Shift)，通過設(shè)置合適的空間和光譜參數(shù)，對NDWI特征參數(shù)進行圖像分割。

(2)計算分割后各個圖斑的NDWI參數(shù)平均值，利用直方圖雙峰分布準(zhǔn)則進行逐圖斑判斷，利用動態(tài)閾值算法自動化實現(xiàn)水體提取的最佳分割閾值。

(3)對提取結(jié)果進行人工目視檢驗并修正，最終得到湖泊面積提取結(jié)果。

3.2.3 VIC模型模擬方法

VIC模型進行水量平衡計算后，輸出研究區(qū)各網(wǎng)格的徑流和蒸發(fā)結(jié)果，同Lohmann等[20]開發(fā)的匯流模型相耦合，將網(wǎng)格內(nèi)的產(chǎn)流轉(zhuǎn)化為流域出口斷面的流量[10]?？稍诤闯隹跀嗝嬖O(shè)置出水口，通過設(shè)置匯流參數(shù)，可以將出口流量分為陸面凈補給量與凍土融水和地下水入湖補給量。

3.2.4 冰川水量估算方法

將第二次冰川編目數(shù)據(jù)集與2017—2018年中國西部冰川編目數(shù)據(jù)集中研究區(qū)內(nèi)冰川面積進行統(tǒng)計，利用式(1)[11]和式(2)[11]計算，取平均值，最終計算出2011年與2018年鹽湖流域冰川體積(單位km3)，將2011—2018年冰川體積變化量做線性插值為2015—2018年冰川體積變化量，并利用冰與水的體積比0.9將其換算成水量，參與湖泊水量平衡計算。

V=0.036 5A1.375；

(1)

V=0.043 3A1.29。

(2)

3.3 湖泊水量計算公式

鹽湖來水量主要包括湖面降雨補給、陸面凈補給入、冰川融水補給、凍土融水和地下水入湖補給4部分，損失量主要是湖面蒸發(fā)，湖面降雨根據(jù)湖面面積與降雨量計算，湖面蒸發(fā)根據(jù)湖面面積與水面蒸發(fā)計算。建立水量平衡方程為

ΔV=P+Rs+Rg+Ri-Ew。

(3)

式中：ΔV為湖泊蓄變量；P為湖面降雨量；Rs為陸面凈補給量；Rg為凍土融水和地下水入湖補給量；Ri為冰川融水補給量；Ew為湖面蒸發(fā)量。

由于VIC模型不包含冰川模塊，冰川融水會不斷補給到河道徑流中[21]，認(rèn)為冰川融水匯入陸面凈補給中，因此模擬結(jié)果中陸面凈補給量需扣除冰川融水量。

3.4 相關(guān)性分析方法

將1989—2010年、2011—2018年湖泊總面積與氣象要素累積距平[14]進行Pearson相關(guān)性計算，分析氣象要素對湖泊面積變化的影響，并分析湖泊面積與湖泊水量間關(guān)系。

利用2015—2018年VIC模型模擬結(jié)果得到的湖泊水量變化，定量分析氣象要素對湖泊水量變化的影響，計算影響湖泊水量變化的主要因素。

4 研究結(jié)果

4.1 鹽湖流域氣候變化特征

由圖4可知：①可可西里湖泊群受季風(fēng)影響，1989—1995年降雨量總體呈下降趨勢，1995—2018年降雨量總體增加趨勢顯著，多年平均變化速率為36.2 mm/(10 a)。②1989—1997年平均氣溫總體呈下降趨勢，1997—2018年平均氣溫總體增加趨勢顯著，多年平均變化速率為0.59 ℃/(10 a)。③1989—2018年間年蒸發(fā)量總體呈下降趨勢，呈波動變化，下降速率為14.5 mm/(10 a)，且自2015年后年蒸發(fā)量迅速下降。

圖4 年降雨量、年平均氣溫、年蒸發(fā)量線性趨勢分析Fig.4 Linear trend analysis of annual rainfall，annual average temperature and annual evaporation

利用M-K檢驗對年降雨量、年平均氣溫、年蒸發(fā)量進行突變檢驗(圖5)，結(jié)果表明降雨量在2006年存在突變點，降雨量不斷累積，且在2010年后UFK統(tǒng)計量超過U0.05特征值，上升趨勢顯著，為后期卓乃湖外溢提供了基礎(chǔ)。年平均氣溫在1989—2018年間沒有明顯突變點，在2003年以后UFK統(tǒng)計量超過U0.05特征值，年平均氣溫呈顯著上升趨勢。年蒸發(fā)UFK統(tǒng)計量與UBK統(tǒng)計量均沒有超過特征值線，下降變化趨勢不顯著。

圖5 年降雨量、年平均氣溫、年蒸發(fā)量M-K檢驗結(jié)果Fig.5 M-K test of annual rainfall，annual average temperature and annual evaporation

4.2 湖泊面積變化特征

1989—2010年4個湖泊面積先減小后增大，且均在1995年發(fā)生轉(zhuǎn)折。2011年卓乃湖面積較前年驟減105.48 km2，此后平均下降速率為2.02 km2/a。2011年庫賽湖面積較前年增加54.58 km2，此后平均下降速率為2.06 km2/a。2011年海丁諾爾湖面積較前年增加29.13 km2，此后平均下降速率為0.58 km2/a。2011年鹽湖接收上游來水，面積迅速擴張，此后平均上升速率為17.46 km2/a，變化最為劇烈。4個湖泊部分水體面積計算結(jié)果、面積變化遙感圖像及變化曲線見表2、圖6、圖7。

表2 部分年份湖泊面積計算結(jié)果

圖6 鹽湖流域部分年份面積變化遙感圖像Fig.6 Remotely sensed images of area change of Salt Lake Basin

圖7 1989—2018年鹽湖流域湖泊面積變化Fig.7 Lake area change in Salt Lake Basin from 1989 to 2018

4.3 氣象要素對湖泊面積變化的影響

利用SPSS軟件分析1989—2010年、2011—2018年兩個時段湖泊面積與氣象要素的相關(guān)性。其中將2011—2018年細(xì)分為2011—2014年和2015—2018年兩個時段。將2015年前后分開是為了與VIC模型定量分析做對比驗證。2011年前研究單獨湖泊面積與其相關(guān)性，2011年后研究四湖總面積與其相關(guān)性。分析中只考慮年降雨量、年平均氣溫和年蒸發(fā)量3個要素，結(jié)果見表3、表4。

表3 鹽湖流域湖泊面積與氣象要素的Pearson相關(guān)系數(shù)

表4 鹽湖流域四湖總面積與氣象要素的Pearson相關(guān)系數(shù)

1989—2010年卓乃湖與鹽湖湖泊面積與年降雨量、年平均氣溫呈顯著正相關(guān)，與年蒸發(fā)量呈顯著負(fù)相關(guān)(P<0.01)，年降雨量相關(guān)性系數(shù)最大。

1989—2010年庫賽湖與海丁諾爾湖泊面積與年降雨量、年平均氣溫呈顯著正相關(guān)，與年蒸發(fā)量呈顯著負(fù)相關(guān)(P<0.01)，年平均氣溫相關(guān)性系數(shù)最大。

2011—2018年四湖總面積與3個氣象要素呈顯著正相關(guān)，年降雨量相關(guān)性系數(shù)最大，年蒸發(fā)量次之，年平均氣溫最小。細(xì)化時段分析可知2011—2014年四湖總面積與3個氣象要素均呈正相關(guān)，年降雨量相關(guān)性系數(shù)最大。2015—2018年四湖總面積與年降雨量、年平均氣溫呈顯著正相關(guān)(P<0.05)，年平均氣溫相關(guān)性系數(shù)最大。

綜上所述有：

(1)2011年前，卓乃湖和鹽湖面積變化可能受降雨量影響較大，庫賽湖與海丁諾爾湖面積變化可能受氣溫影響較大。

(2)2011—2018年4個湖泊面積變化可能受年降雨量影響最大。

(3)細(xì)化時間段，2011—2014年可能受年降雨量影響較大，年平均氣溫影響次之，2015—2018年可能受年平均氣溫影響較大，年降雨量影響次之。

相關(guān)性計算結(jié)果不同，可能原因為2011年4個湖泊受卓乃湖外溢影響，面積變化存在突變，而氣候變化對于湖泊水量變化是一個前期累積的過程，沒有達(dá)到與面積變化一致的突變程度，二者數(shù)據(jù)統(tǒng)計的相關(guān)性受到影響。此外還可能與選取計算的時間序列的長度等因素有關(guān)。

4.4 湖泊面積與湖泊水量間關(guān)系

湖泊水量的變化表現(xiàn)在面積變化，其變化是所在補給流域水量平衡的結(jié)果[22]。借助水位-面積、水位-容積曲線關(guān)系得到湖泊面積-容積間關(guān)系。

利用得到的庫賽湖、鹽湖面積-容積關(guān)系曲線，擬合曲線推求出相應(yīng)水量，最終面積-水量關(guān)系見圖8。用Pearson相關(guān)性分析庫賽湖、鹽湖面積與容積變化之間關(guān)系，相關(guān)系數(shù)分別為0.967、0.960，呈顯著正相關(guān)(P<0.01)。湖泊面積變化與湖泊水量變化有較好的相關(guān)性，湖泊面積的變化間接體現(xiàn)湖泊水量變化，氣象要素與湖泊面積間的變化可以間接體現(xiàn)出氣象要素對湖泊水量變化的影響。

圖8 庫賽湖和鹽湖2011—2018年面積-水量關(guān)系Fig.8 Area-volume relations in 2011-2018 for Kusai Lake and Salt Lake

Pearson相關(guān)性一定程度上可以體現(xiàn)年尺度上氣象要素與湖泊面積變化之間有一定相關(guān)性，雖然可通過面積-容積曲線關(guān)系間接反映氣象要素對水量的影響，但Pearson相關(guān)性分析是對兩個變量之間線性相關(guān)程度的衡量[23]，其結(jié)果可能與實際情況不同，需引入水文模型進行定量計算進行對比驗證。

4.5 定量分析氣象要素對湖泊水量的影響

2011年后湖泊總面積變化與降雨量、年平均氣溫、年蒸發(fā)量相關(guān)性很高，可說明氣象要素與流域水量變化之間存在關(guān)系。對于統(tǒng)計結(jié)果的差異性，需要通過VIC模型定量計算氣象要素對湖泊水量的影響，由于用于模型驗證的數(shù)據(jù)目前只有2015—2018年，因此選用2015—2018年時間段定量分析相關(guān)性。

經(jīng)VIC模型模擬和匯流模型匯流后得到各出水口水量變化結(jié)果。2015年和2018年計算結(jié)果見圖9、表5，結(jié)合計算的流域內(nèi)冰川補給量，最終2015—2018年鹽湖流域水量變化見表6。

圖9 2015年和2018年鹽湖流域水量分析結(jié)果Fig.9 Analysis results of water quantity in the Salt Lake Basin in 2015 and 2018

表5 VIC模型模擬結(jié)果

表6 鹽湖流域水量變化

將結(jié)果進行整理，陸面凈補給變化量扣除冰川補給量后算入降雨量貢獻中，可得降雨增加量、凍土釋水和地下水補給增加量、冰川融水量對鹽湖流域擴張的貢獻約為34.48%、57.66%、7.86%。氣溫變化成為影響湖泊水量變化的主要因素，降雨量影響次之。與2015—2018年P(guān)earson相關(guān)分析結(jié)果一致。

綜合統(tǒng)計結(jié)果及水文模型定量分析可知：

(1)1989—2010年間，卓乃湖、庫賽湖、海丁諾爾、鹽湖面積變化與年降雨量、年平均氣溫呈顯著正相關(guān)，與年蒸發(fā)量均呈顯著負(fù)相關(guān)。卓乃湖、鹽湖水量變化主要受降雨量的影響；庫賽湖、海丁諾爾湖水量變化主要受氣溫的影響。

(2)2011年后四湖連通，研究總面積與水量變化之間的關(guān)系，2011—2014年由統(tǒng)計結(jié)果得知水量變化主要受降雨量的影響，氣溫影響次之。2015—2018年統(tǒng)計結(jié)果結(jié)合水文模型模擬計算得知氣溫變?yōu)橹饕蛩?，降雨量影響次之?/p>

(3)2011年前4個湖泊各自水量平衡，并沒有上游湖泊下溢湖水的補給，而2011年后卓乃湖外溢改變了原有的來水組成，上游來水占補給中較大比例，上游補給水量的變化主要受氣候變化的影響。2015—2018年氣溫升高的速度遠(yuǎn)高于多年平均溫度，隨著氣溫升高，加快了區(qū)域冰川、凍土和積雪的消融速度，氣溫變化成為影響湖泊補給變化的主要因素。

5 結(jié)論與展望

(1)1989—2018年降雨量總體增加趨勢顯著，增長速率為36.2 mm/(10 a)，在2006年突變且在2010年后增加趨勢顯著，為后期卓乃湖外溢提供了基礎(chǔ)。年平均氣溫增長速率為0.59 ℃/(10 a)，2003年后升高趨勢顯著。年蒸發(fā)量下降速率為14.5 mm/(10 a)，并波動變化，且變化趨勢不顯著。

(2)1989—2010年湖泊面積先減小后增大，且均在1995年發(fā)生轉(zhuǎn)折。2011—2018年四湖總面積由664.09 km2擴大至753.63 km2，鹽湖面積擴大是四湖面積增加的原因。年降雨量、年平均氣溫、年蒸發(fā)量與湖泊面積密切相關(guān)，湖泊面積變化與湖泊水量變化相關(guān)性高，氣象要素與湖泊面積間的變化可以間接體現(xiàn)出氣象要素對湖泊水量變化的影響。

(3)結(jié)合統(tǒng)計分析及VIC模型定量分析結(jié)果可知2011年前卓乃湖和鹽湖水量變化主要受降雨量的影響，庫賽湖和海丁諾爾湖水量變化主要受氣溫的影響。2011—2014年4個湖泊水量主要受降雨量影響，2015—2018年4個湖泊水量變化主要受氣溫影響，降雨量影響次之。

(4)青藏高原的湖泊變化對氣候變化具有非常敏感的響應(yīng)，若未來氣候的暖濕化更加強烈，湖泊可能會出現(xiàn)更強烈的擴張。鹽湖流域水量持續(xù)增加可能對于下游的青藏公路、鐵路和保護區(qū)內(nèi)其他設(shè)施造成危害。未來要繼續(xù)細(xì)化水量研究，分析水量時空變化，完善統(tǒng)計分析中的不確定性，更充分評價湖泊水量變化與氣候變化之間的關(guān)系。