一种基于半监督字典学习的高光谱图像分类方法

摘要

本发明公开了一种基于半监督字典学习的高光谱图像分类方法，主要解决高光谱图像高维和在相对小样本情况下分类精度低的问题。其步骤包括：将高光谱图像的像素点用光谱特征向量表示；选取标记样本集、无标记样本集和测试样本集；构造有标记样本的类标矩阵；构造无标记样本的拉普拉斯矩阵；使用交替优化策略和梯度下降法求解半监督字典学习模型；使用学习得到的字典编码有标记样本、无标记样本和测试样本；使用学习得到的稀疏编码作为特征分类高光谱图像。本发明采用半监督思想，相比监督学习方法能够获得较高的分类正确率，可用于精细农业，植被调查，军事侦察等领域。

说明书

技术领域

本发明属于图像处理技术领域，涉及半监督学习和稀疏表示方法，用于相对小样本情景下的高光谱图像分类问题。 

背景技术

高光谱遥感技术产生于20世纪80年代，结合了成像技术和光谱技术，可以在电磁波的紫外到近红外的几十至几百个很窄的连续波段上获取感兴趣地面物体的辐射特性，是对地观测的重要前沿技术。与传统光谱成像技术相比，高光谱遥感不仅波段数和谱分辨率更高，而且波段几乎是连续的，可以为每个像元差生一条连续的光谱曲线，其获取的图像包含了空间、辐射，和光谱三重信息，具有图谱合一的特点。 

目前，许多国家研制了高光谱遥感系统，比如美国国家宇航局(NASA)的AVIRIS、EO-I HYPERION、加拿大的荧光线成像光谱仪FLI、德国的ROSIS-10、ROSIS-22、澳大利亚的HyMap、加拿大ITRES公司的CASI、SASI、我国的OMIS、PHI。常用的高光谱图像数据包括由美国国家宇航局的AVIRIS获取的Indian Pines数据集、Kennedy Space Center(KSC)数据集，以及EO-I HYPERION光谱仪获取的Botswana数据集等。高光谱遥感已经广泛应用于天气预报、环境监测、灾害评估、精细农业、地质调查、军事侦察等许多领域。 

地物分类的任务是确定感兴趣地面物体的所属地物类别，它是高光谱遥感最重要的应用之一，是许多相关应用的基础。对于特定的波长，不同的物质具有不同的电磁辐射特性，高光谱能够捕获从可见光到近红外光谱的连续光谱信息，因而提供了分类不同地物的重要区分信息。高光谱遥感 为分类提供丰富信息的同时，也带来了巨大的挑战：1)很高的波段数(数十到数百)；2)相对较少的标记样本(总样本数大，而标记样本的代价高昂)。在高维和相对小样本条件下，分类器可以任意好地拟合有限的训练数据，而未必能有效地预测测试数据，也即面临着过拟合的风险。对于理想的高光谱分类算法，应该在高维和少量标记样本的条件下仍能给出较高的分类精度。 

为了处理高光谱图像分类的相对小样本问题，已经有许多方法被提出。基于核的方法，如：支持向量机(SVM)，对高维数据不太敏感，在高光谱图像分类中表现出良好的性能。半监督的分类方法，如：基于图的半监督分类算法RLS，同时考虑在标记样本上的分类误差和无标记样本上的预测类标的平滑性。稀疏表示也被成功用于高光谱图像分类，如：稀疏表示分类方法SRC。 

虽然SVM对于高维数据具有一定的鲁棒性，但是研究表明，高维仍然对SVM的分类性能产生较大影响；基于图的半监督方法，相比于有监督算法可以提高分类性能，然而它是直接处理高维的谱段特征，并不能有效的提高分类精度。本发明的半监督字典学习方法，利用高光谱数据的稀疏特性和半监督学习策略，结合少量标记样本的判别信息和大量无标记样本的结构信息改进高光谱图像分类性能。 

发明内容

本发明的目的在于同时利用少量的有标记样本以及大量的无标记样本，提出一种基于半监督字典学习的高光谱图像分类方法，在有标记样本较少的情况下提高分类性能。 

为此，本发明提供了一种基于半监督字典学习的高光谱图像分类方法， 其技术方案如下： 

一种基于半监督字典学习的高光谱图像分类方法，包括如下步骤： 

(1)输入一幅高光谱图像I，包含c类地物共n个像素点，每一个像素点为一个样本，每个样本用光谱特征向量表示，样本的特征维数为d； 

(2)从图像I中选取n_l个有标记样本构成标记样本集其中，表示标记样本集的第i个样本；标记样本集对应的类标集为 其中，是标记样本集中第i个样本的类标；选取n_u个无标记样本构成无标记样本集其中，表示无标记样本集的第i个样本；其余样本构成测试样本集其中，表示测试样本集的第i个样本，n_t表示测试样本个数；R^d表示d维向量空间； 

(3)构造无标记样本的相似度矩阵S： 

其中，S_ij表示矩阵S的第i行、第j列元素，表示样本的k近邻集合，σ为控制高斯核平滑度的参数，e^(·)为指数函数； 

(4)计算无标记样本的拉普拉斯矩阵L： 

L＝D-S 

其中，D为对角矩阵，第i个对角线元素为

(5)计算有标记样本的类标矩阵H： 

(6)交替优化关于分类器参数矩阵W和字典B的半监督字典学习目标 函数： 

$\left(\begin{array}{c}f(W,B)=\underset{W,B}{\arg \min }\underset{i=1}{\overset{n_l}{\Sigma }}{\left|\right|H_{:,i}-W^T{z}_i^{\left(l\right)}\left|\right|}^2+\lambda {\left|\right|W\left|\right|}^2+\mathrm{\mu Tr}\left({\left(W^T{Z}_U\right)}^{T}L\left(W^T{Z}_U\right)\right)\\ \underset{Z_L}{\arg \min }\underset{i=1}{\overset{n_l}{\Sigma }}{\left|\right|x_i^{\left(l\right)}-{\mathrm{Bz}}_i^{\left(l\right)}\left|\right|}^2+\gamma {\left|\right|z_i^{\left(l\right)}\left|\right|}_1\\ \underset{Z_U}{\arg \min }\underset{i=1}{\overset{n_u}{\Sigma }}{\left|\right|x_i^{\left(u\right)}-{\mathrm{Bz}}_i^{\left(u\right)}\left|\right|}^2+\gamma {\left|\right|z_i^{\left(u\right)}\left|\right|}_1\end{array}\right)$

其中，B∈R^d×r表示字典，每一列表示一个字典原子，r表示字典原子的个数； 表示标记样本集中第i个样本的稀疏编码向量，为标记样本集的稀疏编码矩阵，每一列对应一个样本的编码向量；表示无标记样本集中第i个样本的稀疏编码向量，为无标记样本集的稀疏编码矩阵，每一列对应一个样本的编码向量；H_:,i∈R^c为类标矩阵H的第i列列向量，表示标记样本集中第i个样本的类标向量；W∈R^r×c是以稀疏编码为输入的线性分类器参数矩阵；f(·)表示损失函数，Tr(·)表示迹函数，||·||₁表示向量的l₁范数，γ、μ、λ为正则项权重参数； 

(7)对测试样本进行稀疏编码： 

$\underset{Z_T}{\arg \min }\underset{i=1}{\overset{n_t}{\Sigma }}{\left|\right|x_i^{\left(t\right)-{\mathrm{Bz}}_i^{\left(t\right)}\left|\right|}}^2+\gamma {\left|\right|z_i^{\left(t\right)}\left|\right|}_1$

其中，表示测试样本集中第i个样本的稀疏编码向量，为测试样本集的稀疏编码矩阵，每一列对应一个样本的编码向量； 

(8)预测测试集样本类标

$y_i^{\left(t\right)}=\underset{j}{\arg \max }{w}_j^T{z}_i^{\left(t\right)},j=\mathrm{1,2},...,c$

其中，w_j是线性分类器参数矩阵W的第j列。 

上述步骤(6)所述的交替优化关于分类器参数W和字典B的半监督字典学习目标函数，按如下步骤实施： 

6a)从标记样本集和无标记样本集中随机选取样本构建初始字典B，采 用随机矩阵初始化线性分类器参数矩阵W； 

6b)求解如下稀疏编码问题，更新有标记样本和无标记样本的稀疏编码矩阵： 

$\left(\begin{array}{c}\underset{Z_L}{\arg \min }\underset{i=1}{\overset{n_l}{\Sigma }}{\left|\right|x_i^{\left(l\right)}-{\mathrm{Bz}}_i^{\left(l\right)}\left|\right|}^2+\gamma {\left|\right|z_i^{\left(l\right)}\left|\right|}_1\\ \underset{Z_U}{\arg \min }\underset{i=1}{\overset{n_u}{\Sigma }}{\left|\right|x_i^{\left(u\right)}-B{z}_i^{\left(u\right)}\left|\right|}^2+\gamma {\left|\right|z_i^{\left(u\right)}\left|\right|}_1\end{array}\right);$

6c)采用梯度法更新线性分类器参数矩阵W： 

W＝(Z_L(Z_L)^T+λI+μZ_UL(Z_U)^T)^-1Z_LH^T

其中，I为单位矩阵； 

6d)采用梯度下降法更新字典B： 

$\left(\begin{array}{c}\frac{\partial f(W,B)}{\partial B}=\underset{m=1}{\overset{n_l}{\Sigma }}\frac{\partial f(W,B)}{{\partial z}_m^{\left(l\right)}}\frac{\partial z_m^{\left(l\right)}}{\partial B}+\underset{m=1}{\overset{n_u}{\Sigma }}\frac{\partial f(W,B)\partial z_m^{\left(u\right)}}{\partial B}\\ \frac{\partial {\tilde{z}}_m^{\left(l\right)}}{\partial B_{\mathrm{ij}}}={\left({\tilde{B}}_m^{\left(l\right)T}{\tilde{B}}_m^{\left(l\right)}\right)}^{-1}(\frac{\partial {\tilde{B}}_m^{\left(l\right)}{x}_m^{\left(l\right)}}{\partial B_{\mathrm{ij}}}-\frac{\partial {\tilde{B}}_m^{\left(l\right)T}{\tilde{B}}_m^{\left(l\right)}}{\partial B_{\mathrm{ij}}}{\tilde{z}}_m^{\left(l\right)}),i=\mathrm{1,2},...,d;j=\mathrm{1,2},...,r\\ \frac{\partial {\tilde{z}}_m^{\left(u\right)}}{\partial B_{\mathrm{ij}}}={\left({\tilde{B}}_m^{\left(u\right)T}{\tilde{B}}_m^{\left(u\right)}\right)}^{-1}(\frac{\partial {\tilde{B}}_m^{\left(u\right)T}{x}_m^{\left(u\right)}}{\partial B_{\mathrm{ij}}}-\frac{\partial {\tilde{B}}_m^{\left(u\right)T}{\tilde{B}}_m^{\left(u\right)}}{\partial B_{\mathrm{ij}}}{\tilde{z}}_m^{\left(u\right)}),i=\mathrm{1,2},...,d;j=\mathrm{1,2},...,r\\ B=B-\xi \frac{\partial f(W,B)}{\partial B}\end{array}\right)$

其中，和分别是稀疏编码向量和的非零元素构成的向量，和 分别是与和对应选择的原子构成的子字典，B_ij表示字典B的第i行第j列元素，表示目标函数关于变量的导数，0≤ξ≤1是优化步长因子； 

6e)执行步骤6b)-6d)，直至满足最大迭代次数，然后停止。 

本发明的有益效果：本发明结合字典学习和半监督学习策略，利用少量有标记样本提供的判别信息，和大量无标记样本包含的几何结构信息，学习既具有判别性又具有良好的泛化能力的字典，然后利用学习到的字典对样本稀疏编码，用得到的稀疏编码系数作为输入特征构造分类器对高光谱图像进 行分类。本发明与现有的技术相比具有以下优点： 

1、本发明采用稀疏表示方法，能够较好的处理高光谱数据的高维问题。 

2、本发明同时利用了有标记样本和无标记样本，能够充分利用标记样本的判别信息和无标记样本的结构信息，提高分类精度。 

对比实验表明，本发明能较好地处理高光谱数据的高维和小样本问题，提高了高光谱遥感图像的分类准确率。 

附图说明

图1是本发明的流程示意图； 

图2是本发明仿真采用的Indian Pines图像； 

图3是本发明与现有方法对Indian Pines图像的分类结果图，图3(a)-3(d)分别为RLS、SRC、SVM和本发明方法的分类结果图。 

具体实施方式

参照图1，本发明的具体实施步骤包括： 

步骤1，输入一幅高光谱图像I，包含c类地物共n个像素点，每个像素点作为一个样本，每个样本用光谱特征向量表示，样本的特征维数为d。 

步骤2，构建标记样本集X_L、类标集Y_L、无标记样本集X_U和测试样本集X_T。 

2a)从每一类地物像素点中随机选取等量的样本作为有标记样本，共nl个有标记样本，构成标记样本集其中，表示标记样本集的第i个样本；标记样本集对应的类标集为其中，是标记样本集中第i个样本的类标；R^d表示d维向量空间； 

2b)从标记样本集以外的样本中随机选取n_u个样本作为无标记样本，构成无标记样本集其中，表示无标记样本集的第i个 样本； 

2c)将标记样本集和无标记样本集以外的样本作为测试样本，构成测试样本集其中，表示测试样本集的第i个样本，n_t表示测试样本个数。 

步骤3，构造无标记样本的相似度矩阵S： 

其中，S_ij表示矩阵S的第i行、第j列元素，表示样本的k近邻集合，σ为控制高斯核平滑度的参数，e^(·)为指数函数。 

步骤4，计算无标记样本的拉普拉斯矩阵L： 

L＝D-S 

其中，D为对角矩阵，第i个对角线元素为

步骤5，计算有标记样本的类标矩阵H： 

步骤6，交替优化关于分类器参数W和字典B的半监督字典学习目标函数： 

$\left(\begin{array}{c}f(W,B)=\underset{W,B}{\arg \min }\underset{i=1}{\overset{n_l}{\Sigma }}{\left|\right|H_{:,i}-W^T{z}_i^{\left(l\right)}\left|\right|}^2+\lambda {\left|\right|W\left|\right|}^2+\mathrm{\mu Tr}\left({\left(W^T{Z}_U\right)}^{T}L\left(W^T{Z}_U\right)\right)\\ \underset{Z_L}{\arg \min }\underset{i=1}{\overset{n_l}{\Sigma }}{\left|\right|x_i^{\left(l\right)}-{\mathrm{Bz}}_i^{\left(l\right)}\left|\right|}^2+\gamma {\left|\right|z_i^{\left(l\right)}\left|\right|}_1\\ \underset{Z_U}{\arg \min }\underset{i=1}{\overset{n_u}{\Sigma }}{\left|\right|x_i^{\left(u\right)}-{\mathrm{Bz}}_i^{\left(u\right)}\left|\right|}^2+\gamma {\left|\right|z_i^{\left(u\right)}\left|\right|}_1\end{array}\right)$

其中，B∈R^d×r表示字典，每一列表示一个字典原子，r表示字典原子的个数； 表示标记样本集中第i个样本的稀疏编码向量，为标记样本集的稀疏编码矩阵，每一列对应一个样本的编码向量；表示无标记 样本集中第i个样本的稀疏编码向量，为无标记样本集的稀疏编码矩阵，每一列对应一个样本的编码向量；H_:,i∈R^c为类标矩阵H的第i列列向量，表示标记样本集中第i个样本的类标向量；W∈R^r×c是以稀疏编码为输入的线性分类器参数矩阵；f(·)表示损失函数，Tr(·)表示迹函数，||·||₁表示向量的l₁范数，γ、μ、λ为正则项权重参数。 

6a)从标记样本集和无标记样本集中随机选取样本构建初始字典B，采用随机矩阵初始化线性分类器参数矩阵W； 

6b)求解如下稀疏编码问题，更新有标记样本和无标记样本的稀疏编码矩阵： 

$\left(\begin{array}{c}\underset{Z_L}{\arg \min }\underset{i=1}{\overset{n_l}{\Sigma }}{\left|\right|x_i^{\left(l\right)}-{\mathrm{Bz}}_i^{\left(l\right)}\left|\right|}^2+\gamma {\left|\right|z_i^{\left(l\right)}\left|\right|}_1\\ \underset{Z_U}{\arg \min }\underset{i=1}{\overset{n_u}{\Sigma }}{\left|\right|x_i^{\left(u\right)}-B{z}_i^{\left(u\right)}\left|\right|}^2+\gamma {\left|\right|z_i^{\left(u\right)}\left|\right|}_1\end{array}\right);$

6c)采用梯度法更新线性分类器参数矩阵W： 

W＝(Z_L(Z_L)^T+λI+μZ_UL(Z_U)T)^-1Z_LH^T

其中，I为单位矩阵； 

6d)采用梯度下降法更新字典B： 

$\left(\begin{array}{c}\frac{\partial f(W,B)}{\partial B}=\underset{m=1}{\overset{n_l}{\Sigma }}\frac{\partial f(W,B)}{{\partial z}_m^{\left(l\right)}}\frac{\partial z_m^{\left(l\right)}}{\partial B}+\underset{m=1}{\overset{n_u}{\Sigma }}\frac{\partial f(W,B)\partial z_m^{\left(u\right)}}{\partial B}\\ \frac{\partial {\tilde{z}}_m^{\left(l\right)}}{\partial B_{\mathrm{ij}}}={\left({\tilde{B}}_m^{\left(l\right)T}{\tilde{B}}_m^{\left(l\right)}\right)}^{-1}(\frac{\partial {\tilde{B}}_m^{\left(l\right)}{x}_m^{\left(l\right)}}{\partial B_{\mathrm{ij}}}-\frac{\partial {\tilde{B}}_m^{\left(l\right)T}{\tilde{B}}_m^{\left(l\right)}}{\partial B_{\mathrm{ij}}}{\tilde{z}}_m^{\left(l\right)}),i=\mathrm{1,2},...,d;j=\mathrm{1,2},...,r\\ \frac{\partial {\tilde{z}}_m^{\left(u\right)}}{\partial B_{\mathrm{ij}}}={\left({\tilde{B}}_m^{\left(u\right)T}{\tilde{B}}_m^{\left(u\right)}\right)}^{-1}(\frac{\partial {\tilde{B}}_m^{\left(u\right)T}{x}_m^{\left(u\right)}}{\partial B_{\mathrm{ij}}}-\frac{\partial {\tilde{B}}_m^{\left(u\right)T}{\tilde{B}}_m^{\left(u\right)}}{\partial B_{\mathrm{ij}}}{\tilde{z}}_m^{\left(u\right)}),i=\mathrm{1,2},...,d;j=\mathrm{1,2},...,r\\ B=B-\xi \frac{\partial f(W,B)}{\partial B}\end{array}\right)$

其中，和分别是稀疏编码向量和的非零元素构成的向量，和 分别是与和对应选择的原子构成的子字典，B_ij表示字典B的第i行 第j列元素，表示目标函数关于变量的导数，0≤ξ≤1是优化步长因子； 

6e)执行步骤6b)-6d)，直至满足迭代停止条件，即：最大迭代次数Iter。 

步骤7，对测试样本进行稀疏编码： 

$\underset{Z_T}{\arg \min }\underset{i=1}{\overset{n_t}{\Sigma }}{\left|\right|x_i^{\left(t\right)-{\mathrm{Bz}}_i^{\left(t\right)}\left|\right|}}^2+\gamma {\left|\right|z_i^{\left(t\right)}\left|\right|}_1$

其中，表示测试样本集中第i个样本的稀疏编码向量，为测试样本集的稀疏编码矩阵，每一列对应一个样本的编码向量。 

步骤8，预测测试集样本类标

$y_i^{\left(t\right)}=\underset{j}{\arg \max }{w}_j^T{z}_i^{\left(t\right)},j=\mathrm{1,2},...,c$

其中，w_j是线性分类器参数矩阵W的第j列。 

本发明的效果可以通过以下仿真实验说明： 

1.仿真条件: 

仿真实验采用美国国家宇航局(NASA)的AVIRIS于1992年6月在印第安纳西北部获取的Indian Pines图像，如图2所示，图像大小为145x145，共16类地物，220个波段，移除被水域吸收的20个波段，实验中移除样本数较少的7类，只考虑表1所示的9类数据。 

仿真实验在Intel Core(TM)2Duo CPU、主频2.33GHz，内存2G，MATLAB 7.14上进行。 

表1 Indian Pines图像中的9类数据 

类别编号 类别名称 样本个数 1 Corn-notill 1434

[0086] 

2 Corn-min 834 3 Grass/Pasture 497 4 Grass/Trees 747 5 Hay-windrowed 489 6 Soybeans-notill 968 7 Soybeans-min 2468 8 Soybean-clean 614 9 Woods 1294

2.仿真内容及分析： 

使用本发明的半监督字典学习方法SSDL与现有三种方法RLS、SRC、SVM对Indian Pines高光谱图像进行分类。RLS的正则化参数均设置为最优参数，SRC的稀疏度参数设置为0.1，分类方法SVM的核参数和惩罚因子通过5倍交叉验证寻找最优参数，本发明方法的稀疏度参数γ设置为0.1，正则参数λ、μ分别设置为0.3和0.5，无标记样本的k近邻设置为5近邻，算法最大迭代次数Iter设置为20。 

从表1所示的9类数据中每类选取固定数目的像素点作为有标记样本，剩余样本中每类选取20％的样本作为无标记样本，其余样本用于测试。用本发明方法与现有三种方法对9类数据进行20次实验，取平均结果作为最终分类精度。 

表2给出了四种方法在每类有标记样本个数为2、5、8、12、15、20时的分类精度，可以看出本发明获得了比其它三种方法更高的精度，尤其在有标记样本数很少的情况下具有较为明显的优势。 

表2 Indian Pines 9类地物分类精度 

样本数 分类精度(％)

[0093] 

   LapRLS SRC SVM SSDL 2 40.08±5.21 49.83±3.87 40.19±4.66 53.44±3.11 5 46.67±4.28 57.97±2.90 52.61±3.34 58.02±2.73 8 55.09±2.89 65.62±2.66 63.85±2.56 69.92±2.06 12 57.16±2.96 69.56±1.96 66.76±2.43 71.34±1.75 15 59.81±2.58 70.65±1.88 68.35±2.24 74.31±1.45 20 62.48±3.02 73.57±1.56 74.36±1.47 75.18±1.28

图3给出了四种方法在每类选取20个标记样本时的分类结果图。图3(a)-3(d)分别为RLS、SRC、SVM、SSDL的分类结果，可以看出本发明的分类结果图最为清晰，区域致性比现有方法更好。 

综上，本发明方法基于半监督字典学习，充分利用有标记样本的判别信息与无标记的结构信息，能够较好地处理高光谱数据高维和小样本的问题，与现有的方法相比具有一定的优势。 

以上例举仅仅是对本发明的举例说明，并不构成对本发明的保护范围的限制，凡是与本发明相同或相似的设计均属于本发明的保护范围之内。