高分一号汶川极震区滑坡提取研究

黄 汀,白仙富,庄齐枫,徐敬海

(1. 南京工业大学测绘科学与技术学院，江苏 南京 210009； 2. 云南省地震局，云南 昆明 650000)

遥感技术凭借覆盖范围广、信息量大、可连续观测的特点成为地震灾情获取的首选技术。近年来，高分辨率数据越来越多，传统基于像元的提取方法无法充分利用高分辨率影像的信息，而面向对象提取方法不仅考虑影像的光谱信息，而且能结合纹理、空间信息[1]，提高分类精度，逐渐成为了遥感信息提取的主要方法，在灾后信息获取方面也得到了广泛应用。如李强[2]利用影像提取地震灾后的道路受损情况；安立强[3]利用面向对象的分析方法完成了对地震后产生的次生灾害的提取工作；OTHMAN等[4]完成了对伊拉克地区的滑坡提取；MICHAEL等[5]通过光学影像与DEM判别滑坡。但是由于滑坡成因复杂、形态多样的特点，滑坡信息的提取仍存在一定的困难，如分割参数的选取没有统一的标准、分类时影像信息利用不充分、使用国外卫星数据源时间不匹配等。高分一号(GF-1)卫星的投入使用，为滑坡提取提供了新的前景。高分一号卫星具有分辨率高、幅宽较宽的优势，解决了空间分辨率与时间分辨率统一的难题，对地物信息展示更清楚，精度高，更稳定。但目前利用高分一号的滑坡提取研究尚处于探索阶段，文献[6—7]利用高分影像提取地震后滑坡信息。但是上述研究都是对于中小地震后的滑坡提取，对于特大地震如汶川地震的提取效果尚未可知，其所采用的最优尺度、特征参数、分类规则等均无法适用于汶川极震区。本文拟利用高分一号影像并结合eCognition软件，分析高分影像与震后滑坡特点优化参数选取模式，建立分类规则，实现对汶川地震的极震区域的滑坡提取。

1 数据源及预处理

本文采用高分一号高分辨率影像提取滑坡信息，影像时间为2013年9月17日。选择“5.12”汶川地震极震区为研究区(如图1所示)，地理坐标为102°33′E—106°17′E、30°40′N—33°7′N，在四川42个受灾县(市)共有滑坡3000余处，汶川地震后，经过专业人员调查，发生地质灾害约15 000处，其中滑坡1700处[8-9]。

图1 研究区的GF-1遥感影像

将原始影像进行辐射定标与大气校正后，利用影像红外与近红外波段计算NDVI，为后续提取作准备。

2 滑坡提取方法

2.1 多尺度分割

影像分割决定着分类的好坏，其原理是将影像分割为多个对象，对象之间具有较大的差异，对象内部具有相似的灰度、纹理等特征，其最终目的是使得每一个分割出来的对象都与实际地物符合。本文基于eCoginition Developer 8.6软件平台，采用多尺度分割方法对影像进行分割，以任意像元为初始对象，再计算其与相邻像元的异质性值，若异质性在规定范围内则将两者合并，成为新的对象，反之停止合并。多尺度分割结果受分割参数影响，参数包括光谱权重、均质因子、分割尺度，目前对于参数的选取都是基于经验结合目视判断并没有特定标准，主观性大，分割效果不理想。下文将说明各参数的选取过程。

2.1.1 光谱权重

虽然多尺度分割考虑了形状纹理等因素，但是影像的光谱信息仍在分割中作为主导。由于每个波段的信息量不同，在分割时应当根据信息的多少赋予每个波段不同的权重，权重越高，代表该波段的信息量越大，分割过程中该波段信息会更多被参考；如果某波段权重低，代表该波段包含目标信息少，对于提取没有太多贡献。高分一号影像包含4个波段，信息势必存在冗余，本文应用相关系数和协方差建立矩阵来衡量相关性与信息量，并以此作为依据设置权重。

协方差矩阵、相关矩阵能反映多个变量之间的相关程度。矩阵中的每一个元素表示对应两个波段之间的相关系数，两个波段之间的相关系数越大说明这两个波段相关性越大，所包含信息存在冗余。在确定光谱权重时应提高波段间相关性小的权重，提高分割影像时的准确性。表1、表2分别为对试验数据计算的各波段统计协方差矩阵和各波段相关性统计矩阵。

表1 各波段统计协方差矩阵

表2 各波段相关性统计矩阵

在表1中可以看出，方差从大到小排列依次是波段4、波段2、波段1、波段3。在表2中可以看出波段1与波段2、波段3这3个可见光波段相关性很大，可认为它们对分割的贡献程度相似，而波段4信息量大，且与波段1、波段2、波段3相关性小。考虑到本次研究提取目标为滑坡，植被的变化是判断滑坡的重要依据，判断对象是否为植被是提取滑坡的重要环节。可见光3个波段相关性大，且植被在可见光波段并没有较大的响应，应当赋予波段1、波段2、波段3较低的权重；同时，植被在近红外波段有强反射，而且波段4具有的信息量最大，与其他波段的相关性小，因此应赋予波段4较高的权重。根据研究,权重的取值一般在0～2之间[10-12]，故将波段1、波段2、波段3的权重设为1，波段4权重设为2。

2.1.2 均质因子

对象的异质性h表示对象内部像元与其相邻像元的相似程度，异质性越低说明对象越纯净，即像元与像元之间的差异小，可以认为这些像元都代表同一类地物，并且具有相似特征(如方差、亮度等)。反之说明对象内部差异较大，含有多种地物，影响分类精度，不利于提取。异质性由均质因子控制，其包括颜色和形状两个标准因子，表达式如下

h=wcolor×hcolor+wshape×hshape

(1)

式中,h为对象的总体异质性；hcolor为光谱异质性；hshape为形状异质性；wcolor为光谱异质性的权重；wshape为形状异质性的权重。且wcolor与wshape之和为1。

光谱异质性hcolor取决于两个因素：一是组成对象的像元数目n，二是各波段标准差σ，如下

(2)

式中，wn表示参与分割及合并波段的权重；σn为区域的标准差；n为波段序号；m为波段数。

形状异质性公式如下

hshape=wcomhcom+wsmoothhsmooth

(3)

式中，wcom指用户自定义的紧凑度权重；wsmooth指用户自定义的光滑度权重。且wcom+wsmooth=1。

参数设置时要保证光谱信息的充分利用，光谱因子不能过低，过高的形状因子会使对象内异质性增大，进而阻碍信息的正确提取。因此多尺度分割需要遵循两个原则：一是尽可能地保障光谱信息，确保光谱因子在分割中起主导作用；二是提高形状因子的权重来调节分割对象形状，使对象轮廓尽可能地与实际地物符合[13]。根据多次试验，将光谱因子设为0.7，即形状因子为0.3，形状因子中将紧致度因子设为0.7，则光滑度因子设为0.3。

2.1.3 最优分割尺度

分割后影像中包含纯对象与混合对象，就分类而言，纯对象的数量越多分类越准确，进而影响地物信息提取的精度。分割尺度决定了影像中纯对象与混合对象的比例，因此最优分割尺度选取尤为重要。不同地物的最优分割尺度各不相同，因此最优分割尺度只能是一个参考值，也没有严格的标准。尺度的选择可以参照一定的规则使分割尺度尽量接近最优尺度：当某一地物在最优尺度分割时，分割后的对象的轮廓应与实际地物符合，地物可由一个或多个对象组成，即不能过分割使对象太多，也不能欠分割使一个对象含有混合对象。并且不同类对象之间差异大，容易与其他物类别区分开。本文利用局部方差法选择最优分割尺度，并借助ESP(estimation of scale parameters)工具选择最佳分割尺度[14]。

局部方差法将每个分割对象的平均局部方差作为评价参数，假设对象与背景存在着较大的光谱差异，随着分割尺度的增大，分割出的对象面积也随之增大，其内部的方差也会增大，当对象与实际地物范围符合时，分割对象会停止变化，方差也不变，分割出的对象边界会在较大的分割尺度下保留，直到这些对象被合并为更高的类别(如树木合并为森林)，代表这种地物与实际范围最符合，此时的分割尺度为该地物的最优分割尺度。

如果某种地物被最优分割出时，那么这类地物方差的增大与停滞就会对整个图像的平均方差产生影响。本文希望得到的是一条随着尺度增大平均方差上升的曲线，当曲线的上升趋势减小或趋于平坦时，此时的分割尺度为最优分割尺度，但是曲线的变化对应的尺度在图像上不易得到，因此ESP工具引入了ROC(rate of change)

(4)

式中，L为当前分割尺度层所对应的方差；L1为上一个分割尺度层对应的方差。

从图2中看出，ROC曲线整体呈下降趋势，定义在曲线持续下降或突然衰减后的第一个峰值作为最优分割尺度。图像在95、115、125出现了相对的峰值，考虑到滑坡为本次的提取目标，较大的分割尺度会出现欠分割的情况，不利于后续的分类，故选择95为最优分割尺度。

图2 ROC曲线

2.2 提取规则建立

影像分割将图像分割成许多对象，它们都具有相似的特征，这些特征是决定对象归属的关键。但这些对象仍有一定比例的混合像元，可以采用模糊分类方法对分割对象进行分类。它将“是”与“非”的判断转换成0到1之间连续变化的值来描述目标对象属于多种地物的概率，这样给出的分类信息更加全面，使得对混合地物的分类结果更具有客观性[15]。

提取地物特征，并建立一定的分类规则是面向对象分类的关键，规则的好坏对于分类的精度有着巨大影响。对象特征主要包括光谱特征、几何特征和纹理特征。

光谱特征是判别对象的主要特征，研究区滑坡的土层整体偏暗，而非滑坡土层具有较高的亮度，滑坡体与周围也有着显著的差异。

植被的变化也是滑坡重要的解译标志之一，大多数滑坡都主要发生在斜坡上，通常这些斜坡都有植被覆盖，滑坡体受到重力的作用沿着斜坡向下滑动时，导致植被受损、折断、掩埋。

本文利用植被覆盖度来描述植被的变化，估算公式为

Fcover=(NDVI-NDVImin)/(NDVImax+NDVImin)

(5)

式中，NDVImin为影像中NDVI最小值0.05；NDVImax为NDVI最大值0.8。

滑坡一旦发生就会形成一些特殊的地貌形状，对于大部分滑坡，根据它们独特的滑坡地貌是较容易辨识的。根据文献[16]，汶川地震造成的滑坡多呈长条状、舌型。因此，可以将长宽比作为几何特征识别滑坡。

将规则整理后，建立了如图3所示的提取流程。首先判断植被覆盖度，将覆盖度较小的非植被区域剔除，再以长度及负植被覆盖度作为依据剔除水体，保留滑坡疑似区域，最后通过坡度与长宽比筛选滑坡。

图3 滑坡提取流程

3 结果与分析

3.1 提取结果

依据上述所建规则，分别对研究区的遥感影像进行地震型滑坡信息提取，提取结果如图4所示。

图4 滑坡提取结果

从图4可看出，整体提取效果较好，汶川地震型滑坡有沿发震断裂呈条带状分布特点。但是研究区中有部分小滑坡未被提取出来：一是因为滑坡不明显区域光谱特征与植被相似，产生误分割，导致分类时特征不明显，无法正确分类；二是因为高分辨率影像的特征丰富，同类地物之间的差异性增加了区分难度。

3.2 精度评价

本文采用混淆矩阵法对分类结果进行评价，首先选取个数为N的样本，将每个样本的地表真实像元的分类与分类图像中的分类相比较计算建立矩阵即

(6)

矩阵中行向量代表一类地物的实际分类情况，列向量代表等于地面真实地物在分类图像中对应于相应类别的数量，对角线上的数值越大说明分类越准确。主要评价指标有总分类精度、各类别生产精度、用户精度及Kappa系数等。

总体分类精度为正确分类的像素个数占总数的比例，是对分类结果的总体评估，计算公式为

(7)

用户精度是指正确分到某类的像元总数与分类器将整个影像的像元分为像元总数的比值，计算公式为

(8)

生产精度是指将整个影像的像元正确分为某类的像元数与某类真实参考总数的比值，计算公式为

(9)

总体Kappa系数计算公式如下

(10)

本研究滑坡信息提取结果的混淆矩阵见表3。从表中可以看出植被的分类精度高，这是由于植被覆盖度的加入将植被与非植被很好地区分，滑坡与其他地物(包括薄云和裸地)存在混分的现象。影像分割后对象不可能与实际地物完全符合，过分割、欠分割不可避免，混合像元的分类不可能完全正确，从而影响分类精度。总体而言，滑坡生产精度为87%，用户精度为81%，总体Kappa精度为82.2%，说明面向对象提取滑坡具有较高精确性，表明基于多尺度分割技术的滑坡灾害提取能够满足震后评估的要求，该方法快速有效地提取了高分辨率遥感影像上地物信息。

表3 面向对象分类的混淆矩阵

同时，利用研究区域野外滑坡点的资料[17]，经统计得出研究区滑坡56处，如图5所示。并与提取结果进行对比，发现有47处被识别并提取，9处未提取，精度达到83.9%，进一步说明该提取方法具有可行性。

图5 实测滑坡点对比

4 结 语

本文尝试利用GF-1高分辨率卫星，选取汶川地震极震区域影像，提出了一种快速提取滑坡信息的方法。该方法首先确立适合极震区的分割参数，对影像进行分割，分析滑坡特征，建立面向对象遥感分类规则，再通过模糊分类方法提取滑坡灾害信息。分类完成后选取随机样本与实际滑坡点进行对比，证明分类精度较高，验证了方法有效性。试验结果证明了遥感方法尤其是利用高分辨率影像在地质灾害提取方面的优势，能够为震灾调查和快速预评估提供辅助数据，为今后地质灾害的防治和快速评估提供参考。在完成影像对象分类后的处理中，虽然完成了滑坡的提取，但是对于滑坡特征与空间信息的利用并不充分，造成提取目标与相似地物的混淆。因此，在后续的研究中加入道路、河网等矢量数据，可进一步提高提取精度。

[1] BAATZ M, SCHPE A. An Optimization Approach for High Quality Multi-scale Image Segmentation[C]∥Proceedings of Beiträge Zum AGIT-Symposium. [S.l.]: [s.n.], 2000:12-23.

[2] 李强，张景发. GIS与面向对象结合的遥感影像震后损毁道路快速提取[J]. 测绘通报,2015(4):78-81.

[3] 安立强，张景发，赵福军.汶川地震次生灾害提取——面向对象影像分类技术的应用[J]. 自然灾害学报,2011,20(2):160-168.

[4] OTHMAN A A, GLOAGUEN R. Automatic Extraction and Size Distribution of Landslides in Kurdistan Region, NE Iraq[J]. Remote Sensing, 2013, 5(5):2389-2410.

[5] RAU J Y, JHAN J P, RAU R J. Semiautomatic Object-oriented Landslide Recognition Scheme from Multisensor Optical Imagery and DEM[J]. IEEE Transactions on Geoscience & Remote Sensing, 2013, 52(2):1336-1349.

[6] 朱娇.基于高分一号影像的宝兴县滑坡信息提取研究[D].成都：成都理工大学,2015.

[7] 李启源.高分一号遥感影像地质灾害信息提取方法研究[J]. 测绘与空间地理信息,2016,39(2):17-20.

[8] 殷跃平.汶川八级地震地质灾害研究[J].工程地质学报,2008,66(4):433-444.

[9] 邓书斌. ENVI遥感图像处理方法[M].北京：科学出版社，2008．

[10] HÖLBLING D, FRIEDL B. An Object-based Approach for Semi-automated Landslide Change Detection and Attribution of Changes to Landslide Classes in Northern Taiwan[J]. Earth Science Informatics, 2015, 8(2):327-335.

[11] 刘建平，赵英时，孙淑玲.高光谱遥感数据最佳波段选择方法试验研究[J]. 遥感技术与应用,2001,16(1):7-13.

[12] 冯文卿. 利用多尺度融合进行面向对象的遥感影像变化检测[J]. 测绘学报, 2015, 44(10):1142-1151.

[13] 王赛. 基于多源遥感数据的汶川地震型滑坡信息提取研究[D].北京：中国地质大学(北京),2015.

[15] 韩凝. 空间信息在面向对象分类方法中的应用[D].杭州：浙江大学,2011.

[16] 范建容，李秀珍，张怀珍，等.汶川地震灾区崩塌滑坡体几何特征信息遥感定量提取与分析[J].水土保持通报,2012,32(2):118-121.

[17] 秦绪文，杨金中，张志，等.汶川地震灾区航天遥感应急调查[M]，北京：科学出版社,2008.