李 海，楊成生，惠文華，朱賽楠，張 勤

（1.長安大學地質(zhì)工程與測繪學院，陜西 西安 710054；2.中國地質(zhì)環(huán)境監(jiān)測院（自然資源部地質(zhì)災害技術指導中心），北京 100081）

0 引言

山地冰川作為冰凍圈的重要組成部分，其變化與氣候變化、區(qū)域水循環(huán)有著密切的關系[1?2]。對于屬于干旱與半干旱區(qū)的冰川發(fā)育區(qū)，冰川融水作為其徑流的主要補給源，是區(qū)域內(nèi)可靠的水源保證和生態(tài)保障[3?4]。同時冰川的持續(xù)消融大大提高了冰湖潰決的風險，危及下游百姓的生命財產(chǎn)安全[5?7]。中國冰川集中分布在青藏高原及其周圍的喜馬拉雅山脈、天山、喀喇昆侖等山脈，各地冰川分布見表1[8]。研究表明，除喀喇昆侖山地區(qū)少量冰川以外，中國西部冰川普遍處于退縮、變薄狀態(tài)，但是受氣候和地形差異的影響，不同區(qū)域存在差別[9?12]。其中青藏高原地區(qū)冰川退縮速率總體呈現(xiàn)由邊緣向內(nèi)陸腹地、東部向西部變慢的規(guī)律[13]。冰湖是由冰川挖蝕、融化形成洼坑或者冰磧物堵塞致冰川槽谷積水而成的一類湖泊[14?15]。王欣等[16]分析了近30年來我國喜馬拉雅山區(qū)冰湖變化的特征，剔除在氣候變化的作用下，冰川變化是冰湖變化的主要影響因素。近年來喜馬拉雅地區(qū)冰川退縮嚴重，一些以冰川融水為主要補給的湖泊，面積擴張較快，冰湖潰決危險性急劇升高，引起了研究者的廣泛關注。

表1 中國各地冰川面積分布Table 1 Glacier area distribution in China

喜馬拉雅地區(qū)是中國重要的冰川作用區(qū)，也是冰湖潰決等地質(zhì)災害的多發(fā)區(qū)[17?18]。在已有的15次冰湖潰決事件記錄中，12次發(fā)生在本區(qū)。其中藏南的希夏邦瑪峰地區(qū)就曾多次發(fā)生冰湖潰決災害[19?20]，如1968年的阿亞錯、1981年波曲的章臧布錯和2002年聶拉木縣北側(cè)的嘉龍錯，并且上述冰湖潰決引發(fā)的洪水和泥石流均造成多處人工設施損壞和人員傷亡[21?22]。鑒于希夏邦瑪峰地區(qū)冰川廣布、作用強烈、冰湖潰決災害多的特點，文章利用陸地資源衛(wèi)星（Landsat系列）影像，通過面向?qū)ο蟮拈撝捣诸惙椒ㄌ崛×讼Ｏ陌瞵敺宓貐^(qū)1994—2018年共9期冰川冰湖的面積，并探究了冰川冰湖的變化規(guī)律。對指導該地區(qū)的水資源調(diào)控和冰川冰湖災害的防治有重要的參考意義。

1 研究區(qū)與數(shù)據(jù)

1.1 研究區(qū)概況

研究區(qū)位于西藏南部喜馬拉雅山脈中段，內(nèi)有希夏邦瑪峰。該地區(qū)主要受到大陸西北季風和南亞西南季風環(huán)流影響，由于地勢阻隔，山脈南北兩側(cè)氣候相差較大。南側(cè)夏季高溫多雨，北側(cè)相對寒冷干旱[23?24]。研究區(qū)山脈南北兩側(cè)皆以沉積構(gòu)造帶為主，地層疏松、板塊活躍，在頻繁的地震和冰川作用的影響下，滑坡、冰湖潰決等地質(zhì)災害頻繁發(fā)生[25]。此外，研究區(qū)內(nèi)有中國和尼泊爾重要邊境通商口岸城市吉隆鎮(zhèn)和聶拉木縣，政治經(jīng)濟地位獨特。

1.2 數(shù)據(jù)介紹

本次冰川冰湖的監(jiān)測數(shù)據(jù)包括遙感數(shù)據(jù)、DEM數(shù)據(jù)、氣象數(shù)據(jù)、冰川編目數(shù)據(jù)4種。其中遙感影像主要來源于美國陸地資源衛(wèi)星（Landsat系列）TM5、ETM、ETM+、OLI8傳感器影像，皆由美國地質(zhì)調(diào)查局免費提供（https://earthexplorer.usgs.gov/）。首先，經(jīng)過影像篩選發(fā)現(xiàn)該地區(qū)每年10月份前后的影像具有云量少、積雪少的特點。其次，以往研究表明1992年后該地區(qū)氣溫上升速率明顯加快[26]，綜合影像資源、云量、積雪等多方面因素，選取1994—2018年跨度24年9期，每期大致10月份的影像開展研究（表2）。冰川參考數(shù)據(jù)選用2013年倫道夫冰川編目數(shù)據(jù)(RGI, Version 3.2，2013)，基本反映2010年前的冰川分布情況。地形數(shù)據(jù)選用SRTM DEM，高程精度為16 m（95%置信區(qū)間）[27]。

表2 光學影像數(shù)據(jù)詳情Table 2 Details of optical image data

對于氣象數(shù)據(jù)的選取，由于該地區(qū)氣象站點稀疏，僅聶拉木縣設有氣象站，但是聶拉木縣位于山脈南側(cè)峽谷之間，受南亞暖濕季風影響顯著，不能準確反映山脈北側(cè)冰川區(qū)域的氣候[28?29]。全球網(wǎng)格氣象再分析數(shù)據(jù)集，采用特殊的插值方法，充分考慮了地形、距離、海拔等因素對氣候影響。因此采用來源于美國國家海洋和大氣管理局（https://psl.noaa.gov/data/gridded/）提供的NOAANCEP陸地表面0.5°×0.5°網(wǎng)格氣溫再分析數(shù)據(jù)和CMP 0.5°×0.5°網(wǎng)格降水再分析數(shù)據(jù)。其中氣溫數(shù)據(jù)均方根誤差為0.18 ℃，日降水量偏差小于1 mm，精度較高[30?31]。

2 數(shù)據(jù)處理

文章在面向?qū)ο蟮幕A上采用閾值分類方法提取冰川冰湖面積。首先在ENVI5.3軟件中對Landsat系列影像進行大氣校正、影像融合、拼接等影像預處理工作。為了更準確的提取冰川邊界，運用面向?qū)ο蟮拈撝捣诸惙椒▽τ跋襁M行了分割與分類工作，全部處理流程如圖1所示。為了方便后期的交叉驗證和人工修正，采用倫道夫冰川編目數(shù)據(jù)或前一期冰川提取矢量結(jié)果參與的多尺度分割方法。通過多次試驗，最終確定分割因子為60、形狀因子為0.4、緊致度因子為0.5。

圖1 數(shù)據(jù)處理流程Fig.1 Data processing flow

對山體陰影區(qū)的分類主要采用目視解譯進行。為了避免山體陰影對冰川分類結(jié)果的影響，首先對山體陰影進行提取，并在結(jié)合Google-Earth多期歷史影像交叉驗證陰影區(qū)冰川范圍。利用山體陰影與冰湖具有相似的光譜響應特征，以及其在近紅外波段的超低反射，并經(jīng)多次試驗，陰影提取采用歸一化水指數(shù)NDWI>0.18、近紅外光譜NIR<1 000。隨后，采用NDWI>0.18對冰湖進行了提取?？紤]到提取精度和實際意義，冰湖提取排除了面積小于3個像素冰湖。

采用波段比值法提取凈冰川。根據(jù)冰川在紅波段的高反射和短波紅外高吸收的光譜特性，實驗采用R/SWIR1>2.5進行冰川提取。然而研究區(qū)地面常存在的少量殘留積雪，影響冰川提取的準確度。而以往的研究往往只采用2～3期影像對冰川進行長時間序列監(jiān)測，時間間隔較長，無法發(fā)現(xiàn)積雪的變化以降低積雪的影響或者需要大量的人工目視判讀（表2）。而本文收集了研究區(qū)1994—2018年共9期影像，時間間隔短，可以通過前后對比發(fā)現(xiàn)積雪的變化。在冰川持續(xù)退縮的大前提下，自動去除前時相冰川邊界以外的積雪。其中1994年的冰川提取則參考RGI數(shù)據(jù)。由于RGI數(shù)據(jù)反映的是2000—2010年間的冰川狀況，故僅以此作為參考并輔以人工目視判讀，剔除積雪。2000年的冰川提取參考1994年冰川提取結(jié)果。以此類推提取了9期的冰川冰湖面積，結(jié)果如表3所示。冰川冰湖整體提取結(jié)果如圖2所示。

圖2 1994年冰川冰湖輪廓提取圖（TM影像5、4、3合成）和單個冰川、冰湖9期變化結(jié)果Fig.2 Overall extraction of glacial and glaciers lake in 1994 (TM images 5, 4, 3 are synthesized),the results of 9-stage change of single glacial and glaciers lake

由于本次實驗沒有地面測量資料，因此冰川冰湖誤差評估參照第二次冰川編目誤差評估辦法，即將空間分辨率的一半作為緩沖區(qū)域確定面積不確定度[11]。因此冰川冰湖面積變化誤差е的公式（1）為：

式中：ε1、ε2——前后時相冰川、冰湖面積的誤差。

由誤差傳播定律可知，面積變化率誤差ε的計算公式為：

式中：x、y——分別指冰川、冰湖前后兩期的面積。

表磧覆蓋型冰川采用較為可靠的目視解譯。研究區(qū)南側(cè)受到南亞暖濕季風影響，屬于溫帶海洋性冰川區(qū)，表磧覆蓋型冰川廣布。然而由于表磧冰川在光譜特性上和裸地相似，難以實現(xiàn)自動化提取。近年來，不少學者針對表磧冰川的自動化提取提出了許多方法，但是都存在準確度不高且需要專家目視修正的問題，不適用于變化分析[32?33]。鑒于該地區(qū)表磧覆蓋型冰川較集中、面積大且數(shù)量少、退縮緩慢等特點，參考倫道夫冰川編目數(shù)據(jù)，以目視解譯方式提取1994、2018年首尾兩期的表磧冰川輪廓。1994年表磧冰川面積為（168.41±17.6) km2。至2018年，在部分凈冰川退縮裸露出表磧冰川的形勢下，表磧冰川面積增加至(223.97±19.48 ) km2。而表磧冰川的退縮主要表現(xiàn)為表磧冰川向末端冰湖演變，面積為11.03 km2。

3 分析與討論

3.1 冰川冰湖總體變化分析

文章共提取了1994—2018年間共9期冰川冰湖的輪廓，面積變化如圖3所示。結(jié)果顯示24年間冰川總體上持續(xù)退縮，面積從(1 012.66±68.42) km2退縮至(699.89±61.50 ) km2、退縮比例為(30.83±7.66) %、年平均退縮速率為(1.28±0.32) %，冰川的退縮速率經(jīng)歷了加速減速的輪換。該地區(qū)的冰湖在24年間面積從(36.15±4.32 ) km2擴張至(50.65±7.09 ) km2、擴張比率為(47.23±25.79) %、年平均擴張速率為(1.889±1.07)%，冰湖擴張速率高于冰川退縮速率，并且其擴張的快慢變化與冰川退縮的快慢變化基本吻合，相關系數(shù)為?0.98。較高的相關性說明該地區(qū)冰湖擴張水源基本來源于冰川融水，間接證明了冰川、冰湖的提取方法可靠。

圖3 24年冰川冰湖面積變化Fig.3 Changes of glaciers and glacial lakes area in 24 years

3.2 冰川變化與氣候變化的響應

文章探究了該地區(qū)冰川冰湖變化與氣候變化的響應關系。研究表明，冰川冰湖變化和當?shù)氐臍夂蜃兓嬖陲@著的關聯(lián)，其中氣溫和降水影響最大[34]。

用1992—2018年近27年的夏季氣溫和全年氣溫變化分別與冰川面積變化作比較。為了防止因部分異常氣溫影響擬合結(jié)果，采用魯棒性更好的穩(wěn)健估計法擬合氣溫變化。從圖4（a）可以看出全年平均氣溫和夏季平均氣溫均呈上升趨勢，與冰川面積變化呈負相關，這可能是該地區(qū)冰川退縮速率加快的重要原因。夏季平均氣溫以每年0.035 ℃的速率上升，并從 2010年以來上升了1.8 ℃左右，為冰川在夏季的消融創(chuàng)造了良好的條件，這也導致2010年后冰川退縮速率加快。

圖4 研究區(qū)氣溫和降水變化Fig.4 Variation of temperature and precipitation in the study area

同時，分析了該地區(qū)1990—2019年的近30年的降水變化，將其細分為夏季、冬季、全年降水總量三組。如圖4(b)所示，該地區(qū)降水主要集中在夏季，且全年降水和夏季降水均呈現(xiàn)不斷上升趨勢。與此同時，夏季降水以15.8 mm/a的速率上升，趨勢明顯，但冬季降水較少且無明顯上升。30年間降水量增加了約500 mm，年增量為17.2 mm，且大部分來自夏季降水，對應冰川在消融季的退縮。而春冬冰川物質(zhì)積累季，降水并無明顯增多。綜上所述，夏季氣溫和降水的升高是本地區(qū)冰川的退縮加快的重要原因。與此同時降水增加和冰川的加速退縮，共同促進冰湖的快速擴張，提高了該地區(qū)冰湖潰決的風險。

3.3 影響冰川變化的其他因素

為了進一步研究該地區(qū)冰川變化的規(guī)律，本文進一步討論了冰川面積的大小、冰川經(jīng)緯度位置、流向多個因素對冰川變化的影響。首先利用SRTM DEM提取研究區(qū)域內(nèi)的山脊線作為分冰嶺，將冰川與冰川之間按山脊線分開，劃分成獨立的冰川。山脊線主要采用水文分析法和表面分析法相結(jié)合的方法自動化提取[11]。由于DEM精度問題，部分山脊線會偏離山脊位置，需要少量人工修定。整個冰川區(qū)經(jīng)過分冰嶺分割之后，大約有2 000個冰川，其中細碎冰川較多，分別在1994、2018年的結(jié)果中選擇其中100個冰川作為統(tǒng)計樣本，從面積大小、冰川緯度、經(jīng)度位置、流向四個角度進行統(tǒng)計，統(tǒng)計結(jié)果如表3所示。

表3 9期冰川冰湖面積表Table 3 Area of glaciers and glacial lakes in Phase 9

通過統(tǒng)計不同大小冰川的變化率，探究冰川大小對冰川變化的影響。根據(jù)冰川的數(shù)量和面積特征，將100個冰川按面積大小分為6組，分別是0～1 km2、1～2 km2、2～3 km2、3～4 km2、4～7 km2、7～20 km2。統(tǒng)計結(jié)果如表4所示，隨著面積的增大，冰川的退縮速率變小，說明：①面積較大的冰川應對外界氣候變化更具有穩(wěn)定性，同時面積較大的冰川對氣候具有較大的調(diào)節(jié)作用；②小面積的冰川冰量損失極為嚴重，因此隨著冰川的持續(xù)退縮，相同的氣候的變化將帶來更嚴重的退縮。

表4 不同角度下冰川的退縮率Table 4 Retreat rates of glaciers with multipectives

通過統(tǒng)計不同經(jīng)緯度位置的冰川的變化率，探究位置對冰川變化的影響。從經(jīng)緯度位置來看，研究區(qū)北側(cè)冰川和西側(cè)冰川退縮速率較慢，這與其他研究者獲得的兩大規(guī)律基本吻合：①喜馬拉雅山脈冰川退縮速率整體呈現(xiàn)由東側(cè)向西側(cè)逐步降低；②喜馬拉雅山脈冰川退縮速率外側(cè)高于內(nèi)側(cè)[35?36]。該結(jié)果與地形差異而導致的氣候差異有著密切的關系。研究區(qū)受南亞印度洋溫濕季風和大陸西北干冷季風雙重影響[37]，由于山體的阻隔印度洋季風由南向北迅速減弱，而由西北而來的大陸性季風因希夏邦瑪山脈的阻隔向東減弱。致使冰川退縮速率出現(xiàn)了南快北慢、東快西慢的現(xiàn)象。

通過統(tǒng)計不同流向冰川的變化率，探究流向?qū)Ρㄗ兓挠绊?。?00個冰川大致按北、東北、東南等6個方向分類。冰川變化結(jié)果如表3所示。結(jié)合圖1，北向冰川多位于研究區(qū)南側(cè)，冰川面積小，南亞季風影響嚴重，退縮較大。東北流向冰川主要分布于希夏邦瑪峰東北側(cè)，該地區(qū)冰川退縮最慢，其次是西北流向冰川，主要位于右側(cè)的錯朗瑪冰川群北側(cè)，主要原因是希夏邦瑪冰川群相對緯度高，氣溫相對較低。東南流向冰川主要分布在錯朗瑪冰川群的南側(cè)地區(qū)，西南流向冰川為希夏邦瑪峰西南側(cè)冰川，兩者皆為南亞季風的迎風坡，但是西南流向冰川較慢可能是左側(cè)存在山脈阻隔，南亞季風較弱所致。綜上所述，經(jīng)緯度位置和流向?qū)Ρㄗ兓挠绊懀饕蛟谟诰植繗夂虻牟町悺?/p>

因此根據(jù)本文研究結(jié)果，提出以下建議：①持續(xù)關注氣溫和降水對冰川的影響；②對面積擴張嚴重的高危冰湖進行更細致的監(jiān)測；③走健康綠色發(fā)展道路，減緩全球變暖的速率。

4 結(jié)論

本文利用Landsat系列影像對冰川變化顯著的希夏邦瑪峰區(qū)域進行了長時間序列的冰川冰湖變化規(guī)律研究。結(jié)果表明：①1994—2018年該地區(qū)冰川面積退縮比例為30.83%，冰湖面積擴張比率為47.23%，冰湖擴張嚴重；②該地區(qū)冰川退縮速率加快的主要原因在于夏季氣溫的升高；③隨著冰川面積的逐漸減小，相同的氣溫變化將帶來更加嚴重的冰川退縮；④南亞季風、西北大陸季風、地形等多重因素致使該地區(qū)冰川退縮速率整體呈現(xiàn)南快北慢、東快西慢的現(xiàn)象。