基于张量半监督标度切维数约减的高光谱图像分类方法

技术领域

本发明属于图像处理技术领域，涉及高维数据的降维，用于对高光谱遥感图像的分类。

背景技术

近几年新兴的高光谱遥感图像处理是遥感的前沿技术。高光谱遥感图像处理技术利用成像 光谱仪以纳米级的光谱分辨率，以几十或几百个波段同时对地表物成像，从而获得地物的连 续光谱信息，使得同步获取地物的空间信息、辐射信息、光谱信息成为可能，其最显著的特 性就是“图谱合一”，这些技术的发展使人类在对地观测和信息获取能力方面取得了质的飞 跃。

研究中常用的高光谱图像数据有美国宇航局NASA喷气推进实验室的空载可见光/红外成 像光谱仪AVIRIS获得的Indian Pines数据集和ROSIS光谱仪获得的University of Pavia数 据集。

高光谱图像分类就是结合了高光谱图像的特点将高光谱图像中的每个像素分到对应的 类别中去的过程。虽然高光谱遥感技术获得的图像包含的丰富的谱信息、空间信息、辐射信 息使分类成为更容易的事情，但是高光谱图像的分类仍然面对着巨大的困难和挑战：(1)数 据量大，存储、显示等面临着困难，波段数较高，计算量大；(2)维数灾难,过高的维数带来 的冗余信息会使分类精度降低；(3)波段数多，相关性高，从而对训练样本数的需求增大， 因为没有足够的样本，在一定程度上会降低分类模型参数的可靠性。因此，要想得到较高的 分类精度，就要对高光谱的维数进行约减，减少数据量。

目前已经存在很多的维数约简方法，根据对有标记样本的利用情况可以分为有监督维数约 减方法、无监督维数约简方法和半监督维数约简方法，例如主成分分析PCA和局部线性嵌入 LLE方法属于无监督方法，线性判别分析LDA属于有监督方法，半监督判别分析方法SDA属于 半监督方法。有监督维数约减方法就是利用有标记数据进行训练，根据类别信息求得低维空 间。与有监督相反，无监督维数约减方法没有类标信息，而是通过发现数据的内在结构特征， 选择最能体现数据结构的低维特征向量。半监督方法结合了有监督和无监督这两种方法的特 点，考虑类别信息的同时挖掘数据本身的结构信息，将资源的利用最大化，从而得到更好的 低维空间。

而这些基于向量的方法需要将图像向量化，因此，他们仅仅依赖了谱特性而忽略了空间分 布。为了解决这些缺点，学术界引入一种基于张量的高光谱图像表示，将空间和谱间结构同 时进行分析，得到了良好的结果。

Dacheng Tao等人在论文“General Tensor Discriminant Analysis and Gabor Features  for Gait Recognition”(PAMI 2007)中提出了张量LDA的方法用于步态识别，该方法主要 是将LDA的方法推广到张量的计算上。该方法虽然利用了有标签的信息，但是其不能处理异 方差和多态数据。

发明内容

本发明的目的是提出一种新的基于张量半监督标度切维数约减的高光谱图像分类方法， 以利用少量的有标记样本以及大量的无标记数据，更好地挖掘数据的内在结构，解决现有技 术需要将图像向量化、空间信息丢失和不能处理异方差和多态数据的问题，提高分类精度。

本发明的技术思路是：通过对图像中的特征进行学习，恰当的表征出不同数据之间的共性 与异性，利用张量相关计算方法实现对高光谱图像的降维，克服“维数灾难”，通过发现嵌 入在高维数据中的有价值的本征低维结构信息，寻求事物的本质规律，通过映射投影到一个 低维特征空间，获得了原始数据更精简的表示。其实现步骤包括如下：

(1)输入高光谱数据集A∈R^m×n×D，该数据集包含c类地物，其中m×n表示图像空间大小 即像素点的个数，D表示数据集的总波段数；

(2)以A的每个像素点为中心取5*5的邻域块，得到Q个具有全波段的子数据立方体， 每个子数据立方体作为一个样本用一个三阶张量来表示，得到样本集合其中 c_a表示第a个样本，Q表示样本总个数，Q＝m×n；

(3)从样本集合中随机选择N个有标记样本构成有标记训练样本集对 应的类标向量记为：其余Q-N个无标记样本构成测试样本集 ,y_u∈R^5×5×D，其中，χ_t表示有标记训练集的第t个样本，l_t表示第t个有标记训练 样本所属的类别标号，y_u表示测试集的第u个样本；

(4)从Q-N个无标记样本中选出h个无标记样本，与N个有标记样本一起构成总训练 集其中，s_k表示总训练集的第k个训练样本，N+η为总训练集的样本个数， 1≤η≤Q-N；

(5)构造有标记训练集X的类间不相似性矩阵B：

$B=\underset{p=1}{\overset{c}{\Sigma }}\underset{i\in V_p}{\Sigma }\underset{j\in V_P^{\prime }}{\Sigma }\frac{1}{n_p{n}_{c\left(j\right)}}{\mathrm{mat}}_3\left(\left(({\chi }_i-{\chi }_j){\times }_1{U}_1^T{\times }_2{U}_2^T\right)\right){\mathrm{mat}}_3^T\left(\left(({\chi }_i-{\chi }_j){\times }_1{U}_1^T{\times }_2{U}_2^T\right)\right)$

其中，V_p表示第p类有标记训练样本的集合，V_p'表示有标记训练样本集合中除了第p类样本 外的所有样本组成的集合，n_p表示第p类有标记训练集中的样本个数，n_c(j)表示第j个有标 记训练样本所属类的样本个数，χ_i表示第p类的第i个有标记训练样本，χ_j表示V_p'中的第j个 有标记训练样本，U₁表示全波段子数据立方体水平方向的投影矩阵，U₁＝1_5×5，U₂表示全波 段子数据立方体垂直方向的投影矩阵，U₂＝1_5×5，U₁和U₂的每个元素的值都为1，表示张量与矩阵进行模1、模2积，T表示转置，(·)∈R^5×5×D表示一个张量， 表示对两个张量的积的第一阶和第二阶进行缩并，缩并后得到一个大小为 D×D的矩阵,(·)∈R^5×5×D表示一个三阶张量，(·)₍₃₎∈R^D×(5×5)表 示张量(·)按模式3展开成一个矩阵；

(6)构造有标记训练集X的类内不相似性矩阵W：

$W=\underset{p=1}{\overset{c}{\Sigma }}\underset{i\in V_p}{\Sigma }\underset{h\in V_p}{\Sigma }\frac{1}{n_p{n}_p}{\mathrm{mat}}_3\left(\left(({\chi }_i-{\chi }_h){\times }_1{U}_1^T{\times }_2{U}_2^T\right)\right){\mathrm{mat}}_3^T\left(\left(({\chi }_i-{\chi }_h){\times }_1{U}_1^T{\times }_2{U}_2^T\right)\right)$

其中，χ_h表示V_p训练集中的第h个有标记训练样本；

(7)构造总训练集S中所有样本的相似性矩阵M：

$M=\frac{1}{2}\underset{i^{\prime },j^{\prime }}{\overset{N+\eta }{\Sigma }}{m}_{i^{\prime }{j}^{\prime }}\left({\mathrm{mat}}_3\left(\left(({\chi }_{i^{\prime }}-{\chi }_{j^{\prime }}){\times }_1{U}_1^T{\times }_2{U}_2^T\right)\right){\mathrm{mat}}_3^T\left(\left(({\chi }_{i^{\prime }}-{\chi }_{j^{\prime }}){\times }_1{U}_1^T{\times }_2{U}_2^T\right)\right)\right)$

其中，m_i'j'表示样本c_i'和c_j'之间的相似度，χ_i'表示S中的第i'个训练样本，χ_j'表示S中的第j' 个训练样本；

(8)由有标记训练集的类间不相似性矩阵B、有标记训练集的类内不相似性矩阵W和相 似性矩阵M构建张量半监督标度切目标函数：

$U_3^{*}=\underset{U_3}{\arg \max }\frac{\mathrm{tr}\left(U_3^T{\mathrm{BU}}_3\right)}{\mathrm{tr}\left(U_3^T(W+\mathrm{\beta M})U_3\right)}$

其中，参数β是一个微调参数，其值由人为指定为0.001，U₃为所要求的特征维方向上的投 影矩阵，tr表示矩阵的迹；

(9)对张量半监督标度切目标函数求解，得到特征维方向的投影矩阵

(10)分别将有标记训练集和测试集投影到由所张成的低维空间， 得到投影后新的有标记训练集和新的测试集其中为第t个 有标记训练样本的新特征张量，为第u个测试样本的新特征张量，表示 张量与矩阵进行模3积；

(11)将新的有标记训练集类别标记集和新的测试集输 入到支撑矢量机SVM中进行分类，得到测试集的分类结果其中，l_u'表示第u个测试 样本所属的类别标号。本发明与现有技术相比，具有以下优点：

第一，本发明由于采用维数约减算法对高光谱图像数据进行降维，再进行分类，大大减少 了计算量，提高了分类的速度。

第二，本发明基于在一个小范围的空间区域，高光谱图像对应的地物类别单一，像素之 间有很大相似性的特点，将每个样本表示成一个具有全波段的子数据立方体，利用张量计算， 避免了将图像进行向量化处理，可以最大限度地利用区域空间相关性和光谱间相关性；

第三，本发明与现有的张量LDA方法相比，不需要每一类数据满足高斯等方差分布，而 是通过计算样本与样本之间的不相似性来构造不相似性矩阵，消除了类中心的影响。

第四，该发明充分利用了有标记样本提供的信息，寻找能够更好地保持类的可分性的投 影空间，同时由于利用了无标记样本挖掘数据的几何结构信息，故可反应出数据本质的几何 特征；

对比实验表明，本发明有效降低了计算的复杂度，提高了高光谱遥感图像的分类准确率。

附图说明

图1是本发明的实现流程图；

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

图3是本发明与现有方法对Indian Pines图像的分类结果图。

具体实施措施

下面结合附图对发明的技术方案和效果做进一步描述。

参照图1，对本发明实现步骤如下：

步骤1，输入高光谱数据集A∈R^m×n×D，该数据集包含c＝16类地物，其中m×n表示图像 空间大小即像素点的个数，D表示数据集的总波段数，R表示实数域；

步骤2，选取有标记训练集X、测试集Y以及总训练集S。

2a)以A的每个像素点为中心取5*5的邻域块，得到Q个具有全波段的子数据立方体，每 个子数据立方体作为一个样本用一个三阶张量来表示，得到样本集合其中c_a表示第a个样本，Q表示样本总个数，Q＝m×n；

2b)从样本集合中随机选择N个有标记样本构成有标记训练样本集对 应的类标向量记为：其余Q-N个无标记样本构成测试样本集 ,y_u∈R^5×5×D，其中，χ_t表示有标记训练集的第t个样本，l_t表示第t个有标记训练 样本所属的类别标号，y_u表示测试集的第u个样本；

2c)从Q-N个无标记样本中选出h个无标记样本，与N个有标记样本一起构成总训练集 其中，s_k表示总训练集的第k个训练样本，N+η为总训练集的样本个数， 1≤η≤Q-N。

步骤3，构造有标记训练集X的类间不相似性矩阵B和类内不相似性矩阵W。

3a)将有标记训练集X中第p类有标记样本构成同类样本集合p＝1,2,...,c，其中χ_i表示第p类的第i个有标记训练样本，n_p表示集合V_p中有标记训练样本 个数；

3b)将有标记训练样本集合中除了第p类样本外的所有样本构成不同类样本集合 其中χ_j表示V_p'中的第j个有标记训练样本，n_c(j)表示第j个有标记训练样本所 属类的样本个数；

3c)计算同类样本集合V_p中的有标记样本和不同类样本集合V_p'中的有标记样本之间的不 相似性得到每一类的类间不相似性矩阵B_p：

$B_p=\underset{i\in V_p}{\Sigma }\underset{j\in V_P^{\prime }}{\Sigma }\frac{1}{n_p{n}_{c\left(j\right)}}{\mathrm{mat}}_3\left(\left(({\chi }_i-{\chi }_j){\times }_1{U}_1^T{\times }_2{U}_2^T\right)\right){\mathrm{mat}}_3^T\left(\left(({\chi }_i-{\chi }_j){\times }_1{U}_1^T{\times }_2{U}_2^T\right)\right)$

U₁表示全波段子数据立方体水平方向的投影矩阵，U₁＝1_5×5，U₂表示全波段子数据立方体垂 直方向的投影矩阵，U₂＝1_5×5，U₁和U₂的每个元素的值都为1，表示张量与矩阵 进行模1、模2积，T表示转置，表示对两个张量的积的第 一阶和第二阶进行缩并，缩并后得到一个大小为D×D的矩阵，(·)∈R^5×5×D表示一个三阶张量， (·)₍₃₎∈R^D×(5×5)表示张量(·)按模式3展开成一个矩阵_，T表示转置；

3d)计算同类样本集合V_p中的有标记样本之间的不相似性，得到每一类的类内不相似性 矩阵W_p：

$W_p=\underset{i\in V_p}{\Sigma }\underset{h\in V_p}{\Sigma }\frac{1}{n_p{n}_p}{\mathrm{mat}}_3\left(\left(({\chi }_i-{\chi }_h){\times }_1{U}_1^T{\times }_2{U}_2^T\right)\right){\mathrm{mat}}_3^T\left(\left(({\chi }_i-{\chi }_h){\times }_1{U}_1^T{\times }_2{U}_2^T\right)\right)$

χ_h表示V_p训练集中的第h个有标记训练样本；

3e)对步骤3c)每一类有标记训练样本的类间不相似性矩阵B_p进行求和，得到有标记训 练集的类间不相似性矩阵B：

$B=\underset{p=1}{\overset{c}{\Sigma }}{B}_P=\underset{p=1}{\overset{c}{\Sigma }}\underset{i\in V_p}{\Sigma }\underset{j\in V_P^{\prime }}{\Sigma }\frac{1}{n_p{n}_{c\left(j\right)}}{\mathrm{mat}}_3\left(\left(({\chi }_i-{\chi }_h){\times }_1{U}_1^T{\times }_2{U}_2^T\right)\right){\mathrm{mat}}_3^T\left(\left(({\chi }_i-{\chi }_h){\times }_1{U}_1^T{\times }_2{U}_2^T\right)\right)$

3f)对步骤3d)每一类有标记训练样本的类内不相似性矩阵W_p进行求和，得到有标记训 练集的类间不相似性矩阵W：

$W=\underset{p=1}{\overset{c}{\Sigma }}{W}_P=\underset{p=1}{\overset{c}{\Sigma }}\underset{i\in V_p}{\Sigma }\underset{h\in V_p}{\Sigma }\frac{1}{n_p{n}_p}{\mathrm{mat}}_3\left(\left(({\chi }_i-{\chi }_h){\times }_1{U}_1^T{\times }_2{U}_2^T\right)\right){\mathrm{mat}}_3^T\left(\left(({\chi }_i-{\chi }_h){\times }_1{U}_1^T{\times }_2{U}_2^T\right)\right)$

步骤4，根据总训练集S构造无监督样本相似性矩阵M。

4a)计算总训练集S中任意两个样本之间的相似度：

其中，m_i'j'表示样本c_i'和c_j'之间的相似度，χ_i'表示S中的第i'个训练样本，χ_j'表示S中的第j' 个训练样本，δ为核参数；

4b)根据总训练集S计算无监督样本相似性矩阵M：

$M=\frac{1}{2}\underset{i^{\prime },j^{\prime }}{\overset{N+\eta }{\Sigma }}{m}_{i^{\prime }{j}^{\prime }}\left({\mathrm{mat}}_3\left(\left(({\chi }_{i^{\prime }}-{\chi }_{j^{\prime }}){\times }_1{U}_1^T{\times }_2{U}_2^T\right)\right){\mathrm{mat}}_3^T\left(\left(({\chi }_{i^{\prime }}-{\chi }_{j^{\prime }}){\times }_1{U}_1^T{\times }_2{U}_2^T\right)\right)\right)$

步骤5，由有标记训练集的类间不相似性矩阵B、有标记训练集的类内不相似性矩阵W和 相似性矩阵M构建张量半监督标度切目标函数：

$U_3^{*}=\underset{U_3}{\arg \max }\frac{\mathrm{tr}\left(U_3^T{\mathrm{BU}}_3\right)}{\mathrm{tr}\left(U_3^T(W+\mathrm{\beta M})U_3\right)}$

其中，参数β是一个微调参数，其值由人为指定为0.001，U₃为所要求的特征维方向上的投 影矩阵，tr表示矩阵的迹；

步骤6，对张量半监督标度切目标函数求解，得到特征维方向的投影矩阵U₃。

6a)将张量半监督标度切目标函数变换成如下形式：

其中是一个调整参数，为(B+W+β×M)^-B对应的最大特征值，(B+W+β×M)^-表示对 B+W+β×M求逆；

6b)设置维数约减后的特征维数d的值，并通过对上述目标函数中的项进行奇异值分解，得出d个最大特征值和这d个最大特征值所对应的特征向量u₁,u₂,...,u_d， 其中d的值为整数，且0＜d≤D；

6c)用特征向量u₁,u₂,...,u_d构成特征维方向上的投影矩阵

步骤7，分别将有标记训练集X和测试集Y投影到由投影矩阵U^*₃所张成的低维空间，得 到投影后的新的有标记训练集X'和新的测试集Y'。

7a)将原始的有标记训练集投影到由投影矩阵张成的空间中，得到新的有标 记训练集其中，为第t个有标记训练样本的新的特征张量；

7b)将原始的测试集投影到由投影矩阵张成的空间中，得到新的测试集 为第u个测试样本的新的特征张量；

步骤8，将新的有标记训练集类别标记集和新的测试集输 入到支撑矢量机SVM中进行分类，得到测试集的分类结果其中，l_u'表示第u个测试 样本所属的类别标号。

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

1.仿真条件:

仿真实验中所采用的图像是美国宇航局NASA喷气推进实验室的空载可见光/红外成像光 谱仪AVIRIS于1992年6月在印第安纳西北部获取的Indian Pines图像，如图2所示。图2中 共有16类地物信息，每类地物信息的类名和样本个数如表1所示。

表1 Indian Pines数据集类别情况

图2中图像大小为145×145，共220个波段，去除噪声以及大气和水域吸收的波段还有 200个波段。本发明的仿真实验是在AMD(TM)A8CPU、主频1.90GHz，内存8G，Windows 7(64bit)平台上的MATLAB 2011b上实现的。

2.仿真内容及分析

使用本发明与现有两种方法对高光谱图像Indian Pines进行维数约减，现有的两种方法分 别是：标度切降维SC，张量线性判别分析TLDA。

对用本发明和现有的SC、TLDA三种维数约减方法得到的三个降维图像进行分类，其中 分类器SVM的核参数γ采用五重交叉验证方法得到，惩罚因子C设置为100，。相似度矩阵 M的核参数σ设置为1，权重参数β设置为0.001，无标记训练样本个数η固定为2000。

仿真1，从表1所示16类数据中的每一类选出10个样本作为有标记样本，将这16类数 据中的剩余样本作为无标记样本，用本发明与现有两种方法对16类数据进行20次降维分类 实验，取20次分类结果的平均值，作为最终分类正确率，结果如表2所示。

表2 不同方法在Indian Pines数据集上的总体分类精度

从表2可以看出，本发明较基于向量的两种现有方法有很大的优越性；当特征维数大于 10时，本发明方法的分类正确率达到60％以上，明显高于现有方法；

从表2还可以看出，在维数大于25后，本发明的结果趋于稳定，因此只需要采用25维 特征，即可以得到较高的识别率，从而大大减少了计算量。

仿真2，从表1所示的16类数据的每一类像素点中选出10个作为有标记像素点，整幅 Indian Pines图像的剩余像素点作为无标记像素点，用上述三种方法对整幅Indian Pines图像的 所有像素点进行分类，各方法中降维后的特征维数设置为25，结果如图3所示，其中图(3a) 是本发明的分类结果图，图(3b)是现有SC+SVM的分类结果图，图(3c)为现有TLDA+SVM 的分类结果图，

从图(3a)、图(3b)和图(3c)这三个图可以看出，本发明的结果图比由现有的两种 方法得到的更平滑，分类结果更好。

综上，本发明对高光谱图像进行维数约减后使用SVM分类，一方面使用张量的相关运算 避免了将图像向量化，充分利用了空间信息；另一方面，利用有标记与无标记信息，充分挖 掘数据的几何结构信息，提高了分类精度，与现有的方法相比具有优势。