基于SPCA和域变换递归滤波的高光谱图像分类

于 多，黄永东 ，2

1.北方民族大学 图像处理与理解研究所，银川 750021

2.大连民族大学 数学与信息科学研究中心，辽宁 大连 116600

高光谱遥感图像是由几十个到数百个连续波段图像组成的三维数据立方体，其不仅反映了地物的光谱信息，而且包含了地物空间分布信息，因而具有图谱合一的特点，具有较强的地物分类和监测能力[1]，被广泛地应用在地质勘探、考古发现、森林火灾检测、军事侦察、环境检测和目标识别等领域[2-5]。像元分类是高光谱遥感图像研究的重要内容之一。一方面，与全色和多光谱遥感图像相比，高光谱图像有更高的光谱分辨率和相对较低的空间分辨率，单纯依靠光谱信息的高光谱分类方法，往往不能达到用户对地物分类不断提高的类别精细程度和精度的要求，越来越多的学者开始尝试综合光谱与空间信息进行高光谱分类研究[2]。另一方面，由于高光谱遥感图像维数高，训练样本数目往往相对不足，无法有效统计分类模型的相关参数，传统分类方法运算速度慢，分类精度低，在分类过程中可能会导致Houghes现象[6]。为了提高分类精度，一些新的分类方法被提出，如支持向量机[7]、决策树[8]、神经网络[9]、稀疏表示[10]等分类方法。Toksoz等人在2016年基于稀疏表示理论设计出一个新的高光谱图像分类器，即基本阈值分类器（Basic Thresholding Classifier，BTC）[11]，其具有分类精度高、易实现、无需调节参数、快速有效且能够联合空间信息等特点。本文选择BTC作为分类器。

在高光谱遥感图像分类算法中，光谱信息和空间信息的特征提取是关键环节和核心研究内容之一，它们直接决定着空谱特征的维数和影响分类精度。许多学者一直致力于研究光谱特征提取方法，以期达到降低特征空间维数和提高分类精度的目的[12]。常用的特征提取方法有主成分分析方法（Principal Component Analysis，PCA）[13]、独立成分分析方法（Independent Component Analysis，ICA）[14]、局部保护判别分析方法[15]和分割的主成分分析方法（Segmented Principal Component Analysis，SPCA）[16]等。与PCA 相比，SPCA 不仅具有时间复杂度低的特点，而且具有能够有效保留光谱差异性信息的优点[16]，这样有利于分类精度的提高。

高光谱遥感图像的空谱分类方法受到广泛关注并取得了大量研究成果[17-21]。文献[17]通过引入本质特征分解来提取高光谱图像的本质特征，并用于分类，该方法能够有效地移除与目标地物无关的信息。文献[18]通过联合超像素和本质特征分解来进一步改善分类精度。近几年，边缘保持滤波已经成功地应用到图像分类中，边缘保持滤波器在平滑图像的同时又能保留图像边缘信息以及整体轮廓[19]。文献[19]提出了基于边缘保持滤波的空谱分类方法，取得了很好的分类效果。文献[20]提出基于递归滤波的空谱分类方法，有效改善分类精度。文献[21]将边缘保持滤波中的滚动引导滤波引入到高光谱图像分类中，有效地改善了分类精度。边缘保持滤波提取特征的效果表明，其在表达高光谱图像的主要空谱结构方面是有效的。然而，单参数边缘保持平滑操作的局限性在于它有降低不同类地物之间的光谱差异的倾向，进而影响了分类效果。为了克服这个不足，文中采用多层（多参数）边缘保持滤波器进行滤波，并利用它们之间信息的互补性进行特征融合，从而有效地提高了分类精度。

基于以上分析，本文提出一种基于分割的主成分分析（SPCA）和域变换递归滤波（Domain Transform Recursive Filtering，DTRF）的高光谱图像分类算法。主要步骤如下：首先，使用SPCA降低数据维数，即将高光谱图像的波段分成一些相邻波段的子集，PCA作用在每一个子集上；然后，利用DTRF 获得高光谱图像的堆叠的边缘保持滤波图；接着，PCA 作用在堆叠的边缘保持滤波图上进行特征融合；最后，用BTC对上述融合图进行分类。

1 理论简介

1.1 域变换递归滤波

域变换递归滤波（DTRF）[22]是一种实时的边缘保留滤波，其在改善高光谱图像分类性能方面是非常有效的。滤波过程可分为两步：第一步进行域变换；第二步是递归滤波。在实际滤波过程中，为了控制滤波器的尺寸和模糊度，常常将域变换定义为应用近似距离变换。即对于给定的一维信号I，域变换定义为：

这里，Ui是域变换信号，δs和δr是用来调节滤波器平滑程度的两个参数。然后，输入信号通过递归滤波器进行处理，即：

其中，Ji是第i个像元的滤波输出，a∈[0,1)代表反馈系数，b表示变换域中两个相邻样本的距离。当图像接近边缘时，b不断增大，ab趋于0，迭代运算终止，从而达到边缘保留的作用。

对二维信号，域变换递归滤波通过沿着图像的每一个维度执行上述所提的一维操作，进而处理图像。本文用DTRF(I,δs,δr)表示域变换递归滤波。

1.2 BTC算法

BTC 是Toksoz 等人基于稀疏表示理论提出的一种高光谱图像分类算法，其具有分类精度高、易实现、无需调节参数、快速有效且能够联合空间信息等特点[11]。它是一个逐像素的分类器，仅仅利用光谱特征来分类每一个像素。最终对每一个像素，产生一个带有类别标签的残差向量。为了改善BTC 分类精度，将算法扩展到空谱邻域，首先利用SPCA 提取光谱特征，然后利用多层域变换递归滤波提取边缘特征，再利用PCA 融合这些边缘特征，最后利用BTC 对融合后的光谱波段（记作H∈ℝn1×n2×B，共有n1×n2个像素，B为波段数）进行逐像素的分类，产生带有类标签的残差向量ε∈ℝC（C是类别数），残差向量形成一个堆叠的图像ℝn1×n2×C，残差图表示为mapi(mapi∈ ℝn1×n2，其中i∈{1,2,…,C}) ，最后基于最小残差来决定每一个像素的类标签。具体BTC算法大致如下：

第一步是计算测试像素y和字典A中的每一个像素的相关性。第二步是计算索引集Λ，由与y的相关性最大的前M个像素所对应的索引组成。第三步是构造矩阵D，矩阵D是由索引在Λ中的训练像素组成。第四步在求解稀疏编码值x^(Λ)中，利用了Tikhonov 正则化，因为矩阵D所包含的像素可能属于同一类，所以D的列变得线性相关，在算法中的逆运算将无法进行，而Tikhonov 正则化次定方程组是一个有效消除这种限制的方法[23]。最后基于最小残差来决定测试像素的类别，其中残差对应于第j类稀疏编码部分。BTC算法详细过程可参见文献[11]。

算法1BTC算法

1.3 PCA融合

选择PCA作为融合工具，是因为它易于实现，没有参数，可以有效提取高光谱图像中的相关信息[24]。具体原因如下：首先，由多参数的域变换递归滤波获得的滤波图不仅仅包含了互补信息，而且还有大量的冗余信息，而PCA可以很好地去除这些冗余信息。再者，滤波器的平滑操作尽管可以有效移除噪声和小尺度的细节，但是也降低了属于不同类地物的像素的差异值（见图1）。为解决这些问题，PCA 是一个很好的工具，它能够提取在滤波图中最相关的信息，也因此能够有效增强属于不同类地物像素的光谱差异[25]，这提高了下一步BTC分类中判断测试像素与训练像素之间的相关性强弱的准确性。

2 基于SPCA-DTRF高光谱遥感图像分类

本文提出一种基于SPCA-DTRF的高光谱遥感图像分类算法，主要由SPCA 方法来降维和特征提取，并利用不同参数设置的DTRF 滤波器得到堆叠的边缘保持滤波图；接着对上述滤波图进行PCA特征融合，实现光谱信息和空间信息的有效利用；最后利用BTC 进行分类。具体步骤如下：

步骤1对高光谱图像进行波段子集划分。把高光谱遥感图像的光谱波段划分为L个相邻波段的子集，每个子集大小（子集含有的波段数）记作k。第ℓ个子集(ℓ≤L)可以表示为：

实验中通过设定子集的大小k，进而确定波段数L，即当k整除M时，当k不整除M时，L=其中X={X1,X2,…,XM}表示高光谱图像包含M个光谱波段，表示取整。

步骤2提取每一个波段子集的第一主成分。对于高光谱遥感图像，它的相邻光谱波段常常有很强的相关性，同时也包含了冗余信息[26]。为了降低后续特征提取与融合和BTC 分类的计算复杂度，更好地保留光谱差异信息，将PCA作用于波段子集Sℓ上，并提取第一主成分，即有：

步骤3把不同波段子集的第一主成分堆叠在一起，得到高光谱数据，即：

其中，Y表示维数降低后的高光谱数据。

步骤4对各波段子集提取的第一主成分进行域变换递归滤波处理，得到堆叠的边缘保持滤波图。对Yℓ用多层（多参数）域变换递归滤波图进行特征提取，并将它们叠加在一起，即：

图1 参数δr 和δs 对滤波效果的影响

图2 算法流程

步骤5对堆叠的边缘保持滤波图用PCA进行特征融合。如果直接将堆叠的边缘保持滤波图输入到分类器，尽管可利用它们之间的互补信息，但是可能导致Houghes 现象和额外计算负担。另一方面，尽管DTRF滤波在平滑图像的同时也能保留边缘，但是平滑过程中不可避免地降低了不同类地物像素的差异性，从而影响分类精度。为了解决上述问题，本文利用PCA 进行特征融合，具体操作如下：

这里，p是保留在L中的主成分数量，F是堆叠的边缘保持滤波图。换句话说，利用PCA 提取堆叠的边缘保持滤波图中前p个主成分。

步骤6用BTC 来确定每个像素的类别，像素值最小的分类残差图对应的类别即为该像素类别。

本文所提算法流程图见图2所示。

3 实验仿真

为了对本文所提算法的分类效果进行评估，使用了两组众所周知的高光谱遥感图像数据集Indian Pines和Salinas进行实验。

3.1 实验数据

（1）Indian Pines 数据集是由机载可见红外光成像光谱仪（AVIRIS）获得的，是在印第安那州西北部的印度松树测试场采集到的，由波长范围是0.4～2.5 nm的光谱波段组成，图像大小为145×145 像素，空间分辨率为20 m，包含16个类别，去除水汽吸收波段后，从204个波段中选取200 个波段作为研究对象。该数据的伪彩色图、参考图和类别图见图3。

（2）Salinas 数据集是由AVIRIS 在加利福尼亚州的农业区域Salinas Vally采集到的，由波长范围是0.4～2.5 nm的224个光谱波段组成，图像大小为512×217像素，空间分辨率是3.7 m，包含16个类别，移除水汽吸收严重的波段，从224 个波段中选200 个波段作为研究对象。该数据的伪彩色图、参考图和类别图见图3。

图3 Indian Pines和Salinas实验数据

3.2 实验分析

为了验证算法的有效性，将本文算法（SPCA-DTRF）与基本阈值分类方法（BTC）[11]、支持向量机方法（Support Vector Machine，SVM）[7]、正交匹配追踪法（Orthogonal Matching Pursuit，OMP）[27]、稀疏正交匹配追踪法（Simultaneous OMP，SOMP）[27]、融合相关系数和联合稀疏表示方法（Correlation Coefficient and Joint Sparse Representation，CCJSR）[28]和边缘保持滤波方法（Edge-Preserving Filtering，EPF）[19]等7 种方法进行比较，其中BTC、SVM、OMP是光谱分类方法，SOMP、CCJSR、EPF、BTC-WLS是空谱分类方法，OMP和SOMP利用了稀疏表示分类器，EPF所用分类器是支持向量机。采纳这些方法，是因为在高光谱遥感领域它们不仅被大量地引用，而且在一些真实的高光谱数据中有很好的表现。支持向量机是通过五层交叉验证的高斯核来实现的，其他方法均采用相应文献的默认参数来实现。为使算法具有可比性，实验中各类样本被随机分成训练样本和测试样本，所有实验都按这种方法进行10次，并对分类结果取平均值。在实验中，本文采用3 种常用分类指标，即总体分类精度（OA）、平均分类精度（AA）和Kappa 系数（Kappa）。3 种指标所得数值越大，表示分类效果越好。所有的实验在Matlab2017a 上进行，电脑配置是3.5 GHz CPU，8 GB安装内存。

3.2.1 Indian Pines数据集

对于Indian Pines 数据集，10%的数据被随机挑选作为训练样本，剩下90%的数据作为测试样本。本文算法与另外7 种分类方法所得分类精度见表1，分类图见图4。从表1 可以清楚地看到，本文算法OA、AA、Kappa 系数均是最优的，在16 个类别中10 个类别的分类精度最高，2 个类别分类精度排第2，4 个类别分类精度排第3。整体上，本文算法分类精度是最优的，分类效果最好。同BTC 方法相比，总体分类精度显著地增加了32%；与SVM 相比，分类精度高出15%左右；与SOMP、CCJSR、EPF、BTC-WLS 相比，分类精度高出了2%～4%，充分说明了算法的优越性。从图4 可以看出，由本文算法获得的分类图与参考图更加相符。下面分析参数k和p，参数和分类器选取对本文算法的影响。

（1）分析波段子集大小k和融合后的主成分个数p对本文算法SPCA-DTRF 的影响。为了便于分析，将每个波段子集所含波段数k取值分成5种情况（k=5，10，15，20，25），并且在5 种情况下观察参数对分类精度的影响，具体结果见图5。从图5可以得到如下结论：一是影响最优分类精度的最主要因素是波段子集大小。二是在波段子集大小k固定情况下，随着主成分个数p的增加，分类精度均呈现出递增，接着趋于平稳，然后略有降低的趋势。三是波段子集大小k=5 时整体分类精度是最优的。因此，实验中取k=5，p=120。

图4 Indian Pines数据集采用不同方法的分类结果图

表1 Indian Pines数据集采用不同方法的分类精度

图5 Indian Pines数据集参数k 和 p 对分类精度的影响

（3）说明本文算法选择BTC 作为分类器的有效性和合理性。与SVM 作为分类器的方法（记为SPCADTRF-SVM）进行比较，该方法与本文算法的前5 步完全相同，差异仅在分类器选择上，3 种分类指标结果见图6。从图6 可以发现，一定训练样本下本文算法优于SPCA-DTRF-SVM 方法，且时间复杂度大大降低（如10%的训练样本下，本文算法用时27.18 s，而SPCADTRF-SVM 用时是140.74 s）。这是因为在所提取的边缘特征经过PCA 融合之后，有效地提取了图像中的相关信息，增强了属于不同类别地物像素的光谱差异，因而增强了BTC中对各测试像素与训练像素之间相关性强弱判断的准确性，提升了其分类性能，而且BTC算法低耗时的处理过程也极大地增强了分类的准确性[11]。但对SVM 而言，融合产生波段较多，训练样本较少，导致参数估计值不够理想，因而分类精度较低[29]。

（4）说明本文算法在时间复杂度上的优越性。由表3可以看出本文算法在时间复杂度上远远低于BTC、SVM、EPF、BTC-WLS，虽然略微高于SOMP，但分类精度却大大高于其结果。此外，由于CCJSR 方法的运行时间已远超过2小时，故没有放在列表中。

表2 Indian Pines数据集DTRF参数对本文算法分类精度的影响

图6 Indian Pines数据集不同训练样本下选择不同分类器对分类的影响

表3 Indian Pines数据集采用不同方法的时间对比 s

（5）为说明所用PCA 融合算法的合理性，本文与LDA、ICA、KPCA融合方法进行比较。如表4所示，PCA融合算法在OA、AA、Kappa 和时间复杂度上都优于其他融合算法。

表4 PCA方法与其他融合方法的比较

3.2.2 Salinas数据集

对于Salinas 数据集，训练样本占5%，测试样本占95%。本文算法与另外7 种方法所得分类结果见表5。从表5 易见，本文算法在OA、AA、Kappa 系数上均是最高的。在16 个类别中，除4 个类别外，本文算法所得分类精度要么最高，要么与最高分类精度相差不大，至多低0.02%～0.80%。从整体上看，本文算法的分类精度是最优的，分类效果是最好的。和BTC、OMP 算法相比，本文算法的OA显著地增加了10%左右；与SVM算法相比，本文算法OA 增加了6.5%；与SOMP、CCJSR、EPF、BTC-WLS相比，本文算法OA高出1.1%～1.7%。其他的分类指标AA 和Kappa 有类似的结果，这说明本文算法具有优势。从图7可以看出，本文算法获得的分类图与参考图更加相符，在葡萄园和酿酒庄园表现得尤为明显。同理，在Salinas数据集继续分析参数k和p，参数以及分类器对本文算法的影响。

图7 Salinas数据集采用不同方法的分类结果图

（1）分析参数k和p对本文算法的影响。将波段子集大小k取值分成5种情况（k=5，10，15，20，25），并且在5种情况下观察参数p对分类精度的影响，具体结果见图8。从图8 可以得出以下结论：一是影响最优分类精度的最主要因素是波段子集大小。二是在波段子集大小固定情况下，随着主成分个数p的增加，分类精度均呈现出先增长接着趋于平稳的趋势。三是波段子集所含波段数在k=5 时整体分类精度是最优的。这个结论与Indian Pines 数据集是一致的。因此，实验中取k=5,p=120。

表5 Salinas数据集采用不同方法的分类精度

图8 Salinas数据集参数k 和 p 对分类精度的影响

表6 Salinas数据集DTRF参数对本文算法分类精度的影响

图9 滤波堆叠对OA、AA、Kappa的影响

（3）说明本文算法选择BTC 作为分类器的有效性和合理性。与SVM 作为分类器的方法（记为SPCADTRF-SVM）进行比较，该方法与本文算法的前5 步完全相同，差异仅在分类器选择上，3 种分类指标结果见图10。从图10可以发现，一定训练样本下，本文算法优于SPCA-DTRF-SVM，且时间复杂度大大降低（如在6%的训练样本下，本文算法用时158.48 s，SPCA-DTRFSVM用时1 103.6 s）。

图10 Salinas数据集不同训练样本下选择不同分类器对分类的影响

（4）分析本文算法的时间复杂度。由表7可知本文算法的时间复杂度远远低于BTC、SVM、OMP、SOMP、EPF、BTC-WLS。此外，由于CCJSR方法的运行时间已远远超过2小时，故没有放在列表中。

表7 Salinas数据集采用不同方法的时间对比 s

4 总结

针对高光谱遥感图像分类问题，本文提出基于SPCA和DTRF 的图像分类算法。首先对高光谱图像进行波段子集划分，并且提取每个波段子集的主成分用作该子集的特征，这样能在保留原始图像结构的基础上，减少波段数量，降低运算量。再对由这些特征组成的数据集利用多层DTRF得到堆叠的边缘保持滤波图，实现了多尺度空间特征的提取。然后对其进行PCA 特征融合，达到降低冗余信息，增强类与类之间分离性的效果，进而提高了BTC中对训练像素和测试像素之间相关性强弱判断的准确性，因此极大地提升了分类器的性能，提高了分类精度。通过对两组数据集Indian Pines和Salinas进行仿真实验，表明本文算法与其他算法相比，在OA、AA、Kappa 这3 种分类指标上有一定的优势。此外，本文算法利用了BTC快速高效、易于实现的特性，极大地提高了计算效率。不足之处在于多层域变换递归滤波的参数是通过大量实验获得的，因此如何自动地选择域变换递归滤波的参数将会是未来的研究重点。