CN105069478A - 基于超像素张量稀疏编码的高光谱遥感地物分类方法 - Google Patents

Abstract

本发明公开了一种基于超像素张量稀疏编码的高光谱遥感地物分类方法，克服了现有技术中不能充分利用高光谱图像的空间信息进行分类和分类速度慢的缺点。本发明实现的步骤是：(1)输入待分类的高光谱图像；(2)构建层次空间相似度矩阵；(3)获得超像素集；(4)构造标记样本字典；(5)求解稀疏系数矩阵；(6)超像素分类；(7)输出待分类高光谱图像的分类结果。本发明具有保持高光谱图像同质区域的空间一致性和分类速度快的优点，可用于高光谱图像的快速分类。

Description

基于超像素张量稀疏编码的高光谱遥感地物分类方法

技术领域

本发明属于图像处理技术领域，更进一步涉及高光谱图像分类技术领域中的一种基于超像素张量稀疏编码的高光谱遥感地物分类方法。本发明可用于对高光谱遥感图像进行快速的地物分类。

背景技术

目前高光谱数据的研究重心主要是降维和分类两方面。其中，分类主要是根据不同类别的数据具有不同的特征信息对数据进行处理，为每一个像素赋予一个所属的标记，从而实现对地物进行分类和识别。根据是否有标记样本参与处理过程，分类方法可分为无监督分类方法，有监督分类方法和半监督分类方法。在众多的高光谱分类方法中，稀疏编码分类方法是研究的重点方向之一。该分类方法的核心思想是：高光谱数据中的一个样本，可以由一些训练样本进行线性表示，对于一个未知样本，可以对其进行稀疏表示得到一个稀疏向量，该向量对应于训练样本所占的权重，通过求得该向量，可以直接确定未知样本的类别标签。

JianjunLiu，ZebinWu，ZhihuiWei等人在其发表的论文“Spatial-SpectralKernelSparseRepresentationforHyperspectralImageClassification”(AppliedEarthObservationsandRemoteSensing,IEEE,2013,6(6):2462-2471)中提出一种基于空谱核稀疏表示的高光谱图像分类方法。该方法首先定义一种空谱核，并使用该空谱核获取空间信息和谱信息，然后将稀疏表示分类器扩展到可空间，得到核稀疏表示分类器，最后求解稀疏系数，并得到样本的类别。该方法存在的不足之处是，在求解稀疏系数的过程中需要对每个映射后的样本求解其在映射后字典上的稀疏系数，对于数据量巨大的高光谱图像而言，单个像素进行计算会产生很大的计算负担。

西安电子科技大学在其申请的专利文件“基于稀疏核编码SKR的高光谱图像地物识别方法”(申请号：201110260639.X申请日2011-09-05公告号：CN102324047B)中公开了一种基于稀疏核编码SKR的高光谱图像地物识别方法。该方法首先利用高光谱图像中一些已知标签的光谱向量按类排列作为稀疏编码的字典，所有未知标签的光谱向量样本构成测试样本集合，其次，利用近邻方法构造中心样本矩阵，通过构造稀疏核函数将测试样本和字典分别映射至特征空间得到映射后的字典和测试样本，并对映射后的字典进行列归一化，再次，利用归一化后的字典对映射后的测试样本进行稀疏编码，通过误差判别公式判断测试样本的类别。该方法存在的不足之处是，没有考虑到相邻样本间的空间信息，会造成同质区域的空间一致性不好，影响高光谱图像的识别准确率。

发明内容

本发明的目的在于克服上述已有技术的不足，提出一种基于超像素张量稀疏编码的高光谱遥感地物分类方法。本发明可以充分利用样本间的空间和谱间信息，保证分类结果的空间一致性，加快稀疏编码分类过程的速度，实现对高光谱图像快速准确的分类。

本发明实现的具体步骤如下：

(1)输入待分类的高光谱图像：

输入待分类的高光谱图像，将输入的高光谱图像中的每一个像素点设定为一个样本；

(2)构建层次空间相似度矩阵：

(2a)选取待分类的高光谱图像中任意一个样本周围的空间近邻样本，按照该样本与其空间近邻样本的距离，将距离该样本最近的空间近邻样本作为第一层空间近邻样本，将距离该样本较近的空间近邻样本作为第二层空间近邻样本，将距离该样本最远的空间近邻样本作为第三层空间近邻样本，将三层空间近邻样本组成该样本的层次空间近邻集，对待分类的高光谱图像中的所有样本进行相同处理，得到每个样本的层次空间近邻集；

(2b)将待分类的高光谱图像中的一个样本作为中心样本，采用下式，计算中心样本与其层次空间近邻集中的一个样本的相似度：

u_h＝exp(-||x-x_h||²/2(∑||x-x_h||²/N_g))s.t.x_h∈H_g

其中，u_h表示中心样本x与层次空间近邻集中样本x_h的相似度，h表示层次空间近邻集中样本x_h的索引值，h的取值范围为{1,2,…,n}，n表示待分类的高光谱图像中样本的总数，exp(·)表示指数操作，||·||表示范式操作，x表示中心样本，x_h表示层次空间近邻集中的样本，∑表示求和操作，N_g表示层次空间近邻集中第g层近邻样本的样本个数，g表示层数，g的取值范围为{1,2,3}，s.t.表示约束操作，∈表示属于符号，H_g表示中心样本x的层次空间近邻集中第g层近邻样本；

(2c)按照步骤(2b)，计算中心样本与层次空间近邻集中所有样本的相似度；

(2d)在层次空间近邻集的所有样本中，找到与中心样本具有最大相似度的样本，并记录该样本的位置；

(2e)根据步骤(2d)得到的最大相似度样本和该样本的位置，构造中心样本的层次空间相似度向量，其中该向量中与最大相似度样本对应的位置等于最大相似度；

(2f)将待分类的高光谱图像中所有样本依次选为中心样本，执行步骤(2b)、步骤(2c)、步骤(2d)和步骤(2e)，得到每个中心样本的层次空间相似度向量；

(2g)将所有的层次空间相似度向量按行排列，得到层次空间相似度矩阵；

(3)获得超像素集：

(3a)创建可靠度矩阵和责任度矩阵，将可靠度矩阵和责任度矩阵中的元素初始化为零；

(3b)采用责任度矩阵元素更新公式，得到更新后的责任度矩阵元素；

(3c)采用可靠度矩阵非对角线元素更新公式，得到更新后的可靠度矩阵非对角线元素；

(3d)采用可靠度矩阵对角线元素更新公式，得到更新后的可靠度矩阵对角线元素；

(3e)判断可靠度矩阵对角线的元素是否更新了20次，若是，执行步骤(3f)，否则，执行步骤(3b)；

(3f)采用聚类中心计算公式，得到待分类的高光谱图像中所有样本所属的聚类中心；

(3g)将待分类的高光谱图像中具有相同聚类中心的样本聚成一簇，得到多个具有不同聚类中心的簇，将具有不同聚类中心的簇组成簇集；

(3h)将簇集中的单个簇当作超像素，得到一系列超像素，并组成超像素集；

(4)构造标记样本字典：

(4a)在待分类的高光谱图像中，选取L个样本作为标记样本，获取每个标记样本的类别，其中，L表示标记样本的数目，L的取值范围为{50,100,150,200,250,300,350}；

(4b)将标记样本按照类别进行排列，构成标记样本字典；

(5)求解稀疏系数矩阵：

(5a)选取超像素集中的一个超像素作为待处理的超像素，将该超像素与其相邻的超像素组成超像素张量；

(5b)创建残差和标记样本字典的原子列集合，将残差初始化为步骤(5a)中得到的超像素张量，将标记样本字典的原子列集合初始化为空集合；

(5c)采用下式，计算标记样本字典中最接近残差的原子列数：

$d=\mathrm{argmax}\underset{b}{\Σ}\underset{c}{\Σ}\left|\right|R{\×}_1{D}_a^T{\×}_2{P}_b{\×}_3{Q}_c\left|\right|$

其中，d表示标记样本字典中最接近残差的原子列数，arg表示求取满足条件变量值的操作，max表示最大值操作，∑表示求和操作，b表示单位矩阵P的列数，c表示单位矩阵Q的列数，||·||表示范式操作，R表示残差，×₁表示张量的模-1乘积操作，表示转置后标记样本字典的第a列，T表示转置操作，a表示转置后标记样本字典的列数，D表示标记样本字典，×₂表示张量的模-2乘积操作，P_b表示单位矩阵P的第b列，P表示总列数为p的单位矩阵，p表示待处理超像素中样本的总数，×₃表示张量的模-3乘积操作，Q_c表示单位矩阵Q的第c列，Q表示总列数为q的单位矩阵，q表示超像素张量中超像素的总数；

(5d)将标记样本字典中最接近残差的原子列数加入标记样本字典的原子列集合；

(5e)采用下式，计算稀疏系数张量：

$B=\arg min\left|\right|A{\×}_1{D}_{\mathrm{\Γ}}{\×}_{2}P{\×}_{3}Q-Y|{|}_2^2$

其中，B表示稀疏系数张量，arg表示求取满足条件变量值的操作，min表示最小值操作，||·||₂表示2-范式操作，A表示待求的稀疏系数张量，×₁表示张量的模-1乘积操作，D_Γ表示标记样本字典中对应于Γ中的列组成的字典，Γ表示标记样本字典的原子列集合，D表示标记样本字典，×₂表示张量的模-2乘积操作，P表示总列数为p的单位矩阵，p表示待处理超像素中样本的总数，×₃表示张量的模-3乘积操作，Q表示总列数为q的单位矩阵，q表示超像素张量中超像素的总数，Y表示超像素张量；

(5f)采用下式，更新残差：

R′＝Y-B×₁D_Γ×₂P×₃Q

其中，R′表示更新后的残差，Y表示超像素张量，B表示稀疏系数张量，×₁表示张量的模-1乘积操作，D_Γ表示标记样本字典中对应于Γ中的列组成的字典，D表示标记样本字典，Γ表示标记样本的原子列集合，×₂表示张量的模-2乘积操作，P表示总列数为p的单位矩阵，p表示待处理超像素中样本的总数，×₃表示张量的模-3乘积操作，Q表示总列数为q的单位矩阵，q表示超像素张量中超像素的总数；

(5g)判断是否满足下列条件之一，若是，执行(5h)，否则，执行步骤(5c)：

M≥Q,||R′||_F＜ε

其中，M表示更新后标记样本字典的原子列集合中的元素数目，≥表示大于等于符号，Q表示标记样本字典的原子列集合的最大容量，Q的取值范围为{30,31,…,60}，||·||_F表示Frobenius范数操作，R′表示更新后的残差，＜表示小于符号，ε表示残差阈值，ε的取值范围为{10^-7,10^-6,10^-5,10^-4,10^-3}；

(5h)将稀疏系数张量中第三维第一个位置对应的矩阵，作为待处理超像素的稀疏系数矩阵；

(5i)判断是否超像素集中所有的超像素都得到稀疏系数矩阵，若是，则执行步骤(6)，否则，执行步骤(5a)；

(6)超像素分类：

(6a)选取超像素集中的一个超像素作为待分类的超像素，采用超像素分类公式，对该超像素进行分类；

(6b)判断是否对超像素集中所有的超像素都进行了分类，若是，则执行步骤(7)，否则，执行步骤(6a)；

(7)输出待分类高光谱图像的分类结果。

与现有技术相比，本发明具有以下优点：

第一，由于本发明采用了层次空间相似度矩阵获取超像素集，克服了现有技术因没有考虑到相邻样本间的空间信息而无法精确地对高光谱遥感影像分类的缺点，使得本发明具有在同质区域分类效果更好的优点。

第二，由于本发明采用了稀疏系数矩阵求解方法获取超像素的稀疏系数矩阵，克服了现有技术因求解稀疏系数过程中以单个像素为处理单元而增加计算负担的缺点，使得本发明具有对高光谱图像进行快速准确分类的优点。

附图说明

图1是本发明的流程图；

图2是本发明层次空间近邻集的示意图；

图3是本发明采用的高光谱图像IndianPines及其真实标记图；

图4是本发明与现有技术在选取350个标记样本时，对高光谱图像IndianPines分类结果对比图。

具体实施方式

下面结合附图对本发明做进一步的描述。

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

步骤1，输入待分类的高光谱图像。

输入待分类的高光谱图像，将输入的高光谱图像中的每一个像素点设定为一个样本。

步骤2，构建层次空间相似度矩阵。

第1步，选取待分类的高光谱图像中任意一个样本周围的空间近邻样本，按照该样本与其空间近邻样本的距离，将距离该样本最近的空间近邻样本作为第一层空间近邻样本，将距离该样本较近的空间近邻样本作为第二层空间近邻样本，将距离该样本最远的空间近邻样本作为第三层空间近邻样本，将三层空间近邻样本组成该样本的层次空间近邻集，对待分类的高光谱图像中的所有样本进行相同处理，得到每个样本的层次空间近邻集。

附图2为本发明层次空间近邻集的示意图。图2中的黑色区域表示第一层空间近邻样本，浅灰色区域表示第二层空间近邻样本，深灰色区域表示第三层空间近邻样本，三层空间近邻的样本组成层次空间近邻集。

第2步，将待分类的高光谱图像中的一个样本作为中心样本，采用下式，计算中心样本与其层次空间近邻集中的一个样本的相似度：

u_h＝exp(-||x-x_h||²/2(∑||x-x_h||²/N_g))s.t.x_h∈H_g

其中，u_h表示中心样本x与层次空间近邻集中样本x_h的相似度，h表示层次空间近邻集中样本x_h的索引值，h的取值范围为{1,2,…,n}，n表示待分类的高光谱图像中样本的总数，exp(·)表示指数操作，||·||表示范式操作，x表示中心样本，x_h表示层次空间近邻集中的样本，∑表示求和操作，N_g表示层次空间近邻集中第g层近邻样本的样本个数，g表示层数，g的取值范围为{1,2,3}，s.t.表示约束操作，∈表示属于符号，H_g表示中心样本x的层次空间近邻集中第g层近邻样本。

第3步，按照第2步，计算中心样本与层次空间近邻集中所有样本的相似度。

第4步，在层次空间近邻集的所有样本中，找到与中心样本具有最大相似度的样本，并记录该样本的位置。

第5步，根据第4步得到的最大相似度样本和该样本的位置，构造中心样本的层次空间相似度向量，其中该向量中与最大相似度样本对应的位置等于最大相似度。

第6步，将待分类的高光谱图像中所有样本依次选为中心样本，执行第2步、第3步、第4步和第5步，得到每个中心样本的层次空间相似度向量。

第7步，将所有的层次空间相似度向量按行排列，得到层次空间相似度矩阵。

步骤3，获得超像素集。

第1步，创建可靠度矩阵和责任度矩阵，将可靠度矩阵和责任度矩阵中的元素初始化为零。

第2步，采用责任度矩阵元素更新公式，得到更新后的责任度矩阵元素，责任度矩阵元素更新公式如下：

$\begin{array}{ccc}P(u,v)=S(u,v)-max\left(E\right(u,g)+S(u,g\left)\right)& s.t.& g\≠v,\∀u,v,g\∈\{1,...,n\}\end{array}$

其中，P(u,v)表示责任度矩阵中的第u行第v列元素，u表示责任度矩阵的行数，v表示责任度矩阵的列数，S(u,v)表示层次空间相似度矩阵中的第u行第v列元素，u表示层次空间相似度矩阵的行数，v表示层次空间相似度矩阵的列数，max表示取最大值操作，E(u,g)表示可靠度矩阵中的第u行第g列元素，u表示可靠度矩阵的行数，g表示可靠度矩阵的列数，S(u,g)表示层次空间相似度矩阵中的第u行第g列元素，u表示层次空间相似度矩阵的行数，g表示层次空间相似度矩阵的列数，s.t.表示约束操作，≠表示不等于符号，表示取所有值操作，∈表示属于符号，n表示待分类的高光谱图像中样本的总数。

第3步，采用可靠度矩阵非对角线元素更新公式，得到更新后的可靠度矩阵非对角线元素，可靠度矩阵非对角线元素更新公式如下：

$\begin{array}{ccc}E(u,v)=min\{0,P(v,v)+\mathrm{\Σ}max(0,P(k,v\left)\right)\}& s.t.& k\∉\{u,v\},\∀u,v,k\∈\{1,...,n\}\end{array}$

其中，E(u,v)表示可靠度矩阵中的第u行第v列元素，u表示可靠度矩阵的行数，v表示可靠度矩阵的列数，min表示取最小值操作，P(v,v)表示责任度矩阵中的第v行第v列元素，v表示责任度矩阵的行数，∑表示求和操作，max表示取最大值操作，P(k,v)表示责任度矩阵中的第k行第v列元素，k表示责任度矩阵的行数，v表示责任度矩阵的列数，s.t.表示约束操作，表示不属于符号，表示取所有值操作，∈表示属于符号，n表示待分类的高光谱图像中样本的总数。

第4步，采用可靠度矩阵对角线元素更新公式，得到更新后的可靠度矩阵对角线元素可靠度矩阵对角线元素更新公式如下：

$\begin{array}{ccc}E(v,v)=\mathrm{\Σ}max(0,P(z,v\left)\right)& s.t.& z\≠v,\∀v,z\∈\{1,...,n\}\end{array}$

其中，E(v,v)表示可靠度矩阵中的第v行第v列的对角元素，v表示可靠度矩阵的行数，∑表示求和操作，max表示取最大值操作，P(i,v)表示责任度矩阵中的第i行第v列元素，i表示责任度矩阵的行数，v表示责任度矩阵的列数，s.t.表示约束操作，≠表示不等于符号，表示取所有值操作，∈表示属于符号，n表示待分类的高光谱图像中样本的总数。

第5步，判断可靠度矩阵对角线的元素是否更新了20次，若是，执行第6步，否则，执行第2步。

第6步，采用聚类中心计算公式，得到待分类的高光谱图像中所有样本所属的聚类中心，聚类中心计算公式如下：

$w_s=\mathrm{argmax}\left(E\right(s,v)+P(s,v\left)\right),\∀s,v\∈\{1,...,n\}$

其中，w_s表示待分类的高光谱图像中第s个样本所属的聚类中心，s表示待分类的高光谱图像中样本的索引值，s的取值范围为{1,2,…,n}，arg表示求取满足条件变量值的操作，max表示最大值操作，E(s,v)表示可靠度矩阵中的第s行第v列的对角元素，s表示可靠度矩阵的行数，v表示可靠度矩阵的列数，P(s,v)表示责任度矩阵中的第s行第v列元素，s表示可靠度矩阵的行数，v表示可靠度矩阵的列数，表示取所有值操作，∈表示属于符号，n表示待分类的高光谱图像中样本的总数。

第7步，将待分类的高光谱图像中具有相同聚类中心的样本聚成一簇，得到多个具有不同聚类中心的簇，将具有不同聚类中心的簇组成簇集。

第8步，将簇集中的单个簇当作超像素，得到一系列超像素，并组成超像素集。

步骤4，构造标记样本字典。

第1步，在待分类的高光谱图像中，选取L个样本作为标记样本，获取每个标记样本的类别，其中，L表示标记样本的数目，L的取值范围为{50,100,150,200,250,300,350}。

第2步，将L个标记样本按照类别进行排列，构成标记样本字典。

步骤5，求解稀疏系数矩阵。

第1步，选取超像素集中的一个超像素作为待处理的超像素，将该超像素与其相邻的超像素组成超像素张量。

在本发明的实施实例中，选取5个相邻的超像素与待处理的超像素组成超像素张量。

第2步，创建残差和标记样本字典的原子列集合，将残差初始化为第1步中得到的超像素张量，将标记样本字典的原子列集合初始化为空集合。

第3步，采用下式，计算标记样本字典中最接近残差的原子列数：

$d=\mathrm{argmax}\underset{b}{\Σ}\underset{c}{\Σ}\left|\right|R{\×}_1{D}_a^T{\×}_2{P}_b{\×}_3{Q}_c\left|\right|$

其中，d表示标记样本字典中最接近残差的原子列数，arg表示求取满足条件变量值的操作，max表示最大值操作，∑表示求和操作，b表示单位矩阵P的列数，c表示单位矩阵Q的列数，||·||表示范式操作，R表示残差，×₁表示张量的模-1乘积操作，表示转置后标记样本字典的第a列，T表示转置操作，a表示转置后标记样本字典的列数，D表示标记样本字典，×₂表示张量的模-2乘积操作，P_b表示单位矩阵P的第b列，P表示总列数为p的单位矩阵，p表示待处理超像素中样本的总数，×₃表示张量的模-3乘积操作，Q_c表示单位矩阵Q的第c列，Q表示总列数为q的单位矩阵，q表示超像素张量中超像素的总数。

第4步，将标记样本字典中最接近残差的原子列数加入标记样本字典的原子列集合。

第5步，采用下式，计算稀疏系数张量：

$B=\arg min\left|\right|A{\×}_1{D}_{\mathrm{\Γ}}{\×}_{2}P{\×}_{3}Q-Y|{|}_2^2$

其中，B表示稀疏系数张量，arg表示求取满足条件变量值的操作，min表示最小值操作，||·||₂表示2-范式操作，A表示待求的稀疏系数张量，×₁表示张量的模-1乘积操作，D_Γ表示标记样本字典中对应于Γ中的列组成的字典，Γ表示标记样本字典的原子列集合，D表示标记样本字典，×₂表示张量的模-2乘积操作，P表示总列数为p的单位矩阵，p表示待处理超像素中样本的总数，×₃表示张量的模-3乘积操作，Q表示总列数为q的单位矩阵，q表示超像素张量中超像素的总数，Y表示超像素张量。

第6步，采用下式，更新残差：

R′＝Y-B×₁D_Γ×₂P×₃Q

其中，R′表示更新后的残差，Y表示超像素张量，B表示稀疏系数张量，×₁表示张量的模-1乘积操作，D_Γ表示标记样本字典中对应于Γ中的列组成的字典，D表示标记样本字典，Γ表示标记样本的原子列集合，×₂表示张量的模-2乘积操作，P表示总列数为p的单位矩阵，p表示待处理超像素中样本的总数，×₃表示张量的模-3乘积操作，Q表示总列数为q的单位矩阵，q表示超像素张量中超像素的总数。

第7步，判断是否满足下列条件之一，若是，执行第8步，否则，执行第3步：

M≥Q,||R′||_F＜ε

其中，M表示更新后标记样本字典的原子列集合中的元素数目，≥表示大于等于符号，Q表示标记样本字典的原子列集合的最大容量，Q的取值范围为{30,31,…,60}，||·||_F表示Frobenius范数操作，R′表示更新后的残差，＜表示小于符号，ε表示残差阈值，ε的取值范围为{10^-7,10^-6,10^-5,10^-4,10^-3}。

在本发明的实施实例中，标记样本字典的原子列集合的最大容量Q的值为45，残差阈值ε的值为10^-6

第8步，将稀疏系数张量中第三维第一个位置对应的矩阵，作为待处理超像素的稀疏系数矩阵。

第9步，判断是否超像素集中所有的超像素都得到稀疏系数矩阵，若是，则执行步骤6，否则，执行第1步。

步骤6，超像素分类。

第1步，选取超像素集中的一个超像素作为待分类的超像素，采用超像素分类公式，对该超像素进行分类，超像素分类公式如下：

$C=\mathrm{argmin}\underset{r}{\Σ}\left|\right|y_r-D^m{Z}_r^m|{|}_F^2,\∀m\∈\{1,...,K\}$

其中，C表示待分类超像素所属的类别，arg表示求取满足条件变量值的操作，min表示最小值操作，∑表示求和操作，r表示待分类超像素中样本的索引值，r的取值范围为{1,2,…,P}，P表示待分类超像素中样本的总数，||·||_F表示Frobenius范数操作，y_r表示待分类超像素中第r个样本，y表示待分类超像素，D^m表示由标记样本字典中与第m类对应的列组成的字典，m表示标记样本字典中的类别数，m的取值范围为{1,2,…,K}，D表示标记样本字典，表示矩阵Z^m的第r列，Z^m表示待分类超像素的稀疏系数矩阵中与第m类对应的行组成的矩阵，Z表示待分类超像素的稀疏系数矩阵，表示取所有值操作，∈表示属于符号，K表示标记样本字典中的类别总数。

在本发明的实施实例中，标记样本字典中的类别总数K为16。

第2步，判断是否对超像素集中所有的超像素都进行了分类，若是，则执行步骤7，否则，执行第1步。

步骤7，输出待分类高光谱图像的分类结果。

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

1.仿真实验条件：

本仿真实验的硬件测试平台是：处理器为IntelCorei5CPU，主频为2.40GHz，内存4GB，软件平台为：Windows7操作系统和MatlabR2014a。本发明的输入图像为高光谱图像印第安松树IndianPines，图像大小为145×145×220，图像共包含220个波段和16类地物，图像格式为TIF。

2.仿真实验内容：

本发明用到的四个现有技术对比分类方法分别如下：

Melgani等人在“Classificationofhyperspectralremotesensingimageswithsupportvectormachines,IEEETrans.Geosci.RemoteSens.,vol.42,no.8,pp.1778–1790,Aug.2004”中提出的高光谱图像分类方法，简称支撑向量机SVM分类方法。

Chen等人在“Hyperspectralimageclassificationusingdictionary-basedsparserepresentation,IEEETrans.Geosci.RemoteSens.,vol.49,no.10,pp.3973–3985,oct2011”中提出的高光谱图像分类方法，简称稀疏编码SCC分类方法。

Bor-Chen等人在“Kernel-BasedKNNandGaussianClassifiersforHyperspectralImageClassification,IEEEInternationalGeosci.RemoteSens.Symposium,vol.2,pp.II-1006,II-1008,July2008”中提出的高光谱图像分类方法，简称K近邻KNN分类方法。

Fang等人在“Spectral–SpatialHyperspectralImageClassificationviaMultiscaleAdaptiveSparseRepresentation,IEEETrans.Geosci.RemoteSens.,vol.52,no.12，pp.7738-7749,Dec.2014”中提出的高光谱图像分类方法，简称多尺度自适应稀疏表示MASR分类方法。

在仿真实验中，采用以下两个现有指标来评价本发明方法的性能：

第一个评价指标是总精度(OA)，表示正确分类的样本占所有样本的比例，值越大，说明分类效果越好。

第二个评价指标是卡方系数(Kappa)，表示混淆矩阵中不同的权值，值越大，说明分类效果越好。

附图3为本发明输入的高光谱图像印第安松树IndianPines的真实地物分布图。其中，图3(a)为本发明输入的高光谱图像印第安松树IndianPines，该高光谱图像由美国宇航局的航空可见光/红外光成像光谱仪获取，大小为145×145，图3(b)为本发明输入的高光谱图像印第安松树IndianPines的真实标记图，其大小为145×145。

附图4为仿真实验中本发明与现有技术在选取350个标记样本时，对高光谱图像印第安松树IndianPines分类结果对比图。其中，图4(a)为采用本发明方法对高光谱图像印第安松树IndianPines进行分类的结果图，图4(b)为采用现有技术的多尺度自适应稀疏表示MASR分类方法对高光谱图像印第安松树IndianPines进行分类的结果图，图4(c)为采用现有技术的稀疏编码SCC分类方法对高光谱图像印第安松树IndianPines进行分类的结果图，图4(d)为采用现有技术的K近邻KNN分类方法对高光谱图像印第安松树IndianPines进行分类的结果图，图4(e)为采用现有技术的支撑向量机SVM分类方法对高光谱图像印第安松树IndianPines进行分类的结果图。

3.仿真实验结果分析：

本仿真实验采用本发明和现有技术的四种分类方法对仿真实验数据分类。

对本发明和现有技术的四种分类方法选取不同的标记样本的总数，并对本发明和现有技术的四种分类方法进行仿真实验，得到表1。表1为标记样本总数分别取{50,100,150,200,250,300,350}时，从客观评价指标上对本发明和现有技术的四种分类方法的分类结果进行评价。

表1.本发明与现有技术的四种分类方法在标记样本总数不同时的定量分析

从表1中可以看出，现有技术的支撑向量机SVM分类方法、稀疏编码SCC分类方法和K近邻KNN分类方法的总精度(OA)和卡方系数(Kappa)都较低，主要因为缺少高光谱图像的空间信息，使得分类准确度低；同时，现有技术的多尺度自适应稀疏表示MASR分类方法的总精度(OA)和卡方系数(Kappa)有所提高，但是空间信息的提取不够充分，所以也存在明显的错分现象；本发明在定量分析方面优于前四种现有技术分类方法，总精度(OA)和卡方系数(Kappa)有明显的提高，说明本方法能实现准确的高光谱图像分类。

从附图4中可以看出，现有技术的支撑向量机SVM分类方法、稀疏编码SCC分类方法和K近邻KNN分类方法的分类结果比较差，主要因为这三种方法都只利用了高光谱图像的光谱信息，忽略了空间信息，使得分类结果在同质区域都很低；同时，现有技术的多尺度自适应稀疏表示MASR分类方法的分类结果有很大改善，错分率有所下降，但是现有技术的多尺度自适应稀疏表示MASR分类方法只是利用了空域的四个近邻信息，不能充分挖掘图像的近邻空间信息，所以依然存在图像在同质区域的错分的现象；本发明在视觉方面均优于前四种现有技术分类方法，具有在同质区域的空间一致性好的优点。

以上仿真实验表明：本发明方法能够充分利用高光谱图像的光谱域和邻域信息，在高光谱图像的同质区域能取得很好的分类结果，同时，以超像素为基本处理单元，加快分类的速度，能够解决现有技术方法中存在的忽略高光谱图像的邻域信息、分类速度慢等问题，是一种非常实用的高光谱图像分类方法。

Claims (6)

1.一种基于超像素张量稀疏编码的高光谱遥感地物分类方法，包括如下步骤：

(1)输入待分类的高光谱图像：

输入待分类的高光谱图像，将输入的高光谱图像中的每一个像素点设定为一个样本；

(2)构建层次空间相似度矩阵：

(2a)选取待分类的高光谱图像中任意一个样本周围的空间近邻样本，按照该样本与其空间近邻样本的距离，将距离该样本最近的空间近邻样本作为第一层空间近邻样本，将距离该样本较近的空间近邻样本作为第二层空间近邻样本，将距离该样本最远的空间近邻样本作为第三层空间近邻样本，将三层空间近邻样本组成该样本的层次空间近邻集，对待分类的高光谱图像中的所有样本进行相同处理，得到每个样本的层次空间近邻集；

(2b)将待分类的高光谱图像中的一个样本作为中心样本，采用下式，计算中心样本与其层次空间近邻集中的一个样本的相似度：

u_h＝exp(-||x-x_h||²/2(∑||x-x_h||²/N_g))s.t.x_h∈H_g

其中，u_h表示中心样本x与层次空间近邻集中样本x_h的相似度，h表示层次空间近邻集中样本x_h的索引值，h的取值范围为{1,2,…,n}，n表示待分类的高光谱图像中样本的总数，exp(·)表示指数操作，||·||表示范式操作，x表示中心样本，x_h表示层次空间近邻集中的样本，∑表示求和操作，N_g表示层次空间近邻集中第g层近邻样本的样本个数，g表示层数，g的取值范围为{1,2,3}，s.t.表示约束操作，∈表示属于符号，H_g表示中心样本x的层次空间近邻集中第g层近邻样本；

(2c)按照步骤(2b)，计算中心样本与层次空间近邻集中所有样本的相似度；

(2d)在层次空间近邻集的所有样本中，找到与中心样本具有最大相似度的样本，并记录该样本的位置；

(2e)根据步骤(2d)得到的最大相似度样本和该样本的位置，构造中心样本的层次空间相似度向量，其中该向量中与最大相似度样本对应的位置等于最大相似度；

(2f)将待分类的高光谱图像中所有样本依次选为中心样本，执行步骤(2b)、步骤(2c)、步骤(2d)和步骤(2e)，得到每个中心样本的层次空间相似度向量；

(2g)将所有的层次空间相似度向量按行排列，得到层次空间相似度矩阵；

(3)获得超像素集：

(3a)创建可靠度矩阵和责任度矩阵，将可靠度矩阵和责任度矩阵中的元素初始化为零；

(3b)采用责任度矩阵元素更新公式，得到更新后的责任度矩阵元素；

(3c)采用可靠度矩阵非对角线元素更新公式，得到更新后的可靠度矩阵非对角线元素；

(3d)采用可靠度矩阵对角线元素更新公式，得到更新后的可靠度矩阵对角线元素；

(3e)判断可靠度矩阵对角线的元素是否更新了20次，若是，执行步骤(3f)，否则，执行步骤(3b)；

(3f)采用聚类中心计算公式，得到待分类的高光谱图像中所有样本所属的聚类中心；

(3g)将待分类的高光谱图像中具有相同聚类中心的样本聚成一簇，得到多个具有不同聚类中心的簇，将具有不同聚类中心的簇组成簇集；

(3h)将簇集中的单个簇当作超像素，得到一系列超像素，并组成超像素集；

(4)构造标记样本字典：

(4a)在待分类的高光谱图像中，选取L个样本作为标记样本，获取每个标记样本的类别，其中，L表示标记样本的数目，L的取值范围为{50,100,150,200,250,300,350}；

(4b)将标记样本按照类别进行排列，构成标记样本字典；

(5)求解稀疏系数矩阵：

(5a)选取超像素集中的一个超像素作为待处理的超像素，将该超像素与其相邻的超像素组成超像素张量；

(5b)创建残差和标记样本字典的原子列集合，将残差初始化为步骤(5a)中得到的超像素张量，将标记样本字典的原子列集合初始化为空集合；

(5c)采用下式，计算标记样本字典中最接近残差的原子列数：

$d=\mathrm{argmax}\underset{b}{\Σ}\underset{c}{\Σ}\left|\right|R{\×}_1{D}_a^T{\×}_2{P}_b{\×}_3{Q}_c\left|\right|$

其中，d表示标记样本字典中最接近残差的原子列数，arg表示求取满足条件变量值的操作，max表示最大值操作，∑表示求和操作，b表示单位矩阵P的列数，c表示单位矩阵Q的列数，||·||表示范式操作，R表示残差，×₁表示张量的模-1乘积操作，表示转置后标记样本字典的第a列，T表示转置操作，a表示转置后标记样本字典的列数，D表示标记样本字典，×₂表示张量的模-2乘积操作，P_b表示单位矩阵P的第b列，P表示总列数为p的单位矩阵，p表示待处理超像素中样本的总数，×₃表示张量的模-3乘积操作，Q_c表示单位矩阵Q的第c列，Q表示总列数为q的单位矩阵，q表示超像素张量中超像素的总数；

(5d)将标记样本字典中最接近残差的原子列数加入标记样本字典的原子列集合；

(5e)采用下式，计算稀疏系数张量：

$B=\arg min\left|\right|A{\×}_1{D}_{\mathrm{\Γ}}{\×}_{2}P{\×}_{3}Q-Y|{|}_2^2$

其中，B表示稀疏系数张量，arg表示求取满足条件变量值的操作，min表示最小值操作，||·||₂表示2-范式操作，A表示待求的稀疏系数张量，×₁表示张量的模-1乘积操作，D_Γ表示标记样本字典中对应于Γ中的列组成的字典，Γ表示标记样本字典的原子列集合，D表示标记样本字典，×₂表示张量的模-2乘积操作，P表示总列数为p的单位矩阵，p表示待处理超像素中样本的总数，×₃表示张量的模-3乘积操作，Q表示总列数为q的单位矩阵，q表示超像素张量中超像素的总数，Y表示超像素张量；

(5f)采用下式，更新残差：

R′＝Y-B×₁D_Γ×₂P×₃Q

其中，R′表示更新后的残差，Y表示超像素张量，B表示稀疏系数张量，×₁表示张量的模-1乘积操作，D_Γ表示标记样本字典中对应于Γ中的列组成的字典，D表示标记样本字典，Γ表示标记样本的原子列集合，×₂表示张量的模-2乘积操作，P表示总列数为p的单位矩阵，p表示待处理超像素中样本的总数，×₃表示张量的模-3乘积操作，Q表示总列数为q的单位矩阵，q表示超像素张量中超像素的总数；

(5g)判断是否满足下列条件之一，若是，执行(5h)，否则，执行步骤(5c)：

M≥Q,||R′||_F<ε

其中，M表示更新后标记样本字典的原子列集合中的元素数目，≥表示大于等于符号，Q表示标记样本字典的原子列集合的最大容量，Q的取值范围为{30,31,…,60}，||·||_F表示Frobenius范数操作，R′表示更新后的残差，<表示小于符号，ε表示残差阈值，ε的取值范围为{10^-7,10^-6,10^-5,10^-4,10^-3}；

(5h)将稀疏系数张量中第三维第一个位置对应的矩阵，作为待处理超像素的稀疏系数矩阵；

(5i)判断是否超像素集中所有的超像素都得到稀疏系数矩阵，若是，则执行步骤(6)，否则，执行步骤(5a)；

(6)超像素分类：

(6a)选取超像素集中的一个超像素作为待分类的超像素，采用超像素分类公式，对该超像素进行分类；

(6b)判断是否对超像素集中所有的超像素都进行了分类，若是，则执行步骤(7)，否则，执行步骤(6a)；

(7)输出待分类高光谱图像的分类结果。

2.根据权利要求1所述的基于超像素张量稀疏编码的高光谱遥感地物分类方法，其特征在于，步骤(3b)所述的责任度矩阵元素更新公式如下：

$\begin{array}{ccc}P(u,v)=S(u,v)-max\left(E\right(u,g)+S(u,g\left)\right)& s.t.& g\&NotEqual;v,\&ForAll;u,v,g\&Element;\{1,...,n\}\end{array}$

其中，P(u,v)表示责任度矩阵中的第u行第v列元素，u表示责任度矩阵的行数，v表示责任度矩阵的列数，S(u,v)表示层次空间相似度矩阵中的第u行第v列元素，u表示层次空间相似度矩阵的行数，v表示层次空间相似度矩阵的列数，max表示取最大值操作，E(u,g)表示可靠度矩阵中的第u行第g列元素，u表示可靠度矩阵的行数，g表示可靠度矩阵的列数，S(u,g)表示层次空间相似度矩阵中的第u行第g列元素，u表示层次空间相似度矩阵的行数，g表示层次空间相似度矩阵的列数，s.t.表示约束操作，≠表示不等于符号，表示取所有值操作，∈表示属于符号，n表示待分类的高光谱图像中样本的总数。

3.根据权利要求1所述的基于超像素张量稀疏编码的高光谱遥感地物分类方法，其特征在于，步骤(3c)所述的可靠度矩阵非对角线元素更新公式如下：

$\begin{array}{ccc}E(u,v)=min\{0,P(v,v)+\mathrm{\&Sigma;}max(0,P(k,v\left)\right)\}& s.t.& k\&NotElement;\{u,v\},\&ForAll;u,v,k\&Element;\{1,...,n\}\end{array}$

其中，E(u,v)表示可靠度矩阵中的第u行第v列元素，u表示可靠度矩阵的行数，v表示可靠度矩阵的列数，min表示取最小值操作，P(v,v)表示责任度矩阵中的第v行第v列元素，v表示责任度矩阵的行数，∑表示求和操作，max表示取最大值操作，P(k,v)表示责任度矩阵中的第k行第v列元素，k表示责任度矩阵的行数，v表示责任度矩阵的列数，s.t.表示约束操作，表示不属于符号，表示取所有值操作，∈表示属于符号，n表示待分类的高光谱图像中样本的总数。

4.根据权利要求1所述的基于超像素张量稀疏编码的高光谱遥感地物分类方法，其特征在于，步骤(3d)所述的可靠度矩阵对角线元素更新公式如下：

$\begin{array}{ccc}E(v,v)=\mathrm{\&Sigma;}max(0,P(i,v\left)\right)& s.t.& i\&NotEqual;v,\&ForAll;v,i\&Element;\{1,...,n\}\end{array}$

其中，E(v,v)表示可靠度矩阵中的第v行第v列的对角元素，v表示可靠度矩阵的行数，∑表示求和操作，max表示取最大值操作，P(i,v)表示责任度矩阵中的第i行第v列元素，i表示责任度矩阵的行数，v表示责任度矩阵的列数，s.t.表示约束操作，≠表示不等于符号，表示取所有值操作，∈表示属于符号，n表示待分类的高光谱图像中样本的总数。

5.根据权利要求1所述的基于超像素张量稀疏编码的高光谱遥感地物分类方法，其特征在于，步骤(3f)所述的聚类中心计算公式如下：

$w_s=\mathrm{argmax}\left(E\right(s,v)+P(s,v\left)\right),\&ForAll;s,v\&Element;\{1,...,n\}$

其中，w_s表示待分类的高光谱图像中第s个样本所属的聚类中心，s表示待分类的高光谱图像中样本的索引值，s的取值范围为{1,2,…,n}，arg表示求取满足条件变量值的操作，max表示最大值操作，E(s,v)表示可靠度矩阵中的第s行第v列的对角元素，s表示可靠度矩阵的行数，v表示可靠度矩阵的列数，P(s,v)表示责任度矩阵中的第s行第v列元素，s表示可靠度矩阵的行数，v表示可靠度矩阵的列数，表示取所有值操作，∈表示属于符号，n表示待分类的高光谱图像中样本的总数。

6.根据权利要求1所述的基于超像素张量稀疏编码的高光谱遥感地物分类方法，其特征在于，步骤(6a)所述的超像素分类公式如下：

$C=\mathrm{argmin}\underset{r}{\&Sigma;}\left|\right|y_r-D^m{Z}_r^m|{|}_F^2,\&ForAll;m\&Element;\{1,...,K\}$

其中，C表示待分类超像素所属的类别，arg表示求取满足条件变量值的操作，min表示最小值操作，∑表示求和操作，r表示待分类超像素中样本的索引值，r的取值范围为{1,2,…,P}，P表示待分类超像素中样本的总数，||·||_F表示Frobenius范数操作，y_r表示待分类超像素中第r个样本，y表示待分类超像素，D^m表示由标记样本字典中与第m类对应的列组成的字典，m表示标记样本字典中的类别数，m的取值范围为{1,2,…,K}，D表示标记样本字典，表示矩阵Z^m的第r列，Z^m表示待分类超像素的稀疏系数矩阵中与第m类对应的行组成的矩阵，Z表示待分类超像素的稀疏系数矩阵，表示取所有值操作，∈表示属于符号，K表示标记样本字典中的类别总数。