基于多传感器遥感数据和改进的单类随机森林的滑坡提取

刘君茹，吕悦琪，李培军

（1.山西省住房和城乡建设厅，山西 太原 030013；2.北京大学地球与空间科学学院遥感所，北京 100871）

滑坡是土壤、岩石和有机物质在重力作用下的下坡运动，可由自然或人为触发因素引起，如地震、降雨和工程活动。滑坡是一种世界性重大地质灾害，常导致人员伤亡、巨大财产损失和生态环境破坏[1-2]。因此，在抢险救灾和灾害评估中，需及时、准确的滑坡空间分布信息。遥感为提取滑坡或进行滑坡提取提供了强有力的手段。传统滑坡提取主要是航空相片的目视解译[3]。目视解译费时费力，不适合快速、大面积的滑坡提取。来自卫星和航空遥感的数据，特别是具详细空间细节的高分辨率数据的大量获取，使准确地监测滑坡成为可能[2,4-6]。近年来，有不少研究提出一些利用高分辨率遥感数据自动或半自动提取滑坡的方法[4,7-12]。这些方法中，基于对象的图像分析方法优于基于像元的方法[3,13-15]。

从高分辨率图像提取滑坡的主要问题是滑坡的大小、形状和地表特征等的异质性，及滑坡和其它土地覆盖类型间的光谱相似性[16]。为解决这些问题，现有研究主要利用不同数据源的多种特征。如，Li等使用光谱和空间上下文信息从航空相片提取滑坡[10-11]；Martha 等综合利用光谱特征、植被指数、形状、纹理和地形特征提取印度喜马拉雅地区的滑坡[17]；Stumpf 和Kerle 提出一种利用面向对象的图像分析和随机森林算法从高分辨率卫星图像和航空照片提取滑坡的方法[6]。滑坡提取中使用多种特征，一般包括光谱特征、植被指数特征、空间特征和地形特征。常用空间特征包括纹理、形状和上下文信息[5,6,18]，以区分滑坡和具有相似光谱特征的土地覆盖类型。由于滑坡区土地覆盖类型的复杂性，目前还没有一种被广泛接受的单一空间特征可用于滑坡提取。因此，需探索新的、更有效的空间特征。考虑到滑坡前后地表空间同质性和变异性可能发生显著变化，在滑坡提取中采用量化空间同质性和变异性的空间特征是有意义的，这些空间特征在其它应用中被成功应用[19-20]。同时，数字高程模型(digital elevation model,DEM)被广泛应用于滑坡提取中。从DEM 数据中提取的地形特征，如高程、坡度、坡向和曲率，已被用于描述滑坡的地形特征。值得一提的是，大多数现有研究中使用的DEM数据（1 m到10 m高分辨率），通常来自于3 个数据源：机载Lidar 数据[12,21]、立体像对图像和等高线[2,17,21]。这些高分辨率DEM 数据在滑坡区并不总是可获得的。相比之下，可免费获得的30 m空间分辨率的全球DEM 数据有多种，如高级星载热发射和反射辐射计(ASTER)全球DEM(Global DEM,GDEM)数据和航天飞机雷达地形测绘任务(SRTM)数据。目前还不清楚这些中分辨率DEM数据是否能为滑坡提取提供有效的判别信息。

研究发现，变化检测和分类两种主要技术在滑坡提取中得到广泛应用。一些研究中，通过对同一地理区域在不同日期获取的DEM 或图像计算差值来提取滑坡[10-11,22]。现有研究大多利用多类分类方法提取滑坡，如最大似然法、支持向量机和随机森林等[6,15,23-25]。这些分类方法虽取得很好效果，但利用多类分类方法提取滑坡的一个问题是在分类训练中需利用非滑坡类(非目标类)的样本。实际滑坡提取中，由于非滑坡类别的复杂性，很难从所有这些非滑坡类别中采集足够和完整的样本。此外，在分类中使用来自不同数据来源的多个特征，即数据具高维度。因此，在这种情况下，只需用目标类(即滑坡)的训练样本、且可有效地处理高维数据的单类分类方法[26-27]，可能更适合提取滑坡。因此，本研究提出一种基于多传感器数据的多特征和能有效处理高维数据的单类分类方法的滑坡提取新方法。

1 方法

本研究提出一种基于对象多特征和单类分类的滑坡提取方法。所使用数据包括灾后获取的高分辨率多光谱(或彩色)图像、灾前高分辨率多光谱图像和中分辨率DEM 数据。方法的建立主要考虑以下方面：①即使发生滑坡事件，研究区地形基本保持不变；②滑坡是在处于移动边缘的斜坡上开始的[3,21]；③滑坡改变了地表的外观，这种变化主要反映在光谱特征、植被覆盖和空间变异性上。本文提取方法中，首先通过图像分割生成图像对象，将其作为后续步骤的基本处理单元；再将光谱、植被指数、空间特征和地形特征等对象的多特征相结合，采用改进的单类随机森林(iOCRF)对这些特征图像进行分类[28]，提取滑坡。此外，还计算和分析了不同特征对滑坡提取中的重要性。方法主要包括4 个步骤，即图像分割、对象特征计算、iOCRF分类和特征重要性分析。

1.1 图像分割

为实现面向对象的滑坡提取，首先对图像进行分割，生成合适的图像对象。本研究中对两时相多光谱图像进行分割，产生一致的图像对象。采用一种多层次图像分割方法[29]，该方法已成功用于不同应用中[19,29-31]。

多层次图像分割方法包括3个步骤[29]：多光谱梯度提取、经典分水岭变换和分水线动态值计算。首先使用多通道形态学方法从多光谱图像中提取梯度图像，然后利用经典分水岭变换得到初始的分割结果。为从分水岭变换中产生多层次分割结果，利用分水线动态值评估每条分水线的显著程度。通过计算峰值(即分水线中的鞍点)与周围极小值间的灰度差（即二者的对比度）建立分水线动态值。通过对分水线动态值设置不同阈值，得到多层次分割结果。该方法关键参数是分割层次数(即分水线动态值)[29]。本研究通过试错法确定了最佳分割层次数。

1.2 对象特征计算

为全面反映滑坡引起的地表变化和滑坡发生地的地形特征，将光谱特征、植被指数、空间特征和地形特征综合应用于滑坡提取。从两时相高分辨率多光谱图像中提取光谱、植被指数和空间特征，从分辨率DEM数据中提取地形特征。

1.2.1 光谱特征

本研究中，对每个光谱波段的每个对象内的像元值计算平均值，作为对象光谱特征。

1.2.2 植被指数特征

在植被覆盖密集地区，滑坡通常导致植被覆盖的显著变化，如由植被覆盖区变为非植被覆盖区，因此，采用植被指数表达滑坡引起的植被覆盖变化。有许多不同的植被指数，每一种都有其特定用途的适用性、优点和局限性。利用红光和近红外波段组合计算得到的归一化植被指数(NDVI)是最常用的指数之一。但对于只有可见光波段的彩色图像，可使用可见光波段的植被指数，如绿叶指数(green leaf index,GLI)、归一化绿红差异指数(normalized greenred difference index,NGRDI)和植被颜色指数(color index of vegetation,CIVE)[32-33]。GLI可从土壤和非生物物质中区别出植物[32]。据所用数据是否有近红外波段，选择NDVI 或GLI，计算每个对象内像元的植被指数值平均值，作为对象的植被指数值。

1.2.3 空间特征

本研究中，所采用空间特征包括多变量纹理、对象同质性指数和角点特征[19，20，34]。

多变量纹理是利用多变量变异函数从多光谱图像的多个波段计算得到的[34-35]。多变量纹理计算方法可参考Li等人的文章[34]。多变量纹理有两个关键参数：窗口大小和滞后距离(包括大小和方向)。大多数情况下，使用滞后距离为一个像元计算纹理值。为减少纹理的边缘效应，采用4 个方向（N-S、E-W、NE-SW、NW-SE）最小纹理值(MTmin)为最终纹理值。本研究还使用两个特定方向性纹理，即沿坡向的纹理(MTaspect)和垂直于坡向的纹理(MT⊥aspect)。计算MTaspect或MT⊥aspect时，每个像元纹理方向从DEM 数据提取的坡向确定。首先利用n×n移动窗口从DEM数据中计算出每个像元的坡向值，并将其重采样到与纹理图像一致的像元大小。然后将0°到360°范围内的坡向值等分成4 个方向(即N-S、E-W、NE-SW、NW-SE)，得到最终坡向图。最后计算对象内平均的纹理值MTaspect和MT⊥aspect。

对象同质性指数是一种图像对象的同质性指标[19]，表达滑坡引起的空间同质性差异。计算概括如下：①计算像元的空间同质性指数[36]。对原始多光谱图像计算4个方向(N-S、E-W、NE-SW、NW-SE)的高斯滤波值，得到原始图像与滤波后图像之间的均方根值，取最小值为每个像元的空间同质性指数值；②计算每个对象内像元空间同质性指数值的平均值为对象同质性指数[19]。对象同质性指数值越低，对象的空间同质性越高。

角点特征用来表征滑坡引起的空间变异性变化。角点特征代表所有方向上强度显著变化的图像特征[37]，已被用于提取人造建筑物或建筑物损毁[18，20]。本研究中，加速分段测试特征(Features from Accelerated Segment Test,FAST)角点检测器被用于角点检测[39]。由于角点检测的输出是一组不能直接量化为图像对象的点，因此，本研究计算每个对象内的角点密度[20]，即对象中被检测为角点的像元的比例，并作为对象的角点特征。

1.2.4 地形特征

本研究采用的地形特征包括高程和坡度。由于采用的地形特征和图像的像元大小不同，将高程和坡度图像重采样到与高分辨率图像相同的像元大小，然后，计算相应的对象高程和坡度值平均值。

1.3 iOCRF分类

本研究采用一种最近提出的iOCRF分类方法[28]，它具有精度高、对过拟合训练数据不敏感、能有效地处理高维数据和能够度量特征的重要性等特性。iOCRF 是单类随机森林的改进版本[26]。单类随机森林是一种基于随机森林的集成学习方法[26]，它利用人工离群值(即非目标类的样本)生成过程，将单类分类任务转换为两类分类问题。因此，单类随机森林的关键是准确地生成或估计非目标类的样本。iOCRF中[28]，采用一种从未标记样本生成非目标类样本的新方法。iOCRF 分类概括：①研究区图像中手动选取目标类样本，从不包含选取的目标样本的剩余图像中随机选择未标记的样本，将其视为非目标样本，这些样本既包含非目标类样本，也包含目标类样本；②将这些随机生成的未标记样本和人工选取的目标样本用于随机森林分类中，得到初始分类结果；③得到初始分类结果后，选择分类后验概率高于指定阈值的非目标类样本作为非目标样本，即更新后的非目标样本，将更新后的非目标样本和目标样本用于随机森林，得到最终单类分类结果。对iOCRF 分类，应确定与随机森林训练相关的两个参数，即用于随机特征选择的决策树的数量和特征的数量。决策树的数量通常设置为200，这个值被认为可确保算法在统计学意义上收敛[26]。特征数量为所用数据的特征维度的算术平方根加1[26]。

1.4 特征重要性分析

为分析不同特征对提取滑坡的贡献，使用基于随机森林的特征重要性度量方法计算所使用特征的重要性。特征重要性计算如下[6]，对随机森林中每一棵树，使用训练集中留下的对应的袋外(out-of-bag,OOB)数据计算OOB 误差。通过重排OOB 数据中每个特征的值，计算扰动前后OOB误差的差值，估计每个特征对每棵树的重要性。差值越大，特征越重要。通过平均每个特征在所有树上的重要性排名,得到每个特征最终的重要性。

1.5 结果评价

本研究采用总体精度、生产者精度、用户精度等精度指标评估滑坡提取结果[40]。此外，还采用F1分数[41]，F1分数已被广泛用于单类分类的精度评估[42-43]。为全面地评价方法的性能，从特征集和分类方法两个方面进行对比实验。iOCRF 使用了4 个不同特征集评估特征在滑坡提取中的性能，即光谱特征分别加上其它类型特征(植被指数、空间特征、地形特征)和所有特征集(光谱、植被指数、空间和地形特征)。此外，还与经典的多类分类方法——随机森林方法进行了比较。

2 研究区和数据

利用研究区两时相高分辨率多光谱图像和中分辨率DEM 数据对本研究提出的方法在滑坡提取应用中进行评价。

研究区位于甘肃省舟曲地区，该地区在2010 年8月8日发生了严重的泥石流灾害，淹没了舟曲县中心和周边地区。这是七十年来破坏最严重、死亡人数最多、救灾难度最大的一次泥石流灾害。舟曲北部沟谷坡度大、高差大、植被覆盖度低，为泥石流灾害的发生提供了有利条件。在强烈的季风雨之后，舟曲地区发生了两场大规模的泥石流，流出物从东北部向西南方向滑下山谷，摧毁了大片农田、建筑物和道路。

本研究使用的数据由灾前航空彩色图像、灾后QuickBird 图像和ASTER GDEM 数据组成。航空彩色图像于2008 年7 月采集，有3 个可见光波段(蓝、绿、红)，图像空间分辨率约1 m。QuickBird 图像获取于2010 年8 月15 日，是滑坡灾害发生8 天后拍摄的。QuickBird 图像包括2.4 m 分辨率的4 个多光谱波段(蓝、绿、红和近红外)和0.6 m 分辨率的全色波段。利用Gram-Schmidt 全色锐化方法(Laben and Brower，2000)对QuickBird图像的多光谱和全色图像进行融合，得到像元大小为0.6 m的融合后4波段多光谱图像。采用融合后的QuickBird 多光谱图像(0.6 m)作为参考图像，对航空彩色图像（1 m 分辨率）配准，并重采样到0.6 m的像元大小。最终使用的两类图像（不同时相）的大小为3 000×3 000像元(图1-(a)，(b))。所采用的ASTER GDEM 数据（V2版）是2011年发布的。GDEM数据覆盖了北纬83°至南纬83°间的地球表面，空间分辨率为1角秒(约30 m)[44]。将研究区GDEM 数据从WGS-84 投影到UTM 投影（图1-(c)）。研究区滑坡通过目视解译得到，为参考数据(图1-(d))。

图1 舟曲地区灾前航空彩色图像(a)、灾后融合的QuickBird多光谱图像(R、G、B分别为波段3、2、1)(b)、ASTER GDEM数据(c)和参考图(d)Fig.1 Pre-disaster aerial color image of Zhouqu area(a),post-disaster fusion QuickBird multispectral image(R,G,B are bands 3,2,1)(b),ASTER GDEM data(c)and reference map(d)

3 结果和讨论

3.1 对象特征的生成和分析

首先利用多层次图像分割方法对舟曲地区灾前航空彩色图像和灾后融合的QuickBird 多光谱图像进行分割得到分割结果[29]，使用分割层级为35。然后分别从灾前航空彩色图像和灾后多光谱图像提取光谱、植被指数和空间特征，从ASTER GDEM 数据提取地形特征。选择灾前航空彩色图像的植被指数为GLI，灾后多光谱图像的植被指数为NDVI。使用3×3像元的窗口大小进行多变量纹理计算。对象同质性指数采用窗口大小为5×5、标准差为1.5、0.25（两个方向）的方向性高斯滤波器生成。使用FAST 角点检测器提取角点特征[39]。使用3×3 像元的窗口大小计算坡度和坡向图像后，将高程(即ASTER GDEM数据)、派生的坡度和坡向图像重采样到0.6 m像元大小，最后生成21 个对象特征，包括每个光谱波段的对象平均值、两时相植被指数、MTmin，MTaspect，MT⊥aspect，对象同质性指数和角点密度、高程和坡度。

图2 展示了一个两时相高分辨率图像局部和相应的对象同质性指数图像。通过比较灾后图像和灾前图像(图2-(a)，(b))，发现滑坡区空间同质性增加，即建筑物区变为滑坡，滑坡区对象同质性指数值显著降低(图2-(d)与图2-(c)对比)。因此，对象同质性指数为区分滑坡和非滑坡地区提供了有用的空间特征。另一个两时相高分辨率图像局部和相应的角点密度图像见图3。图中显示滑坡区的空间变异性显著增加(图3-(b)与图3-(a)对比)。相应地，角点密度值也增加(图3-(d)与图3-(c)对比)。因此，所采用的角点密度有效地量化了滑坡引起的空间变异性的变化。图4 展示了整个研究区滑坡和非滑坡地区的高程和坡度的相对频率分布。从图4-(a)可看出，滑坡区高程主要集中在1 320～1 480 m，非滑坡区高程分布是平缓的，高程值的分布非常宽。从图4-(b)可看出，滑坡区坡度分布呈尖峰形，峰值在8 度左右，非滑坡区坡度分布很宽。因此，由ASTER GDEM数据提取的高程和坡度为区分滑坡区和非滑坡区提供了有用信息。

图2 灾前(a)和灾后(b)高分辨率图像局部与相应的灾前(c)和灾后(d)对象同质性指数图像Fig.2 Pre-disaster(a)and post-disaster(b)high-resolution image local and corresponding pre-disaster(c)and post-disaster(d)object homogeneity index image

图3 灾前(a)和灾后(b)高分辨率图像局部与相应的灾前(c)和灾后(d)角点密度(DCP)图像Fig.3 Pre-disaster(a)and post-disaster(b)high-resolution image local and corresponding pre-disaster(c)and post-disaster(d)corner density(DCP)images

图4 舟曲地区滑坡区与非滑坡区高程(a)和坡度(b)的相对频率分布Fig.4 Relative frequency distribution of elevation(a)and slope(b)between landslide area and non-landslide area in Zhouqu area

3.2 滑坡提取结果

通过使用iOCRF 分类对所有21个对象特征图像进行分类，所用参数设置如下：决策树的数量为200个，随机特征选择的属性数量为5个，随机选取的未标注样本的数量是滑坡训练样本(目标类)数量的两倍。选择滑坡类训练样本为82 841 个像元，参考数据的所有剩余像元都被用作检验样本。此外，还利用4种滑坡提取方法（即利用不同的特征和随机森林分类方法）进行滑坡提取，作为对比。对随机森林分类，选取82 841个滑坡像元点和97 176个非滑坡像元点为训练样本。

表1 展示了使用不同方法进行滑坡提取的精度。首先分析使用不同特征集的基于iOCRF 的滑坡地图的精度(表1-a-d)。从表中可看出，使用光谱特征和附加特征的滑坡提取的OA 和F1 分数相似且相对较高(F1 分数在84.79%～86.42%) (方法a-c)，利用光谱特征加空间特征和光谱特征加地形特征滑坡提取精度略高于光谱特征加植被指数的精度。将植被指数、空间和地形特征与光谱特征相结合的方法（本文方法）(方法d)的F1分数比3种使用光谱特征和附加特征的方法的F1分数高2.32%～3.95%。本文方法的生产者精度和用户精度不仅更平衡，且比3 种对比方法的生产者精度和用户精度更高或更具可比性。因此，综合利用光谱、植被指数、空间和地形特征进行滑坡提取，显著提高了提取的精度。

表1 不同方法进行滑坡提取的精度Table 1 Accuracy of landslide extraction by different methods 单位：%

从表1 中发现，所提出的使用iOCRF 分类的方法(本文方法，方法d)比使用随机森林分类的对比方法(方法e)得到的精度高得多。特别是该方法的F1分数远高于基于随机森林的滑坡提取方法，提高了9.49%。此外，与基于随机森林的方法相比，该方法的生产者精度和用户精度分别提高了6.58%和11.95%。

图5展示了使用不同方法的滑坡提取结果。从图中可看出，当使用不同特征集时，大多数滑坡被正确提取出来，但仍存在不同程度低估(即滑坡被误判为非滑坡)和高估(即非滑坡被误判为滑坡)。当使用光谱、植被指数、空间和地形特征时(本文方法，图5-(d))，低估和高估最少。通过与研究区滑坡灾害前后图像的对比(图1)，发现许多发生在建筑区的滑坡在基于随机森林的提取结果(图5-(e))中被明显遗漏。在基于随机森林的提取结果(图5-(e))中，许多非滑坡被误认为滑坡，如裸土。与使用随机森林分类(图5-(e))相比，使用iOCRF分类(本文方法，图5-(d))的提取结果的低估和高估显著减少。

图5 使用不同方法的研究区滑坡提取结果Fig.5 Landslide extraction results in the study area using different methods

为进一步验证所使用特征集的有效性，图6 展示了使用不同特征集的滑坡提取结果的两个局部。从图6中发现，当使用光谱和植被指数特征时，一些建筑物和裸土被错误地识别为滑坡(图6-(b))；当使用光谱和空间特征时，高估主要分布在裸土区(图6-(c))；当使用光谱和地形特征时，高估主要分布在完好的建筑区(图6-(d))。当植被指数、空间和地形特征与光谱特征相结合用于滑坡提取时(本文方法，图6-(e))，在建筑区和裸土区，高估都很少。

图6 两个灾后图像的局部(a)、基于iOCRF使用不同的特征集的滑坡提取结果：光谱和植被指数特征(b)、光谱和空间特征(c)、光谱和地形特征(d)、光谱、植被指数、空间和地形特征(e)Fig.6 Local(a)of two post-disaster images,landslide extraction results using different feature sets based on iOCRF:spectral and vegetation index features(b),spectral and spatial features(c),spectral and topographic features(d),spectral,vegetation index,spatial and topographic features(e)

3.3 特征重要性分析

本研究使用的21个对象特征的重要性如图7所示。从图7可看出，光谱特征是最重要的特征，特征重要性为43.06%；植被指数、空间特征和地形特征的重要性分别为14.81%、29.78%和12.35%，表明空间特征比植被指数特征和地形特征对滑坡提取的贡献更大。从图8 中还看出，灾后图像特征比灾前图像具更高的特征重要性，可能是因为灾后图像特征与滑坡目标直接相关。据特征重要性排序，前3 位特征分别是灾后红波段平均值(12.99%)、灾后植被指数平均值(11.35%)和灾后蓝波段平均值(9.38%)。地形特征(高程和坡度)具较高的特征重要性(分别为6.40%和5.95%)，表明它们可有效减少滑坡和非滑坡区域的混淆。虽各空间特征的重要性相对较低(2.61%～3.57%)，但它们是滑坡提取中不可缺少的特征，为提取滑坡提供了有用信息。

图7 本文方法中使用的21个对象特征的重要性Fig.7 The importance of the 21 object features used in this method

3.4 讨论

（1）不同来源的数据得到的多种特征(光谱、植被指数、空间和地形特征)被综合应用于滑坡提取。每种类型特征在滑坡提取中起到不同作用。由于滑坡和其它土地覆盖类型间的光谱相似性，单独使用光谱特征精度有限。植被指数可有效突出滑坡引起的植被覆盖变化。本研究首次将多变量纹理、对象同质性指数和角点密度等新的空间特征应用于滑坡提取，有效地刻画了滑坡引起的空间变异性的变化。考虑到高分辨率DEM 数据在滑坡区并不总是可获得的，本研究使用从中分辨率DEM 数据(ASTER GDEM)提取的地形特征，为滑坡提取提供了有用的地形信息。其它可免费获得的30 m分辨率的全球DEM 数据(如SRTM 和ALOS World 3D数据)也可用于滑坡提取。

（2）采用iOCRF 分类对滑坡提取。基于iOCRF分类方法只需滑坡类样本，这种简捷样本选择提高了滑坡提取的效率。与传统同时需滑坡和非滑坡类样本的多类随机森林分类方法相比，使用iOCRF 分类方法得到的滑坡提取结果更准确。此外，本文采用iOCRF 还提供了一种特征的重要性度量，用于量化和分析不同特征对滑坡提取的贡献。

4 结论

本文提出一种利用两时相高分辨率多光谱图像和中分辨率DEM 数据的面向对象的滑坡提取方法。该方法综合利用光谱、植被指数、空间和地形特征等多种对象特征来表达滑坡引起的地表变化和滑坡发生地的地形特征。通过使用iOCRF分类对所有特征图像进行分类来提取滑坡。结果表明，该方法在研究区优于对比的方法，为提取滑坡提供了一种有效方法。

致谢：本研究得到国家自然科学基金(41371329)资助。感谢北京师范大学的唐宏教授提供了本研究中使用的Quick-Bird图像。