閆斌，賈洪果，任文靜，吳仁哲，黃心茹

西南交通大學(xué) 地球科學(xué)與環(huán)境工程學(xué)院,成都 611756

1 引 言

冰川是重要的水資源，由于其對(duì)河川徑流起著自然地調(diào)節(jié)作用，被譽(yù)為“高山固體水庫(kù)”；冰川的運(yùn)動(dòng)消融也極易引發(fā)自然災(zāi)害，如冰川躍動(dòng)、冰湖潰決、冰川泥石流等，因此冰川研究對(duì)人類(lèi)社會(huì)具有重要意義（鄔光劍等，2019；張九天等，2012；陳寧生等，2019）。冰川作為一種動(dòng)態(tài)資源，與氣候變化密切相關(guān)（劉潮海等，2002；趙曉艷等，2022）。相關(guān)研究（陳虹舉等，2017；安國(guó)英等，2019）表明：中國(guó)冰川變化趨勢(shì)主要呈現(xiàn)對(duì)夏季平均氣溫變化表現(xiàn)為正響應(yīng)，而對(duì)年降水量變化表現(xiàn)為負(fù)響應(yīng)，盡管冰川分布區(qū)年降水量的增加導(dǎo)致冰川積累量增多，但這并不足以抵消因溫度升高而增加的消融量，因此溫度的升高現(xiàn)已成為中國(guó)西部冰川快速退縮的主導(dǎo)性因素。

冰湖是冰川運(yùn)動(dòng)的產(chǎn)物，通常由冰川挖蝕成的洼坑和冰磧物堵塞冰川槽谷積水而形成。冰湖個(gè)數(shù)增多、面積增大的同時(shí)，冰川則往往處于退縮狀態(tài)，因此冰湖的形成演變與冰川的變化密切相關(guān)（姚曉軍等，2017）。在全球氣候變暖的背景下，由冰湖潰決引發(fā)的洪水和泥石流災(zāi)害呈現(xiàn)數(shù)量增多、危害程度加劇特征（Wang 和Zhou，2017）。由于冰湖多位于高海拔地區(qū)，其潰決所產(chǎn)生的地質(zhì)災(zāi)害具有突發(fā)性強(qiáng)、破壞性大、持續(xù)時(shí)間短、波及范圍廣的特點(diǎn)（Bolch 等，2008；韓金良等，2007；劉建康等，2019；楊瑞敏等，2019）。據(jù)《昌都地區(qū)志》記載，1974 年7 月6日，丁青縣的波戈冰川冰湖潰決，大量湖水泄入河谷，沿河木橋全部沖走，黑（河鎮(zhèn)）昌（都）公路部分受損（姚曉軍等，2014）。因此對(duì)冰湖潰決災(zāi)害研究具有重要的社會(huì)意義。

本文所選研究區(qū)為位于藏東南地區(qū)唐古拉山東段的布加崗日冰川，該冰川的冰湖周?chē)匚镏饕ū嗪完懙?。目前，在冰湖區(qū)基于歸一化水指數(shù)NDWI（Normalized Difference Water Index）進(jìn)行冰湖提取已有先例借鑒（李均力等，2011；Li 和Sheng，2012；駱劍承 等，2009；Watson 等，2018；Wendleder 等，2018），此種方法適用于目標(biāo)冰湖與周?chē)匚镄畔⒃贜DWI的直方圖中呈雙峰特性的情況；并實(shí)驗(yàn)發(fā)現(xiàn)，若冰湖與周?chē)匚镄畔⒃贜DWI直方圖上呈現(xiàn)三峰，如果直接使用閾值分割法極易造成錯(cuò)分現(xiàn)象。針對(duì)上述問(wèn)題，本文提出基于NDWI?NDSI 組合閾值的冰湖面積提取方法。首先利用歸一化雪被指數(shù)NDSI（Normalized Difference Snow Index）具有分離陸地（目標(biāo)）與冰水混合地物（背景）的特性制作去陸地掩膜，再對(duì)掩膜后的NDWI進(jìn)行閾值分割，以準(zhǔn)確提取冰湖面積信息。為驗(yàn)證方法的可行性，首先基于冰川區(qū)域分布的主要地物的NDWI?NDSI 散點(diǎn)圖和冰湖區(qū)的NDSI 直方圖進(jìn)行陸地與冰水混合地物分離，接著分別基于NDWI 和NDWI?NDSI 兩種方法對(duì)冰湖提取進(jìn)行對(duì)比分析。

2 研究區(qū)概況及數(shù)據(jù)介紹

如圖1所示，布加崗日冰川位于距離川藏北線(xiàn)G317 國(guó)道不遠(yuǎn)的西藏自治區(qū)那曲市索縣、巴青縣及昌都市丁青縣的交界處，與唐古拉山脈一致為東西走向。其海拔較高，有16 座海拔超過(guò)6000 m的高峰，最高峰海拔為6328 m。受來(lái)自印度洋的西南季風(fēng)作用，年降水量達(dá)到600—700 mm，高聳的地勢(shì)和豐富的降水使這里成為較大的山岳冰川發(fā)育集中地（Liu等，2016；王寧練和丁良福，2002；王聰強(qiáng)等，2016）。據(jù)光學(xué)影像目視解譯，該冰川發(fā)育面積較為明顯的冰湖約為12 個(gè)，它們分別對(duì)應(yīng)編號(hào)L?1，L?2，…，L?12，其位置如圖2所示。

圖1 布加崗日冰川地理位置及主要冰川分布Fig.1 Geographical location of Bujiagangri glacier and distribution of major glaciers

圖2 布加崗日冰川冰湖分布圖（2016?08?19）Fig.2 Distribution map of Bujiagangri glacial and glacial lakes（2016?08?19）

本文所用數(shù)據(jù)源為美國(guó)NASA 的陸地衛(wèi)星（Landsat）TM/OLI，衛(wèi)星波段介紹如表1 所示。此次實(shí)驗(yàn)所使用的Landsat5/8 衛(wèi)星波段有：Green 波段、NIR波段、SWIR1波段，所選多光譜波段分辨率均為30 m?；诠鈱W(xué)影像成功提取冰湖面積的前提條件是影像獲取時(shí)間段內(nèi)研究區(qū)上空無(wú)云或云霧較少、且冰湖處于非結(jié)冰期，而布加崗日所在地區(qū)的夏季多數(shù)影像數(shù)據(jù)被云霧覆蓋嚴(yán)重?zé)o法使用，同時(shí)部分冬季影像由于冰湖結(jié)冰或被雪覆蓋也導(dǎo)致數(shù)據(jù)不符合要求?？紤]以上因素，本文篩選出1988 年—2018 年間共16 景質(zhì)量較好且覆蓋布加崗日冰川的Landsat（TM/OLI）影像。對(duì)所選影像首先裁剪處理，并進(jìn)行輻射定標(biāo)和大氣校正，以消除或減弱由傳感器和大氣帶來(lái)的各種誤差，避免其對(duì)后期處理結(jié)果的精確性產(chǎn)生影響。

表1 Landsat 5/8衛(wèi)星波段介紹Table 1 Introduction of Landsat 5/8 satellite bands

3 研究方法

NDWI（McFeeters，1996）又稱(chēng)為歸一化水指數(shù)，是基于水體的反射從可見(jiàn)光到中紅外波段逐漸減弱的原理，利用水體的高反射率主要集中在可見(jiàn)光中的藍(lán)綠光波段，而在可見(jiàn)光其它波段及近、中紅外的反射率很低的特性，從而使水體的亮度信息得到最大程度增強(qiáng)，非水體受到普遍的抑制，進(jìn)而達(dá)到突出水體的目的（Benoudjit 和Guida，2019；徐涵秋，2005），其表達(dá)式為

由Hall 等（1995） 提出的歸一化雪被指數(shù)NDSI 是目前光學(xué)遙感提取積雪的通用方法。實(shí)驗(yàn)發(fā)現(xiàn)，該指數(shù)除了能突出積雪信息外，對(duì)水體和冰同樣具有高敏感性，同時(shí)又對(duì)陸地信息具有很強(qiáng)的抑制作用，本文正是利用這一特性實(shí)現(xiàn)對(duì)冰湖區(qū)陸地信息的識(shí)別。如圖3所示，結(jié)合冰湖區(qū)影像的NDWI?NDSI 散點(diǎn)圖中各地物信息表現(xiàn)出的差異性與其N(xiāo)DSI 灰度直方圖呈現(xiàn)雙峰的條件，可由OTSU 法（Otsu，1979）計(jì)算得到分割閾值（紅線(xiàn)所示），據(jù)此可以分離出陸地信息。計(jì)算公式為

圖3 基于NDSI分離陸地與冰水的方法驗(yàn)證Fig.3 Method verification of separation of land，ice and water based on NDSI

基于NDWI?NDSI 組合法實(shí)施冰湖提取的過(guò)程如圖4 所示，與通常基于NDWI 提取冰湖不同的是，本文首先利用OTSU 對(duì)NDSI 進(jìn)行閾值分割得到閾值T1，將NDSI 中小于T1 的部分用于制作掩膜以剔除陸地信息，避免了NDWI在灰度直方圖呈現(xiàn)三峰的可能，再利用OTSU 對(duì)掩膜后的NDWI 所包含的地物信息（積雪、冰舌和冰湖）直方圖進(jìn)行分割得到閾值T2。

圖4 冰湖提取流程Fig.4 Flow chart of glacial lake extraction

由于局部影像中積雪量較少，不會(huì)在NDWI直方圖中產(chǎn)生三峰現(xiàn)象；其次，積雪與冰舌在可見(jiàn)光的藍(lán)綠波段反射率較低，而水體則相反，在NDWI 直方圖中能很好地利用閾值分割法進(jìn)行提取。因此，可繼續(xù)使用閾值分割法對(duì)閾值小于T2 的積雪和冰舌信息進(jìn)行剔除，并統(tǒng)計(jì)大于T2的像素個(gè)數(shù)，最后結(jié)合30 m 的影像分辨率成功實(shí)現(xiàn)冰湖面積的提取。為了驗(yàn)證該方法的正確性，本文分別選取冰湖L?10 和L?12 用于分析（圖5）。L?10 與冰舌和陸地相連，其N(xiāo)DWI 的直方圖呈現(xiàn)三峰，直接基于NDWI 提取的冰湖與RGB 影像對(duì)比可知，提取結(jié)果出現(xiàn)非冰湖區(qū)域；而基于NDWI?NDSI 組合法提取的冰湖則與目視判斷的基本一致，表明了該方法的正確性。與L?10 不同的是，L?12 周?chē)匚镏挥嘘懙兀銷(xiāo)DWI 直方圖顯示雙峰，分別基于NDWI 和NDWI?NDSI 兩種方法實(shí)施冰湖提取結(jié)果基本相同，說(shuō)明該組合法仍適用于雙峰情況；此外，將兩種方法提取的冰湖面積與目視解譯結(jié)果進(jìn)行對(duì)比可知（表2），NDWI?NDSI 組合法在上述兩種情況下提取冰湖的精度明顯優(yōu)于NDWI法。

表2 冰湖提取方法精度評(píng)價(jià)Table 2 Accuracy evaluation of glacial lake extraction method

圖5 基于NDWI與NDWI?NDSI的冰湖提取方法對(duì)比驗(yàn)證Fig.5 Comparative verification of two glacial lake extraction methods based on NDWI and NDWI?NDSI

4 結(jié)果分析

本文利用NDWI?NDSI 組合法成功提取了1988 年—2018 年的布加崗日冰川的冰湖面積，并從單個(gè)和整體兩方面做了時(shí)序研究和統(tǒng)計(jì)分析。

4.1 單個(gè)冰湖變化分析

由圖6 可知，30 年來(lái)布加崗日冰川的多數(shù)冰湖都有不同程度變化。除L?07、L?09、L?12這3個(gè)冰湖面積變化不大外，L?01、L?03、L?11 面積有一定增加，L?02、L?04、L?05、L?06、L?08、L?10面積增加趨勢(shì)明顯，其中L?04和L?06面積從無(wú)到有，這說(shuō)明該冰川冰湖個(gè)數(shù)在增加、面積在擴(kuò)大，也表明其融水逐年增多，融化加速。結(jié)合圖2冰川冰湖分布可知，變化不大的3個(gè)冰湖均未與冰舌相連，但面積增加明顯的6 個(gè)冰湖中有5 個(gè)與冰舌相連，而冰舌前端亦即冰川末端，其融水用于補(bǔ)充冰湖，這說(shuō)明冰川消融是引起冰湖面積變化的主要因素。此外，有必要指出的是冰湖面積出現(xiàn)有→無(wú)→有的變化是由于冰湖結(jié)冰引起的，如L?01 在1990 年、1992 年、1993 年、2002 年、2018 年均處于結(jié)冰狀態(tài)。

圖6 冰湖面積年際變化Fig.6 Interannual change of glacial lakes area

為了進(jìn)一步分析冰湖變化，對(duì)該冰川的12 個(gè)冰湖均做面積分級(jí)及其變化量分級(jí)處理。首先將其逐一按30 年中擴(kuò)張至最大年份的面積進(jìn)行分級(jí)得到表3，其次以1988 年冰湖面積為基準(zhǔn)，用面積最大年份的冰湖面積最減去1988年的冰湖面積，可得到近30 年中每個(gè)冰湖面積的最大變化量，再將變化量進(jìn)行分級(jí)得到表4，最后將12 個(gè)冰湖的面積及其變化在空間的分布情況進(jìn)行了分級(jí)展示，見(jiàn)圖7。由表3 可知，面積超過(guò)0.5 km2的冰湖有6個(gè)，面積在0.1—0.5 km2之間的冰湖有3個(gè)，其余均小于0.1 km2；如表4所示，面積變化最大的冰湖分別為L(zhǎng)?04、L?05、L?10，三者變化量均超過(guò)0.5 km2；其次為L(zhǎng)?02、L?06、L?08，變化量在0.1—0.5 km2之間；L?07、L?09、L?12 這3 個(gè)冰湖面積變化均小于0.01 km2。另外對(duì)冰湖面積及其變化量取交集可得，面積最大且變化量最大的3個(gè)冰湖分別為L(zhǎng)?04、L?05、L?10，這3 個(gè)冰湖極易發(fā)生潰堤且由此引發(fā)的地質(zhì)災(zāi)害會(huì)更為嚴(yán)重。

圖7 冰湖面積及其變化空間分布Fig 7 Spatial distributions of the area of glacial lakes and their changes

表3 冰湖面積分級(jí)表Table 3 Classification of glacial lakes area

表4 冰湖面積變化分級(jí)表Table 4 Classification of glacial lakes area changes

4.2 整體面積變化分析

在4.1 節(jié)中針對(duì)單個(gè)冰湖面積的時(shí)序演變情況做了分析，本節(jié)將對(duì)12 個(gè)冰湖面積累加以分析整個(gè)冰川的冰湖變化情況。但由于影像數(shù)據(jù)受結(jié)冰、落雪及云霧遮擋等影響導(dǎo)致部分較大冰湖面積未被正確提取，因此從16 景影像中篩選出1988 年、1994 年、1998 年、2013 年、2016 年 這5 組質(zhì) 量較好（即冰湖不受上述幾類(lèi)影響因素過(guò)多干擾而被正確提取）的冰湖數(shù)據(jù)進(jìn)行分析。

如圖8（a）與圖8（b）所示，布加崗日冰川冰湖總面積在1988年約為2.7666 km2，1994年增至約2.9799 km2，1988年—1994年期間年平均變化速率約為0.0356 km2/a；1998年增至約3.1959 km2，1994年—1998 年期間年平均變化速率約為0.0540 km2/a；2013 年增至約5.0409 km2，1998 年—2013 年期間年平均變化速率約為0.1230 km2/a；2016 年增至約5.2308 km2，2013年—2016年期間年平均變化速率約為0.0633 km2/a?？梢钥闯鼋?0 年布加崗日冰川的冰湖面積變化呈慢→較快→很快→較快趨勢(shì)，雖然面積變化不能準(zhǔn)確反映冰川的消融量及湖水變化量，但其結(jié)果仍然可以表明布加崗日冰川退縮明顯。結(jié)合圖9可知，冰湖總面積變化與溫度呈現(xiàn)正相關(guān)特性（該地區(qū)溫度數(shù)據(jù)為中國(guó)氣象局提供），表明該地區(qū)逐年氣溫的不斷升高正是冰川消融過(guò)快的主要原因；其次，經(jīng)統(tǒng)計(jì)發(fā)現(xiàn)，布加崗日冰川地區(qū)近30 年來(lái)平均氣溫在0 ℃以上的年份共7 個(gè)，而1998 年—2013 年（紅色虛線(xiàn)框內(nèi)）冰湖面積增長(zhǎng)速率之所以最快，與這段時(shí)間內(nèi)年均氣溫在0 ℃以上的年數(shù)較多緊密相關(guān)（5 個(gè)年份）；2013 年—2016 年（橙色虛線(xiàn)框內(nèi)）增長(zhǎng)速率變化次之，其原因?yàn)槟昃鶜鉁卦?℃以上的年數(shù)相應(yīng)較少（2個(gè)年份）。

圖8 冰湖總面積年際變化及其變化速率Fig 8 Interannual change and rate of change of total area of glacial lakes

圖9 布加崗日地區(qū)年均氣溫變化Fig 9 Annual average temperature change in Bujiagangri area

5 結(jié) 論

本文利用改進(jìn)的NDWI?NDSI 的組合閾值法，對(duì)1988 年—2018 年研究區(qū)內(nèi)的16 景Landsat 影像進(jìn)行了冰湖提取試驗(yàn)，研究的主要結(jié)論如下：

（1）歸一化積雪指數(shù)（NDSI）除了能夠提取積雪外還可以用于分離試驗(yàn)區(qū)內(nèi)陸地信息和冰雪水混合信息，據(jù)此可制作陸地掩膜；

（2）基于NDWI的閾值分割法在冰湖區(qū)提取含有冰舌、陸地兩種地物信息的冰湖時(shí)會(huì)出現(xiàn)錯(cuò)分現(xiàn)象，改進(jìn)的NDWI?NDSI 的組合閾值法成功解決了這一問(wèn)題，增加了閾值分割法提取冰湖的適用范圍；

（3）通過(guò)對(duì)試驗(yàn)區(qū)提取的冰湖做時(shí)序研究和統(tǒng)計(jì)分析可知，30 年來(lái)布加崗日冰川的冰湖面積從約2.7666 km2增加至約5.2308 km2，增加最快時(shí)的年際變化率約0.1230 km2/a。12 個(gè)冰湖中多數(shù)冰湖面積都有不同程度增加，其中L?04、L?05 和L?10 這3 個(gè)冰湖的面積（0.5—1.0 km2）及其擴(kuò)張面積（0.5—1.0 km2）最大，這3 個(gè)冰湖受夏季受高溫天氣影響蓄水過(guò)快隨時(shí)有潰堤的可能性，對(duì)下游居民的人身財(cái)產(chǎn)安全及川藏公路（G317）等形成潛在威脅，希望有關(guān)部門(mén)給予重視。此外，基于單個(gè)和整體的冰湖面積變化分析均表明該冰川30年消融過(guò)快，退縮嚴(yán)重。

本文的不足之處為提出的NDWI?NDSI 組合閾值法只能適用于單個(gè)冰湖提取，且該方法準(zhǔn)確性受冰湖周?chē)匚锏膹?fù)雜性影響，是否能在其它冰川用于冰湖提取有待驗(yàn)證；另外值得說(shuō)明的是，由于受實(shí)驗(yàn)區(qū)域在該時(shí)間段內(nèi)所獲取的光學(xué)影像受云霧遮擋較為嚴(yán)重所限，很難于每年相同或相近月份收集到符合要求的數(shù)據(jù)。因此，只能通過(guò)比較年際帶有一定時(shí)間跨度的冰湖面積進(jìn)行后續(xù)分析。在今后的研究中若能獲取到時(shí)間跨度較小的高質(zhì)量的影像數(shù)據(jù)，將會(huì)考慮如季節(jié)性問(wèn)題等對(duì)冰湖提取所帶來(lái)的影響，從而進(jìn)行更加精確細(xì)致的結(jié)果分析。

志 謝此次所用光學(xué)數(shù)據(jù)來(lái)自美國(guó)NASA的陸地衛(wèi)星（Landsat 5/8），溫度數(shù)據(jù)為中國(guó)氣象局提供，在此表示衷心的感謝！