波密县冰川泥石流遥感判识与提取

刘彦培

(河南理工大学测绘与国土信息工程学院，河南 焦作 454000)

1 概述

冰川泥石流是常发生在山区沟谷的一种自然灾害现象，是由强降雨、冰湖溃决、冰川积雪融水或者冰崩雪崩引发的一种特殊洪流，夹带有大量块石和泥沙。我国泥石流地质灾害分布广泛，以青藏高原、云贵高原、四川盆地和台湾等地的泥石流灾害最为严重[1]。波密县自然条件复杂，生态环境极其脆弱，对外界的干扰十分敏感，易引发多类地质灾害或隐患。因此，解决波密县突出的冰川泥石流灾害问题是国家基础建设的重要需求，对于开展地质灾害监测、评估、预警工作具有重要的现实意义。

冰川泥石流所处地势险要，人迹罕至，常规的手段难以开展相关工作，而3S技术[2]以及合成孔径雷达技术为冰川泥石流研究提供了新方法。目前西藏地区冰川泥石流研究主要集中在藏东南地区、川藏公路沿线及帕隆藏布江流域等。冰川泥石流研究始于施雅风等[3]及杜榕桓等[4]，此后，国内外学者针对冰川泥石流[5-6]的形成机理[7-8]、流域特征等[9-10]方面展开了深入研究，发现冰川泥石流的发生与冰川变化有关[11-12]，但冰川泥石流起动能否用暴雨泥石流理论[13-14]解释仍需探讨。关于冰川泥石流提取模型与方法，国内外运用遥感技术从光学遥感到雷达遥感，从低分辨率到高分辨率，从单一数据源到多源数据，结合地形特征[15]、气象条件等提取冰川泥石流灾害信息[16-23]及分类提取物源信息[24-25]，取得了巨大进展。但大部分仍停留在目视解译、监督与非监督等分类方法，自动提取的模型方法仍然薄弱。

本文选取波密县2015年3月、2015年7月、2015年11月三期Landsat 8 OLI影像为数据源，引入空间判识度与时空判识度，建立基于不同时相遥感影像的冰川泥石流范围的判识方法，解决因冰川面积季节性变化及山体阴影影响的冰川泥石流范围难判定的问题；提取冰川泥石流流域内各地物纹理变化特征，采用二元逻辑回归分析算法建立冰川泥石流初步判识模型及二次判定规则，剔除其他类型泥石流，得到冰川型泥石流流域边界；通过空间叠置分析，实现波密县冰川泥石流范围与类型的完整提取。该研究可以为冰川泥石流的遥感判识与信息提取研究提供方法或技术借鉴，为波密县冰川泥石流灾害的防灾减灾等工作提供数据支撑。

2 研究区与数据

2.1 研究区概况

波密县位于西藏东南部、念青唐古拉山与喜马拉雅山交界处，是318国道川藏线上重要的交通枢纽和商贸重镇。全县面积1.67 km2，人口3.8万人，平均海拔3 300 m，境内矿产资源40余种，国家保护野生动物80余种，植物资源400余种，分布有帕隆藏布峡谷、贡布省级自然保护区等多个自然保护区，属于特殊的高寒生态敏感区[26]。波密县地质构造活跃，地形起伏较大，地形条件复杂，是我国地质灾害相对严重的地区。印度洋和太平洋暖湿气流作用造就了国内最大面积的海洋性冰川区，区域内松散堆积物特别发育，为滑坡、泥石流等地质灾害提供了丰富的物源，研究区平均气温在8.7 ℃[27]，年均降雨量以+3 mm/(10 a)的趋势增加[28]，而近40 a降雪量呈逐年降低的趋势[29]，造成冰川处于强烈的消融阶段[30]，为冰川泥石流提供了充沛的水源及动力条件。充足的水源、丰富的物源加上波密县的地势、气候条件，使得该区域内冰川泥石流具有突发性强、链生性及破坏力强等特点[31-32]，成为该地区的主要地质灾害类型。

2.2 数据源与预处理

2.2.1 数据源

本文选择波密县2015-03，2015-07，2015-11三期共12景Landsat 8 OLI影像为基础数据，云量覆盖均在20%以下，削弱云的影响。辅助数据为分辨率为30 m的DEM。Landsat数据及DEM均下载自地理空间数据云网站(http://www.gscloud.cn/)，数据具体信息如表1所示。

表1 Landsat 8 OLI数据信息表

2.2.2 数据预处理

Landsat 8 OLI影像L1T级产品已做过几何校正，消除了因系统辐射畸变、地形起伏而引起的成像误差，可以直接使用。因此，利用ENVI 5.3遥感图像处理平台Radiometric Calibration模块的辐射定标工具和FLAASH模块的大气校正模型对Landsat 8 OLI数据进行辐射定标和大气校正，以及影像镶嵌、影像裁剪等预处理。

3 冰川泥石流遥感判识

3.1 随机森林分类

3.1.1 影像特征

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

选取Landsat 8 OLI影像的第5波段和第7波段计算冰川泥石流指数(GDFI)，计算公式见式(1)：

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

(1)

其中，B5为Landsat 8 OLI影像的近红外波段；B7为Landsat 8 OLI影像的热红外波段。

在冰川泥石流指数特征图(如图1所示)中，冰川泥石流流域呈灰黑色，其中形成区冰碛物(椭圆内)呈现黑色，易于判识。冰川泥石流流域特征明显的还有主成分分析(PCA)后的1，2，3波段假彩色(Red：band1，Green：band2，Blue：band3)显示(见图2)，在PCA特征图中，冰碛物(椭圆内)呈黄色，主要分布在形成区，少量分布在流通区沟道及坡度较缓的堆积扇。

2)形成区。

对Landsat 8影像进行最小噪声分离(MNF)处理，得到MNF特征图(见图3(a))，形成区(上圆圈)呈亮白色，与周围植被(下圆圈)的灰黑色区别明显；对Landsat 8影像进行缨帽变换(K-T变换)得到亮度指数(见图3(b))、绿度指数(见图3(c))、湿度指数(见图3(d))，亮度指数特征图中，形成区冰川(椭圆)呈亮白色；湿度指数特征图中，由于形成区松散物质及冰碛物(上圆圈)含水量较高而呈现黑色，与周围植被及裸地(下圆圈)的灰色差异明显；绿度指数特征图中，植被呈亮白色，易于判识。

3)流通区。

冰川泥石流在狭窄的流通区侵蚀沟床、侧壁，发生前后沟道内纹理粗糙、植被损毁严重，与冰川泥石流其他区域不同。归一化植被指数(NDVI)反映植被在近红外波段、红光波段反射率差异，用来监测植被；利用灰度共生矩阵提取影像均值、方差、协同性等8个纹理特征，用来表征纹理。

4)堆积区。

对Landsat 8影像进行独立主成分分析(ICA)得到ICA特征图(见图4)，冰川泥石流堆积区(上椭圆内)整体呈黄色，因少量表碛物覆盖而在色彩上掺杂些许紫色，部分新堆积区因无植被覆盖在ICA特征图中呈现白色，与呈淡紫色的河漫滩(中下圆圈)差异明显。

综上，选取冰川泥石流指数特征、MNF特征波段、PCA特征波段、ICA特征波段组成多特征波段图像用于随机森林地物分类，减少参与分类的波段数量，避免维数灾难。

3.1.2 地物分类

采用随机森林算法将三期多特征波段影像中的地物分为云、阴影、冰川、冰碛物、水体、植被、裸地、居民地、河漫滩及冰川泥石流形成区、流通区、堆积区，分类结果见图5a～图7a。以波密县古乡沟(如图5b～图7b所示)为例，2015-03提取的冰川泥石流形成区冰川覆盖范围最广，流通区、堆积区相对集中；随着冰川融化，形成区内冰碛物逐渐增多且集中在冰舌前端，2015-11影像中提取的形成区冰碛物最多，且形成区至堆积区有明显的破坏迹象，流通区两侧植被受损明显。

3.1.3 精度评价

采用基于样本的混淆矩阵算法对三期影像分类结果进行评价，得到三期影像的总体精度与Kappa系数(见表2)。由表2可知，三期影像分类精度均达到83%以上，分类精度较好。

表2 三期影像分类精度评价表

3.2 冰川泥石流范围判识

空间判识度：单期影像分类结果中，3×3窗口内，冰川泥石流像元个数占窗口内像元总个数的比值。

时空判识度：多期影像中，同一位置像元的空间判识度的平均值。计算时空判识度过程中，若遥感影像中有阴影，则先计算有阴影影像的时空判识度，再与无阴影影像空间判识度求均值。

三期影像空间判识度求平均得到时空判识度，利用时空判识度提取冰川泥石流范围信息，充分利用了遥感影像的时相特征，一定程度上削弱冰雪季节性变化的影响。图8～图10为波密县及县内古乡沟流域三期影像的空间判识度图，图中色调由黑至白表示该处像元属于冰川泥石流的可能性越大。由图8，图10可知：影像中形成区所含的冰川泥石流信息较少，且流通区(方框内)受阴影的影响呈灰黑色，表明该区域属冰川泥石流的可能性较小；2105-07(如图9所示)影像中形成区冰川泥石流信息最为丰富，且流通区(方框内)受阴影影响小，色调呈灰白色，冰川泥石流整体特征明显。利用三期空间判识度计算时空判识度(见图11)，原空间判识度图中有阴影覆盖的流通区部分区域色调由灰黑色转而呈灰白色，表明通过计算时空判识度削弱了阴影对冰川泥石流范围判识的影响，使得冰川泥石流信息表现更丰富。

3.3 冰川泥石流沟道判识

3.3.1 纹理变化特征

本文以2015-07与2015-11两期影像(如图12a，图12b所示)为例计算纹理变化特征，图12c，图12d为单期影像的方差纹理特征，图12e为两期影像方差的变化特征。与单期影像纹理特征相比，纹理变化特征图更突出冰川泥石流形成区的冰碛物等物源(上方框)、流通区沟道(左下方框)、沟道两侧植被、裸地等地物(右下方框)；而单期纹理特征图中，除形成区冰川、冰雪易于判识外，形成区、流通区、堆积区及各区域内植被、裸地等地物特征不明显，难以区分。因此，两期影像的纹理变化特征图(如图13所示)比单期影像纹理特征图更易于地物判识及分类提取。

3.3.2 沟道判识模型

1)初步判识模型。基于随机森林地物分类结果，选取冰川泥石流流域样本33个及其他类型泥石流流域样本11个，利用流域单元统计冰川泥石流前后纹理变化最为明显的冰川、冰碛物、植被、裸地及冰川泥石流形成区、流通区、堆积区7类地物的纹理变化均值，采用二元逻辑回归分析建立冰川泥石流沟道初步判识模型，见式(2)：

Logistic(Y)=86.130×m(冰川)-
23.943×m(冰碛物)-638.006×m(裸地)-
136.183×m(植被)+53.142×m(形成区)+
39.646×m(流通区)+133.022×m(堆积区)

(2)

其中，m()为某地物的纹理变化均值；Y为0或1，1表示该流域是冰川型泥石流，0表示该流域是其他类型泥石流。

2)二次判定规则。部分冰川泥石流因形成区常年被冰川、冰雪覆盖，冰川泥石流发生前后形成区内冰川纹理变化不大，几乎为0，而被错分成同样没有冰川纹理变化特征的其他类型泥石流。因此，需要对流通区和堆积区进行二次判定。

根据不同类型泥石流沟道内各地物纹理变化规律(如表3所示)可知冰川泥石流及其各类地物纹理变化特征值范围在-1.36～10.55，而其他类型泥石流及其各类地物纹理变化特征值范围在-0.67～6.37，冰川泥石流纹理变化特征值范围明显大于其他类型泥石流的纹理变化特征值范围，尤其冰碛物的纹理变化最为明显；另外，非冰川型泥石流虽没有冰川纹理变化特征，但同样没有明显冰川纹理变化特征的冰湖溃决型泥石流易与之混淆；此外，植被、冰川泥石流形成区的纹理变化特征明显区别于其他类型泥石流。依据表3，针对冰川纹理变化为0的流域单元集，以冰碛物、形成区、冰川的纹理变化特征建立冰川泥石流二次判定规则(见图14)。

表3 不同类型泥石流流域内各地物纹理变化范围

4 冰川泥石流提取

4.1 冰川泥石流范围提取

对冰川泥石流时空判识度图进行密度分割，通过不断尝试及目视提取效果，获取最佳提取阈值7.54提取冰川泥石流范围(见图15)。

4.2 冰川泥石流沟道提取

利用式(2)冰川泥石流沟道初步判识模型对研究区内流域单元类型进行初步判识(见图16)，利用冰川泥石流二次判定规则进行二次判定，得到最终的冰川泥石流沟道判识结果，与提取的冰川泥石流范围信息进行叠置，得到最终的冰川泥石流信息(见图17)。

4.3 精度评价

任意选取波密县三条冰川泥石流(见图18)，目视解译其流域范围，分别计算随机森林分类得到的冰川泥石流及本文方法提取的冰川泥石流面积(见表4)。

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

由表4可知，2015-03期影像分类提取的冰川泥石流面积与目视解译的面积相差最大，误差最大为9.9 km2；而2015-07期夏季影像中冰川泥石流受冰雪覆盖影响较小，提取较为完整，与目视解译面积最为接近；但由于影像中仍存在部分阴影区域，导致单期影像分类所获取的冰川泥石流面积与目视解译面积存在差异。使用本文方法提取的冰川泥石流面积，比单期影像分类结果更接近于目视解译的面积，且误差较小。

5 结语

本文针对影像中山体阴影、冰川冰雪季节变化以及其他类型泥石流沟道干扰的问题，提出冰川泥石流范围与冰川泥石流类型判识与提取方法，主要结论有：

1)利用特征波段组合的多波段图像进行随机森林地物分类，减少波段冗余，分类结果表明单期影像精度均在83%以上，分类效果较好，为后续冰川泥石流遥感判识与信息提取奠定良好的基础。

2)基于2015-03，2015-07，2015-11三期影像地物分类结果，充分利用波密县三期影像的时相特征，分别计算冰川泥石流空间判识度与时空判识度，通过阈值分割及不断尝试分割阈值，得到最佳分割阈值7.54，提取冰川泥石流范围信息。

3)基于冰川泥石流前(2015-07)、冰川泥石流后(2015-11)两期影像纹理变化特征，结合地物分类结果，统计不同类型泥石流流域内各地物的纹理变化规律，据此采用二元逻辑回归方法建立波密县冰川泥石流初步判识模型及冰川泥石流二次判定规则，剔除其他类型泥石流的干扰，准确获取冰川型泥石流流域边界信息。

4)将冰川泥石流范围信息与冰川泥石流边界信息进行叠置分析，实现冰川泥石流信息完整提取。对比单期影像提取冰川泥石流范围、本文方法提取冰川泥石流范围与目视解译范围，发现本文方法提取的冰川泥石流比单期影像提取结果更接近于目视解译的面积，且误差相较最小。