劉彥培

(河南理工大學(xué)測繪與國土信息工程學(xué)院，河南 焦作 454000)

1 概述

冰川泥石流是常發(fā)生在山區(qū)溝谷的一種自然災(zāi)害現(xiàn)象，是由強降雨、冰湖潰決、冰川積雪融水或者冰崩雪崩引發(fā)的一種特殊洪流，夾帶有大量塊石和泥沙。我國泥石流地質(zhì)災(zāi)害分布廣泛，以青藏高原、云貴高原、四川盆地和臺灣等地的泥石流災(zāi)害最為嚴(yán)重[1]。波密縣自然條件復(fù)雜，生態(tài)環(huán)境極其脆弱，對外界的干擾十分敏感，易引發(fā)多類地質(zhì)災(zāi)害或隱患。因此，解決波密縣突出的冰川泥石流災(zāi)害問題是國家基礎(chǔ)建設(shè)的重要需求，對于開展地質(zhì)災(zāi)害監(jiān)測、評估、預(yù)警工作具有重要的現(xiàn)實意義。

冰川泥石流所處地勢險要，人跡罕至，常規(guī)的手段難以開展相關(guān)工作，而3S技術(shù)[2]以及合成孔徑雷達技術(shù)為冰川泥石流研究提供了新方法。目前西藏地區(qū)冰川泥石流研究主要集中在藏東南地區(qū)、川藏公路沿線及帕隆藏布江流域等。冰川泥石流研究始于施雅風(fēng)等[3]及杜榕桓等[4]，此后，國內(nèi)外學(xué)者針對冰川泥石流[5-6]的形成機理[7-8]、流域特征等[9-10]方面展開了深入研究，發(fā)現(xiàn)冰川泥石流的發(fā)生與冰川變化有關(guān)[11-12]，但冰川泥石流起動能否用暴雨泥石流理論[13-14]解釋仍需探討。關(guān)于冰川泥石流提取模型與方法，國內(nèi)外運用遙感技術(shù)從光學(xué)遙感到雷達遙感，從低分辨率到高分辨率，從單一數(shù)據(jù)源到多源數(shù)據(jù)，結(jié)合地形特征[15]、氣象條件等提取冰川泥石流災(zāi)害信息[16-23]及分類提取物源信息[24-25]，取得了巨大進展。但大部分仍停留在目視解譯、監(jiān)督與非監(jiān)督等分類方法，自動提取的模型方法仍然薄弱。

本文選取波密縣2015年3月、2015年7月、2015年11月三期Landsat 8 OLI影像為數(shù)據(jù)源，引入空間判識度與時空判識度，建立基于不同時相遙感影像的冰川泥石流范圍的判識方法，解決因冰川面積季節(jié)性變化及山體陰影影響的冰川泥石流范圍難判定的問題；提取冰川泥石流流域內(nèi)各地物紋理變化特征，采用二元邏輯回歸分析算法建立冰川泥石流初步判識模型及二次判定規(guī)則，剔除其他類型泥石流，得到冰川型泥石流流域邊界；通過空間疊置分析，實現(xiàn)波密縣冰川泥石流范圍與類型的完整提取。該研究可以為冰川泥石流的遙感判識與信息提取研究提供方法或技術(shù)借鑒，為波密縣冰川泥石流災(zāi)害的防災(zāi)減災(zāi)等工作提供數(shù)據(jù)支撐。

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

2.1 研究區(qū)概況

波密縣位于西藏東南部、念青唐古拉山與喜馬拉雅山交界處，是318國道川藏線上重要的交通樞紐和商貿(mào)重鎮(zhèn)。全縣面積1.67 km2，人口3.8萬人，平均海拔3 300 m，境內(nèi)礦產(chǎn)資源40余種，國家保護野生動物80余種，植物資源400余種，分布有帕隆藏布峽谷、貢布省級自然保護區(qū)等多個自然保護區(qū)，屬于特殊的高寒生態(tài)敏感區(qū)[26]。波密縣地質(zhì)構(gòu)造活躍，地形起伏較大，地形條件復(fù)雜，是我國地質(zhì)災(zāi)害相對嚴(yán)重的地區(qū)。印度洋和太平洋暖濕氣流作用造就了國內(nèi)最大面積的海洋性冰川區(qū)，區(qū)域內(nèi)松散堆積物特別發(fā)育，為滑坡、泥石流等地質(zhì)災(zāi)害提供了豐富的物源，研究區(qū)平均氣溫在8.7 ℃[27]，年均降雨量以+3 mm/(10 a)的趨勢增加[28]，而近40 a降雪量呈逐年降低的趨勢[29]，造成冰川處于強烈的消融階段[30]，為冰川泥石流提供了充沛的水源及動力條件。充足的水源、豐富的物源加上波密縣的地勢、氣候條件，使得該區(qū)域內(nèi)冰川泥石流具有突發(fā)性強、鏈生性及破壞力強等特點[31-32]，成為該地區(qū)的主要地質(zhì)災(zāi)害類型。

2.2 數(shù)據(jù)源與預(yù)處理

2.2.1 數(shù)據(jù)源

本文選擇波密縣2015-03，2015-07，2015-11三期共12景Landsat 8 OLI影像為基礎(chǔ)數(shù)據(jù)，云量覆蓋均在20%以下，削弱云的影響。輔助數(shù)據(jù)為分辨率為30 m的DEM。Landsat數(shù)據(jù)及DEM均下載自地理空間數(shù)據(jù)云網(wǎng)站(http://www.gscloud.cn/)，數(shù)據(jù)具體信息如表1所示。

表1 Landsat 8 OLI數(shù)據(jù)信息表

2.2.2 數(shù)據(jù)預(yù)處理

Landsat 8 OLI影像L1T級產(chǎn)品已做過幾何校正，消除了因系統(tǒng)輻射畸變、地形起伏而引起的成像誤差，可以直接使用。因此，利用ENVI 5.3遙感圖像處理平臺Radiometric Calibration模塊的輻射定標(biāo)工具和FLAASH模塊的大氣校正模型對Landsat 8 OLI數(shù)據(jù)進行輻射定標(biāo)和大氣校正，以及影像鑲嵌、影像裁剪等預(yù)處理。

3 冰川泥石流遙感判識

3.1 隨機森林分類

3.1.1 影像特征

1)冰川泥石流流域特征。

選取Landsat 8 OLI影像的第5波段和第7波段計算冰川泥石流指數(shù)(GDFI)，計算公式見式(1)：

GDFI=(B5-B7)/(B5+B7)

(1)

其中，B5為Landsat 8 OLI影像的近紅外波段；B7為Landsat 8 OLI影像的熱紅外波段。

在冰川泥石流指數(shù)特征圖(如圖1所示)中，冰川泥石流流域呈灰黑色，其中形成區(qū)冰磧物(橢圓內(nèi))呈現(xiàn)黑色，易于判識。冰川泥石流流域特征明顯的還有主成分分析(PCA)后的1，2，3波段假彩色(Red：band1，Green：band2，Blue：band3)顯示(見圖2)，在PCA特征圖中，冰磧物(橢圓內(nèi))呈黃色，主要分布在形成區(qū)，少量分布在流通區(qū)溝道及坡度較緩的堆積扇。

2)形成區(qū)。

對Landsat 8影像進行最小噪聲分離(MNF)處理，得到MNF特征圖(見圖3(a))，形成區(qū)(上圓圈)呈亮白色，與周圍植被(下圓圈)的灰黑色區(qū)別明顯；對Landsat 8影像進行纓帽變換(K-T變換)得到亮度指數(shù)(見圖3(b))、綠度指數(shù)(見圖3(c))、濕度指數(shù)(見圖3(d))，亮度指數(shù)特征圖中，形成區(qū)冰川(橢圓)呈亮白色；濕度指數(shù)特征圖中，由于形成區(qū)松散物質(zhì)及冰磧物(上圓圈)含水量較高而呈現(xiàn)黑色，與周圍植被及裸地(下圓圈)的灰色差異明顯；綠度指數(shù)特征圖中，植被呈亮白色，易于判識。

3)流通區(qū)。

冰川泥石流在狹窄的流通區(qū)侵蝕溝床、側(cè)壁，發(fā)生前后溝道內(nèi)紋理粗糙、植被損毀嚴(yán)重，與冰川泥石流其他區(qū)域不同。歸一化植被指數(shù)(NDVI)反映植被在近紅外波段、紅光波段反射率差異，用來監(jiān)測植被；利用灰度共生矩陣提取影像均值、方差、協(xié)同性等8個紋理特征，用來表征紋理。

4)堆積區(qū)。

對Landsat 8影像進行獨立主成分分析(ICA)得到ICA特征圖(見圖4)，冰川泥石流堆積區(qū)(上橢圓內(nèi))整體呈黃色，因少量表磧物覆蓋而在色彩上摻雜些許紫色，部分新堆積區(qū)因無植被覆蓋在ICA特征圖中呈現(xiàn)白色，與呈淡紫色的河漫灘(中下圓圈)差異明顯。

綜上，選取冰川泥石流指數(shù)特征、MNF特征波段、PCA特征波段、ICA特征波段組成多特征波段圖像用于隨機森林地物分類，減少參與分類的波段數(shù)量，避免維數(shù)災(zāi)難。

3.1.2 地物分類

采用隨機森林算法將三期多特征波段影像中的地物分為云、陰影、冰川、冰磧物、水體、植被、裸地、居民地、河漫灘及冰川泥石流形成區(qū)、流通區(qū)、堆積區(qū)，分類結(jié)果見圖5a～圖7a。以波密縣古鄉(xiāng)溝(如圖5b～圖7b所示)為例，2015-03提取的冰川泥石流形成區(qū)冰川覆蓋范圍最廣，流通區(qū)、堆積區(qū)相對集中；隨著冰川融化，形成區(qū)內(nèi)冰磧物逐漸增多且集中在冰舌前端，2015-11影像中提取的形成區(qū)冰磧物最多，且形成區(qū)至堆積區(qū)有明顯的破壞跡象，流通區(qū)兩側(cè)植被受損明顯。

3.1.3 精度評價

采用基于樣本的混淆矩陣算法對三期影像分類結(jié)果進行評價，得到三期影像的總體精度與Kappa系數(shù)(見表2)。由表2可知，三期影像分類精度均達到83%以上，分類精度較好。

表2 三期影像分類精度評價表

3.2 冰川泥石流范圍判識

空間判識度：單期影像分類結(jié)果中，3×3窗口內(nèi)，冰川泥石流像元個數(shù)占窗口內(nèi)像元總個數(shù)的比值。

時空判識度：多期影像中，同一位置像元的空間判識度的平均值。計算時空判識度過程中，若遙感影像中有陰影，則先計算有陰影影像的時空判識度，再與無陰影影像空間判識度求均值。

三期影像空間判識度求平均得到時空判識度，利用時空判識度提取冰川泥石流范圍信息，充分利用了遙感影像的時相特征，一定程度上削弱冰雪季節(jié)性變化的影響。圖8～圖10為波密縣及縣內(nèi)古鄉(xiāng)溝流域三期影像的空間判識度圖，圖中色調(diào)由黑至白表示該處像元屬于冰川泥石流的可能性越大。由圖8，圖10可知：影像中形成區(qū)所含的冰川泥石流信息較少，且流通區(qū)(方框內(nèi))受陰影的影響呈灰黑色，表明該區(qū)域?qū)俦嗍鞯目赡苄暂^?。?105-07(如圖9所示)影像中形成區(qū)冰川泥石流信息最為豐富，且流通區(qū)(方框內(nèi))受陰影影響小，色調(diào)呈灰白色，冰川泥石流整體特征明顯。利用三期空間判識度計算時空判識度(見圖11)，原空間判識度圖中有陰影覆蓋的流通區(qū)部分區(qū)域色調(diào)由灰黑色轉(zhuǎn)而呈灰白色，表明通過計算時空判識度削弱了陰影對冰川泥石流范圍判識的影響，使得冰川泥石流信息表現(xiàn)更豐富。

3.3 冰川泥石流溝道判識

3.3.1 紋理變化特征

本文以2015-07與2015-11兩期影像(如圖12a，圖12b所示)為例計算紋理變化特征，圖12c，圖12d為單期影像的方差紋理特征，圖12e為兩期影像方差的變化特征。與單期影像紋理特征相比，紋理變化特征圖更突出冰川泥石流形成區(qū)的冰磧物等物源(上方框)、流通區(qū)溝道(左下方框)、溝道兩側(cè)植被、裸地等地物(右下方框)；而單期紋理特征圖中，除形成區(qū)冰川、冰雪易于判識外，形成區(qū)、流通區(qū)、堆積區(qū)及各區(qū)域內(nèi)植被、裸地等地物特征不明顯，難以區(qū)分。因此，兩期影像的紋理變化特征圖(如圖13所示)比單期影像紋理特征圖更易于地物判識及分類提取。

3.3.2 溝道判識模型

1)初步判識模型?；陔S機森林地物分類結(jié)果，選取冰川泥石流流域樣本33個及其他類型泥石流流域樣本11個，利用流域單元統(tǒng)計冰川泥石流前后紋理變化最為明顯的冰川、冰磧物、植被、裸地及冰川泥石流形成區(qū)、流通區(qū)、堆積區(qū)7類地物的紋理變化均值，采用二元邏輯回歸分析建立冰川泥石流溝道初步判識模型，見式(2)：

Logistic(Y)=86.130×m(冰川)-
23.943×m(冰磧物)-638.006×m(裸地)-
136.183×m(植被)+53.142×m(形成區(qū))+
39.646×m(流通區(qū))+133.022×m(堆積區(qū))

(2)

其中，m()為某地物的紋理變化均值；Y為0或1，1表示該流域是冰川型泥石流，0表示該流域是其他類型泥石流。

2)二次判定規(guī)則。部分冰川泥石流因形成區(qū)常年被冰川、冰雪覆蓋，冰川泥石流發(fā)生前后形成區(qū)內(nèi)冰川紋理變化不大，幾乎為0，而被錯分成同樣沒有冰川紋理變化特征的其他類型泥石流。因此，需要對流通區(qū)和堆積區(qū)進行二次判定。

根據(jù)不同類型泥石流溝道內(nèi)各地物紋理變化規(guī)律(如表3所示)可知冰川泥石流及其各類地物紋理變化特征值范圍在-1.36～10.55，而其他類型泥石流及其各類地物紋理變化特征值范圍在-0.67～6.37，冰川泥石流紋理變化特征值范圍明顯大于其他類型泥石流的紋理變化特征值范圍，尤其冰磧物的紋理變化最為明顯；另外，非冰川型泥石流雖沒有冰川紋理變化特征，但同樣沒有明顯冰川紋理變化特征的冰湖潰決型泥石流易與之混淆；此外，植被、冰川泥石流形成區(qū)的紋理變化特征明顯區(qū)別于其他類型泥石流。依據(jù)表3，針對冰川紋理變化為0的流域單元集，以冰磧物、形成區(qū)、冰川的紋理變化特征建立冰川泥石流二次判定規(guī)則(見圖14)。

表3 不同類型泥石流流域內(nèi)各地物紋理變化范圍

4 冰川泥石流提取

4.1 冰川泥石流范圍提取

對冰川泥石流時空判識度圖進行密度分割，通過不斷嘗試及目視提取效果，獲取最佳提取閾值7.54提取冰川泥石流范圍(見圖15)。

4.2 冰川泥石流溝道提取

利用式(2)冰川泥石流溝道初步判識模型對研究區(qū)內(nèi)流域單元類型進行初步判識(見圖16)，利用冰川泥石流二次判定規(guī)則進行二次判定，得到最終的冰川泥石流溝道判識結(jié)果，與提取的冰川泥石流范圍信息進行疊置，得到最終的冰川泥石流信息(見圖17)。

4.3 精度評價

任意選取波密縣三條冰川泥石流(見圖18)，目視解譯其流域范圍，分別計算隨機森林分類得到的冰川泥石流及本文方法提取的冰川泥石流面積(見表4)。

表4 不同方法提取的冰川泥石流面積

由表4可知，2015-03期影像分類提取的冰川泥石流面積與目視解譯的面積相差最大，誤差最大為9.9 km2；而2015-07期夏季影像中冰川泥石流受冰雪覆蓋影響較小，提取較為完整，與目視解譯面積最為接近；但由于影像中仍存在部分陰影區(qū)域，導(dǎo)致單期影像分類所獲取的冰川泥石流面積與目視解譯面積存在差異。使用本文方法提取的冰川泥石流面積，比單期影像分類結(jié)果更接近于目視解譯的面積，且誤差較小。

5 結(jié)語

本文針對影像中山體陰影、冰川冰雪季節(jié)變化以及其他類型泥石流溝道干擾的問題，提出冰川泥石流范圍與冰川泥石流類型判識與提取方法，主要結(jié)論有：

1)利用特征波段組合的多波段圖像進行隨機森林地物分類，減少波段冗余，分類結(jié)果表明單期影像精度均在83%以上，分類效果較好，為后續(xù)冰川泥石流遙感判識與信息提取奠定良好的基礎(chǔ)。

2)基于2015-03，2015-07，2015-11三期影像地物分類結(jié)果，充分利用波密縣三期影像的時相特征，分別計算冰川泥石流空間判識度與時空判識度，通過閾值分割及不斷嘗試分割閾值，得到最佳分割閾值7.54，提取冰川泥石流范圍信息。

3)基于冰川泥石流前(2015-07)、冰川泥石流后(2015-11)兩期影像紋理變化特征，結(jié)合地物分類結(jié)果，統(tǒng)計不同類型泥石流流域內(nèi)各地物的紋理變化規(guī)律，據(jù)此采用二元邏輯回歸方法建立波密縣冰川泥石流初步判識模型及冰川泥石流二次判定規(guī)則，剔除其他類型泥石流的干擾，準(zhǔn)確獲取冰川型泥石流流域邊界信息。

4)將冰川泥石流范圍信息與冰川泥石流邊界信息進行疊置分析，實現(xiàn)冰川泥石流信息完整提取。對比單期影像提取冰川泥石流范圍、本文方法提取冰川泥石流范圍與目視解譯范圍，發(fā)現(xiàn)本文方法提取的冰川泥石流比單期影像提取結(jié)果更接近于目視解譯的面積，且誤差相較最小。