文广超，张哲玮，肖学军，谢洪波，王文博，杨运航，陈红旗

(1.河南理工大学资源环境学院, 河南 焦作 454000；2.河北省煤田地质局第四地质队，河北 张家口 075100；3.南京师范大学地理科学学院，江苏 南京 210023；4.中国地质环境监测院，北京 100081)

0 引言

滑坡在自然灾害造成的经济损失中位列第二[1]。我国西南地区位于中国地形第一阶梯向第二阶梯过渡的狭长地带，地形起伏大，地质构造复杂，是滑坡灾害发生的极高风险区。我国约80%的大型滑坡都发生在此范围内，给当地人民群众的生命财产造成了巨大损失[2]，如1989年7月四川华蓥山溪口滑坡灾害，2004年7月四川宣汉天台乡滑坡灾害，2017年6月四川茂县新磨村滑坡灾害等。滑坡灾后区域的道路场地狭小、抢险情势复杂、人员搜救艰难，救援善后工作需要动用巨大的社会力量[3]。因此，在灾害发生后快速获得灾区的滑坡规模等灾情信息就显得尤为重要。

目前国内外获取滑体信息主要靠野外调查，无法满足灾后快速查明灾情，指导应急救援力量进场的需要。遥感技术具有获取速度快、覆盖面积大、周期性强等优点[4]，可以很好地弥补野外调查的不足。为了能及时获取灾区信息，为救援工作提供准确资料，采用遥感技术这一先进手段是必要的，这也是空-天-地-现场一体化发展的必然趋势。基于遥感技术提取滑坡信息的方法主要有三大类：目视解译法、面向对像分类法和基于像元分类法。通过影像目视解译识别滑体可充分利用专家经验，具有较高的准确度，是一种重要的技术方法。缺点在于受主观因素影响强，效率低下，无法实现信息的快速提取。结合数字图像处理技术的滑体信息自动、半自动提取，在工作效率和客观性上具有很强的优势，是近年来的研究热点。其中利用多源数据，融合多维特征信息的面向对像分类法被诸多学者研究和应用[5-7]。该方法的实质是基于人的思维模式把对像作为信息的内容在地理空间进行分类[8]。面向对像分类方法对遥感数据分辨率的要求较高，需辅以多源数据来构建分割规则，在多尺度分割后产生的被忽略的细小对象也会影响对影像进一步的应用[9]。相较于此，基于像元的滑体信息提取方法可较好地保留影像细节信息。诸多学者用这种方法进行了研究，如苏凤环等[10]基于TM、ETM+影像，利用缨帽变换、DN值差值增强技术对大规模泥石流及滑坡进行了提取试验，获得了较好地提取效果。NICHOL J等[11]对灾前、灾后SPOT影像利用最大似然法进行变化检测，识别出了研究区70%的滑体。张海涛[12]利用高分辨率的遥感影像，以像素为单位将ICA原理引入，提取影像的独立成分，基于变化检测算法识别滑体，识别精度达75%。基于像元的识别方法一般未考虑像素间空间位置信息而仅考虑光谱值的变化[13]，对数据地利用不够充分，易受影像分辨率、大气干扰等许多因素的影响，产生“椒盐”噪音，导致滑体识别的精度整体不高。

因此，本文基于Landsat 系列遥感数据和数字高程数据(DEM)对 “6·24”四川茂县、“8·28”贵州纳雍等滑坡区的地物特征进行深入分析，综合滑体与非滑体光谱反射率差异和空间位置信息建立滑体提取方法，为区域滑坡调查研究及灾后应急响应提供技术支持。

1 数据源及预处理

Landsat系列遥感数据和Aster GDEM数据为免费开源的数据，获取便捷，信息丰富，适用于滑体快速提取研究(表1)。其中Landsat系列遥感数据的空间分辨率为30 m，已经做过地形参与的几何校正、辐射定标和大气校正等预处理。Aster GDEM数据根据美国航空航天局的新一代对地观测卫星TERRA的观测结果制作完成，空间分辨率为30 m，经预处理生成坡度信息。以上数据均源自美国地质勘探局(USGS，http://glovis.usgs.gov)。

表1 数据类型

2 滑体特征遥感识别分析

2017 年 6 月 24 日，四川省茂县叠溪镇新磨村后山发生高位山体滑坡。滑源区岩体沿岩层层面滑出，解体后沿坡面高速运动，沿途铲刮坡面原有松散崩滑堆积物，直至运动到河谷底部[14]。整个滑坡过程对原始地层的整体性以及地表盖层造成严重破坏，在遥感数据中表现出与背景地物在光谱特征上的巨大差异。本文利用多期遥感影像与Google Earth高清影像对滑坡及背景地物目视解译并进行采样(图1)。其中，滑体的采样区为滑坡主体，植被的采样区为滑坡周边山坡NDVI值大于0.8的区域，水体的采样区为离新磨村较近的叠溪海子内，阴影的采样区为滑坡周边厚云层下部阴影区。因滑坡区附近的植被覆盖度较好，因此裸地的采样区位于茂县与汶川县交界处，该处裸地结合多期影像解译处于一个较稳定的状态。根据以上采样数据，计算各类地物的光谱均值并生成光谱曲线(图2)。

图1 茂县滑坡区域概况及采样示意图Fig.1 Maoxian landslide area overview and sampling diagram

图2 茂县滑坡区地物光谱曲线图Fig.2 Comprehensive spectrum curve of ground objects in the study areas

对茂县滑坡区域地物光谱地综合分析后得到以下规律：(1)在红光波段，植被与阴影的最大值均在500以下，而滑体的最小值基本都在500以上，最大值在2 000以下。(2)从红光波段到近红外波段，土壤的反射率快速增加，但滑体反射率升高并不明显。因为滑体表层岩土结构松散，具有较大的孔隙度和含水率，导致其在近红外波段的反射率较低。因此滑体从红光波段到近红外波段的“斜率”会小于裸地。据统计裸地的“斜率”均在0.3以上，滑体的“斜率”在0.3以下。(3)在对水体的去除试验中，本文使用NDWI等多种水体指数进行了提取试验，最终发现使用NEW[15]水体指数计算后水体光谱值集中在-1～0.5的范围内，与其它水体指数相比具有最大的区分度。

同时在目视解译中发现滑体与相邻地物在色调，本质即是光谱上有明显的差异。借鉴边界检测中寻找图像灰度变化不连续点的思想，提出使用“跨越”滑体边界，同时包含滑体光谱信息与相邻非滑体光谱信息的光谱剖面线方法(图3-a)，对边界附近的信息进行线状连续采样，以多段线的形式展示每个象元在7个波段上的反射率值。对采样结果进行分析发现，七个波段均在滑体边界处发生了突变(图3-b)。其中，可见光波段的光谱亮度值曲线突变规律性较为明显，在第15像元和第28像元处，也即目视解译中滑体的边界处发生了突变。红外波段的光谱亮度值曲线也具有一定的突变，但曲线形态表现为不规则锯刺状，突变规律不明显，突变点难以确定。因此可见光波段对滑坡边界的检测效果强于近红外、短波红外波段。其次，在可见光波段的曲线波动中，红光曲线具有最大的变化值，这主要是因为研究区滑体周围植被盖度较高，而滑体表面为破碎岩土体，突变点记录的是这二者之间的差异，这与设置红光波段以增强植被探测的目的相一致。

图3 滑体剖面线示意图及各波段光谱剖面 (以“6·24”茂县滑坡为例)Fig.3 Schematic diagram of the section line and the spectral profile of each band (take the “6·24” Maoxian landslide as an example)

对6个案例滑坡区域进行相同地采样和光谱分析后发现，地物反射率数值变化规律以及反射率在滑体边界两侧的突变规律依然存在，说明滑坡灾害对地形地貌的破坏在遥感数据上表现出的规律具有一般性，可用以对滑坡进行识别和提取。

同时，滑坡作为斜坡岩土体沿着贯通的剪切破裂面产生的整体滑移现象，坡度在其中扮演着重要角色。 据统计，“6·24”茂县滑坡灾害所在斜坡平均坡度34.5°，最大坡度57°。西南地区的滑坡多发生在25～55°的坡度范围内[16-17]。为了尽量保留完整滑体的同时限定滑体潜在区域，经反复试验最终选择以10°为下限对遥感数据进行阈值分割(图4)。

图4 茂县滑坡区5°、10°、15°坡度分割效果图Fig.4 Maoxian landslide area 5°, 10°, 15° slope segmentation effect diagram

3 滑体信息提取方法

根据前一节对茂县滑坡区地物光谱特征的分析与光谱变化规律曲线图(图2)，确立了识别规则的初步公式与阈值。应用于其它滑坡区域进行验证，结合验证过程中存在的问题对公式和阈值做出完善与调整，最终确立滑体空间位置识别规则(式1)，其中阈值均来自统计分析。

(1)

式中：Red——红光波段地表真实反射率；

Blue——蓝光波段地表真实反射率；

NIR——近红外波段地表真实反射率；

SWIR2——中红外波段地表真实反射率；

Slope——坡度/(°)。

滑体实质上是综合多种地物要素的复杂地质对象，如基岩、土壤和植被。其复杂性给滑体的完整识别造成困难。因此本文公式(1)作为对滑体的初步提取规则，首先确定滑体空间位置。在此基础上，使用光谱剖面线方法，对红光波段的影像数据进行检测，提取滑体边界。

影响滑体边界的主要因素是扫描起点，扫描半径，扫描间隔，需保证扫描线的均匀分布，使检测出的光谱突变点拥有较好的连续性。经过多次试验，发现使用长宽比小于等于2的外接矩形的形心为扫描线起点，以形心至外接矩形顶点长度为初始扫面半径，以3°～5°为间隔进行扫描时，可以在提高效率的同时最大程度上避免扫描线某处过于密集或稀疏，最终获得的滑体边界效果较好(图5)。

图5 不同角度间隔下光谱剖线检测结果对比图(以茂县为例)Fig.5 Comparison of spectral line detection results at different angular intervals (taking Maoxian as an example)

综合以上两个步骤建立滑体信息提取模型(图6)，使用C#语言基于Microsoft Visual Studio 2017平台完成开发，实现了对滑体信息的自动化提取。

图6 滑体信息提取模型Fig.6 Flow chart of landslide body information extraction

4 结论与讨论

为验证本文提出的滑体提取模型的应用效果，分别对6处滑体进行了提取试验，提取结果见图7。从提取结果看，第一步工作基本圈定滑体位置和大致形状，第二步工作更加详细的勾画出滑体边界，两者结合实现了滑体的快速提取。因本研究所使用的数据为空间分辨率30 m的Landsat 8 OLI数据，混合像元的存在导致模型将滑体与周边地物的混合像元归入到了滑坡区当中，且滑体一般不具有明晰的边界，光谱剖面线法监测的滑体边界突变，将一部分滑坡影响区纳入识别范围，导致各滑体提取面积均合理的在一定程度上超出了前人研究的结果，所以结合既往文献资料中的滑体面积与提取的结果进行对比发现，模型提取面积平均约占文献资料中滑坡面积的109%(表2)。滑体提取形态特征完整，清晰的展现出滑体的“舌”“梨”状形态，可清晰辨别滑体的物源区、流通区和堆积区，提高了滑体信息提取的精度和识别效率，可以为地区滑坡灾害调查与灾后应急救援提供支持。

本文研究的滑坡发生在植被发育良好地区，所使用的影像在灾害发生后一年以内，对发生在干旱半干旱地区的滑坡或是形成许久的老滑坡，因其与周围地物的光谱特征趋同，滑坡识别的效果会显著降低。在研究过程中同时发现，中等分辨率遥感影像对小型滑坡(体积<1×105m3)识别提取的响应较差。随着社会的进一步发展，滑坡灾害调查和灾后应急救援工作将会对滑体信息提取提出更加快速和精准的要求，对于老滑坡、干旱半干旱地区的滑坡、复杂环境下的滑坡，如何综合多源数据，进一步优化模型，放大和捕捉滑坡与非滑体在光谱特征上的差异，实现更加普适和高精度的滑坡自动识别，是后续要开展的研究工作。

图7 滑体遥感影像与模型提取结果对比图Fig.7 Comparison of remote sensing images and model extraction results of landslides

表2 案例区模型提取滑坡面积比对表