高光谱和空间分辨率遥感影像近几年受到越来越多的关注,这类高光谱数据被广泛应用于城市地区土地覆盖/土地利用制图上。为了广泛利用高光谱图像的空间结构以及光谱信息,国内外学者已展开了若干研究。其中Fauvel将光谱信息与提取的目标特征信息重构目标特征矢量,并通过分类实验验证了提出方法的有效性[1]。Camps-Valls在分类过程中分别构造谱域信息与空域信息的核函数,并且研究了几种核函数组合方式[2]。Tarabalka利用聚类算法分割高光谱图像,进而结合分割结果以及光谱特征分类结果,总体分类结果有一定提高[3]。
纹理分析在计算机视觉领域有着广泛的应用。研究表明,Gabor纹理信息能够很好地表示有用的空间信息,在纹理分割研究中展现了良好的应用价值[4]。文献[5]将Gabor纹理信息应用于面部图像的表示和识别。文献[6]讨论了在高光谱图像分类技术中引入Gabor纹理信息。
经验模态分解(EMD)是在1998年由N.E.Huang[7]提出的一种对于非平稳、非线性信号处理效果理想的方法,它能够将信号分解成有限个固有模态函数(IMF)和一个残差信号。EMD提取信号本质特征的方式是自适应的,对比小波变换,EMD表现出了更好的时频特性,在高光谱数据特征提取上的潜力和优势是十分明显的[8]。
在变换域上进行高光谱数据分析和处理已在分类、压缩、超分辨等方面广泛使用,具有减小信息冗余、降低噪声干扰等优点,显示出良好的效果[9]。针对现有Gabor纹理提取方法对高光谱图像空间信息提取不足的问题,提出一种基于经验模态分解的Gabor纹理提取算法。
1 支持向量机分类算法支持向量机(SVM)是一种快速,高效的机器学习算法,因其较好的分类性能而广泛应用于高光谱图像分类问题[10]。SVM基本原理介绍如下。
设样本数据为xi∈Rd,yi∈{+1,-1}为相应类别的样本标号,i=1,2,…,n。一般形式的线性判别函数为g(x)=〈ω,x〉+b分类面表示为〈ω,x〉+b=0。式中x=[x1 x2…xd]T是d维特征矢量,权向量表示为ω=[ω1 ω2…ωd]T。b为常数,称为阈权值。
当所要解决的问题是非线性时,SVM引入非线性映射φ(·)的方法,这样,原始空间不能线性分开的数据点就转化成了高维空间上的可分点。此时通过解决式(1)来获得最优分类超平面:
$\min \left\{ {\frac{1}{2}{{\left\| \omega \right\|}^2} + C\sum\limits_{i = 1}^n {{\xi _i}} } \right\}$ | (1) |
式中:ξi为松弛项,常值C>0,它代表对分类结果不正确样本的惩罚程度。约束条件为
$\begin{array}{c} {y_i}\left[{\left\langle {\omega ,\varphi \left( {{x_i}} \right)} \right\rangle + b} \right] \ge 1 - {\xi _i}\\ {\xi _i} \ge 0,i = 1,2,\cdots ,n \end{array}$ | (2) |
为了求解问题(1),建立并求解其对偶问题,即最大化下式函数:
$L\left( {\omega ,b,\alpha } \right) = \sum\limits_{i = 1}^n {{\alpha _i}} - \frac{1}{2}\sum\limits_{i,j = 1}^n {{\alpha _i}{\alpha _j}{y_i}{y_j}K\left( {{x_i},{x_j}} \right)} $ | (3) |
式中:αi为与每个样本对应的Lagrange乘子。在约束条件0≤αi≤C和$\sum\limits_{i = 1}^n {{\alpha _i}{y_i} = 0} $下,最大化目标函数(3)的解是唯一的[11]。设(α*,b*)为该优化问题的最优解。对于测试向量,最终的决策函数为
$f\left( x \right) = {\mathop{\rm sgn}} \left\{ {\sum\limits_{i = 1}^n {\alpha _i^*{y_i}K\left( {{x_i},x} \right) + {b^*}} } \right\}$ | (4) |
SVM是二元分类器,将其扩展到多类别问题通常是组合多个二值分类器。1-a-r(1-against-rest)多类分类器和1-a-1(1-against-1)多类分类器是2种常用的多类分类器结构。实验采用1-a-r型支持向量机来进行高光谱数据分类。
2 空谱结合分类 2.1 经验模态分解经验模态分解(empirical mode decomposition,EMD)法是一种新型自适应信号时频处理方法[12],它能够将复杂数据表示成数量有限并且较少的固有模态函数(intrinsic mode function,IMF)分量的加和形式。固有模态函数满足如下条件[8]:
1) 函数的局部极值点和过零点数在整个时间范围内必须相等,或最多相差为一;
2) 在任意时刻,由局部最大值确定的上包络和由局部最小值确定的下包络的平均必须为零。
EMD获得IMF是通过反复循环的筛选步骤来实现的。Xl,m(n)(i,j)用于显示在第n次迭代用来寻找第l个频带的第m个IMF的当前值。迭代过程从高光谱图像第一维数据开始,筛选过程如下:
1) 找到当前波段上全部局部极大值点和局部极小值点;
2) 使用插值算法把局部极大值拟合成上包络面Emax(i,j),局部极小值拟合成下包络面Emin(i,j);
3) 求得均值包络面Zm(n)(i,j)为
$Z_m^{\left( n \right)}\left( {i,j} \right) = \frac{{\left( {{E_{\max }}\left( {i,j} \right) + {E_{\min }}\left( {i,j} \right)} \right)}}{2}$ | (5) |
4) 原始信号与均值包络面相减得到差值Dm(n)(i,j)为
$D_m^{\left( n \right)}\left( {i,j} \right) = X_m^{\left( n \right)}\left( {i,j} \right) - Z_m^{\left( n \right)}\left( {i,j} \right)$ | (6) |
5) 检查结果是否符合迭代结束条件SD<τ,SD表示为
${\rm{SD}} = \max \left( {{\rm{abs}}\left( {Z_m^{\left( n \right)}\left( {i,j} \right)} \right)} \right)/\max \left( {{\rm{abs}}\left( {{R_m}\left( {i,j} \right)} \right)} \right)$ | (7) |
其中残余量Rm(i,j)以原始波段为首次输入。如果不满足停止条件,(此时为第N次迭代),将以Dm(n)(i,j)为初始信号从步骤1)开始重新迭代,直到满足终止条件,此时获得第1个IMF;
6) 当得到IMF函数后,更新残余量Rm(i,j)为
${R_m}\left( {i,j} \right) = {R_m}\left( {i,j} \right) - {\rm{IM}}{{\rm{F}}_{l,m}}\left( {i,j} \right)$ | (8) |
如果此残余量不含极值,此EMD程序终止。如果残余量含极值,将此残余量Rm(i,j)作为初始输入从步骤1)继续求解。
综上所述,光谱原始波段图像Xl(i,j)被所求得的IMF函数及残余量重构,表达式为
${X_i}\left( {i,j} \right) = \sum\limits_{m = 1}^M {{\rm{IM}}{{\rm{F}}_{l,m}}\left( {i,j} \right) + {R_M}\left( {i,j} \right)} $ | (9) |
Gabor函数形成的二维Gabor滤波器具有在空间域和频率域同时取得最优局部化的特性。因此能够很好地描述对应于空间频率(尺度)、空间位置及方向选择性的局部结构信息[13]。在空域,一个二维的Gabor滤波器是一个正弦平面波调制的高斯核函数,一般表达式为
$\varphi \left( {x,y,f,\theta ,\gamma ,\sigma } \right) = \exp \left( { - \frac{{{x^{'2}} + {\gamma ^2}{y^{'2}}}}{{2{\sigma ^2}}}} \right)\cos \left( {2{\rm{\pi }}fx' + \varphi } \right)$ | (10) |
式中:图像中某点像素位置参数x和y,x′=xcos θ+ysin θ,y′=-xsin θ+ycos θ,f和θ各自表示Gabor函数的频率和方向,φ为Gabor函数相角,γ为空间宽高比,Gabor函数高斯因子的标准差是σ,σ的大小根据正实数带宽b(一般是1)改变。它们之间的关系为
$b = {\rm{lb}}\frac{{\frac{\sigma }{\lambda }{\rm{\pi }} + \sqrt {\frac{{\ln 2}}{2}} }}{{\frac{\sigma }{\lambda }{\rm{\pi }} - \sqrt {\frac{{\ln 2}}{2}} }},\frac{\sigma }{\lambda } = \frac{1}{{\rm{\pi }}}\sqrt {\frac{{\ln 2}}{2}} \cdot \frac{{{2^b} + 1}}{{{2^b} - 1}}$ | (11) |
式中:λ=1/f。
由于图像中各种纹理包含各自的带宽和中心频率,单独的Gabor滤波器并不能提取全部的纹理信息,因而对纹理图像滤波时,需要采用一系列包含不同带宽和频率的Gabor滤波器。实验中设置了5个尺度6个方向的滤波器组:
$\begin{array}{c} {\varphi _{u,v}} = \varphi \left( {x,y,{f_u},{\theta _v},\gamma ,\sigma } \right),\\ {f_u} = {f_{\max }}/{\sqrt 2 ^u},{\theta _v} = \frac{v}{6}{\rm{\pi ,}}\\ {\rm{u = 0,1,}} \cdots {\rm{,4,v = 0,1,}} \cdots {\rm{,5}} \end{array}$ | (12) |
所得Gabor特征集通过输入图像与这30个滤波器组卷积而得到
${O_{u,v}}\left( {x,y} \right) = \left( {I*{\varphi _{u,v}}} \right)\left( {x,y} \right)$ | (13) |
不同于原始Gabor纹理提取方法从原始高光谱图像中直接提取纹理信息的方式,提出一种结合二维经验模态分解的变换域纹理提取方法,并通过分类实验验证提出方法的有效性,实验流程图如图 1所示。
图 1中,以L波段高光谱图像为例,每一波段包含M×N个像素点;分量1、分量2、分量h分别表示第h个EMD变换域分量;Q表示Gabor滤波器组数。
首先,求输入信号X′的局部极值,利用插值算法求出上包络、下包络和均值包络。然后将输入信号减去均值包络,获得第一个分量D1(i,j),判断D1(i,j)是否满足终止迭代条件。如果满足则D1(i,j)是第1个IMF,如果不满足则把D1(i,j)作为初始信号开始重新迭代,直到满足终止条件,此时获得第1个IMF,并求得残余量R1(i,j)。判断残余量是否包含极值,如果此残余量不含极值,EMD程序终止。假如残余量含极值,以此残余量R1(i,j)当做初始信号继续求解下一个IMF和下一个残余量,当残余量不含极值,获得全部IMF,EMD程序终止。
由于固有模态函数的本质是信号在不同频率范围内原信号的本质特征,原信号噪声的集合表示为残差形式,所以提取信号的本质特征可通过保留IMF分量和丢弃残余量来实现。在改进纹理提取算法中,经验模态分解只保留低阶的IMF分量,高阶的IMF分量及残余量不予考虑。这样,就可以保证包含空间关系信息的低阶IMF分量被保留,而较高阶的IMF分量(通常缺乏空间结构信息)将被丢弃。
其次,经验模态分解之后,原始高光谱图像被表示为有限个IMF分量的集合,对各IMF集合利用主成分分析法(PCA)提取各集合数据的第1主成分。
最后,进行空间纹理滤波处理。空间纹理滤波建立在各主成分分析提取的第1主成分上。
针对改进的纹理提取算法,采用两种信息融合方式,并通过SVM分类算法验证结果的有效性。
方法1:特征层融合,实验中所用特征层融合采用复和核函数的方法,各项权重均为1。xisn∈RNs,n=1,2,…,h为样本点空域特征,xiw∈RNw为像元的谱域特征,Ksn、Kw分别为各自的核矩阵。核函数的复合形式表示如下:
$K\left( {{x_i},{x_j}} \right) = {K_w}\left( {x_i^w,x_{jj}^w} \right) + \sum\limits_{n = 1}^h {{K_{sn}}\left( {x_i^{sn},x_j^{sn}} \right)} $ | (14) |
方法2:图像层融合,实验中采用原始高光谱图像经PCA提取前三维,与纹理数据进行重组,利用所得重组数据进行分类。
3 分类实验及结果评价 3.1 实验数据实验选用2组常用来验证高光谱图像分类算法的高光谱数据,分别为印第安针叶林、帕维亚大学数据。
第1数据集为1992年6月由AVIRIS拍摄的印第安纳西北部农业区高光谱航空影像,像素信息144×144,200维波段。该数据光谱范围0.4~2.5 μm,空间分辨率20 m。在原始的16类样本数据中忽略小样本集数据,采用样本数量最多的9类地物进行分类研究。这9类样本类别及样本数量如表 1所示。
地物名称 | 像素个数 |
Corn-no till | 1 434 |
Corn-min till | 834 |
Grass/Pasture | 497 |
Grass/Trees | 747 |
Hay-windrowed | 489 |
Soybean-no till | 968 |
Soybean-min till | 2 468 |
Soybean-clean till | 614 |
Woods | 1294 |
第2数据集是包含340×610像素和103维波段数据。该数据光谱范围是0.4~0.85 μm,空间分辨率为1.5 m。本实验选取数据包含8类地物的样本数据来验证文中所提方法的有效性,这8类样本及数量如表 2所示。
地物名称 | 像素个数 |
Asphalt | 559 |
Meadows | 785 |
Trees | 154 |
Metal sheets | 721 |
Soil | 4 021 |
Bitumen | 785 |
Bricks | 129 |
Shadows | 285 |
实验结果通过不同的训练样本数据率(training data rates,TDR)来进行证明。TDR为5%表示随机抽取每类样本数据的5%为训练样本,剩余样本做测试样本。二维经验模态分解部分,对于每一波段图像,以尺度100将图像镜像延拓;插值方式采用cubic三次多项式插值;τ取0.3。SVM采用高斯径向基核函数,并采用网格搜索法搜索最佳参数(数据集1中惩罚系数为1 000,高斯径向基核函数参数σ=1.8,数据集2中分别为700和1.2)。
3.3 结果与讨论按照前文算法流程对高光谱图像进行分类实验,实验结果如图 2所示。
图 2(a)为原始高光谱图像;图 2(b)为利用SVM对原始高光谱数据9类地物分类结果;图 2(c)为对原始高光谱图像PCA降维后取第一主成分分量利用Gabor滤波器组提取纹理信息,与原始高光谱图像特征层融合经SVM分类结果;图 2(d)为改进方法。从图 2可以看出,改进方法进一步消除了“麻点”现象。
图 3为印第安针叶林9类地物各类别精度对比,可以看出,改进方法对各类地物分类精度较Gabor-PCA方法有不同程度的提高。
由表 3分析得出,对于只应用SVM分类的结果来说,结合空间信息后分类精度已经有很大程度的提高。然而,由于原始高光谱图像存在噪声干扰、信息冗余等因素,直接对原始图像提取纹理信息并不能达到高精度分类的目的。改进方法在对原始图像经过自适应的二维经验模态分解之后,对于分解所得IMF分量分别提取纹理信息,去除掉部分冗余信息的同时也增加了可融合的纹理信息。在TDR为5%的情况下可使整体精度从91.05%提高到94.79%。总体分类精度得到明显提高。
方案 | 5%训练样本 | 10%训练样本 | ||
精度/% | Kappa系数 | 精度/% | Kappa系数 | |
SVM | 81.26 | 0.777 8 | 84.45 | 0.816 6 |
Gabor-PCA | 91.05 | 0.894 7 | 95.08 | 0.942 2 |
改进算法 | 94.79 | 0.938 8 | 97.75 | 0.973 6 |
此外,从表 4可以看出,随着融合信息的增加,精度的提高已经不是很明显了。特别的,在训练样本数量为总体5%时,融合前4个IMF对比融合前3个IMF的精度反而有下降。这是纹理信息冗余所导致的,因此,综合计算量及分类精度的折中考虑,只利用前3个IMF即可获得较好效果。
IMF数量 | 5%训练样本 | 10%训练样本 | ||
总体精度/% | Kappa系数 | 总体精度/% | Kappa系数 | |
1IMF | 89.17 | 0.872 3 | 93.22 | 0.920 3 |
2IMFs | 93.04 | 0.918 2 | 96.45 | 0.958 3 |
3IMFs | 94.79 | 0.938 8 | 97.75 | 0.973 6 |
4IMFs | 94.76 | 0.938 3 | 97.97 | 0.976 1 |
为了证明改进方法的有效性,利用帕维亚大学数据进行验证。此实验的信息融合方式采用图像层融合,TDR为5%,分类结果如图 4所示。
精度评价结果如表 5所示。由表 5可得改进方法对于图像层融合方式同样有效。
表 6给出针对帕维亚大学数据不同IMF分量数对分类结果的影响,可以看出在保留3个IMF分量时结果最优。
IMF数量 | 5%训练样本 | |
总体精度/% | Kappa系数 | |
1IMF | 97.43 | 0.961 3 |
2IMFs | 97.63 | 0.964 4 |
3IMFs | 98.01 | 0.970 1 |
4IMFs | 97.96 | 0.969 2 |
此外,为了验证改进方法对含噪声高光谱图像分类处理的有效性,对两数据集所有波段添加方差为σ2的白噪声,并利用所得噪声数据进行分类实验。数据集1中TDR为5%,数据集2中TDR为2.5%,采用整体精度OA和Kappa系数表征分类效果,精度评价结果如表 7、8所示。
指标 | 方法 | 方差0 | 方差20 | 方差40 | 方差60 | 方差80 | 方差100 | 方差120 |
总体精度/% | SVM(含噪) | 81.26 | 80.98 | 80.68 | 80.21 | 79.42 | 79.24 | 78.91 |
Gabor-PCA | 89.26 | 89.16 | 89.30 | 89.34 | 89.31 | 89.18 | 89.12 | |
改进算法 | 94.29 | 94.06 | 94.30 | 93.42 | 93.90 | 93.80 | 93.98 | |
Kappa系数 | SVM(含噪) | 0.777 8 | 0.774 5 | 0.771 1 | 0.765 2 | 0.756 0 | 0.753 7 | 0.749 7 |
Gabor-PCA | 0.873 8 | 0.872 6 | 0.874 3 | 0.874 8 | 0.874 4 | 0.872 8 | 0.872 1 | |
改进算法 | 0.932 9 | 0.930 2 | 0.933 0 | 0.922 6 | 0.928 3 | 0.927 2 | 0.929 2 |
指标 | 方法 | 方差0 | 方差20 | 方差40 | 方差60 | 方差80 | 方差100 | 方差120 |
总体精度/% | SVM(含噪) | 89.06 | 89.00 | 88.52 | 88.02 | 88.17 | 88.25 | 88.01 |
Gabor-PCA | 94.89 | 94.86 | 95.13 | 94.80 | 94.91 | 94.77 | 95.16 | |
改进算法 | 95.56 | 95.59 | 95.55 | 95.47 | 95.58 | 95.75 | 95.87 | |
Kappa系数 | SVM(含噪) | 0.834 1 | 0.832 8 | 0.825 4 | 0.817 4 | 0.819 3 | 0.821 5 | 0.817 5 |
Gabor-PCA | 0.922 3 | 0.921 9 | 0.926 1 | 0.920 8 | 0.922 6 | 0.920 4 | 0.926 6 | |
改进算法 | 0.932 5 | 0.932 9 | 0.932 3 | 0.931 0 | 0.932 7 | 0.935 4 | 0.937 4 |
将噪声数据分类精度评价结果以折线图方式表示,如图 5、图 6所示。从图中可以看到,噪声的加入对普通SVM分类结果的影响是很大的,当噪声方差在20~120变化时,两数据集SVM分类精度损失分别超过2%和1%。随着噪声的增大,此时的含噪图像在SVM分类方法下终将因为分类精度的降低而失去研究意义。而原始Gabor-PCA方法可以缓解噪声带来的影响,原因是PCA方法本身包含去除噪声和冗余的特点,经过空间信息融合之后分类精度有很大提高。改进方法也具有良好的抗噪声性能,经验模态分解根据信号本身的特点滤去包含噪声的高频IMF分量及残余量,这对接下来的处理提供了消噪数据。同时,IMF分量的有限性也保证了最大程度融合纹理信息的同时又不至于纹理信息融合过度。当噪声方差在20~120变化时,整体分类精度并没有随着噪声强度的增加而有明显的下降。在含噪图像测试下,改进方法对分类性能的提升仍优于原有Gabor-PCA方法,再次证明了该方法的有效性。
4 结论1)与传统的Gabor-PCA算法相比,该算法在精度上明显提高很多,而且在抗噪性能上也有优秀的表现。
2)与传统高光谱Gabor纹理提取方法研究思路不同,提出算法将滤波操作转移到变换域上进行,最大程度保留原始数据信息的同时更加高效地挖掘了空间信息,为高光谱图像的纹理提取开辟了新的思路。
3)经实验验证,该算法很大程度消除了分类问题中样本错分现象,而且精度很高,不受噪声影响,是一种优势明显的算法策略。
[1] | FAUVEL M, BENEDIKTSSON J A, CHANUSSOT J, et al. Spectral and spatial classification of hyperspectral data using SVMs and morphological profiles[J]. IEEE transactions on geoscience and remote sensing, 2008, 46(11): 3804-3814. |
[2] | CAMPS-VALLS G, GOMEZ-CHOVA L, MUNOZ-MARI J, et al. Composite kernels for hyperspectral image classification[J]. Geoscience and remote sensing letters, 2006, 3(1): 93-97. |
[3] | TARABALKA Y, BENEDIKTSSON J A, CHANUSSOT J. Spectral-Spatial classification of hyperspectral imagery based on partitional clustering techniques[J]. IEEE transactions on geoscience and remote sensing, 2009, 47(8): 2973-2987. |
[4] | LI Chaorong, DUAN Guiduo, ZHONG Fujin. Rotation invariant texture retrieval considering the scale dependence of Gabor wavelet[J]. IEEE transactions on image processing, 2015, 24(8): 2344-2354. |
[5] |
王科俊, 邹国锋. 基于子模式的Gabor特征融合的单样本人脸识别[J]. 模式识别与人工智能, 2013, 26(1): 50-56. WANG Kejun, ZOU Guofeng. A sub-pattern Gabor features fusion method for single sample face recognition[J]. Pattern recognition and artificial intelligence, 2013, 26(1): 50-56. |
[6] | JIA Sen, SHEN Linlin, DENG Lin. Band selection-based Gabor wavelet feature extraction for hyperspectral imagery classification[C]//Proceedings of the 4th workshop on hyperspectral image and signal processing. Shanghai, China, 2012: 1-4. |
[7] | HUANG N E, SHEN Zheng, LONG S R, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[C]//Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. London, 1998, 454: 903-995. |
[8] |
沈毅, 张敏, 张淼. 基于互信息波段选择和经验模态分解的高精度高光谱数据分类[J]. 激光与光电子学进展, 2011, 48(9): 59-66. SHEN Yi, ZHANG Min, ZHANG Miao. Mutual information bands selection and empirical mode decomposition based support vector machines for hyperspectral data high-accuracy classification[J]. Laser and optoelectronics progress, 2011, 48(9): 59-66. |
[9] |
王立国, 赵妍. 基于MAP的高光谱图像超分辨率方法[J]. 光谱学与光谱分析, 2010, 30(4): 1044-1048. WANG Liguo, ZHAO Yan. MAP based super-resolution method for hyperspectral imagery[J]. Spectroscopy and spectral analysis, 2010, 30(4): 1044-1048. |
[10] | KUO B C, HO H H, LI C H, et al. A kernel-based feature selection method for SVM with RBF kernel for hyperspectral image classification[J]. IEEE journal of selected topics in applied earth observations and remote sensing, 2014, 7(1): 317-326. |
[11] |
张学工. 关于统计学习理论与支持向量机[J]. 自动化学报, 2000, 26(1): 32-42. ZHANG Xuegong. Introduction to statistical learning theory and support vector machines[J]. Acta automatica sinica, 2000, 26(1): 32-42. |
[12] | DEMIR B, ERTURK S. Empirical mode decomposition of hyperspectral images for support vector machine classification[J]. IEEE transactions on geoscience and remote sensing, 2010, 48(11): 4071-4084. |
[13] |
肖倩. 结合空间信息与光谱信息的高光谱图像分类研究[D]. 哈尔滨: 哈尔滨工程大学, 2013: 35-38. XIAO Qian. Combination of spatial information and spectral information for hyperspectral imagery classification[D]. Harbin: Harbin Engineering University, 2013: 35-38. |