基于多源遥感时空谱特征融合的滑坡灾害检测方法

陈善静 向朝参 康 青 吴 韬 刘 凯 冯 亮 邓 涛

1(陆军勤务学院 重庆 401311) 2(重庆大学计算机学院 重庆 400044) 3(国防科技大学电子对抗学院 合肥 230037) 4(78092部队 成都 610031)chengshanjing_11@163.com)

我国约70%的地域为山区，地质灾害发生频率高，每年因各种地质灾害造成的损失数以亿计，给国家和人民带来巨大损失.我国西南地区多山地，地质、地貌稳定性差，在地震或暴雨等诱发因素作用下滑坡等地质灾害可能呈大面积高发状态.在大面积滑坡灾害发生时快速获取滑坡区域分布、数量、规模等灾情信息对救援决策和防灾减灾都有着重要意义.当前针对滑坡的检测识别与灾情信息提取多以现地勘测为主，存在耗时长、工作量大、效费比低，提取灾情信息不直观、不精确、不全面，受人为因素和地理环境影响大，难以对灾害地区实现快速大面积精确检测与综合评估等问题.

在航空遥感滑坡检测识别方面，文献[8]提出了一种基于航空遥感图像的半主动滑坡危险区域制图方法，该方法利用变化向量分析和阈值分割提取候选滑坡危险区域，以形态学方法消除误差，最后以边缘分级进化和区域分级进化对滑坡目标进行快速识别.该方法主要适用于高分辨率航空正射影像，而对多类型、多谱段、低分辨率遥感图像其适用效果不够理想.文献[9]以2013年芦山地震震区发生的滑坡为研究对象，通过构建区域地质特色的地震滑坡样本库，引入迁移学习机制实现基于无人机高分辨率遥感影像地震滑坡信息提取.该方法对地震滑坡识别精度较高，满足地震滑坡灾害环境宏观调查、滑坡灾害体监测等应急需求，但是在训练样本完整性、全面性和特征准确提取方面还有待进一步完善和优化.文献[10]以机载激光雷达(light detection and ranging, LiDAR)提取的DEM数据和QuickBird卫星图像为基础，通过小波变换对2类数据进行融合，并根据目标光谱、纹理和局部空间位置特征对其进行面向对象的目标分类，最终实现了对热带城市地区滑坡的精确检测，但该方法对LiDAR遥感数据的空间分辨率要求较高，而低分辨率遥感图像融合处理效果不够理想.

本文从大面积滑坡发生区域灾情信息快速获取的客观需求出发，以灾前灾后多波段光学遥感图像为基础，在综合分析滑坡区域灾前灾后时空谱特征变化规律的基础上，提出了一种基于多源遥感时空谱特征融合的滑坡灾害检测方法，通过对滑坡发生前后遥感图像进行光谱空间和尺度空间配准、多特征融合遥感影像数据构建、基于SVM(support vector machine)的地物目标识别、典型形状特征分类，实现对滑坡目标的精确识别.

1 基于多维特征融合的遥感影像集构建

1.1 多源遥感影像光谱空间配准

Fig. 1 Results of spectral space registration for remote sensing images(remote sensing images of Sentinel-2)图1 遥感图像进行光谱空间配准后结果(Sentinel-2卫星遥感影像)

在空天遥感成像中各种光学成像设备的光谱辐射分辨率各不相同，获取的遥感图像定标和校正方法也存在一定差异，这就造成了多源光学遥感图像在多个谱段上的辐射亮度分布呈现较大差异.即使同一卫星对相同场景多次成像其获取的光学图像也可能出现较为明显的光谱特征变化.以Sentinel-2卫星遥感成像为例，其对同一区域滑坡前后遥感成像获得的多光谱图像就存在明显的色差，如图1(a)(b)所示.滑坡后遥感图像各地物颜色与真实目标较为接近，但滑坡前的遥感图像整体偏绿.同一地物目标在灾前灾后图像的颜色或光谱特征差异将直接影响后期数据融合处理和目标检测应用的准确性和可靠性，因此，本文在深入分析滑坡灾前灾后遥感图像光谱特征的基础上提出对同一场景多时相多源遥感图像进行光谱空间分布配准.对于第p波段的图像，其直方图分布概率函数可表示为

(1)

其中，sk为图像第k级灰度值，0≤sk≤1，k=0,1,…,L-1，np,k为第p波段图像中灰度值为k的像素的个数，L为图像灰度级，N为该波段图像中像素的总数.直方图提供了特定谱段遥感图像中各种灰度分布情况，是一幅图像所有灰度或亮度值的整体描述[11].将直方图分布转换为积累分布函数：

(2)

Fig. 2 The histogram distribution of unregistrated/ registrated image for landslide area图2 配准前后滑坡区域的遥感图像直方图分布情况

本文在多源遥感影像光谱空间配准中以灾后滑坡区域各波段遥感图像直方图分布函数为基准，对灾前遥感图像进行对应的直方图相似变换.根据单映射规则从小到大依次找到能使式(3)最小的k和l，将灾前图像直方图概率分布Qp,before(sk)对应到灾后图像直方图概率分布Qp,after(sl).

minψ(k,l)=|tp,before(k)-tp,after(l)|,

(3)

其中，tp,before(k)和tp,after(l)为滑坡区域灾前和灾后p波段图像的直方图积累分布函数，k∈[0,Lbefore-1]，l∈[0,Lafter-1].Lbefore和Lafter为灾前和灾后p波段图像的灰度级数.对灾前遥感图像进行变换的结果如图1(c)所示，变换前灾前灾后遥感图像中红色分量的直方图分布如图2(a)(b)所示，通过对比可以发现，灾前图像中红色分量像素灰度级整体偏低，与灾后图像存在红色分量的直方图分布差异明显，造成2幅图像出现显著色差.通过直方图相似变换后，灾前图像的红色分量直方图分布整体与灾后图像红色分量直方图分布趋于相似，如同图2(c)所示，对比图1(b)(c)可以发现，利用直方图变换后获得的滑坡发生前遥感图像与滑坡后遥感图像整体色差基本消除，相同地物目标在不同遥感图像中颜色或光谱特征基本趋于一致.

1.2 多源遥感影像尺度空间配准

多源遥感平台由于成像高度、角度和探测器空间分辨率等因素各不相同，其获得的图像空间尺度上也存在一定差异，需要开展灾前灾后多源遥感影像尺度空间配准.SIFT(scale invariant feature transform)算法在特征点提取方面具有缩放不变性、旋转不变性、仿射不变性，还有一定的抗光照变化和抗视点变换性能[12-13]，因此本文利用SIFT算法对灾前灾后滑坡区域遥感图像进行尺度空间配准，其基本步骤包括：尺度空间极值点提取、提取特征点量化表征和特征点匹配.

对于二维图像I(x,y)的高斯尺度空间定义为L(x,y,t)，该尺度空间可由参数为t的高斯核G(x,y,t)与I(x,y)卷积得到：

(4)

其中

(5)

其中，x,y,t分别为L(x,y,t)的像素位置坐标和尺度参数.尺度参数t为尺度空间因子，它是高斯正态分布的方差，表征图像被平滑的程度.为高效地在尺度空间内检测出稳定的特征点，使用尺度空间中差分高斯(difference of Gaussian, DOG)极值作为判断依据.DOG算子定义为：

D(x,y,t)=(G(x,y,kt)-G(x,y,t))*
I(x,y)=L(x,y,kt)-L(x,y,t).

(6)

利用多级金字塔模型检测DOG局部极值.将每个像素跟同一尺度的周围邻域8个像素和相邻尺度对应位置的9×2个像素总共26个像素点进行比较，寻找极值点并作为特征点保存.根据邻域像素梯度分布特性，以梯度模m(x,y)和梯度方向θ(x,y)对高斯图像中每个点L(x,y)进行量化表征.

m(x,y)={[L(x+1,y)-L(x-1,y)]2+
[L(x,y+1)-L(x,y-1)]2},

(7)

(8)

利用SIFT算法分别提取灾前灾后遥感影像对应特征点的m(x,y)与θ(x,y)，利用二维范数寻找最优匹配特征点，实现滑坡灾前灾后遥感影像的特征点匹配，图1(b)(c)的特征点匹配结果如图3所示：

Fig. 3 Registration result of characteristic points for remote sensing images before and after the landslide occurrence图3 滑坡灾前灾后遥感影像特征点匹配结果

1.3 基于时空谱特征融合的遥感影像数据集构建

以经过光谱空间和尺度空间配准的灾前灾后遥感图像构建新的多维特征融合影像数据集.将滑坡区域灾前灾后时序变化特征转化为多光谱图像的光谱特征，为后期时空谱特征融合应用奠定基础.以Sentinel-2遥感图像为例，在综合分析滑坡地区背景地物光学特征、遥感成像有效波段和图像分辨率的基础上，选择能覆盖滑坡颜色和光谱特征的B(440～538 nm),G(537～582 nm),R(646～684 nm)和NR(760～908 nm)这4个典型特征波段灾前灾后图像构建新的影像数据，如图4所示.新构建的影像数据中，前4个波段为目标区域灾前遥感图像，后4个波段为目标区域灾后遥感图像.由于经过遥感光谱空间和尺度空间配准，图像中相同位置地物目标光谱信息和灾前灾后图像变化信息均整合为新构建的8波段影像数据的光谱特征信息.在新构建的数据中部分地物目标的光谱分布曲线如图5所示.从图5中可以看出，植物和建筑物目标灾前灾后图像中未发生明显变化，因此灾前和灾后B,G,R,NR这4个波段的光谱曲线空间分布基本相似，而新翻土地在灾前灾后图像中发生了变化，灾前和灾后4个波段的光谱曲线空间分布出现差异，这种差异主要为目标区域的变化信息，通过对该信息的有效利用即可实现对特定地物目标的变化检测和识别.

Fig. 4 The data of new eight bands images图4 新构建的8波段遥感影像数据集

Fig. 5 The object spectrum of new eight bands images set图5 新构建的8波段影像数据中部分地物光谱分布曲线

2 基于时空谱特征融合的滑坡检测方法

2.1 基于SVM的滑坡地物目标检测

SVM是一种基于VC(vapnick-chervonenkis)维理论和结构风险最小化原理的经典机器学习分类方法，在图像处理和目标识别等方面都有着重要应用.本文将其用于新构建8波段遥感图像中滑坡目标检测.SVM求解最优超平面问题等价于求解如下方程[14]：

(9)

(10)

其中K(·)是满足Mercer条件的核函数.

2.2 滑坡典型空间形状特征分选

发生滑坡区域山体表面的植被会受到不同程度的破坏，新裸露出的土石在亮度、色调、空间形状方面与周围植被和人工地物有明显差异.本文选取滑坡中较为常见的2种基本形状——簸箕形和三角形，构建基础形状模型，如图6所示：

Fig. 6 Fundamental shape models of landslides图6 滑坡常见基础形状模型

通过轴向长宽比、面积参数和不变矩特征对遥感图像滑坡区域进行分选，实现对滑坡目标的精确识别.

1) 轴向长宽比

通过对滑坡基础形状模型的分析，建立表征滑坡形状特征的轴向长宽比(如图7所示)：

(11)

Fig. 7 The axis length ratio of fundamental shape models of landslides图7 滑坡常见基础形状模型的轴向长宽比

2) 面积参数

当大面积区域内发生滑坡灾害时，小型滑坡数量较多，但通常造成危害不大；而特大型滑坡或者巨型滑坡数量通常较少，前期可能已有一定的监测数据.这2类滑坡在遥感检测与识别中可适度弱化，而对于中、大型滑坡其危害性和未知性均较大，将成为遥感检测与识别的重点.因此本文通过设置合理的滑坡区域面积参数Rsize，实现对典型滑坡的分选与识别.

3) 不变矩特征

考虑到滑坡在遥感图像中可能存在旋转和尺度的不确定性，本文利用不变矩特征作为检测依据，对滑坡遥感检测结果进行特征匹配识别.对于图像I(x,y)的阶几何矩(标准矩)定义为[11]：

(12)

其中，p=0,1,…，W;q=0,1,…，H;W和H为图像的宽高.p+q阶中心矩定义为

(13)

(14)

利用二阶和三阶归一化中心矩构造7个不变矩特征参数M1～M7如下：

(15)

提取滑坡检测结果中的符合面积参数要求的连通区域作为疑似滑坡区域，分别计算基础形状模型与各疑似滑坡区域的不变矩特征参数，以不变矩的绝对差值作为疑似滑坡区域的选择标准.

(16)

其中，Mi-possible为疑似滑坡区域不变矩特征参数M1～M7，i∈{1,2,…,7}，Mi-standard为基础形状模型不变矩特征参数M1～M7，i∈{1,2,…,7}.

2.3 多源遥感图像时空谱特征融合检测流程

本文提出的基于多源遥感时空谱特征融合的滑坡灾害检测方法以滑坡区域灾前灾后多光谱遥感图像为基础，经过多源遥感影像图像光谱空间配准和尺度空间配准后，构建包含地物光谱和时变特征信息的多特征融合遥感影像数据集，在此基础上利用SVM对滑坡区域进行检测识别，获得疑似滑坡区域，再结合滑坡空间形状基础特征模型，依据轴向长宽比、面积参数和不变矩等特征指标对滑坡区域进行精确分选识别，实现对大面积范围内的滑坡灾害快速检测和精确识别.本文提出方法的基本流程如图8所示:

Fig. 8 Basic flow chart of the proposed method in our paper图8 本文滑坡遥感检测方法流程图

Fig. 9 Field investigation and information collection by UAV图9 现地勘察和无人机信息采集

3 实验结果与分析

为验证本文提出方法的有效性，利用Sentinel-2获得的重庆市沙坪坝大学城某地区的滑坡前后多光谱遥感图像(B,G,R,NR这4个波段)进行了滑坡遥感检测实验.实验将目标检测和识别中常用的大津法(Otsu)、Bayes分类器、传统的SVM、传统的SVDD(support vector data description)、神经网络(neural networks, NN)等方法[15-17]与本文方法的检测识别效果进行了对比.算法运行平台：Intel Core i5-4210 2.4 GHz CPU， 4 GB内存，MATLAB2010.通过现地勘察和无人机信息采集后确定训练样本和检验样本区域(如图9和图10所示)，在遥感图上标记训练样本分布如图10所示，训练用滑坡样本25个、训练用背景样本112个，其中，植被样本32个、建筑物样本32个、水体样本16个、训练场样本32个.检测效果量化评估中所用到的检测样本如图11所示，其中，滑坡样本137个、背景地物样本514个.各对比方法参数设置和使用的波段如表1所示.各种检测方法对滑坡及背景地物的识别情况如图12所示.对滑坡目标识别精度如表2所示.

Fig.10 Distribution of training samples for each detection method图10 各检测方法在所用到训练样本分布

Fig. 11 Distribution of test samples图11 检测样本分布情况

Table 1 Used Bands and Parameter of Comparative Methods表1 各对比方法使用的波段和参数设置

Fig. 12 The detection result of comparative methods on Sentinel-2 remote sensing images图12 Sentinel-2遥感图像中各对比方法滑坡识别结果

Table 2 Accuracy of Each Method for Landslide Detection表2 各种方法滑坡目标识别精度表 %

结合图12和表2可以看出，传统的Otsu方法主要根据地物目标的亮度信息对目标进行分割和识别，难以区分建筑物或道路等亮度较高的人工地物目标，整体识别精度较差，虚警和漏警的像素均较多.Bayes分类器对滑坡识别精度较高，但虚警率也较高，检测结果图出现了大量的噪声像素点.利用灾后的R,G,B,NR这4个波段进行基于SVM的滑坡检测，其正确识别率达67.88%，但漏警率较高，并且与滑坡颜色相近的许多训练场也存在识别错误.SVDD方法整体识别效果较好，但虚警率较高，结合图12(d)可以看出，该方法对建筑和道路等目标仍难以区分.NN方法能区分大部分建筑和道路目标，但对滑坡的正确识别率和漏警率均不理想.利用灾前灾后8波段图像进行基于SVM的滑坡检测方法整体识别效果较好，结合图12(f)可以看出，部分新建工程开挖的地面由于其光谱特征和变化特征与滑坡较为相似，因此被错误地识别成了滑坡目标，所以其误警率达10.70%.本文方法在基于SVM的检测方法基础上加入滑坡基础形状特征，实现了对检测结果的精确分选和识别，最终的识别效果如图12(g)所示，结合表2可知，本文方法有效去除了误分类点，滑坡的正确识别率达到95%以上，同时误警率下降到了4.47%.通过多种滑坡检测和识别方法的对比可以看出，本文提出的方法整体识别精确优于其他多种方法.

4 检测方法讨论

本文提出的滑坡灾害检测与识别方法主要利用了遥感图像中目标的光谱特征、变化特征和空间形状特征.具体实现过程如下：

1) 在经光谱-尺度空间配准的灾前和灾后图像基础上构建多维特征融合的遥感影像数据集，将灾前灾后遥感图像中光谱信息和变化信息融合后表征为新建遥感影像数据集的各种光谱曲线.

2) 利用SVM对各种模式的光谱曲线进行识别，提取疑似滑坡区域.此时SVM识别滑坡目标的主要依据为目标的光谱信息和变化信息，结合表2 中SVM(8 Bands)和SVM(4 Bands)两栏数据可以看出，融合变换信息和光谱信息后，滑坡正确识别率提升了27.01%.

3) 结合滑坡典型空间形状特征对SVM的识别结果进行分选和优化.该环节将滑坡的空间形状信息融合到目标识别算法中，最终实现光谱特征、变化特征和形状特征的融合检测应用.结合表2中实验数据可知，本文最终方法的滑坡正确识别率得到进一步提升，同时误警率较未加形状特征的SVM(8 Bands)下降了6.23%.

在滑坡灾害遥感探测和灾情提取方面，灾前灾后遥感图像对比分析提取常常是比较有效的方法，但受卫星重返周期和气象条件的影响，近实时的灾后图像和滑坡前短期内灾前图像往往不容易获取，灾前图像和灾后图像较大的时间差可能引入较多的变化信息，影响本文方法对滑坡识别的精度.

5 结论与展望

本文从遥感图像时序变化特征、光谱特征和空间形状特征融合利用角度出发，提出一种基于多源遥感时空谱特征融合的滑坡灾害检测方法，在多源遥感图像光谱空间配准和尺度空间配准的基础上构建基于多维特征融合的遥感影像数据集，实现目标区域时序变化和光谱特征的有效融合，利用SVM对滑坡进行检测与识别，结合空间形状特征对检测结果进行分选和优化，实现光学遥感影像时、空、谱多维特征信息的滑坡检测应用，为滑坡灾害的遥感检测和识别提供了一种新方法.

下一步工作中我们将研究航空、航天多源遥感信息融合处理和变化特征提取方法，通过多源遥感信息融合和互补提升变化特征提取的有效性和准确性，进而实现对滑坡目标的快速精确探测.