文章快速检索     高级检索
  北京化工大学学报(自然科学版)  2019, Vol. 46 Issue (4): 116-121   DOI: 10.13543/j.bhxbzr.2019.04.017
0

引用本文  

马方, 赵丽娜, 何磊, 杨宏伟. 基于潜在低秩图判别分析的高光谱分类[J]. 北京化工大学学报(自然科学版), 2019, 46(4): 116-121. DOI: 10.13543/j.bhxbzr.2019.04.017.
MA Fang, ZHAO LiNa, HE Lei, YANG HongWei. Hyperspectral image classification based on latent low-rank graphs[J]. Journal of Beijing University of Chemical Technology (Natural Science), 2019, 46(4): 116-121. DOI: 10.13543/j.bhxbzr.2019.04.017.

基金项目

国家自然科学基金(11301021/11571031)

第一作者

马方, 女, 1993年生, 硕士生.

通信联系人

赵丽娜, E-mail: zhaoln@mail.buct.edu.cn

文章历史

收稿日期:2018-12-13
基于潜在低秩图判别分析的高光谱分类
马方 1, 赵丽娜 1, 何磊 1, 杨宏伟 2     
1. 北京化工大学 理学院, 北京 100029;
2. 北京化工大学 信息中心, 北京 100029
摘要:提出一种基于潜在低秩图判别分析(LatLGDA)算法,利用数据的自表示对数据的列表示系数矩阵和行表示系数矩阵同时施加低秩约束,得到保留数据结构的亲和矩阵,再与图嵌入模型相结合实现高光谱图像的流形降维并进行分类。与其他基于稀疏图或稀疏低秩图的高光谱特征提取算法相比,LatLGDA可利用数据的行信息弥补列信息的不足或缺失,对噪音的抗干扰能力更强;在真实数据集上的实验结果表明,LatLGDA算法具有较高的分类精度和运算效率,应用前景广阔。
关键词稀疏图    稀疏低秩图    高光谱分类    
Hyperspectral image classification based on latent low-rank graphs
MA Fang1 , ZHAO LiNa1 , HE Lei1 , YANG HongWei2     
1. Faculty of Science, Beijing University of Chemical Technology, Beijing 100029, China;
2. Center for Information Technology, Beijing University of Chemical Technology, Beijing 100029, China
Abstract: In this paper, latent low-rank graph discrimination analysis (LatLGDA) is proposed. Our algorithm uses self-representation of the data to apply low-rank constraints to the column and row representation coefficient matrix in order to obtain the affinity matrix of the retained data structure. Combined with a graph embedding model, both manifold dimension reduction and classification of hyperspectral images can be realized. Compared with other hyperspectral feature extraction algorithms based on principles such as sparse graphs or sparse and low-rank graph discrimination analysis, LatLGDA can use the row information data to compensate for the lack of column information and has better resistance to interference from noise. Experiments on a real hyperspectral data set from the University of Pavia demonstrate that LatLGDA has the advantages of high classification accuracy, fast operation efficiency and broad application prospects.
Key words: sparse graph    sparse and low-rank graph    hyperspectral image classification    
引言

近年来,随着数据的爆炸式增长,数据量之大、维度之高和结构之复杂,都为数据分析带来了挑战。高光谱遥感是一种多维数据信息获取技术,它将成像技术和光谱技术相结合,在航空、卫星和飞机等多种平台上搭建传感器,通过远距离采集地面物体的电磁波反射数据来辨别地面物体的大小和位置信息。高光谱遥感在探测物体二维空间信息的同时还能探测其光谱信息,得到的高光谱图像具有立方体结构,因此可用x轴和y轴组成的平面表示空间信息,z轴表示光谱信息。由于包含信息丰富,高光谱图像在环境灾害检测、大气污染、农业生产及天文观测等领域都有着广泛的应用[1-2]

高光谱分辨率高、波段数足够多,能够获取地物连续的光谱曲线。由于不同的物体在光谱曲线上存在显著差异,可以依据不同地物的光谱维特征来实现地物类别分类和目标检测。高光谱遥感图像所具有的大量光谱波段提供了物体的丰富信息,一般的高光谱遥感图像都包含几百个波段,数据量为单波段的几百倍。高光谱遥感图像的光谱维数据分布在高达几百维的光谱特征空间中,相邻波段之间的高度相关性所带来的数据冗余增加了计算的复杂度,因此有必要对光谱维进行特征提取以得到低维数据,这无论在数据压缩还是数据分类方面都有着重要的意义[3]

压缩感知和稀疏表示是高光谱图像处理领域的热点研究问题。Ly等[4]最先提出稀疏图判别分析(SGDA)算法,在图嵌入模型上引入稀疏表示,同时保留降维后的低维空间样本间表示的稀疏特性。由于SGDA算法未能在低维流形空间中保留数据的全局特性,Li等[5]提出了基于稀疏和低秩图嵌入(SLRGE)算法以及基于稀疏和低秩图判别分析(SLGDA)算法,通过加入核范数约束来保留具有全局约束的样本间表示的低秩特性。但SGDA和SLGDA算法仅考虑了列空间,还不能有效应对数据样本不足的情形。本文结合潜在低秩表示[6]提出潜在低秩图判别分析(latent low-rank graph discrimination analysis,LatLGDA)算法,首先通过不精确的增广拉格朗日乘子法(IALM)[7]求解出最优的列表示系数矩阵,再与图嵌入模型相结合实现高光谱数据的特征提取,最终利用支持向量机(SVM)得到分类结果。LatLGDA算法从行空间和列空间重构数据,利用行信息弥补列信息的不足,并在低维空间里保留数据的全局低秩特性。特别是在处理被噪音污染的数据或样本数量不足的数据时,LatLGDA算法能够提取出更加丰富的数据信息,从而提升数据子空间表示能力。

1 相关算法 1.1 稀疏低秩表示

基于大多数样本可以被某些带有重要信息的样本稀疏表示的特性,SLRGE算法及SLGDA算法通过对列表示系数矩阵Wl1范数构建稀疏性来保留数据的局部结构,通过核范数构造低秩特性来保留数据的整体低秩结构。

给定数据集XRD×N, D为光谱特征维数,N为样本数目,稀疏和低秩表示模型如式(1)所示

$ \left\{ \begin{array}{l} \arg \mathop {\min }\limits_\mathit{\boldsymbol{W}} {\left\| \mathit{\boldsymbol{W}} \right\|_1} + \lambda {\left\| \mathit{\boldsymbol{W}} \right\|_ * }\\ {\rm{s}}.\;{\rm{t}}.\;\;\;\mathit{\boldsymbol{X}} = \mathit{\boldsymbol{XW}},{\mathit{\boldsymbol{W}}_{ii}} = 0 \end{array} \right. $ (1)

式中,‖·‖ 1是列表示系数矩阵Wl1范数,其值等于矩阵中每个非零元素的绝对值之和;‖·‖*W的核范数,其值为该矩阵所有奇异值之和;λ是正则化参数;Wii=0, i=1, 2, …, N是用于去掉平凡解、防止数据自表示的约束条件。式(1)中求解得到的WN×N维矩阵,其中每个N×1维列向量Wi是其余样本点对第i个点的稀疏和低秩表示的系数。

1.2 图嵌入模型

设数据集合X=[x1, x2, …, xn]∈Rd×nn个数据样本点,每个样本点分布于d维空间中。图嵌入的模型是指找到X在某个低维空间中的表示,X的图模型有n个点,每个点可以写成d维向量,其中矩阵W的每个元素Wij代表数据点xixj之间的权重,则W包含了X内所有的数据样本点之间的关联度大小和数据的流形结构[8]。对于一个线性映射PTRk×d, k < d,将数据集X映射至一个新的低维特征空间,即Y=PTX, yi=PTxiX的特征维可从d维降至k维,则最优的投影方式由目标函数式(2)表示

$ \mathop {\min }\limits_\mathit{\boldsymbol{P}} \sum\limits_{i,j} {{{\left( {{\mathit{\boldsymbol{y}}_i} - {\mathit{\boldsymbol{y}}_j}} \right)}^2}{\mathit{\boldsymbol{W}}_{ij}}} $ (2)

式中,P为要求解的最优投影矩阵。

常见的刻画数据点间关系的方法是令Wij等于两个数据点欧式距离的倒数,这样数据点相距越远,相应权重值就越小;反之亦然。

若数据点xixj之间的权重越大,则在降维时对目标函数式(2)的惩罚越大,经过P投影后的对应低维空间中的点yiyj相距越近;反之若数据点之间的权重越小,则投影后低维空间中两点间的距离越远。所以找到一个合适的矩阵W来较好地表示数据点之间的关系,使得在投影后保持这种流形结构,是图嵌入模型中最关键的步骤。

2 潜在低秩图判别分析算法

将数据本身作为字典,从列空间和行空间共同重构原始数据,同时考虑列信息和行信息,以便恢复出对噪音更加鲁棒的原始数据。对列表示系数矩阵W和行表示系数矩阵G同时施加核范数正则化约束,最终将观测数据X分解为3部分:低秩部分XW、显著特征部分GX和误差项E。目标函数如式(3)

$ \left\{ \begin{array}{l} \mathop {\min }\limits_{\mathit{\boldsymbol{W}},\mathit{\boldsymbol{G}},\mathit{\boldsymbol{E}}} {\left\| \mathit{\boldsymbol{W}} \right\|_ * } + {\left\| \mathit{\boldsymbol{G}} \right\|_ * } + \lambda {\left\| \mathit{\boldsymbol{E}} \right\|_{2,1}}\\ {\rm{s}}.\;{\rm{t}}.\;\;\mathit{\boldsymbol{X}} = \mathit{\boldsymbol{XW}} + \mathit{\boldsymbol{GX}} + \mathit{\boldsymbol{E}} \end{array} \right. $ (3)

式中,${\left\| \mathit{\boldsymbol{W}} \right\|_*} = \sum\limits_i {{\sigma _i}} $W的核范数,σiW的奇异值;${\left\| \mathit{\boldsymbol{E}} \right\|_{2, 1}} = \sum\limits_i {\sqrt {\sum\limits_j {{\bf{E}}_{ij}^2} } } $El2, 1范数,用于刻画数据的误差项,其中Eij是矩阵的非零元素。

式(3)为凸优化问题,可通过IALM算法求解。加入两个辅助矩阵J, K,则可将式(3)转化为如下同等问题

$ \left\{ \begin{array}{l} \mathop {\min }\limits_{\mathit{\boldsymbol{J}},\mathit{\boldsymbol{K}},\mathit{\boldsymbol{W}},\mathit{\boldsymbol{G}},\mathit{\boldsymbol{E}}} {\left\| \mathit{\boldsymbol{J}} \right\|_ * } + {\left\| \mathit{\boldsymbol{K}} \right\|_ * } + \lambda {\left\| \mathit{\boldsymbol{E}} \right\|_{2,1}}\\ {\rm{s}}.\;{\rm{t}}.\;\;\mathit{\boldsymbol{X}} = \mathit{\boldsymbol{XW}} + \mathit{\boldsymbol{GX}} + \mathit{\boldsymbol{E}},\mathit{\boldsymbol{W}} = \mathit{\boldsymbol{J}},\mathit{\boldsymbol{G}} = \mathit{\boldsymbol{K}} \end{array} \right. $ (4)

构造式(4)的增广拉格朗日函数如式(5)所示[7]

$ \begin{array}{l} L\left( {\mathit{\boldsymbol{J}},\mathit{\boldsymbol{K}},\mathit{\boldsymbol{W}},\mathit{\boldsymbol{G}},\mathit{\boldsymbol{E}}} \right) = {\left\| \mathit{\boldsymbol{J}} \right\|_ * } + {\left\| \mathit{\boldsymbol{K}} \right\|_ * } + \lambda \\ {\left\| \mathit{\boldsymbol{E}} \right\|_{2,1}} + \left\langle {{\mathit{\boldsymbol{Y}}_2},\mathit{\boldsymbol{W}} - \mathit{\boldsymbol{J}}} \right\rangle + \left\langle {{\mathit{\boldsymbol{Y}}_3},\mathit{\boldsymbol{G}} - \mathit{\boldsymbol{K}}} \right\rangle + \left\langle {{\mathit{\boldsymbol{Y}}_1},\mathit{\boldsymbol{X}} - } \right.\\ \left. {\mathit{\boldsymbol{XW}} - \mathit{\boldsymbol{GX}} - \mathit{\boldsymbol{E}}} \right\rangle + \frac{\mu }{2}\left\| {\mathit{\boldsymbol{W}} - \mathit{\boldsymbol{J}}} \right\|_F^2 + \frac{\mu }{2}\left\| {\mathit{\boldsymbol{G}} - \mathit{\boldsymbol{K}}} \right\|_F^2 + \frac{\mu }{2}\\ \left\| {\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{XW}} - \mathit{\boldsymbol{GX}} - \mathit{\boldsymbol{E}}} \right\|_F^2 \end{array} $ (5)

式中Y1Y2Y3为拉格朗日乘子,μ>0为惩罚参数。

对式(5)最小化来交替更新变量,即逐一更新每个变量时固定其他变量,则每个变量的更新表达式如下。

固定K, W, G, E,更新J,目标函数转化为式(6)

$ \begin{array}{l} \mathit{\boldsymbol{J}} = \mathop {\arg \min }\limits_\mathit{\boldsymbol{J}} \frac{1}{\mu }{\left\| \mathit{\boldsymbol{J}} \right\|_ * } + \frac{1}{2}\left\| {\mathit{\boldsymbol{J}} - \left( {\mathit{\boldsymbol{W}} + \frac{{{\mathit{\boldsymbol{Y}}_2}}}{\mu }} \right)} \right\|_F^2 = \\ {{\bf{\Theta }}_{\frac{1}{\mu }}}\left( {\mathit{\boldsymbol{W}} + \frac{{{\mathit{\boldsymbol{Y}}_2}}}{\mu }} \right) \end{array} $ (6)

固定J, W, G, E,更新K,目标函数转化为式(7)

$ \begin{array}{l} \mathit{\boldsymbol{K}} = \mathop {\arg \min }\limits_\mathit{\boldsymbol{K}} \frac{1}{\mu }{\left\| \mathit{\boldsymbol{K}} \right\|_ * } + \frac{1}{2}\left\| {\mathit{\boldsymbol{K}} - \left( {\mathit{\boldsymbol{G}} + \frac{{{\mathit{\boldsymbol{Y}}_3}}}{\mu }} \right)} \right\|_F^2 = \\ {{\bf{\Theta }}_{\frac{1}{\mu }}}\left( {\mathit{\boldsymbol{G}} + \frac{{{\mathit{\boldsymbol{Y}}_3}}}{\mu }} \right) \end{array} $ (7)

式中${\Theta _{\frac{1}{\mu }}}(\cdot)$为奇异值阈值算子[9]

固定J, K, G, E,更新W,目标函数转化为式(8)

$ \begin{array}{l} L\left( \mathit{\boldsymbol{W}} \right) = \frac{\mu }{2}\left( {\left\| {\mathit{\boldsymbol{W}} - \mathit{\boldsymbol{J}} + \frac{{{\mathit{\boldsymbol{Y}}_2}}}{\mu }} \right\|_F^2 + } \right.\\ \left. {\left\| {\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{XW}} - \mathit{\boldsymbol{GX}} - \mathit{\boldsymbol{E}} + \frac{{{\mathit{\boldsymbol{Y}}_1}}}{\mu }} \right\|_F^2} \right) \end{array} $ (8)

式(8)中只包含矩阵的Frobenius范数,因此最小化时可直接对W求导,整理可得

$ \begin{array}{l} \mathit{\boldsymbol{W}} = {\left( {{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{X}} + \mathit{\boldsymbol{I}}} \right)^{ - 1}}\left( {{\mathit{\boldsymbol{X}}^{\rm{T}}}\left( {\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{GX}} - \mathit{\boldsymbol{E}}} \right) + \mathit{\boldsymbol{J}} + } \right.\\ \left. {\frac{{{\mathit{\boldsymbol{X}}^{\rm{T}}}{\mathit{\boldsymbol{Y}}_1} - {\mathit{\boldsymbol{Y}}_2}}}{\mu }} \right) \end{array} $ (9)

固定J, K, W, E,更新G,目标函数转化为式(10)

$ \begin{array}{l} L\left( \mathit{\boldsymbol{G}} \right) = \frac{\mu }{2}\left\| {\mathit{\boldsymbol{G}} - \mathit{\boldsymbol{K}} + \frac{{{\mathit{\boldsymbol{Y}}_3}}}{\mu }} \right\|_F^2 + \frac{\mu }{2}\left\| {\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{XW}} - } \right.\\ \left. {\mathit{\boldsymbol{GX}} - \mathit{\boldsymbol{E}} + \frac{{{\mathit{\boldsymbol{Y}}_1}}}{\mu }} \right\|_F^2 \end{array} $ (10)

最小化式(10)时可直接对G求导,整理可得更新表达式如式(11)

$ \begin{array}{l} \mathit{\boldsymbol{G}} = \left( {\left( {\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{XW}} - \mathit{\boldsymbol{E}}} \right){\mathit{\boldsymbol{X}}^{\rm{T}}} + \mathit{\boldsymbol{K}} + \frac{{{\mathit{\boldsymbol{Y}}_1}{\mathit{\boldsymbol{X}}^{\rm{T}}} - \mathit{\boldsymbol{Y}}3}}{\mu }} \right)\\ {\left( {\mathit{\boldsymbol{X}}{\mathit{\boldsymbol{X}}^{\rm{T}}} + \mathit{\boldsymbol{I}}} \right)^{ - 1}} \end{array} $ (11)

固定J, K, W, G,更新E,目标函数转化为式(12)

$ \begin{array}{l} \mathit{\boldsymbol{E}} = \mathop {\arg \min }\limits_\mathit{\boldsymbol{E}} \frac{\lambda }{\mu }{\left\| \mathit{\boldsymbol{E}} \right\|_{2,1}} + \frac{1}{2}\left\| {\mathit{\boldsymbol{E}} - \left( {\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{XW}} - } \right.} \right.\\ \left. {\left. {\mathit{\boldsymbol{GX}} + \frac{{{\mathit{\boldsymbol{Y}}_1}}}{\mu }} \right)} \right\|_F^2 \end{array} $ (12)

式(12)的最小化可由引理1得到。

引理1[10]  当式(13)最优化问题中的H为已知矩阵时,最优解E*的第i列元素值如式(14)所示

$ \mathop {\arg \min }\limits_\mathit{\boldsymbol{E}} \alpha {\left\| \mathit{\boldsymbol{E}} \right\|_{2,1}} + \frac{1}{2}\left\| {\mathit{\boldsymbol{E}} - \mathit{\boldsymbol{H}}} \right\|_F^2 $ (13)
$ {\left[ {{\mathit{\boldsymbol{E}}^ * }} \right]_{:,i}} = \left\{ \begin{array}{l} \frac{{{{\left\| {{\mathit{\boldsymbol{H}}_{:,i}}} \right\|}_2} - \alpha }}{{{{\left\| {{\mathit{\boldsymbol{H}}_{:,i}}} \right\|}_2}}}{\mathit{\boldsymbol{H}}_{:,i}},\;\;\;\;\;\;\;\;\;{\left\| {{\mathit{\boldsymbol{H}}_{:,i}}} \right\|_2} > \alpha \\ 0,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;其他 \end{array} \right. $ (14)

通过式(15)和式(16)更新Y1, Y2, Y3μ

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{Y}}_1} = {\mathit{\boldsymbol{Y}}_1} + \mu \left( {\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{XW}} - \mathit{\boldsymbol{GX}} - \mathit{\boldsymbol{E}}} \right)\\ {\mathit{\boldsymbol{Y}}_2} = {\mathit{\boldsymbol{Y}}_2} + \mu \left( {\mathit{\boldsymbol{W}} - \mathit{\boldsymbol{J}}} \right)\\ {\mathit{\boldsymbol{Y}}_3} = {\mathit{\boldsymbol{Y}}_3} + \mu \left( {\mathit{\boldsymbol{G}} - \mathit{\boldsymbol{K}}} \right) \end{array} \right. $ (15)
$ \mu = \min \left( {\rho \mu ,{{\max }_\mu }} \right) $ (16)

式(16)中ρ>1是标量,用于更新μ的值,根据IALM算法,ρ取1.5。

对目标函数式(3)采用IALM算法求解的具体流程如算法1所示。

算法1  利用IALM求解式(3)

1) 输入  XRd×n, λ>0

2) 初始化  W0, J0, K0, G0, E0=0, μ0=10-6

$ {\max _\mu } = {10^6},\varepsilon = {10^{ - 6}},\rho = 1.5 $

3) 迭代  (t=0, 1, 2, …)

固定其他变量,通过式(6)更新J

固定其他变量,通过式(7)更新K

固定其他变量,通过式(9)更新W

固定其他变量,通过式(11)更新G

固定其他变量,通过式(12)更新E

由式(15)更新Y1, Y2, Y3

由式(16)更新μ

判断迭代结束条件为

$ \left\{ \begin{array}{l} {\left\| {\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{XW}} - \mathit{\boldsymbol{GX}} - \mathit{\boldsymbol{E}}} \right\|_\infty } < \mathit{\boldsymbol{\varepsilon }}\\ {\left\| {\mathit{\boldsymbol{W}} - \mathit{\boldsymbol{J}}} \right\|_\infty } < \varepsilon ,{\left\| {\mathit{\boldsymbol{G}} - \mathit{\boldsymbol{K}}} \right\|_\infty } < \varepsilon \end{array} \right. $

4) 输出  最优解W, G, E

应用算法1求解出W之后,通过基于图嵌入模型找到一个最优的投影矩阵P,经过投影变换后的低维数据表示为Y=PTX。为了在低维空间中保留数据的原有流形结构,根据1.1节的图嵌入模型式(2)得到目标函数

$ \left\{ \begin{array}{l} \begin{array}{*{20}{c}} {\mathop {\arg \min }\limits_\mathit{\boldsymbol{P}} \sum\limits_{i = 1}^n {{{\left\| {{\mathit{\boldsymbol{P}}^{\rm{T}}}{\mathit{\boldsymbol{x}}_i} - \sum\limits_{j = 1}^n {{\mathit{\boldsymbol{W}}_{ij}}{\mathit{\boldsymbol{P}}^{\rm{T}}}{\mathit{\boldsymbol{x}}_j}} } \right\|}^2}} = }\\ {\mathop {\arg \min }\limits_\mathit{\boldsymbol{P}} tr\left( {{\mathit{\boldsymbol{P}}^{\rm{T}}}\mathit{\boldsymbol{XL}}{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{P}}} \right)} \end{array}\\ {\rm{s}}.\;{\rm{t}}.\;\;\;{\mathit{\boldsymbol{P}}^{\rm{T}}}\mathit{\boldsymbol{X}}{\mathit{\boldsymbol{L}}_\mathit{\boldsymbol{P}}}{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{P}} = \mathit{\boldsymbol{I}} \end{array} \right. $ (17)

式中,L=D-W,其中n×n维对角矩阵DW的度矩阵,LW的拉普拉斯矩阵;I=LP为单位矩阵;PTXLPXTP=I为拉格朗日约束条件。

式(17)的求解过程是广义特征值分解问题[11],即XLXTP=VXLPXTP,其中V为对角矩阵,对角线上每一个元素对应一个广义特征值,投影矩阵P由对应的特征值向量组成。

综上,本文提出的LatLGDA算法用于高光谱图像分类的具体步骤如下。

1) 数据预处理。输入高光谱数据集,先对数据进行归一化运算,再将归一化后的三维高光谱数据转化为二维数据矩阵,矩阵的行数等于光谱的波段数,列数等于样本点数,则矩阵每列表示每个样本点的光谱特征。

2) 选择训练集和测试集。从归一化后的数据矩阵的每一类别中随机选取一定数量数据组成训练样本集,剩余作为测试样本集。

3) 通过算法1求解最优列表示系数矩阵W

4) 求解最优投影矩阵P。通过W计算出拉普拉斯矩阵L,再根据广义特征值求解最优投影矩阵P,则训练集降维后的数据变成Y=PTX

5) 高光谱图像分类。利用SVM对降维后的训练集和测试集进行分类。

3 实验分析

为了说明本文提出算法的有效性,选取真实的高光谱图像数据集University of Pavia[12]进行高光谱分类实验,并与不同图判别分析方法的结果进行对比。计算测试样本每一类别的分类准确率和总体分类准确率作为各个算法的评价指标。仿真实验环境为Matlab R2016b版,机器配置为ThinkCentre M8400t-N000,内存4 GB,频率3.4 GHz,操作系统Windows 7。

University of Pavia实验数据集由ROSIS传感器采集,图像场景是位于意大利Pavia城的一所大学。此高光谱图像共有610×340个像元,波段数103,空间分辨率1.3 m,其中有标签的像元为42 776个,分别属于9个类别。在分类实验中,首先选取每一类的5%作为训练集,9个类别的样本信息如表 1所示。为保证实验公平性,不同算法都先通过同样的训练样本数据求解各自的最优投影矩阵进行降维,得到降维后的训练数据集与测试数据集,再通过SVM分类来比较这些算法特征提取的效果。

下载CSV 表 1 University of Pavia数据集选择的训练集和测试集 Table 1 Training and test samples selected from University of Pavia data set

SVM工具选用台湾大学林智仁开发的LIBSVM工具包。惩罚系数c取值范围0.000 1~10 000,c越大对错误例惩罚程度就越大,本文实验中所有算法均取c=10 000;选择核函数类型为高斯径向基核函数(RBF)时的参数σ,分类时σ从{0.01, 0.05, 0.5, 1, 5, 10, 50, 100}中依次循环取值进行计算,最终的分类结果取8个值中最高者。

图 1为降维数变化时,LatLGDA算法与已有算法的分类精确度对比。从图中可以看出,随着最终降维数从10增加到50,本文算法的分类精确度始终高于SGDA、SLRGE和SLGDA。然后统一将光谱特征维度降至22,得到不同方法对每一类别和总体的分类精确度如表 2所示。可以看出本文算法的类别和总体分类精确度均为最高,其中总体分类精确度分别比SGDA、SLRGE和SLGDA方法高8.9%,5.12%和3.58%。

图 1 降维后不同维度对分类精确度的影响 Fig.1 The effect of dimensionality on classification accuracy
下载CSV 表 2 不同算法的分类精确度 Table 2 Classification accuracy of different algorithms

实验中选取5%的训练集时,样本数小于高光谱维度的第5、7和9类别出现小样本问题,这与文献[4]中选取8%训练集时SLGDA算法的分类结果相似,即第7类别的分类精确度是9个类别中最低的,这是由于训练样本少造成的。表 2中本文算法关于第7类的分类精确度最低,为65.64%,但仍比SGDA和SLGDA方法分别高出29.95%和9.47%,说明本文方法利用行信息弥补列信息的不足,能更好地应对小样本问题。

为验证本文算法在训练样本数目很少时的分类效果,分别选取每一类的1%、2%、3%、4%、5%作为训练集,分类结果如图 2所示。可以看出随着训练集的增大,SGDA、SLRGE、SLGDA和LatLGDA的分类精确度都有所提升,其中LatLGDA分类效果最好,SLRGE和SLGDA次之,SGDA最差。

图 2 训练样本数目变化对分类精确度的影响 Fig.2 The influence of the training ratio on the classification

为了从计算时间和效率方面比较不同算法的差异,在训练样本占比1%、2%、3%、4%、5%条件下,记录不同算法的运算时间如图 3所示。可以看出,随着训练比率的提高,各算法的计算时间都随之增加,但比较来看,本文算法计算时间最短、效率最高,更适用于真实高光谱数据分类。

图 3 不同算法的运算时间比较 Fig.3 Computation time of different algorithms
4 结束语

将潜在低秩表示和图嵌入模型相结合提出一种高光谱图像光谱维特征提取算法,该算法同时考虑行信息和列信息,能够恢复出对噪音更加鲁棒的原始数据;与其他基于图判别分析的高光谱分类算法相比,本文算法对噪音更加鲁棒,分类精度更高,计算效率更快,并且在训练样本较少时仍具有较高的分类准确率。

参考文献
[1]
赵春晖, 王立国, 齐滨. 高光谱遥感图像处理方法及应用[M]. 北京: 电子工业出版社, 2016: 369-381.
ZHAO C H, WANG L G, QI B. Hyperspectral remote sensing image processing methods and applications[M]. Beijing: Electronics Industry Press, 2016: 369-381. (in Chinese)
[2]
CHIVASA W, MUTANGA O, BIRADAR C. Application of remote sensing in estimating maize grain yield in heterogeneous African agricultural landscapes:a review[J]. International Journal of Remote Sensing, 2017, 38(23): 6816-6845. DOI:10.1080/01431161.2017.1365390
[3]
YANG H, DU Q, CHEN G S. Particle swarm optimization-based hyperspectral dimensionality reduction for urban land cover classification[J]. IEEE Journal of Selected Topics in Applied Earth Observations & Remote Sensing, 2012, 5(2): 544-554.
[4]
LY N H, DU Q, FOWLER J E. Sparse graph-based discriminant analysis for hyperspectral imagery[J]. IEEE Transactions on Geoscience & Remote Sensing, 2014, 52(7): 3872-3884.
[5]
LI W, LIU J B, DU Q. Sparse and low-rank graph for discriminant analysis of hyperspectral imagery[J]. IEEE Transactions on Geoscience & Remote Sensing, 2016, 54(7): 4094-4105.
[6]
LIU G C, YAN S C. Latent low-rank representation for subspace segmentation and feature extraction[C]//Proceedings of the 2011 IEEE International Conference on Computer Vision. Barcelona, 2011: 1615-1622.
[7]
LIN Z C, CHEN M M, MA Y. The augmented lagrange multiplier method for exact recovery of corrupted low-rank matrices[J/OL]. (2013-10-18)[2018-11-02]. https: //arXiv.org/abs/1009.5055.
[8]
YAN S C, XU D, ZHANG B Y, et al. Graph embedding and extensions:a general framework for dimensionality reduction[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2007, 29(1): 40-51. DOI:10.1109/TPAMI.2007.250598
[9]
CAI J F, CANDES E J, SHEN Z. A singular value thresholding algorithm for matrix completion[J]. SIAM Journal on Optimization, 2010, 20(4): 1956-1982. DOI:10.1137/080738970
[10]
YANG J F, YIN W T, ZHANG Y, et al. A fast algorithm for edge-preserving variational multichannel image restoration[J]. SIAM Journal on Imaging Sciences, 2009, 2(2): 569-592. DOI:10.1137/080730421
[11]
GUATTERY S, MILLER G L. Graph embeddings and laplacian eigenvalues[J]. SIAM Journal on Matrix Analysis and Applications, 2000, 21(3): 703-723. DOI:10.1137/S0895479897329825
[12]
PENG B, LI W, XIE X M, et al. Weighted-fusion-based representation classifiers for hyperspectral imagery[J]. Remote Sensing, 2015, 7(11): 14806-14826. DOI:10.3390/rs71114806