2021574
45号钢
高光谱遥感图像是由几十个到数百个连续波段图像组成的三维数据立方体，其不仅反映了地物的光谱信息，而且包含了地物空间分布信息，因而具有图谱合一的特点，具有较强的地物分类和监测能力[1]，被广泛地应用在地质勘探、考古发现、森林火灾检测、军事侦察、环境检测和目标识别等领域[2-5]。像元分类是高光谱遥感图像研究的重要内容之一。一方面，与全和多光谱遥感图像相比，高光谱图像有更高的光谱分辨率和相对较低的空间分辨率，单纯依靠光谱信息的高光谱分类方法，往往不能达到用户对地物分类不断提高的类别精细程度和精度的要求，越来越多的学者开始尝试综合光谱与空间信息进行高光谱分类研究[2]。另一方面，由于高光谱遥感图像维数高，训练样本数目往往相对不足，无法有效统计分类模型的相关参数，传统分类方法运算速度慢，分类精度低，在分类过程中可能会导致Houghes 现象[6]。为了提高分类精度，一些新的分类方法被提出，如支持向量机[7]、决策树[8]、神经网络[9]、稀疏表示[10]等分类方法。Toksoz等人在2016年基于稀疏表示理论设计出一个新的高光谱图像分类器，即基本阈值分类器（Basic Thresholding Classifier，BTC）[11]，其具有分类精
于多1，黄永东1，2
1.北方民族大学图像处理与理解研究所，银川750021
2.大连民族大学数学与信息科学研究中心，辽宁大连116600
摘要：提出一种基于分割的主成分分析（Segmented Principal Component Analysis，SPCA）和域变换递归滤波（Domain Transform Recursive Filtering，DTRF）的高光谱图像分类算法。利用SPCA方法降低高光谱图像的维数和提取各波段子集的第一主成分。使用不同参数的域变换递归滤波器对各波段子集第一主成分进行滤波，形成堆叠的边缘保持滤波图。采用主成分分析（Principal Component Analysis，PCA）将堆叠的边缘保持滤波图进行特征融合。利用基本阈值分类器（Basic Thresholding Classifier，BTC）对融合后的主成分进行分类。仿真实验表明，所提方法能够提高分类精度，且在总体分类精度、平均分类精度、Kappa系数等方面优于已有方法。 关键词：主成分分析（PCA）；分割的主成分分析（SPCA）；域变换递归滤波（DTRF）；高光谱图像分类；基本阈值分类器文献标志码：A中图分类号：TP753doi：10.3778/j.issn.1002-8331.1912-0240
Hyperspectral Image Classification Based on SPCA and Domain Transform Recursive Filtering YU Duo1,HUANG Yongdong1，2
1.Institute of Image Processing and Understanding,North Minzu University,Yinchuan750021,China
2.Center for Mathematics and Information Science,Dalian Minzu University,Dalian,Liaoning116600,China
Abstract：A new hyperspectral image classification method is proposed based on Segmented Principal Component Analysis （SPCA）and Domain Transform Recursive Filtering（DTRF）.First,the SPCA method is used to reduce the dimension of hyperspectral image and extract the first principal component of each band subset.Then,DTRF with different parameters deals with the first principal component of each band subset to form a stacked edge-preserving filter map.And the Princi-pal Component Analysis（PCA）is used to fuse the features of the stack-preserving filter map.At last,the Basic Thresholding Classifier（BTC）classifies the fused principal components.Simulation experiments show that the proposed method can improve the classification accuracy,and the overall accuracy,average accuracy and Kappa coefficient are higher than the some existing methods.
Key words：Principal Component Analysis（PCA）;Segmented Principal Component Analysis（SPCA）;Domain Trans-form Recursive Filtering（DTRF）;hyperspectral image classification;basic thresholding classifier
基金项目：国家自然科学基金（11761001，11761003）；宁夏科技创新领军人才项目（KJT2016002）。
作者简介：于多（1991—），男，硕士研究生，CCF学生会员，研究领域为高光谱图像分类，E-mail：；黄永东（1974—），通信作者，男，博士，教授，研究领域为模式识别与图像处理。
收稿日期：2019-12-16修回日期：2020-04-22文章编号：1002-8331（2021）04-0199-10
199
2021574
度高、易实现、无需调节参数、快速有效且能够联合空间信息等特点。本文选择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 Transfor
m Recursive Filtering，DTRF）的高光谱图像分类算法。主要步骤如下：首先，使用SPCA降低数据维数，即将高光谱图像的波段分成一些相邻波段的子集，PCA作用在每一个子集上；然后，利用DTRF获得高光谱图像的堆叠的边缘保持滤波图；接着，PCA作用在堆叠的边缘保持滤波图上进行特征融合；最后，用BTC对上述融合图进行分类。
1理论简介
1.1域变换递归滤波
域变换递归滤波（DTRF）[22]是一种实时的边缘保留滤波，其在改善高光谱图像分类性能方面是非常有效的。滤波过程可分为两步：第一步进行域变换；第二步是递归滤波。在实际滤波过程中，为了控制滤波器的尺寸和模糊度，常常将域变换定义为应用近似距离变换。即对于给定的一维信号I，域变换定义为：
U i=I0+∑
j=1
i
(1+
δs
δr|
|I j-I j-1)（1）这里，U i是域变换信号，δs和δr是用来调节滤波器平滑程度的两个参数。然后，输入信号通过递归滤波器进行处理，即：
J i=(1-a b)I i+a b J i-1（2）其中，J i是第i个像元的滤波输出，a∈[0,1)代表反馈系数，b表示变换域中两个相邻样本的距离。当图像接近边缘时，b不断增大，a b趋于0，迭代运算终止，从而达到边缘保留的作用。
对二维信号，域变换递归滤波通过沿着图像的每一个维度执行上述所提的一维操作，进而处理图像。本文用DTRF(I,δs,δr)表示域变换递归滤波。
1.2BTC算法
BTC是Toksoz等人基于稀疏表示理论提出的一种高光谱图像分类算法，其具有分类精度高、易实现、无需调节参数、快速有效且能够联合空间信息等特点[11]。它是一个逐像素的分类器，仅仅利用光谱特征来分类每一个像素。最终对每一个像素，产生一个带有类别标签的残差向量。为了改善BTC分类精度，将算法扩展到空谱邻域，首先利用SPCA提取光谱特征，然后利用多层域变换递归滤波提取边
缘特征，再利用PCA融合这些边缘特征，最后利用BTC对融合后的光谱波段（记作H∈ℝn1×n2×B，共有n1×n2个像素，B为波段数）进行逐像素的分类，产生带有类标签的残差向量ε∈ℝC（C是类别数），残差向量形成一个堆叠的图像ℝn1×n2×C，残差图表示为map i(map i∈ℝn1×n2，其中i∈{1,2,⋯,C})，最后基于最小残差来决定每一个像素的类标签。具体BTC算法大致如下：
第一步是计算测试像素y和字典A中的每一个像素的相关性。第二步是计算索引集Λ，由与y的相关性最大的前M个像素所对应的索引组成。第三步是构造矩阵D，矩阵D是由索引在Λ中的训练像素组成。第四步在求解稀疏编码值x(Λ)中，利用了Tikhonov正则化，因为矩阵D所包含的像素可能属于同一类，所以D的列变得线性相关，在算法中的逆运算将无法进行，而Tikhonov正则化次定方程组是一个有效消除这种限
制的方法[23]。最后基于最小残差arg min
j
ε(j)来决定测试像素的类别，其中残差ε(j)=
y-A j x j
2
,x j对应于x第
j类稀疏编码部分。BTC算法详细过程可参见文献[11]。
算法1BTC算法
输入：HSI dictionary A∈ℝB×N；
200
2021574Te s t pixel y ∈ℝB ；Threshold M ∈ℕ；
Regularization constant α∈(0,1)；Initial sparse coefficients with zeros X ∈ℝN
过程：
1.Class of ν←A T y ；
2.Λ←L M (ν)；
3.D ←A (Λ)；
4.x
(Λ)⇐(D T
D +αI )-1
D T
y ；5.ε(j )⇐ y -A j x
j 2
,∀j ∈1,2,⋯；6.class (y )←arg min j
ε(j )礼仪的功能
输出：
class o f y ；Residual vector ε∈ℝC
1.3PCA 融合
选择PCA 作为融合工具，是因为它易于实现，没有
参数，可以有效提取高光谱图像中的相关信息[24]。具体原因如下：首先，由多参数的域变换递归滤波获得的滤波图不仅仅包含了互补信息，而且还有大量的冗余信息，而PCA 可以很好地去除这些冗余信息。再者，滤波器的平滑操作尽管可以有效移除噪声和小尺度的细节，但是也降低了属于不同类地物的像素的差异值（见图1）。为解决这些问题，PCA 是一个很好的工具，它能够提取在滤波图中最相关的信息，也因此能够有效增强属于不同类地物像素的光谱差异[25]，这提高了下一步BTC 分类中判断测试像素与训练像素之间的相关性强弱的准确性。
2基于SPCA-DTRF 高光谱遥感图像分类
本文提出一种基于SPCA-DTRF 的高光谱遥感图像
分类算法，主要由SPCA 方法来降维和特征提取，并利用不同参数设置的DTRF 滤波器得到堆叠的边缘保持滤波图；接着对上述滤波图进行PCA 特征融合，实现光谱信息和空间信息的有效利用；最后利用BTC 进行分类。具体步骤如下：
步骤1对高光谱图像进行波段子集划分。把高光谱遥感图像的光谱波段划分为L 个相邻波段的子集，每个子集大小（子集含有的波段数）记作k 。第ℓ个子集(ℓ≤L )可以表示为：
S ℓ
=ìíî
(X (ℓ-1)k +1,⋯,X ℓk ),ℓk ≤M
(X (ℓ-1)k +1,⋯,X M
),ℓk >M （3）实验中通过设定子集的大小k ，进而确定波段数
L ，
即当k 整除M 时，L =M
k
；当k 不整除M 时，L =éëêùû
úM k +1。其中X ={X 1,X 2,⋯,X M }表示高光谱图像包含M 个光谱波段，éëêùû
úM k 表示取整。步骤2提取每一个波段子集的第一主成分。对于高光谱遥感图像，它的相邻光谱波段常常有很强的相关性，同时也包含了冗余信息[26]。为了降低后续特征提取与融合和BTC 分类的计算复杂度，更好地保留光谱差异信息，将PCA 作用于波段子集S ℓ上，并提取第一主成分，即有：
Y ℓ=PCA (S ℓ,1),ℓ=1,2,⋯,L
（4）
步骤3把不同波段子集的第一主成分堆叠在一起，得到高光谱数据，即：
Y ={Y 1,Y 2,⋯,Y L }
（5）
其中，
Y 表示维数降低后的高光谱数据。步骤4对各波段子集提取的第一主成分进行域变换递归滤波处理，得到堆叠的边缘保持滤波图。对Y ℓ用多层（多参数）域变换递归滤波图进行特征提取，并将它们叠加在一起，即：
F x ℓ=DTRF (Y ℓ,δx s ,δx
r ),x =1,2,⋯,w,ℓ=1,2,⋯,L （6）F x ={F x 1,F x 2,⋯,F x L },x =1,2,⋯,w
（7）F ={F 1,F 2,⋯,F w }
（8）
其中，δx s 和δx
r 是域变换递归滤波器的第x 个参数，δx s 控制滤波窗口大小，δx r 控制滤波器的平滑程度。具体
来讲，对维数降低后的第ℓ个波段，通过设置不同的参
数(δ1s ,δ1r ),(δ2s ,δ2r ),⋯,(δw s ,δw r )来获得不同的边缘保持滤
波器进而实施多层边缘保持滤波。多层边缘保持滤波的主要目的是去捕捉由边缘保持滤波器获得的多尺度结构特征，如图1所示。通过调整域变换递归滤波器的参数，滤波后的图像在保留边缘或者边界方面可能看起来有很大差别，但是有一定的互补性。一方面，对于有噪声的区域，平滑滤波能够有效地降低噪声，同时增加像素的光谱纯度。然而，当模糊程度相对大的时候，属于不同类别地物像素的光谱可分离性可能降低。这些意味着，在表达这些目标或不同尺度的特征方面，不同模糊程度的滤波所得到的边缘保持图像有各自不同的优势。因此通过整合堆叠的边缘保持滤波图中的这些互补信息，
预期在改善分类精度方面起到很好的作用。
（a ）输入高光谱波段（b ）(δr =0.2,δs =6)（c ）(δr =0.2,δs =15)
（d ）(δr =0.2,δs =40)（e ）(δr =12,δs =50)（f ）(δr =80,δs =100)
图1
参数δr 和δs 对滤波效果的影响
妈妈的奖励201
2021574步骤5对堆叠的边缘保持滤波图用PCA 进行特征融合。如果直接将堆叠的边缘保持滤波图输入到分类器，尽管可利用它们之间的互补信息，但是可能导致Houghes 现象和额外计算负担。另一方面，尽管DTRF 滤波在平滑图像的同时也能保留边缘，但是平滑过程中不可避免地降低了不同类地物像素的差异性，从而影响分类精度。为了解决上述问题，本文利用PCA 进行特征融合，具体操作如下：
L =PCA (F,p )
（9）
这里，p 是保留在L 中的主成分数量，F 是堆叠的边缘保持滤波图。换句话说，利用PCA 提取堆叠的边缘保持滤波图中前p 个主成分。
步骤6用BTC 来确定每个像素的类别，像素值最小的分类残差图对应的类别即为该像素类别。
本文所提算法流程图见图2所示。
3实验仿真
为了对本文所提算法的分类效果进行评估，使用了
两组众所周知的高光谱遥感图像数据集Indian Pines 和Salinas 进行实验。
3.1实验数据
（1）Indian Pines 数据集是由机载可见红外光成像
光谱仪（AVIRIS ）获得的，是在印第安那州西北部的印度松树测试场采集到的，由波长范围是0.4~2.5nm 的光谱波段组成，图像大小为145×145像素，空间分辨率为20m ，包含16个类别，去除水
汽吸收波段后，从204个波段中选取200个波段作为研究对象。该数据的伪彩图、参考图和类别图见图3。
（2）Salinas 数据集是由AVIRIS 在加利福尼亚州的农业区域Salinas Vally 采集到的，由波长范围是0.4~2.5nm 的224个光谱波段组成，图像大小为512×217像素，空间分辨率是3.7m ，包含16个类别，移除水汽吸收严重的波
段，从224个波段中选200个波段作为研究对象。该数据的伪彩图、参考图和类别图见图3。
3.2实验分析
为了验证算法的有效性，将本文算法（SPCA-DTRF ）
与基本阈值分类方法（BTC ）[11]
、支持向量机方法（Sup-port Vector Machine ，SVM ）[7]、正交匹配追踪法（Orthog-onal 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 所用分类器是支持向量机。采纳这些方法，是因为在高光谱遥感领域它们不仅被大量地引用，而且在一些真实的高光谱数据中有很好的表现。支持向量机是通过五层交叉验证的高斯核来实现的，其他
第1波段子集
第2波段子集第k 波段子集
提取第1主成分
提取第1主成分提取第1主成分
︙︙堆叠的主成分堆叠的滤波图融合后的滤波图
波段子集划分
P C A 融合
B T
C 分类
不同参数设置的域变换递归滤波
图2
算法流程
Corn-no till
Corn-min till Corn
Soybeans-no till Soybeans-min till Soybeans-clean till Alfalfa
Grass/pasture
Grass/trees
Grass/pasture-mowed Hay-windrowed Oats Wheat Woods
Bldg-Grass-Tree-Drives Stone-steel
towers
Brocoli_green_weeds 1
Brocoli_green_weeds 2Fallow
Fallow_rough_plow Fallow_smooth Stubble Celery
Grapes_untrained Soil_vinyard_develop Corn_senesced_weeds Lettuce_romaine_4wk Lettuce_romaine_5wk Lettuce_romaine_6wk Lettuce_romaine_7wk Vinyard_untrained
Vinyard_treils
Indian Pines
Salinas
（a ）伪彩图（b ）地面参考数据
（c ）颜编码
图3
Indian Pines 和Salinas 实验数据
202
2021574
方法均采用相应文献的默认参数来实现。为使算法具有可比性，实验中各类样本被随机分成训练样本
和测试样本，所有实验都按这种方法进行10次，并对分类结果取平均值。在实验中，本文采用3种常用分类指标，即总体分类精度（OA）、平均分类精度（AA）和Kappa系数（Kappa）。3种指标所得数值越大，表示分类效果越好。所有的实验在Matlab2017a上进行，电脑配置是3.5GHz CPU，8GB安装内存。
3.2.1Indian 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，参数(δx s,δx r)和分类器选取对本文算法的影响。
（1）分析波段子集大小k和融合后的主成分个数p 对本文算法SPCA-DTRF的影响。为了便于分析，将每个波段子集所含波段数k取值分成5种情况（k=5，10，15，20，25），并且在5种情况下观察参数对分类精度的影响，具体结果见图5。从图5可以得到如下结论：
一是（a）BTC（b）SVM（c）OMP（d）SOMP
（e）EPF（f）CCJSR（g）BTC-WLS（h）SPCA-DTRF
图4Indian Pines数据集采用不同方法的分类结果图
类别
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
训练样本
10
143
83
24
49
73
10
48
10
98
246
60
21
127
39
10
OA/%
AA/%
Kappa/%
测试样本
36
1285
747
213
434
657
18
430
10
874
2209
533
184
1138
347
83
分类精度/%
BTC
85.18
59.42
54.01
36.66
83.73
87.64
74.44
82.11
88.00
60.21
61.76
51.10
95.27
87.69
53.11
59.03
66.97
69.93
62.39
SVM
76.09
79.65
80.04
67.77
89.90
91.42
87.37
97.56
50.89
77.09
81.70
77.71
93.42
87.64
73.72
98.48
83.48
82.16
81.09
OMP
67.53
63.77
58.13
42.02
65.07
64.33
20.44
97.07
59.43
73.30
63.23
64.72naieve
87.64
35.44
70.69
67.05
58.30
62.51
SOMP
93.62
94.02
94.57
96.10
94.44
95.53
77.12
99.94
37.20
93.97
96.44
94.26
97.98
97.25
97.57
92.57
95.39
90.79
94.81
CCJSR
98.51
93.66
94.90
91.79
94.18
95.52
74.16
100.00
84.62
93.59
97.33
89.96
99.55
99.32
95.18
95.86
95.62
93.63
95.01
EPF
99.34
97.56
刘震云单位98.13
92.89
98.13
96.31
100.00
99.77
100.00
93.85
92.66
95.41
100.00
93.88
可发性聚苯乙烯93.41
93.28
95.16
96.54
94.46
BTC-WLS
99.38
93.35
91.58
93.24
94.72
100.00
97.22
100.00
100.00
90.52
99.24
97.59
100.00
99.87
89.42
97.23
95.16
96.46
95.78
本文算法
99.06（3）
98.27
98.25
99.11
98.63
99.30（2）
96.94（3）
99.97（3）
99.09（3）
98.19
99.29
98.41
100.00
99.78（2）
99.13
98.52
98.96
98.88
98.81
表1Indian Pines数据集采用不同方法的分类精度
注：最优结果用粗体表示，括号内的数字表示该类别本文算法分类精度的排名。
203
