基于NDWI-NDSI組合閾值法的布加崗日冰湖提取及其變化分析

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

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

1 引 言

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

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

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

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

如圖1所示，布加崗日冰川位于距離川藏北線G317 國道不遠(yuǎn)的西藏自治區(qū)那曲市索縣、巴青縣及昌都市丁青縣的交界處，與唐古拉山脈一致為東西走向。其海拔較高，有16 座海拔超過6000 m的高峰，最高峰海拔為6328 m。受來自印度洋的西南季風(fēng)作用，年降水量達(dá)到600—700 mm，高聳的地勢和豐富的降水使這里成為較大的山岳冰川發(fā)育集中地（Liu等，2016；王寧練和丁良福，2002；王聰強等，2016）。據(jù)光學(xué)影像目視解譯，該冰川發(fā)育面積較為明顯的冰湖約為12 個，它們分別對應(yīng)編號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ù)源為美國NASA 的陸地衛(wèi)星（Landsat）TM/OLI，衛(wèi)星波段介紹如表1 所示。此次實驗所使用的Landsat5/8 衛(wèi)星波段有：Green 波段、NIR波段、SWIR1波段，所選多光譜波段分辨率均為30 m。基于光學(xué)影像成功提取冰湖面積的前提條件是影像獲取時間段內(nèi)研究區(qū)上空無云或云霧較少、且冰湖處于非結(jié)冰期，而布加崗日所在地區(qū)的夏季多數(shù)影像數(shù)據(jù)被云霧覆蓋嚴(yán)重?zé)o法使用，同時部分冬季影像由于冰湖結(jié)冰或被雪覆蓋也導(dǎo)致數(shù)據(jù)不符合要求?？紤]以上因素，本文篩選出1988 年—2018 年間共16 景質(zhì)量較好且覆蓋布加崗日冰川的Landsat（TM/OLI）影像。對所選影像首先裁剪處理，并進行輻射定標(biāo)和大氣校正，以消除或減弱由傳感器和大氣帶來的各種誤差，避免其對后期處理結(jié)果的精確性產(chǎn)生影響。

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

3 研究方法

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

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

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

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

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

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

表2 冰湖提取方法精度評價Table 2 Accuracy evaluation of glacial lake extraction method

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

4 結(jié)果分析

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

4.1 單個冰湖變化分析

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

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

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

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

表3 冰湖面積分級表Table 3 Classification of glacial lakes area

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

4.2 整體面積變化分析

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

如圖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 年布加崗日冰川的冰湖面積變化呈慢→較快→很快→較快趨勢，雖然面積變化不能準(zhǔn)確反映冰川的消融量及湖水變化量，但其結(jié)果仍然可以表明布加崗日冰川退縮明顯。結(jié)合圖9可知，冰湖總面積變化與溫度呈現(xiàn)正相關(guān)特性（該地區(qū)溫度數(shù)據(jù)為中國氣象局提供），表明該地區(qū)逐年氣溫的不斷升高正是冰川消融過快的主要原因；其次，經(jīng)統(tǒng)計發(fā)現(xiàn)，布加崗日冰川地區(qū)近30 年來平均氣溫在0 ℃以上的年份共7 個，而1998 年—2013 年（紅色虛線框內(nèi)）冰湖面積增長速率之所以最快，與這段時間內(nèi)年均氣溫在0 ℃以上的年數(shù)較多緊密相關(guān)（5 個年份）；2013 年—2016 年（橙色虛線框內(nèi)）增長速率變化次之，其原因為年均氣溫在0℃以上的年數(shù)相應(yīng)較少（2個年份）。

圖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é) 論

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

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

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

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

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

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