结合纹理特征分析与比辐射率估计的震后滑坡提取[PDF全文]
杨明生, 李山山, 冯钟葵
摘要: 滑坡提取是获取滑坡灾情信息的首要环节。针对背景复杂的震后灾区,为了更加有效地增强裸地信息并降低土壤背景差异对NDVI植土分离能力的影响,利用NDVI进行滑坡提取时结合纹理特征分析和比辐射率估计来提高滑坡提取的精度。首先,根据纹理特征差异将主导背景因子不同的区域划分为独立的处理单元,以避免全局滑坡提取所出现的过提取或欠提取现象;然后,利用NDVI估算比辐射率,以增强裸地信息;最后,通过阈值分割实现更高精度的滑坡提取。围绕5·12大地震后的汶川县城及其周边区域,采用Landsat 5 TM和Terra ASTER数据展开实验,该方法的滑坡提取结果与人工提取结果吻合度较好,Kappa系数分别为0.8531和0.9271,优于基于NDVI的全局阈值分割法和其他一些典型的监督分类法。实验结果表明,在背景复杂的灾区环境下,结合纹理特征分析和比辐射率估计的滑坡提取方法能明显降低漏检率和错检率,进而有效提高滑坡提取的精度,且对中等分辨率的不同遥感数据源的适用性较强。
关键词: 滑坡提取     比辐射率     纹理分析     灾害遥感     复杂背景     阈值分割    
DOI: 10.11834/jrs.20187261    
收稿日期: 2017-06-28
中图分类号: TP79/P642.22    文献标识码: A    
作者简介: 杨明生,1992年生,男,硕士研究生,研究方向为遥感数据处理与系统。E-mail:yangmingsheng15@mails.ucas.ac.cn
通信作者: 李山山,1979年生,男,副研究员,研究方向为遥感图像处理与遥感应用。E-mail:liss@radi.ac.cn.
基金项目: 中国科学院遥感与数字地球研究所所长青年基金(A类)资助
Post-seismic landslide extraction by combining texture analysis and emissivity estimating
YANG Mingsheng, LI Shanshan, FENG Zhongkui
1. Institute of Remote Sensing and Digital Earth, Chinese Academy of Sciences, Beijing 100094, China
2. University of Chinese Academy of Sciences, Beijing 100049, China
Abstract: Landslides are the most common geological disasters that result from large seismic activities and heavy rainfall in mountainous areas. They are destructive events that occur suddenly and have a wide distribution range. Landslide extraction is the primary factor in collecting destruction information and plays a key role in disaster prevention and emergency rescue. This paper presents a method that combines texture analysis and emissivity estimation to extract landslides from different backgrounds. Many researchers have focused on landslide extraction and recognition, and they proposed methods based on the ability of the Normalized Difference Vegetation Index (NDVI) to separate vegetation. However, although NDVI easily extracts landslides from high vegetation areas, it is affected by soil interference in sparse vegetation areas. Thus, extracting landslides is difficult. An NDVI-based method is proposed for landslide extraction in a complicated environment and to make NDVI effective in sparse vegetation zones. In accordance with the research area and remote sensing data, this method consists of four steps on the assumption that image preprocessing has been completed. First, the NDVI image is calculated by using near-infrared and red bands; this approach shows that some local regions have a sharp contrast, whereas others have a weak contrast. Second, texture analysis is conducted to divide the research area into several blocks according to this contrast distribution. NDVI mean and near-infrared Angular Second Moment (ASM) based on Gray-Level Co-occurrence Matrix (GLCM) could be selected as the texture feature to segment zones easily by providing simple thresholds. A continual large vegetation area is outstanding in the Mean feature, and a mountainous landslide area is rough in the ASM feature. The remaining area, except the first two, varies in the ASM feature related to soil components. Thus, the combination of Mean and ASM features facilitates texture analysis. Third, emissivity is estimated in different blocks based on the NDVI image. In this section, the NDVIs of pure soil and pure vegetation from statistics of the NDVI image in different blocks are critical parameters in calculating the Percentage of Vegetation (PV), thereby obtaining the emissivity. In a limited area, emissivity contrasts within different objects, especially for soil and vegetation, thereby making it more suitable than NDVI for landslide extraction. Finally, the landslides are extracted through emissivity threshold segmentation technology with an appropriate threshold manually. Landsat 5 TM and Terra ASTER data are tested by this method, and the result is favorable given that it is more consistent with that of artificial landslide extraction than that of the other methods, such as maximum likelihood supervised classification, neural net supervised classification, support vector machine-based supervised classification, and NDVI global threshold segmentation. For objective and quantitative evaluation, the result of artificial extraction is considered as the ground-truth data to calculate the confusion matrix of other results, with the Kappa coefficient being used to demonstrate the performance of each method. As a result, the method described in this article achieves a high Kappa coefficient, namely, 0.8531 (TM) and 0.9271 (ASTER). By contrast, maximum likelihood classification achieves 0.7634 (TM), neural net classification achieves 0.66 (TM), SVM-based classification achieves 0.6896 (TM), NDVI global threshold (0.3) segmentation achieves 0.622 (TM), and NDVI global threshold (0.5) segmentation achieves 0.7487% (TM). Evidently, this method can effectively eliminate omission and misclassification. With the increase in the resolution of remote sensing data, ASTER (15 m) provides a better result compared with TM (30 m), thereby showing that this method does not rely on data sources, and a high resolution contributes to a good result. Comparative analysis between different methods and data sources indicates the following results: The NDVI mean and near-infrared ASM are good texture features for separating different background environment blocks, especially with soil as the leading factor. Emissivity reflects the spatial difference in objects, although it is estimated by NDVI. The estimated emissivity not only separates vegetation from NDVI but also increases the difference between soils, thereby resulting in good landslide extraction. The advantages of texture analysis and emissivity estimation enhance landslide extraction in complicated research areas. Finally, the threshold selection problem in NDVI global threshold segmentation method is solved, and the sample selection and learning process in supervised classification are avoided. This method is designed for a medium resolution that ranges between 10 m to 50 m; thus, this method can be used with other data sources with a similar resolution, such as Landsat 5 TM and Terra ASTER. We proposed an effective method for landslide extraction for medium resolution. To handle high-resolution remote sensing data, this method will be combined with object-oriented methods in follow-up studies for accurate landslide extraction. The proposed method has a strong dependency on manual threshold selection; this dependency is inconvenient in the auto-extraction process. Therefore, self-adaptive threshold extraction is another challenge in future work.
Key Words: landslide extraction    emissivity    texture analysis    disaster remote sensing    diverse background    threshold segmentation    
1、引 言

本文所指的滑坡均为广义滑坡。它包括狭义滑坡、崩塌、泥石流和碎屑流等,是指在重力侵蚀作用下构成斜坡的物质向下或向外移动的现象(王治华,2012)。滑坡是最为常见的自然灾害之一,也是地震等地质构造活动中极易诱发的次生灾害,对人类的生存环境造成极大的威胁,一直是灾害防治工作中的重中之重。

遥感技术具有广覆盖、高时效、低成本、低风险等优越性,国内外将其应用于滑坡灾害调查与研究已有50多年。范一大等人(2016)从实际应用出发,总结了不同灾害要素监测的时空尺度范围,其中针对滑坡灾害的遥感数据最大适宜空间尺度是50 m。因此,中等及更高分辨率的遥感数据是进行滑坡提取的重要数据来源。

对于灾情分布区域较大的滑坡提取,研究者会根据实际应用需求优先考虑分辨率一般介于10—50 m的中等或中等稍高分辨率的遥感数据,如Landsat 5 TM、Landsat 7 ETM+、Landsat 8 OLI、Terra ASTER、SPOT 4 HRVIR等。针对中等或中等稍高分辨率的遥感数据,滑坡提取主要以基于像元的方法为主。傅文杰和洪金益(2006)基于Landsat 7 ETM+数据首次提出用支持向量机SVM(Support Vector Machine)进行滑坡提取,许高程等人(2009)将该方法应用到Landsat 5 TM数据上,李松等人(2015)又将该方法应用到Formosat-2多光谱数据上,一系列的研究证明了该方法的有效性和可行性。许冲(2013)以北川县为例分析了最大似然法阈值选取对滑坡提取的影响。Manfré等人(2016)利用Landsat 5 TM数据对比分析了不同遥感图像分类方法及其组合方法提取滑坡的有效性。除此之外,有学者利用多时相变化检测结果进行滑坡提取,也取得较好的效果(Li和Hua,2009Tang 等,2010Mondini 等,2011许积层 等,2012);但该方法容易引入非滑坡变化噪声,同时需要灾前数据,在灾害应急调查时可能面临无灾前数据可用的情况。

在进行滑坡提取时,目前大多数研究均基于滑坡引起植被变化这一假设。利用滑坡与背景环境之间存在较大反差的特点,使得滑坡提取的关键在于提取裸地(Yang 等,2014)。但是滑坡区域在图像上具有明显的裸地特征,与居民地、道路等在很多方面都难以区分。多数学者从人居环境选址的实际情况出发,利用DEM等辅助数据对居民地、道路等人工地物作滤除或者掩膜处理。

在区分裸地和植被时,多数文献使用到归一化差分植被指数NDVI(Normalized Difference Vegetation Index),利用其在高植被区分离植被的优势,在滑坡提取上取得了一定的效果。杨树文等人(2012)在NDVI的基础上提出了常量化修正的土壤调节植被指数CMSAVI(Constant Modified Soil Adjusted Vegetation Index),有效地减少了土壤背景的影响,并在SPOT 5数据上取得了较好的滑坡提取效果,但该方法涉及的参数设置对数据源的依赖性较强。

在面对灾情分布范围较广且包含多种地物背景的滑坡区域进行滑坡提取时,高植被区、稀疏植被区、裸地、居民地等往往交替分布,许多研究方法容易受到土壤背景因素的影响。特别是在植被覆盖率较低的稀疏植被区,土壤背景的干扰非常强烈(唐怡 等,2006),单纯利用NDVI全局阈值分割分离裸地和植被的效果并不理想,往往导致滑坡提取精度较低。

图像纹理特征与地物背景类型密切相关,利用纹理特征差异将研究区划分为多个背景条件相对单一的局部处理单元,能有效避免全局阈值分割的过提取或欠提取现象。同时,与NDVI图像相比,比辐射率的图像差异更大,尤其表现在稀疏植被区;采用比辐射率可以实现裸地增强,使其更加有利于分离滑坡。

因此,为了能保留NDVI在高植被区分离植被的优点,并在稀疏植被区实现裸地增强以弥补NDVI分离能力被削弱的缺点,本文提出了结合纹理特征分析与比辐射率估计的滑坡提取方法,通过滑坡分区和局部阈值调整实现更高精度的滑坡提取。

2、研究区及数据     (2.1) 研究区概述

2008年5月12日,四川省阿坝藏族羌族自治州汶川县发生了有记录以来最强烈的大地震,震中位于31° N、103.4° E,造成了非常严重的人员伤亡和经济损毁。受强烈的主震及反复多次余震的影响,震后引发了大量山体滑坡,又恰逢雨季的到来,雨水的强烈冲刷使滑坡灾情更为严重。时至今日,仍有较多滑坡裸露表面尚未得到较好的恢复。据统计,地震诱发的地质灾害点有3万—4万处,受灾面积达44万km2(魏永明 等,2014),其中重大滑坡灾害点有300余处(黄润秋,2008),具有数量大、分布广、群发性、影响持久等特点。本文选取的研究区位于汶川县境内,面积约169 km2,以暖温带大陆性半干旱季风气候为主,降雨集中在5—8月,平均海拔高于2000 m,县城及其附近地势较平缓,植被较稀疏,其余地方大都属于植被茂盛的高山峡谷地带,地质构造复杂,地处龙门山断裂带,属于滑坡灾害多发区(苏凤环 等,2008),又是汶川县城所在地,一直是灾害防治的重要监测区域。

研究区独特的地形地貌、气候条件以及地质构造活动,使其在较长的一段时间内成为滑坡常发地,为开展滑坡灾害研究提供理想场所。如图1所示,选取红色方形部分作为研究区。

图 1 研究区地理位置 Figure 1 Geographic location map of the study area
    (2.2) 实验数据说明

本文使用到的实验数据有3类:卫星影像数据、数字高程数据和高分辨率航空影像数据。

         2.2.1. 用于提取滑坡的卫星影像数据

基于当前滑坡提取的研究现状,结合具体的应用需求,考虑到灾情分布较广以及数据的可用性与易用性,本文利用Landsat 5 TM数据开展实验。该数据来源于美国地质调查局(USGS)官网,轨道号为WRS-2 Path130/Row38,空间分辨率为30 m,成像于2008年7月18日,覆盖整个研究区,且在研究区内无云覆盖。数据产品等级为L1T,相当于国内的4级卫星数据产品,经过了辐射校正、系统几何校正、地面控制点几何精校正以及地形校正(冯钟葵 等,2016)。Landsat 5 TM数据有7个波段,如表1所示,本文主要使用到红波段和近红外波段。为了进一步消除大气的影响,利用ENVI的FLAASH模块对该数据进行了大气校正。

表 1 Landsat 5 TM和Terra ASTER的波段设计 Table 1 Band design of Landsat 5 TM and Terra ASTER

为了检验本文所提出方法在不同数据源之间的适用性,本文利用同一地区的同期Terra ASTER数据再次进行实验,具体的波段分布范围如表1所示。该数据来源于美国地质调查局(USGS)官网,轨道号为WRS-2 Path130/Row38,空间分辨率为15 m,成像于2008年5月23日,覆盖整个研究区,且在研究区内无云覆盖,其产品等级为L1T。为了便于比较分析,需将Terra ASTER数据配准Landsat TM数据,均方根误差小于0.5个像元。

         2.2.2. 数字高程数据DEM

人工地物一般分布在坡度较小的平缓地带,DEM能反映地形的起伏变化特点,因此,可根据DEM计算坡度,然后将其用于剔除道路、居民地等人工地物。ASTER GDEM是目前唯一覆盖全球陆地表面的高分辨率高程影像数据,其分辨率为30 m,本文中的地形数据采用ASTER GDEM V1数据。根据Lu等人(2012)的研究成果,本文将坡度小于25°的区域看作非滑坡区,在数据处理时制作掩膜以消除人工地物的干扰。

         2.2.3. 用于目视解译的高分辨率航空影像数据

本文选取与滑坡提取所采用卫星影像数据时相相近的航空影像数据作为人工提取滑坡时目视解译的依据,其分辨率为2 m,成像于2008年5月15日。目视解译主要用于验证阶段辅助人工提取滑坡以实现精度分析。除此以外,该数据还被用于制作掩膜以滤除与滑坡高度混淆的种植作物区。

3、研究方法

本文提出的一种结合纹理特征分析和比辐射率估计的滑坡提取方法流程如图2所示。对经过辐射、几何和大气校正处理后的中等分辨率遥感卫星数据,首先采用纹理分析和阈值分割方法,根据其近红外数据和NDVI数据的纹理特征差异将研究区划分为若干局部分区,使各分区具有相近的背景条件,这些条件主要体现在土壤和植被的成分及其构成比例上。然后,在各分区内基于NDVI估算比辐射率,使各分区的比辐射率具有较为一致的图像差异,从而实现整个研究区单一阈值的滑坡提取。

图 2 本文方法的具体流程图 Figure 2 Flow chart of the proposed method
    (3.1) 滑坡与背景地物的NDVI图像差异

NDVI是常用的遥感指数之一,它利用植被在近红外波段具有高反射这一特性,可将植被从水体和土壤中分离出来的。其表达式(孙家抦,2009)如下:

${\rm{NDVI}} = \frac{{{\rho _ {\rm{nir}}} - {\rho _{{\rm{red}}}}}}{{{\rho _ {\rm{nir}}} + {\rho _{{\rm{red}}}}}}$ (1)

式中, ${\rho _{{\rm{nir}}}}$ ${\rho _{{\rm{red}}}}$ 分别为近红外波段和红光波段的反射率。

滑坡提取的主要目标在于利用两者的图像差异分离滑坡和背景地物。背景地物由土壤、植被、人工地物等构成,成分复杂,是影响NDVI图像差异的关键因素。如图3所示,在汶川城区集聚地带和城外沿河分布地带,植被较稀疏,覆盖率较低,土壤背景突出,导致该环境下滑坡与背景的NDVI图像差异小,其NDVI分离阈值低;在城外高山峡谷地带,植被茂密,覆盖率高,土壤背景弱,导致该环境下滑坡与背景的NDVI图像差异较大,其NDVI分离阈值高。对以土壤为主导背景因子和以植被为主导背景因子的滑坡区域进行取样,基于Landsat 5 TM数据分别统计各样区NDVI的均值和标准差,得到表2。由表2可见,土壤背景滑坡样区的均值和标准差均小于植被背景滑坡样区,说明在不同主导背景因子下局部区域之间图像差异程度不一致,这也反映在2 m高分辨率的局部自然真彩色示意图上。

图 3 近红外波段和NDVI图像的纹理特征 Figure 3 Texture feature of near infrared and NDVI

表 2 不同主导背景因子下滑坡样区的NDVI统计 Table 2 Statistics of NDVI in different landslide samples

这种各局部之间图像差异程度的不一致在近红外影像和NDVI图像相结合的纹理特征上得到较好体现。如表3所示,按照主导背景因子类型和居民地分布差异,将整个研究区划分为城外高山峡谷地带、城外沿河分布地带和汶川城区集聚地带。在近红外影像上,纹理平滑的区域主要分布在城外高山峡谷地带的无滑坡区、城外沿河分布地带和汶川城区集聚地带;纹理粗糙的区域主要分布在城外高山峡谷地带的有滑坡区。虽然在近红外影像纹理特征上无法区分城外高山峡谷地带的无滑坡区与土壤背景占优势的滑坡区域,但在NDVI图像纹理特征上,高山峡谷地带的无滑坡区亮度最高,与其他区域形成鲜明比照,可利用这一特点将其独立分离出来。这为复杂背景条件下根据局部主导背景因子不同进行分区提取滑坡指明了方向。

表 3 不同分区的背景特点及其纹理特征 Table 3 Background characteristics and texture feature of different areas
    (3.2) 纹理特征分析

灰度共生矩阵GLCM(Gray Level Co-occurrence Matrix)是通过统计图像上一定大小 $m$ 的基准窗口内某一 $\theta $ 方向上相距为 ${d}$ 的两个灰度级 $i$ $j$ 同时出现的联合概率分布而获得的对称矩阵 ${P}_{m, \theta, d}(i, j)$ ,能表现基准窗口中心点附近像素灰度的空间相关性。以灰度共生矩阵 ${P}_{m, \theta, d}(i, j)$ 为基础作进一步统计分析,并得到相应的纹理描述子,便是基于灰度共生矩阵的纹理分析,它是有效表达纹理特征的常用方法。由Haralick等人提出的描述子有14种(冯建辉 等,2007),本文用到均值 $Mean$ 和角二阶矩ASM(Angular Second Moment),其计算公式如下:

$Mean = \sum\limits_{i = 0}^k {\sum\limits_{j = 0}^k {i \cdot {P_{m, \theta, d}}(i, j)} } $ (2)
${\rm{ASM}} = \sum\limits_{i = 0}^k {\sum\limits_{j = 0}^k {P_{_{m, \theta, d}}^2(i, j)} } $ (3)

式中, $k$ 为最大灰度级。 $Mean$ 主要反映局部区域纹理的平均亮度信息,灰度平均值越大,其值越大,反之越小。ASM是图像灰度分布均匀性的度量,局部区域内图像灰度分布越均匀,其值越大,反之越小。

结合3.1.1节滑坡与背景的NDVI图像差异分析,可根据图4所示流程进行区域划分。首先,将 $Mean$ 用于分离城外高山峡谷地带的无滑坡区,然后利用ASM分离城外高山峡谷地带的有滑坡区,合并成城外高山峡谷地带作为以植被为主导背景因子的区域,剩下的则是以土壤为主导背景因子的区域。以上根据不同纹理特征,通过简单的阈值分割便可界定不同背景条件的区域分布,本文将这一过程称之为纹理特征分析。

图 4 纹理特征分析流程图 Figure 4 Flow chart of texture analysis

在进行纹理特征分析时,以土壤为主导背景因子的区域还需要根据土壤背景的复杂程度确定适宜的分区数量。土壤背景的复杂程度取决于裸地、居民地、道路等土壤成分的多少。在本文所选研究区内,汶川城区一带居民地、道路等较为集中,土壤背景最复杂;城区以外的岷江两岸的居民地、道路等较少,土壤背景复杂度次之;周边其他区域人烟稀少,植被茂盛,主要以植被背景为主。

    (3.3) 比辐射率估计

比辐射率是热红外波段遥感数据反演地表温度的关键物理参量,已被广泛应用到城市热岛分析。李霞等人(2016)将其应用到裸地提取上,首先通过比辐射率实现热红外亚像元分解,再利用该热红外数据计算不透水面指数NDISI(Normalized Difference Impervious Surface Index),最后结合土壤指数NDSI(Normalized Difference Soil Index)实现城市裸地提取。但在该方法中比辐射率只是用于细化热红外数据,并没有直接用于分离地物。事实上,覃志豪等人(2004)在提出基于NDVI的比辐射率估算方法时就指出在一定范围的局部区域内,比辐射率具有反映地物空间差异的重要作用。因此,在有限的局部区域内可以直接利用估算的比辐射率分离地物。

根据覃志豪等人(2004)的研究成果,比辐射率估算方法主要包括统计NDVI纯净值、计算植被构成比例PV(Proportion of Vegetation)、计算比辐射率3个步骤。

         3.3.1. 统计NDVI纯净值

本文将各分区自然表面看作是由不同比例的植被叶冠和裸土组成的混合像元。在计算混合像元植被构成比例PV时,首先需要给定裸土纯净值NDVIS和植被纯净值NDVIV。这两个参数的确定具有重要意义,将直接影响PV计算结果。在实际应用中因缺少大面积地表实测数据作为参考,很难从影像上获取到纯净像元值,所以通常结合NDVI统计直方图和给定置信区间进行取值,本文采用NDVI累计直方图频率约为5%和95%时的NDVI值分别作为NDVIS和NDVIV的值。

         3.3.2. 计算PV

植被构成比例(PV),也叫植被覆盖度,是反映混合像元中植被成分含量的重要参量。纯净裸土像元PV=0,表示不含植被成分。纯净植被像元PV=1,表示完全的植被成分,不含其他杂质。实际上,纯净像元是理想化模型,在实际应用中,只要 ${\rm{NDVI}} \leqslant {\rm{NDVI}}_{\rm{S}}$ 的像元均被视为纯净裸土, ${\rm{NDVI}} \geqslant {\rm{NDVI}}_{\rm{V}}$ 的像元均被视为纯净植被。根据像元二分模型的原理,PV值与NDVI之间的关系(李苗苗 等,2004)如下:

${\rm{PV}} = \left\{ {\begin{array}{*{20}{l}} 0& {\rm{NDVI}} \leqslant {\rm{NDVI}}_ {\rm{S}} \\ \displaystyle\frac{{{\rm{NDVI}} - {\rm{NDVI}}_ {\rm{S}}}}{{{\rm{NDVI}}_ {\rm{V}} - {\rm{NDVI}}_ {\rm{S}}}}& {\rm{NDVI}}_ {\rm{S}} < {\rm{NDVI}} < {\rm{NDVI}}_ {\rm{V}}\\ 1& {\rm{NDVI}} \geqslant {\rm{NDVI}}_ {\rm{V}} \end{array}} \right.$ (4)

从式(4)可以看出,PV是基于 ${\rm{NDVI}}_ {\rm{S}}$ ${\rm{NDVI}}_ {\rm{V}}$ 这两个参数所做的线性拉伸,两者的反差越小,其拉伸越大,反之越小。其作用是将大气、土壤背景与植被类型等的影响降至最低(郭伟伟 等,2012)。在各局部分区中, ${\rm{NDVI}}_ {\rm{S}}$ ${\rm{NDVI}}_ {\rm{V}}$ 是分区统计的,NDVI图像差异小的分区拉伸大,NDVI图像差异大的分区拉伸小,使整个研究区的PV差异程度趋向一致。

         3.3.3. 计算比辐射率 $\varepsilon $

纯净像元的比辐射率一般选取经验值,而混合像元的比辐射率则与纯净像元、植被构成比例PV、地形因子等密切相关,其关系式(覃志豪 等,2004)如式(5):

$\varepsilon = \left\{ \begin{gathered} {\varepsilon _{\rm{S}}}(1 - PV) + {\varepsilon _{\rm{V}}}PV + 0.003796PV\;\;\;\;\;\;\;\;\;\;, PV \leqslant 0.5 \\ {\varepsilon _{\rm{S}}}(1 - PV) + {\varepsilon _{\rm{V}}}PV + 0.003796(1 - PV)\;\;\;, PV > 0.5 \\ \end{gathered} \right.$ (5)

式中, $\varepsilon $ 为估算的比辐射率, ${\varepsilon _ {\rm{S}}}$ ${\varepsilon _ {\rm{V}}}$ 分别为纯净裸土和纯净植被的比辐射率经验值。此处根据覃志豪等人(2004)的研究成果,选取经验值 ${\varepsilon _ {\rm{S}}} =$ 0.97215、 ${\varepsilon _ {\rm{V}}} = 0.986$

    (3.4) 基于阈值的滑坡提取

比辐射率是在假定研究区由植被、裸土及其混合像元构成的基础上估算得来的,对植被和裸地具有较强的分离能力。滑坡在本文中可看作较为纯净的裸地,因此,可以选择合适的阈值将其从植被或者植被与土壤的混合像元中提取出来。阈值的确定主要依靠滑坡分布情况来判定,带有一定实践经验。

4、结果与讨论

根据上述方法,本文在ENVI平台上基于Landsat 5 TM数据(大小为433像素×433像素)进行了实验探究,通过实验结果分析纹理特征和比辐射率的作用,并将滑坡提取结果与相同数据下其他方法的提取结果进行定性与定量分析,以评价本文方法的提取精度。同时,本文还基于同一地区的同期Terra ASTER数据(大小为866像素×866像素)再次进行了实验探究,通过将同一方法应用到不同数据源(Landsat 5 TM和Terra ASTER)上,以说明本文方法对不同数据源的适用性。

    (4.1) 纹理特征分析与比辐射率估计的优点          4.1.1. 纹理特征分析的优点

纹理特征计算过程需要提供基准窗口大小、移动步长和灰度量化等级这3个参数,移动步长均设为1×1像素,灰度量化等级均设为16,基准窗口大小为奇数,且与遥感数据空间分辨率密切相关,分辨率每提高一倍,基准窗口大小也相应地增大一倍,本文经过反复试验,最终确定30 m空间分辨率的基准窗口大小为55像素×55像素(Landsat 5 TM),相应地,15 m分辨率的基准窗口大小为109像素×109像素(Terra ASTER)。基于纹理特征图像 $Mean$ 和ASM,按照实践经验分别选取不同阈值实现Landsat 5 TM和Terra ASTER数据的纹理特征分析,得到如图5所示的结果。无论是Landsat 5 TM还是Terra ASTER数据,汶川城区集聚地带、城外沿河分布地带等以土壤为主导背景因子的区域大部分被划分到分区Ⅰ和分区Ⅱ上,且分区Ⅰ主要以汶川城区集聚地带为主,分区Ⅱ主要以城外沿河分布地带为主,分区Ⅰ的土壤复杂度较分区Ⅱ高;城外高山峡谷地带等以植被为主导背景因子的区域大部分被划分到分区Ⅲ上。对比3.1.1小节中所得的分区情况,可以看出实验所得纹理分区与研究区的背景类别具有较高的一致性。

图 5 纹理特征分析结果 Figure 5 Result of texture analysis

事实上,通过纹理特征分析,将整个研究区划分为若干个局部分区,其中,以植被为主导背景因子的区域被划分为一个分区,但为了体现土壤因素对NDVI分离能力抑制程度的差异,以土壤为主导背景因子的区域则根据土壤背景的复杂度将其细分为一个或多个分区。这种分区滑坡提取方法充分考虑了局部区域之间的背景条件差异,尤其是土壤背景,以背景条件相近的局部代替全局作为滑坡提取的处理单元,能有效避免全局滑坡提取所出现的过提取或欠提取现象。

         4.1.2. 比辐射率估计的优点

以Landsat 5 TM数据为例,图6(a)是通过近红外和红光波段计算得到的NDVI,图6(b)则是利用NDVI估算得到的比辐射率。通过两者的比较,本文发现比辐射率比NDVI能更加有效地增强滑坡信息。

图 6 比辐射率和NDVI的比较 Figure 6 Compare between emissivity and NDVI

NDVI整体反差较小,其图像明暗对比不强烈,且沿河滑坡区和高山峡谷地带滑坡区的NDVI图像差异程度较大,沿河滑坡的NDVI主要分布在0.3以下,高山峡谷地带滑坡区的NDVI则主要分布在0.5以下,导致整体上采用单一阈值分离滑坡的效果不理想。比辐射率整体反差大,其图像明暗对比强烈,且不同滑坡区域比辐射率分布具有较好的一致性,大都介于0.97215—0.97450,不管是土壤背景占优还是植被背景占优的滑坡区域,它们的比辐射率图像差异较为一致,能在整体上获得较好的滑坡分离效果。

事实上,覃志豪等人(2004)指出,在背景条件相似的局部范围内,比辐射率的空间差异是地表温度反演的决定性参数,它反映了物体的热辐射特性。利用比辐射率提取目标地物实际上是根据热辐射差异分离地物,即便是在土壤背景复杂的条件下,也能较好地区分各种类型的目标对象。由于易受土壤因素影响,NDVI在植被稀疏的汶川城区集聚地带和城外沿河分布地带的分离能力受到不同程度的抑制,相比而言,比辐射率整体反差更大且在不同局部之间图像差异较为一致,能有效地增强滑坡信息,从而有利于提高滑坡提取的精度。

    (4.2) 滑坡提取精度的定性与定量分析

为了检验滑坡提取精度,本文在获得该区域滑坡前后变化检测结果的基础上,结合高分航空遥感数据进行目视解译,修正变化检测结果,并以此作为人工提取结果(图7),用于对比分析其他方法的滑坡提取效果。为了从视觉上突出滑坡区域,本文将提取到的滑坡设置成黄色,并叠加到标准假彩色影像上。

图 7 人工提取结果 Figure 7 Result of artificial extraction
         4.2.1. 相同数据不同方法的比较

以Landsat 5 TM数据为例,将本文方法与基于NDVI全局阈值分割法、多种监督分类法的提取结果进行比较分析。基于NDVI全局阈值分割法分别选择阈值为0.3和0.5进行全局阈值分割。监督分类法选择NDVI、band1—5、band7共7个数据作为输入特征量,训练样本分为滑坡和非滑坡两类;其中,参考许冲(2013)的研究成果,基于最大似然监督分类法的似然度阈值(probability threshold)设为none;参考Manfré等人(2016)的研究成果,基于神经网络监督分类法选择对数活化函数、训练贡献阈值(training threshold contribution)设置为0.9、权重调节速率(training rate)设置为0.2、训练终止RMS误差标准(training RMS exit criteria)设置为0.1、隐藏层数量(number of hidden layer)设置为1、迭代次数设置为1000;参考许高程等人(2009)的研究成果,基于SVM监督分类法选择二次多项式核函数模型、偏置(Bias)设为1、因子(Grama)设为输入特征量个数的倒数、惩罚因子设为100。

图8的对比中可以定性地看到不同提取方法的提取效果。本文方法的提取结果(图8(a))与人工提取结果(图7)最为吻合,尤其表现在土壤背景最为复杂的汶川县中心城区沿河一带。基于最大似然监督分类方法也取得了较好的结果(图8(b)),但中心城区沿河一带与人工提取结果存在较大的差异。基于神经网络的监督分类方法(图8(c))和基于SVM的监督分类方法(图8(d))提取效果并不理想,中心城区沿河一带遗漏现象较为严重。基于NDVI全局阈值法的提取效果也不理想,当阈值取0.3时(图8(e)),能较好地提取到中心城区沿河一带的滑坡区域,但在提取高山峡谷地带的滑坡时,遗漏现象严重;当阈值为0.5时(图8(f)),能较好地提取到高山峡谷地带的滑坡区域,但对于中心城区沿河一带,其提取结果包含了大量的非滑坡区域。

图 8 不同方法的提取结果(TM) Figure 8 Results of different methods(TM)

为了能够定量分析滑坡提取的精度以及有效性,本文将人工提取的结果看作地面真实数据,然后根据其他方法的提取结果,分别计算混淆矩阵,作为评价各方法性能的指标来源。如表4所示,本文采用Kappa系数进行客观评价本文方法的性能,并与其他方法进行比较。相对于基于最大似然的监督分类法、基于神经网络的监督分类法、基于SVM的监督分类法和基于NDVI的全局阈值分割法,本文方法的提取精度有了明显提高。在这些方法中,基于NDVI的全局阈值分割法容易出现过提取和欠提取现象,导致提取精度不高,其平均Kappa系数为0.6854。基于神经网络的监督分类法和基于SVM的监督分类法不仅耗时长,而且欠提取现象较为严重,其Kappa系数分别为0.6600和0.6896。基于最大似然监督分类法的整体提取效果较好,但局部过提取现象仍较明显,其Kappa系数为0.7634。本文方法在充分考虑背景复杂性的基础上,以估算的比辐射率代替NDVI进行阈值分割,其精度大大提高,其Kappa系数达到0.8531。这说明本文方法能在背景复杂的环境下降低错分率与漏分率,从而实现更为有效的滑坡提取。

表 4 不同方法的提取精度(TM) Table 4 Accuracy of different methods (TM)
         4.2.2. 不同数据同一方法的比较

本文方法是针对中等或中等稍高分辨率遥感数据而设计的,理论上能推广到类似分辨率的其他数据源上。为了证明该方法的适用性,本文将该方法应用到同一地区时相相近的分辨率为15 m的Terra ASTER数据上,并取得了较好的提取结果,如图9所示。对比图7图8(a),基于Terra ASTER数据的滑坡提取结果与人工提取结果、同期Landsat 5 TM提取结果吻合度都比较高,且从Terra ASTER和Landsat 5 TM提取结果的局部放大图的对比中可以看出,Terra ASTER数据提取到的滑坡更加精细,细节更加丰富。基于Terra ASTER数据的Kappa系数为0.9271,相比Landsat 5 TM数据的0.8531,精度更高。同时,本文也将该数据应用到基于NDVI全局阈值分割法以及多种监督分类方法上以便于比较分析,其滑坡提取精度如表5所示。相比其他方法,本文所提方法在Terra ASTER数据上的滑坡提取精度有较大提升,优于其他方法。由此可见,本文方法具有较好的数据适用性,且分辨率的提高有利于滑坡提取精度的提高。

图 9 Terra ASTER数据滑坡提取结果 Figure 9 Extraction result of Terra ASTER

表 5 不同方法的提取精度(ASTER) Table 5 Accuracy of different methods (ASTER)
5、结 论

复杂多变的地物对象是影响滑坡提取精度的重要因素。滑坡具有明显的裸地特征,针对滑坡提取,必须较好地消除背景环境的干扰,尤其是以土壤为主导背景因子的稀疏植被区。本文在基于NDVI全局阈值分割方法的基础上,考虑到局部主导背景因子的差异性而引入局部分区的思想,了解到NDVI在植被稀疏区分离能力受抑制而引入比辐射率以增强裸地信息,最终顾及这两种考虑提出了结合纹理特征分析和比辐射率估计的滑坡提取方法,并在稀疏植被区和高植被区共存的复杂研究区开展实验,结果表明:

(1)结合近红外波段和NDVI的纹理特征,能较好地划分各种类型的背景环境,尤其是植被稀疏的滑坡地区,这些地区主要以土壤为主导背景因子。

(2)比辐射率在一定的局部区域内具有反映地物空间差异的重要作用,它能较好地表征地物之间的对比差异,尤其是土壤和植被。估算的比辐射率既能有效保留NDVI分离植被的优越性能,又能增强裸地信息,从而提高滑坡提取能力。

(3)复杂背景下结合纹理特征分析和比辐射率估算提取滑坡区域,既能明显改善基于NDVI全局阈值分割法整体提取效果不理想的现状,也能避免监督分类样本选取、训练等耗时过程。这是一种高效的基于像元的滑坡提取方法,其精度优于基于NDVI全局阈值分割法、基于最大似然监督分类法、基于神经网络监督分类法、基于SVM监督分类法等常见滑坡提取算法,且对数据源的依耐性弱,适用于设有近红外波段和红光波段的中等或中等稍高分辨率的不同来源的遥感数据。

本文提出了一种结合纹理特征分析和比辐射率估计的滑坡提取方法,并以汶川县城及其周边一带山区震后滑坡为例展开实验并取得良好效果,较好地解决了稀疏植被区和高植被区共存的复杂灾区环境下的滑坡提取。该方法主要针对背景复杂的滑坡区和中等或中等稍高分辨率的遥感数据。对于背景单一的研究区,并没有使用纹理特征分析的必要;对于更高分辨率的遥感数据,在以后的工作中可以考虑将本文方法和面向对象的方法相结合应用到更高分辨率的遥感数据中。除此以外,本文方法在纹理特征分析和基于比辐射率的阈值提取过程均依赖实践经验选择阈值,能否实现阈值的自适应是实现滑坡自动化提取的关键,这需要在后续的研究中作深入的挖掘与探讨。

参考文献
[1] 范一大, 吴玮, 王薇, 刘明, 温奇. 中国灾害遥感研究进展[J]. 遥感学报, 2016, 20 (5) : 1170 –1184. Fan Y D, Wu W, Wang W, Liu M and Wen Q. Research progress of disaster remote sensing in China[J]. Journal of Remote Sensing, 2016, 20 (5) : 1170 –1184. DOI: 10.11834/jrs.20166171
[2] 冯建辉, 杨玉静. 基于灰度共生矩阵提取纹理特征图像的研究[J]. 北京测绘, 2007, 2007 (3) : 19 –22. Feng J H and Yang Y J. Study of texture images extraction based on gray level co-occurence Matrix[J]. Beijing Surveying and Mapping, 2007, 2007 (3) : 19 –22. DOI: 10.3969/j.issn.1007-3000.2007.03.006
[3] 冯钟葵, 葛小青, 张洪群, 李安. 2016. 遥感数据接收与处理技术. 北京: 北京航空航天大学出版社 Feng Z K, Ge X Q, Zhang H Q and Li A. 2016. Remote Sensing Data Receiving and Processing Technology. Beijing: Beihang University Press
[4] 傅文杰, 洪金益. 基于支持向量机的滑坡灾害信息遥感图像提取研究[J]. 水土保持研究, 2006, 13 (4) : 120 –121, 124. Fu W J and Hong J Y. Discussion on application of support vector machine technique in extraction of information on landslide hazard from remote sensing images[J]. Research of Soil and Water Conservation, 2006, 13 (4) : 120 –121, 124. DOI: 10.3969/j.issn.1005-3409.2006.04.038
[5] 郭伟伟, 王秀兰, 冯仲科, 姜磊. 基于NDVI的植被覆盖度变化的研究与分析——以河北省张家口市为例 [J]. 测绘与空间地理信息, 2012, 35 (7) : 63 –66. Guo W W, Wang X L, Feng Z K and Jiang L. The research on vegetation coverage changes based on NDVI: a case study of Zhangjiakou [J]. Geomatics and Spatial Information Technology, 2012, 35 (7) : 63 –66. DOI: 10.3969/j.issn.1672-5867.2012.07.018
[6] 黄润秋. " 5•12”汶川大地震地质灾害的基本特征及其对灾后重建影响的建议[J]. 中国地质教育, 2008, 17 (2) : 21 –24. Huang R Q. Characteristics of geological disasters of 5•12 Wenchuan earthquake and recommendation on its impact on reconstruction[J]. Chinese Geological Education, 2008, 17 (2) : 21 –24. DOI: 10.16244/j.cnki.1006-9372.2008.02.007
[7] 李苗苗, 吴炳方, 颜长珍, 周为峰. 密云水库上游植被覆盖度的遥感估算[J]. 资源科学, 2004, 26 (4) : 153 –159. Li M M, Wu B F, Yan C Z and Zhou W F. Estimation of vegetation fraction in the upper basin of Miyun reservoir by remote sensing[J]. Resources Science, 2004, 26 (4) : 153 –159. DOI: 10.3321/j.issn:1007-7588.2004.04.022
[8] 李松, 邓宝昆, 徐红勤, 王治福. 地震型滑坡灾害遥感快速识别方法研究[J]. 遥感信息, 2015, 30 (4) : 25 –28. Li S, Deng B K, Xu H Q and Wang Z F. Fast interpretation methods of landslides triggered by earthquake using remote sensing imagery[J]. Remote Sensing Information, 2015, 30 (4) : 25 –28. DOI: 10.3969/j.issn.1000-3177.2015.04.005
[9] Li S and Hua H Q. 2009. Automatic recognition of landslides based on change detection//Proceedings of SPIE 7348, International Symposium on Photoelectronic Detection and Imaging 2009: Advances in Imaging Detectors and Applications. Beijing, China: SPIE: 73842E [DOI: 10.1117/12.836109]
[10] 李霞, 徐涵秋, 李晶, 郭燕滨. 基于NDSI和NDISI指数的SPOT-5影像裸土信息提取[J]. 地球信息科学学报, 2016, 18 (1) : 117 –123. Li X, Xu H Q, Li J and Guo Y B. Extraction of bare soil features from SPOT-5 imagery based on NDSI and NDISI[J]. Journal of Geo-Information Science, 2016, 18 (1) : 117 –123. DOI: 10.3724/SP.J.1047.2016.00117
[11] Lu T, Zeng H C, Luo Y, Wang Q, Shi F S, Sun G, Wu Y and Wu N. Monitoring vegetation recovery after China’s May 2008 Wenchuan earthquake using Landsat TM time-series data: a case study in Mao County[J]. Ecological Research, 2012, 27 (5) : 955 –966. DOI: 10.1007/s11284-012-0976-y
[12] Manfré L A, de Albuquerque Nóbrega R A and Quintanilha J A. Evaluation of multiple classifier systems for landslide identification in Landsat Thematic Mapper (TM) images[J]. ISPRS International Journal of Geo-Information, 2016, 5 (9) : 164 . DOI: 10.3390/ijgi5090164
[13] Mondini A C, Chang K T and Yin H Y. Combining multiple change detection indices for mapping landslides triggered by typhoons[J]. Geomorphology, 2011, 134 (3/4) : 440 –451.
[14] 覃志豪, 李文娟, 徐斌, 陈仲新, 刘佳. 陆地卫星TM6波段范围内地表比辐射率的估计[J]. 国土资源遥感, 2004, 16 (3) : 28 –32, 36, 41. Qin Z H, Li W J, Xu B, Chen Z X and Liu J. The estimation of land surface emissivity for Landsat TM6[J]. Remote Sensing for Land and Resources, 2004, 16 (3) : 28 –32, 36, 41. DOI: 10.3969/j.issn.1001-070X.2004.03.007
[15] 苏凤环, 刘洪江, 韩用顺. 汶川地震山地灾害遥感快速提取及其分布特点分析[J]. 遥感学报, 2008, 12 (6) : 956 –963. Su F H, Liu H J and Han Y S. The extraction of mountain hazard induced by Wenchuan earthquake and analysis of its distributing characteristic[J]. Journal of Remote Sensing, 2008, 12 (6) : 956 –963. DOI: 10.11834/jrs.200806128
[16] 孙家抦. 2009. 遥感原理与应用. 武汉: 武汉大学出版社 Sun J B. 2009. Principles and Applications of Remote Sensing. Wuhan: Wuhan University Press
[17] Tang L L, Hu D Y, Li X J and Lian J. 2010. Change detection of landslides and debris in south Taiwan after " Morakot” typhoon based on HJ-1-B Satellite images//Proceedings of 2010 IEEE International Geoscience and Remote Sensing Symposium (IGARSS). Honolulu, HI: IEEE: 3186–3189
[18] 唐怡, 刘良云, 黄文江, 王纪华. 土壤背景对冠层NDVI的影响分析[J]. 遥感技术与应用, 2006, 21 (2) : 142 –148. Tang Y, Liu L Y, Huang W J and Wang J H. Analysis on the influence of soil backgrounds on canopy NDVI[J]. Remote Sensing Technology and Application, 2006, 21 (2) : 142 –148. DOI: 10.3969/j.issn.1004-0323.2006.02.011
[19] 王治华. 2012. 滑坡遥感. 北京: 科学出版社 Wang Z H. 2012. Remote Sensing for Landslide. Beijing: Science Press
[20] 魏永明, 魏显虎, 陈玉. 岷江流域映秀—茂县段地震次生地质灾害分布规律及发展趋势分析[J]. 国土资源遥感, 2014, 26 (4) : 179 –186. Wei Y M, Wei X H and Chen Y. Analysis of distribution regularity and development tendency of earthquake secondary geohazards in Yingxiu-Maoxian section along the Minjiang River[J]. Remote Sensing for Land and Resources, 2014, 26 (4) : 179 –186. DOI: 10.6046/gtzyyg.2014.04.28
[21] 许冲. 基于最大似然法的地震滑坡信息自动提取及其可靠性检验[J]. 中国地质灾害与防治学报, 2013, 24 (3) : 19 –25. Xu C. Automatic extraction of earthquake-triggered landslides based on maximum likelihood method and its validation[J]. The Chinese Journal of Geological Hazard and Control, 2013, 24 (3) : 19 –25.
[22] 许高程, 张文君, 王卫红. 支持向量机技术在遥感影像滑坡体提取中的应用[J]. 安徽农业科学, 2009, 37 (6) : 2781 –2782. Xu G C, Zhang W J and Wang W H. Application of supporting vector machine technology in extraction of remote sensing image of landslide[J]. Journal of Anhui Agricultural Sciences, 2009, 37 (6) : 2781 –2782. DOI: 10.3969/j.issn.0517-6611.2009.06.179
[23] 许积层, 卢涛, 石福孙, 唐斌, 慕楠. 基于NDVI监测5•12震后岷江河谷映秀汶川段滑坡体植被恢复[J]. 植物研究, 2012, 32 (6) : 750 –755. Xu J C, Lu T, Shi F S, Tang B and Mu N. Monitoring the vegetation recovery at landslides along the Minjiang river valley after 5•12 earthquake using NDVI: a case study of the Yingxiu-Wenchuan section[J]. Bulletin of Botanical Research, 2012, 32 (6) : 750 –755.
[24] Yang S W, Li Y K, Feng G S and Zhang L F. A method aimed at automatic landslide extraction based on background values of satellite imagery[J]. International Journal of Remote Sensing, 2014, 35 (6) : 2247 –2266. DOI: 10.1080/01431161.2014.890760
[25] 杨树文, 谢飞, 韩惠, 冯光胜. 基于SPOT5遥感影像的浅层滑坡体自动提取方法[J]. 测绘科学, 2012, 37 (1) : 71 –73, 88. Yang S W, Xie F, Han H and Feng G S. Automatic extraction of shallow landslides based on SPOT-5 remote sensing images[J]. Science of Surveying and Mapping, 2012, 37 (1) : 71 –73, 88. DOI: 10.3969/j.issn.1673-6338.2012.01.017