高分辨率遥感影像滑坡纹理特征提取方法

王庆雅，张 锦

(太原理工大学 矿业工程学院，太原 030006)

滑坡是我国发生数量最多、分布最广的地质灾害，每年因滑坡而造成直接或间接的经济损失不可估量。因此，在滑坡发生后，利用高分辨率遥感影像进行滑坡快速提取对应急救援、灾损评估具有重要的意义。

由于地形地貌、地质条件等因素的不同，滑坡体在遥感影像上的特征有较大差别[1]，因此，很难建立起通用的滑坡识别特征体系。目前，滑坡自动提取主要依靠基于高分辨率遥感影像结合多特征体系的滑坡自动识别方法。MARTHA et al[2]利用光谱、纹理、形态等特征实现了对不同类型滑坡的半自动提取。KURTZ et al[3]利用光谱、纹理、几何、上下文等特征采用自上而下的基于区域的层次框架实现了滑坡的分层提取。刘辰等[4]使用不同分辨率的DEM进行面向对象的滑坡提取，指出空间分辨率越高的DEM数据，越能准确地识别出滑坡。彭令等[5]结合光谱、纹理、几何、地形等特征属性建立规则集进行区域滑坡信息提取。充分利用高分影像的纹理、几何和上下文等空间特征及光谱特征，结合地形特征，建立多特征滑坡识别体系，虽然是提高滑坡识别与提取能力的有效途径，但其中存在特征冗余、滑坡识别规则难以统一等问题，而且在目前滑坡提取中用到的纹理特征主要是灰度共生矩阵这种统计特征，对于单个滑坡在遥感影像上具体的纹理表现方式及规律性研究还不够充分。

纹理反映了物体表面颜色和灰度的某种变化规律，是区分具有相似光谱、不同空间分布结构特性地物的有效信息，被广泛应用在影像信息的提取中[6]。宋荣杰等[7]分别用灰度共生矩阵、分形和空间自相关3种纹理特征方法进行苹果园提取并对比了识别精度，表明灰度共生矩阵的效果较好。闫利等[8]利用一种新的结构指数特征表达种植园的纹理结构，实现依靠单一或者少数特征的种植园自动提取。这些方法要求地物具有较为明显规则的纹理结构，对于没有明显规律的地物适应性较弱。徐望明等[9]利用多尺度多方向Gabor变换提取GIST特征实现人脸识别。游永发等[10]利用多尺度多方向Gabor小波变换结果提取建筑物特征点并实现了建筑物的分级提取。但是这些方法都存在特征冗余、计算过程较为繁琐等问题。

上述纹理特征提取方法主要是针对有规则的人工地物具有较好的效果，而滑坡的纹理结构往往不规则，其规律性不够明显，因此上述方法对于滑坡的纹理提取效果不明显。本文探索了单体滑坡的纹理特征，提出了一种多尺度多方向的滑坡纹理表达方法，能够有效减少特征冗余，简化提取过程，实现遥感影像上滑坡多尺度多方向的纹理提取。

1 研究区及数据

黑方台位于甘肃省永靖县，地处湟水河与黄河交汇口上游，台塬分布广泛，其独特的地形地貌地质环境及人类耕地灌溉活动，使得区域滑坡灾害频发。实验选用从谷歌地球上获取的2015年9月14日黑方台党川段滑坡数据，空间分辨率为0.29 m，面积为1.05 km2，如图1所示。1#滑坡和2#滑坡均由北向南滑动，分布在台塬边缘，滑动距离较远，滑坡区域较亮与周围的农田有较大的光谱差异；形态上看呈两头大中间长的哑铃状，滑坡后缘呈圈椅状；滑坡体上的溯流状明显，有明显的波纹状起伏，如图2所示。

图1 实验区遥感影像Fig.1 Remote sensing image of the experimental area

图2 滑坡体上的波浪状起伏[11]Fig.2 Wavy ups and downs on the landslide

2 滑坡纹理特征分析

2.1 滑坡区域不同地物影像纹理特征

纹理作为一种重要的视觉表达，广泛存在于物体的表面，在遥感影像上则表现为纹理特征值反复出现的局部模式和它们的排列规则[12]。周期性和方向性作为描述纹理图像的视觉感知，所以利用周期和方向进行纹理分析是十分重要的。

研究区滑坡特殊的运动现象，使得其在遥感影像上具有明显的方向性纹理，形成与周围地物有明显不同的纹理特征，如图3所示。从几种典型地物的遥感影像来看，滑坡体上具有明显的垂直于滑动方向的波纹状纹理，人工建筑区多呈排列规则的网格状纹理，耕地多呈规则块状并具有沿一定方向排列的颗粒状纹理，斜坡的纹理方向杂乱不规则。

(a)泥流型黄土滑坡；(b)人工建筑区；(c)耕地；(d)斜坡图3 不同地物遥感影像Fig.3 Remote sensing images of different features

2.2 滑坡区域不同地物频谱能量表征

为对上述不同地物的纹理方向进行定量分析，利用二维傅里叶变换将影像从空域转换到频域，计算不同地物的频谱能量并用文献[13]提出的楔特征(角向分布特征)计算方法，在[0°，180°]方向上以1°为采样间隔，对从原点出发的扇面求能量和，分析不同地物的纹理主方向。

傅里叶变换能量谱的自配准性质，使得同一方向不同分布地物特征线的能量叠加在一起，其所对应的谱线经过频谱中心并垂直于原地物的特征线方向[14]。即若原地物的纹理主方向为θ，则该地物影像对应的角向分布曲线的峰值出现在θ+π/2上；若原地物没有明显的方向特点，则其对应的角向分布曲线较为平稳，没有明显的峰值。图3中4种地物相对应的角向分布曲线如图4所示。

图4 不同地物角向分布曲线Fig.4 Angular distribution curves of different features

分析上述不同地物的角向分布曲线可以得出，滑坡的能量集中在91°的方向上，纹理主方向为1°；人工建筑区的能量集中在42°和133°两个方向上，纹理特征方向为132°和43°；耕地的能量集中在104°的方向上，纹理特征方向为14°；斜坡的能量分布较为分散，没有明显的方向性纹理。由此可知，一定区域内，滑坡与周围地物具有不同的方向性纹理。

3 研究方法

考虑到在一定区域内滑坡与周围地物这种不同的方向性纹理表现，并利用此种纹理差异进行滑坡的提取，本文设计多尺度多方向Gabor滤波器组进行滑坡纹理特征的提取，首先对滑坡图像进行多方向多尺度的Gabor变换，得到多方向多尺度的Gabor特征；其次根据熵值对同一方向不同尺度的特征图进行融合，得到多方向Gabor特征；最后对多方向Gabor特征进行最大值过滤构成多尺度多方向Gabor纹理特征。为了分析此方法对滑坡提取的有效性，分别计算了多尺度多方向Gabor纹理特征、局部二值模式(LBP)和局部空间自相关三种纹理特征，并将纹理特征与光谱特征融合，利用支持向量机(SVM)提取滑坡，根据混淆矩阵分析不同方法滑坡的提取精度。

3.1 多尺度多方向Gabor纹理特征

3.1.1Gabor滤波

Gabor变换是一种基于信号处理的纹理描述方法，自适应Gabor滤波器能够从不同尺度和方向充分对纹理特征进行表征。在空域上，一个二维Gabor函数实质上是一个二维高斯函数与二维复数正弦函数的叠加，表达如下：

(1)

其中

xr=xcosθ+ysinθ.

(2)

yr=-xsinθ+ycosθ.

(3)

式中：σx，σy分别代表高斯函数在x、y方向上的标准差，决定滤波器的形状；θ为高斯函数的旋转方向，代表滤波器的方向；f为滤波器的中心频率，决定空间尺度因子的选择，中心频率越小空间尺度越大；(x,y)为图像像素坐标。

从Gabor函数可以看出，Gabor滤波器为复数滤波器，实部为偶对称滤波器，虚部为奇对称滤波器，文中考虑复数响应，方向取值在[0，2π]范围内8等分：

(4)

中心频率依据JAIN et al[15]提出的滤波器同一方向上不同中心频率的确定方法，取不同尺度为：

(5)

式中：m为中心频率的个数。文中所用数据大小为1 024×1 024，m为8.但是实验结果发现当尺度过大时，滤波结果不能有效描述图像纹理变化，如图5所示。在固定90°的方向上，当m≥4时，滤波图像上强弱响应变化太大，无法解释纹理变化，因此文中m=0，1，2，3.

图5 同一方向不同尺度滤波图像Fig.5 Filtered images of different scales in the same direction

故本实验中有8个方向(u=0,1,2…7)，4个尺度(v=0,1,2,3)共32个Gabor滤波器，用这32个滤波器分别对原始图像进行卷积运算。

Fu,v(x,y)=I(x,y)*Gu,v(x,y).

(6)

式中：*为卷积操作符；I(x,y)为原始滑坡图像；Gu,v(x,y)为Gabor滤波器；Fu,v(x,y)为不同方向和尺度的Gabor特征图。

文中使用复数Gabor核函数，卷积结果也会产生实部和虚部的复数响应，相关研究发现，滤波器的实部对纹理比较敏感，大多数研究在计算中只取实部作为纹理计算结果；本文研究发现，在滤波后的滑坡影像中，综合实部与虚部的幅值图能够突出影像局部的能量特征，并能够较好地体现滑坡与周围地物的差异。因此，本文根据滤波后的幅值图反映滑坡影像的纹理信息。

Mu,v(x,y)=Amplitude[Fu,v(x,y)]=Ru,v2+Iu,v2.

(7)

式中：Ru,v和Iu,v分别为实部和虚部图像，Mu,v(x,y)为幅值图像。

图6是u=2，v=1时的滑坡影像经滤波处理后的结果。

从图6可以看出，滑坡与非滑坡的滤波响应是不同的，滑坡区的滤波响应整体上小于非滑坡区，这是因为，在非滑坡区有大量排列规则的人工建筑及耕地，纹理结构简单，在单方向上会有较大的响应，而滑坡的能量虽然集中在大约90°的方向上，但是滑坡的纹理结构较为复杂，拥有更多方向的纹理信息，局部纹理也更加丰富。Gabor滤波后的实部和虚部图像都能够在一定程度上反映滑坡的特征信息，但与综合实部和虚部的幅值图像相比，幅值图像中的信息量较多，能显著减少影像信息的缺失。

图6 滑坡影像与Gabor滤波结果Fig.6 Landslide image and Gabor filtering result

8个方向4个尺度的滤波结果如图7所示。

在图7中，每行方向一致，尺度从左至右为0，1，2，3，每列的尺度固定，方向从下向上分别为0，1，2，3，4，5，6，7.从图中可以看出，在方向固定的条件下，小尺度下图像的边缘更细腻，大尺度下图像的轮廓更清晰，这是因为尺度越大频率越小，高频凸显图像的边缘特征，低频凸显图像的轮廓特征。在尺度固定的情况下，不同方向的滤波器允许通过地物的方向是不一样的，这在规则排列的人工建筑区有很好的体现，而滑坡的方向特征相对其它地物不明显，在滤波结果图上具有较小的响应。

图7 不同尺度和方向的滤波结果Fig.7 Filter results of different scales and directions

3.1.2多尺度Gabor特征图融合

如上图7所示，同一方向上的单一尺度往往不足以表征地物的纹理信息，为充分表征地物的轮廓和边缘信息，根据熵值将同一方向不同尺度的特征图进行融合。由信息论可知，熵是对不确定性的测量，信息熵值越大其包含的信息也就越多。本文根据熵值对同一方向不同尺度的特征图进行融合，熵值越大对应特征图被赋予的权重也就越多，得到多方向Gabor特征Su(x,y)(u=0,1,2…7)，计算方式如下:

(8)

其中权重计算方法为：

(9)

式中：Eu,v为滤波后图像的信息熵；Pu,v为对应特征图的权重。

如图8所示，融合后的图像能够在保留全局轮廓的基础上有效增加图像的边缘信息。

图8 同一方向不同尺度融合结果Fig.8 Fusion results of different scales in the same directio

3.1.3多方向最大值滤波

如图7所示，在固定像素处Gabor纹理特征在不同方向上响应不同，具有八维特征，存在一定程度的信息冗余，将8个不同方向特征图中各像素的最大响应值作为Gabor纹理特征，构成多尺度多方向Gabor纹理特征G(x,y).

G(x,y)=max(Su(x,y)) .

(10)

如图9所示，经过最大值过滤后，具有多种方向的地物在对应方向上的纹理特征都得到了增强，如人工建筑区已经呈现出明显的建筑轮廓，滑坡体虽然有一主方向，但这种方向相对于图像并不明显，因此滑坡的纹理值在图像上是低于其它地物的。并且，经过最大值过滤之后的图像，相同地物具有相似的纹理特征表现出同类的聚集现象；不同方向的地物纹理特征值相差较大，表现为不同类的分区现象。

图9 多尺度多方向Gabor纹理特征Fig.9 Multiscale and multidirectional Gabor texture features

3.2 LBP纹理特征

LBP(local binary pattern，LBP)是根据中心像素及其邻域系统的联合分布来进行量化以获得纹理特征[12]，具有灰度不变性。在窗口内，将窗口中心像素与邻域像素的差值与0进行比较，大于等于0量化为1；否则，量化为0，最后将8位二进制数转化为十进制构成图像LBP纹理特征。本文用3×3窗口计算图像的LBP特征，表达式如下：

(11)

其中

(12)

式中：hi为中心像素值，hn为邻域像素值。

3.3 局部空间自相关

空间自相关是指空间上邻近位置属性信息的相似性，局部空间统计量可以度量局部的空间相关性，揭示局部空间关联模式，能够较好地理解和度量区域内部空间关联的结构模式与变异性[16]。本文用基于距离的Getis指数来表示影像的局部空间自相关，其表达如下：

(13)

其中

(14)

3.4 SVM图像分类

SVM(support vector machine，SVM)是建立在统计学习理论中的VC维理论和结构风险最小化准则基础上的一种监督分类算法，在核函数确定的特征空间上寻找一个能够正确划分训练数据集，并且几何间隔最大的分离超平面，即求解最优分类决策函数：

(15)

核函数类型及参数设置会影响SVM分类结果精度，常见核函数有线性函数、多项式函数、径向基函数和Sigmoid函数。本文在ENVI5.3平台上选择常用的径向基函数作为核函数，试错法对比分析获得最优核参数g和惩罚函数C.

4 数据结果处理与分析

4.1 纹理特征提取分析

图10 不同方法纹理特征提取结果Fig.10 Texture feature extraction results by different methods

4.2 滑坡提取结果分析

为验证本文提出的纹理计算方法对滑坡提取的有效性，并对不同纹理特征方法进行定量比较，本文设计了如表1所示的4种滑坡提取方案，并通过滑坡提取结果来分析3种纹理计算方法的优劣性。

表1 不同滑坡提取方案Table 1 Different landslide extraction methods

各方案经SVM分类后的结果有大量的小图斑，对于小图斑的处理，分别用以下步骤对分类后的图像进行处理：首先利用5×5的变换核函数对待分类滑坡图像进行Majority分析，将大量小图斑归到邻域内占主导地位的类别中；其次用形态学闭运算进一步填补滑坡内部空洞，使得图像更加平滑；最后对滑坡图像中仍存在的较大独立图斑进行过滤处理，过滤阈值设为150，并将过滤掉的图斑归为背景类。

经分类后处理得到的滑坡提取结果如图11(b)-(d)所示。目视解译结果如图11(a)所示。从图中可以看出，方案Ⅰ提取结果中错分的滑坡较少，但是有漏分现象，其中1#滑坡的滑源区缺失严重；方案Ⅱ和方案Ⅲ提取结果中错分的滑坡较多，都包含了大量的斜坡，方案Ⅲ提取结果中错分滑坡最多。

图11 目视解译与不同方案滑坡提取结果Fig.11 Visual interpretation and different methods of landslide extraction results

混淆矩阵是目前分析分类结果精度最常用的方法，为对提取结果进行定量评价，通过目视解译结果与不同方案的提取结果计算混淆矩阵的总体分类精度(overall accuracy，OA)、Kappa系数、错分误差(commission error，CE)和漏分误差(omission error，OE)等指标来评价提取结果精度和可靠性。3种方案的提取精度如表2所示。

对比3种方案提取结果并结合表2可知，多尺度多方向Gabor纹理特征结合光谱特征的提取效果最好，总体分类精度达到了95.89%，Kappa系数为0.87；其次是光谱结合LBP纹理特征，总体分类精度为90.1%，Kappa系数为0.73；光谱结合空间自相关的总体分类精度最低为89.13%，Kappa系数为0.70.

表2 不同方案滑坡提取精度Table 2 Extraction accuracy of landslides in different methods

同时，多尺度多方向Gabor纹理特征提取结果的错分误差最小，仅为8.10%；而经LBP和局部空间自相关纹理特征提取的滑坡，错分误差均达到了30%以上，尤其局部空间自相关提取结果中的错分误差最大，为33.22%.从提取结果来看，错分误差较大是由于错提了大量的斜坡而引起的，表明多尺度多方向Gabor纹理特征能够较好地对滑坡和斜坡做出区分，而LBP和局部空间自相关纹理特征对滑坡和斜坡的区分能力较差。

从漏分误差来看，多尺度多方向Gabor纹理特征遗漏提取的滑坡像元最多，漏分误差最大为12.57%，而LBP和局部空间自相关纹理特征的漏分误差都在8%以下，分析提取结果可知，这是由于1#滑坡的滑源区与整体滑坡方向不一致，没有较好地提取出滑源区。整体提取结果表明，本文提出的多尺度多方向Gabor纹理特征对具有不同方向纹理特征的地物具有较好的可区分性，在能够较好地提取出滑坡的基础上有效减少斜坡基岩对滑坡提取的影响，提高滑坡提取精度。

5 结论

本文利用谷歌地球影像和SVM图像分类技术，综合运用光谱信息和纹理信息，在对泥流型黄土滑坡纹理特征分析的基础上，设计了一种能够有效描述泥流型黄土滑坡纹理结构的特征，多尺度多方向Gabor特征，与LBP和局部空间自相关纹理提取方法进行了对比，结果表明：

1) 与原始多维Gabor特征相比，此特征经过融合与过滤，能有效降低特征维数，在保留原有重要信息的基础上有效减少冗余特征。

2) 通过对比可知，此特征相比较LBP、局部空间自相关特征对滑坡具有较好的提取结果，提取精度达到了95%以上。

3) 本文不足之处在于实验中相关参数的设置包括LBP纹理窗口、局部空间自相关距离设置均为经验性值，这些与影像分辨率、滑坡尺度大小有关，缺乏有效的依据，需要进一步研究，同时本文提出的纹理特征方法是否能进一步将单体滑坡划分为滑源区、滑移区和堆积区是接下来需要继续研究的内容。