赵志鹏,张 锦
(太原理工大学 矿业工程学院,山西 太原 030024)
我国是发生滑坡灾害较多的国家之一,滑坡的发生常造成严重的人员伤亡和财产损失[1]。随着空间遥感技术的发展,遥感图像已经成为人们获取滑坡信息的重要数据源。遥感影像滑坡的轮廓特征具有直观性强,且不受光谱、纹理、地形和种类不同等因素的影响,因此准确快速地获取滑坡轮廓特征,对应急救援和灾后重建等工作具有重要意义。
国内外研究者对滑坡灾害进行了多方面的研究,在滑坡特征提取方面,文献[2]对滑坡影像利用Gabor变换来提取滑坡的纹理特征。文献[3]从滑坡目标和背景颜色差异建立特征模型来提取滑坡。文献[4-5]利用滑坡影像多种特征实现了滑坡自动识别提取。文献[6-7]根据滑坡特征建立规则集实现滑坡的分层提取。显然研究者对滑坡的轮廓特征大都是结合光谱和纹理等特征来共同识别提取滑坡,对于滑坡单一的轮廓特征在遥感影像上的表现形式和规律性研究不是很充分。
目标轮廓是图像中用于识别的非常重要的一种特征,文献[8]采用傅里叶描述子来提取光斑轮廓特征。文献[9]利用傅里叶描述子来表达高分辨率遥感影像4种不同地物的形状特征。文献[10]利用人工神经网络和傅里叶形状描述子来识别造船部件。文献[11]使用傅里叶描述子描述服装轮廓特征。以上研究证明了傅里叶描述子提取目标轮廓特征的优越性,但针对提取滑坡遥感影像轮廓特征,傅里叶描述子是否具有实用性还有待研究。
本文提出基于傅里叶描述子的遥感影像滑坡轮廓特征提取方法,可以通过少量描述子重构滑坡轮廓。实验结果表明,傅里叶描述子可以有效提取滑坡轮廓特征,计算量小,具有良好的普适性。
本文研究区位于滑坡灾害频发的四川省绵竹市北部山区,影像来源于谷歌地球。滑坡影像数据获取日期为2014年12月23日,空间分辨率为0.25 m。
滑坡平面有舌形、椭圆形、长椅形、倒犁形、牛角形和平行四边形等11种[12]形状,本文实验数据选取倒犁形、舌形和长条形3种不同形状的滑坡高分辨率遥感影像。本次实验采用MatlabR2018b编程实现,运行环境为Intel(R)-Core(TM)-i5-1135G7,2.00 GHz,内存8 GB,64位Window10操作系统。
在对滑坡轮廓特征提取之前,需要对滑坡图像进行预处理。预处理的目的是消除图像中与滑坡目标无关的信息,去除噪声,最大限度地简化数据,增强滑坡信息的可检测性,提高滑坡轮廓提取精度。预处理的一般过程有灰度化、去噪、图像平滑、形态学处理和边缘检测等。预处理过程如图1所示。
图1 滑坡遥感影像预处理流程Fig.1 Flow chart of landslide remote sensing image preprocessing
获取滑坡轮廓的过程如图2所示,图2(a)为滑坡原始遥感影像,可以看出滑坡在影像上的色调与周围环境有明显的差别,呈灰白色,滑坡体表面几乎没有植被覆盖。图2(b)将滑坡遥感影像转化为灰度图像,利用中值滤波来去除图像中的噪声。图2(c)选择OTSU(最大类间方差)阈值分割法对滑坡图像进行二值化,选择最佳阈值使滑坡目标与背景进行分割,但同时受周围地物的干扰,使滑坡目标与无关信息分离的不够充分。图2(d)中首先对分割后的图像删除小面积像素信息的干扰,其次对图像进行孔洞填充,然后构造结构元素,进行形态学开运算,使滑坡目标与背景完全分离开。图2(e)利用Canny边缘检测,获得滑坡轮廓。
(a) 滑坡原始遥感影像
(a) 100个描述子
2.1.1 傅里叶描述子的原理
傅里叶描述子就是将图像的轮廓从空间域转换到频域,提取频域信息来表示目标形状特征的方法。首先提取目标的边界函数,然后将其展开后的系数幅度值作为目标轮廓的特征向量[13]。本文将遥感影像中滑坡的边界信息通过傅里叶变换将滑坡的轮廓特征从空间域变换到频域内,提取滑坡的频域信息作为特征向量,将滑坡轮廓数字化,从而区分滑坡的不同轮廓,进而达到识别不同形状滑坡的目的。
假设由K个点构成图形边界,从起始点(x0,y0)开始,沿逆时针方向追踪至边界点,将这些点的坐标当作复数来处理,可表示为:
z(k)=x(k)+y(k),k=0,1,2,…,K-1。
(1)
通过一维离散傅里叶变换,变换系数z(k)可表示为:
(2)
Z(u)称为边界的傅里叶描述子,通过傅里叶反变换可以重构z(k),表示为:
(3)
2.1.2 傅里叶描述子的归一化
通过上述方法得到的傅里叶描述子为了使其具有平移、旋转和尺度不变性,需要对其进行归一化处理。根据傅里叶变换的性质,将目标轮廓边界起始点位置平移α长度,放大r倍,旋转角度φ和平移(x0,y0)后,新的傅里叶变换系数Z′(u)为:
Z′(u)=F[(x′+iy′)rejφ+(x0+iy0)]=
rejφF(x(l+α)+iy(l+α))+F(x0+iy0)=
(4)
u=0,1,2,…,K-1,
(5)
(6)
d(u)对于旋转、位移、尺度和起始点位置不存在依赖关系,即得到归一化傅里叶描述子,且具有平移、旋转和尺度不变性,是一个鲁棒性较好的图像特征[14]。
2.1.3 傅里叶描述子提取滑坡轮廓特征
利用傅里叶变换得到滑坡边界的傅里叶描述子,在计算傅里叶反变换进行图像还原时,高频分量决定滑坡细节部分,低频分量决定滑坡总体轮廓。如果描述子过多就会造成特征冗余,增加计算量。随着特征向量的减少,如果描述子过少则无法表达滑坡的总体轮廓特征。因此,需要选择合适的描述子数量来提取滑坡的轮廓特征。
从图3可以看出,用100个描述子重构得到的图3(a)产生的轮廓与原始滑坡轮廓十分相似。图3(b)显示了50个描述子描述的结果,可以看出丢失了少量的细节信息,存在一些稍平滑的部分。由于高频分量决定滑坡细节信息,低频分量决定滑坡总体轮廓。随着描述子的减少,描述的滑坡轮廓会丢失许多细节信息。图3(c)显示了30个描述子描述的结果,丢失了大量的细节信息,但保留了滑坡基本的细节特征。图3(d)显示了20个描述子描述的结果,保留了滑坡总体轮廓特征。图3(e)是10个描述子描述的结果,已经产生了无法接受的失真,产生的边界越来越平滑,丢失了关键的轮廓信息。因此,使用30个描述子可以基本描述滑坡轮廓的细节特征,使用20个描述子可以基本描述滑坡总体轮廓特征。
(a) 原始图像
(a) 漏勺形滑坡
矩特征主要表征了图像区域的几何特征,又称为几何矩。Hu几何矩具有平移、旋转和尺度不变性,7个不变矩M1~M7是由二阶和三阶归一化中心矩组成的[15-16]。在本文中,通过Hu不变矩对滑坡图像计算得到的7个不变矩构成一组特征向量,可以验证滑坡轮廓的平移、旋转和尺度的不变性。其公式定义如下:
M1=η20+η02,
M2=(η20-η02)2+4η112,
M3=(η30-3η12)2+(3η21-η03)2,
M4=(η30+3η12)2+(3η21+η03)2,
M5=(η30-3η12)(η30+η12)[(η30+η12)2-3(η21+η03)2]+
(3η21-η03)(η21+η03)[3(η30+η12)2-(η21+η03)2],
M6=(η20-η02)[(η30+η12)2-(η21+η03)2]+
4η11(η30+3η12)(η21+η03),
M7=(3η21-η03)(η30+η12)[(η30+η12)2-3(η21+η03)2]+
(3η12-η30)(η03+η12)[3(η30+η12)2-(η21+η03)2]。
(7)
选取长条形滑坡为原始图像,对滑坡进行平移、旋转、镜像、放大和缩放,验证Hu不变矩描述滑坡轮廓具有平移、旋转和缩放不变性,图像变换如图4所示。
通过Hu不变矩计算图像变换后和另外2种形状的滑坡轮廓形状的矩不变量,其结果如表1所示。从表1可以看出,通过对滑坡图像平移、旋转和缩放变换后的Hu不变矩特征值与原始滑坡图像的Hu不变矩特征值基本一致,同时也容易看出,原始滑坡图像与不同形状的滑坡图像的Hu不变矩特征值具有较大差异,特别是M2,M4,M5,M6和M7可以明显区别出不同形状的滑坡。
表1 滑坡轮廓的矩不变量Tab.1 Moment invariants of landslide contour
通过傅里叶描述子和Hu不变矩提取滑坡轮廓特征后,选取支持向量机作为滑坡轮廓特征的分类器。本文采用LibSVM工具箱[17-18],对滑坡轮廓进行支持向量机多分类。首先,将描述滑坡轮廓特征的傅里叶描述子、Hu不变矩和滑坡不同形状类别标签输入到支持向量机分类器中。然后,选择线性核函数,其中Svmstrain参数选用“惩罚因子c=2,Gamma值g=0.07”。最后,将测试集的特征矩阵输入到支持向量机分类器中得到滑坡形状类别的分类准确率。
基于研究区的滑坡遥感影像,创建了一个小型滑坡遥感影像样本库,其中包含倒犁形、舌形和长条形滑坡影像各30幅。为提高运算速度,将影像尺寸压缩为600 pixel×500 pixel,随机抽取样本库中40%作为训练集,剩余60%作为测试集,滑坡不同形状的典型样本如图5所示。
为选取合适数量的傅里叶描述子描述滑坡轮廓,分别采用不同数量的傅里叶描述子(n取10,20,30,50,100)作为特征向量。利用LibSVM工具箱对样本集进行分类和测试,其分类准确率如图6所示。20个傅里叶描述子识别准确率最高,因此选择特征向量长度为20来提取滑坡轮廓特征。
图6 不同数量傅里叶描述子滑坡分类准确率Fig.6 Classification accuracy of sub-landslides described by different number of Fourier descriptors
由上文可知,选择20个傅里叶描述子作为支持向量机分类器特征向量的输入。融合特征是将傅里叶描述子的20个特征向量和Hu不变矩的7个特征向量采用串联的方式进行融合。由于2种特征取值范围存在较大的差异,且不同特征表达的意义也不同,因此融合之后对特征向量进行归一化,将每一个特征向量归一化至[0,1][19]。为了验证傅里叶描述子可以很好地提取滑坡轮廓特征,对比分析了Hu不变矩和其二者的融合特征的识别效果。
对样本集分别提取傅里叶描述子、Hu不变矩和二者的融合特征,利用支持向量机对样本集进行分类训练与测试,其识别结果如表2所示。结果显示,3种方法消耗时间相近,Hu不变矩的分类准确率最低(87.04%),融合特征(傅里叶描述子和Hu不变矩)的分类准确率(90.74%)低于仅使用傅里叶描述子的识别效果(92.59%)。因此,傅里叶描述子可以有效地描述滑坡影像的轮廓特征,并具有较好的性能,故傅里叶描述子更适用于滑坡遥感影像轮廓特征提取。
表2 不同特征提取方法实验对比Tab.2 Experimental comparison of different feature extraction methods
本文针对遥感影像不同形状的滑坡,提出了一种基于傅里叶描述子的滑坡影像轮廓特征提取方法,实验结果表明:
① 运用少量描述子就可以很好地描述滑坡遥感影像的轮廓特征,减少了数据冗余,提高了计算效率。
② 利用Hu不变矩验证了滑坡轮廓具有平移、旋转和尺度不变性,并且不同形状滑坡之间的不变矩具有较大差异。
③ 采用多分类支持向量机对3种不同形状的滑坡进行分类,通过对90幅滑坡影像的样本库进行实验测试发现:傅里叶描述子比Hu不变矩和融合特征能更好地描述滑坡的轮廓特征。因此,傅里叶描述子是一种稳定有效、适用于滑坡遥感影像滑坡轮廓特征提取的方法。
④ 本文不足之处在于利用部分傅里叶描述子代替滑坡轮廓存在误差,选择多少个傅里叶描述子能够表达滑坡的轮廓特征为经验性值,缺乏有效的依据。优化傅里叶描述子数量重构滑坡轮廓和增加样本库容量提高分类精度是接下来需要继续研究的内容。