基于递归自动编码的高光谱特征学习方法

摘要

本发明属于图像处理技术领域，具体提供了一种基于递归自动编码的高光谱特征学习方法，其步骤包括：将高光谱图像中每个像素点用特征向量表示，进行归一化处理，选出训练集和测试集；对每个特征向量构建邻域窗口块；在训练集中，每个特征向量及其邻域根据特征向量间重构误差最小化准则，合并特征向量对，再根据整个网络的重构误差对整个网络进行训练；分别将训练集和测试集输入到训练好的网络中，得到新的训练集和测试集；将新的训练集和测试机输入到支撑矢量机进行分类，得到分类结果。本发明采用了无监督的思想，克服了高光谱数据获取类标难的问题，获得较高的分类正确率，可用于矿物勘探，环境管理，军事防御领域。

说明书

技术领域

本发明属于图像处理技术领域，具体涉及遥感图像技术领域中的一种 基于递归自动编码的高光谱特征学习方法。

背景技术

高光谱图像技术在近几年来有很快的发展，其研究主要致力于寻找使 计算机智能地学习和识别高光谱图像真实地物类的技术方法。高光谱图像 在城市规划、环境检测、植被分类、军事目标探测以及矿物地质识别等诸 多方面都有着巨大的应用前景。一种普遍的高光谱图像识别方法通常是： 首先从高光谱图像中得到每个像素的光谱特征，通过对光谱特征的提取得 到高级特征，并在此基础上对高光谱图像进行分类。其中的一个关键问题 就是如何对不同的高光谱图像进行恰当的表示，因为表示的合理与否决定 了后续分类的性能上限。另外，由于高光谱图像具有数据量大、冗余信息 多、含有噪声等不利因素，因此要求在对光谱的特征学习时用到的技术方 法高效、简单且有一定抗噪声干扰能力。

重庆大学申请的专利“一种基于相对梯度的ISA模型的行为识别方法” (申请号：200910250850.6，公开号：CN101770584A，公布日：2010.07.07) 公开了一种高光谱特征提取的方法，该方法的具体步骤为：第1步：部分 样本数据点标注类标；第2步：相似图和相异图；第3步：计算投影矩阵； 第4步：利用投影矩阵将高光谱数据投影到低维数据。该方法虽然能够对 高光谱图像进行有效的特征提取，但是，仍然存在的不足之处是：该方法 是基于部分类标信息的，而对于高光谱来说，类标信息比较困难得到，会 消耗大量的时间，这将会加大了时间复杂度。

PCA(主成分分析)与LDA(线性判别分析)是最常用的两种特征提 取方法。PCA以最大方差为准则，能够在均方误差最小意义下最优表达数 据，却不适合于分类。LDA以最大类间离散度矩阵于类内离散矩阵之比为 准则，通过求解广义瑞利商问题获得投影矩阵，但LDA算法在小样本问题 会使算法失效。虽然PCA和LDA方法提取特征可以减少算法复杂度，但 是存在的不足之处是：这两种方法都是线性特征提取方法，对于非线性问 题性能可能会下降。

Richard Socher等人在论文“Recursive Autoencoders”(EMNLP，2011) 中提出Semi-unsupervised Recursive Autoencoders方法，该方法的具体步骤 为第1步：提取样本的初级特征，就是样本的表示；第2步：构建半监督 递归自编码网络；第3步：训练网络权值；第4步：把样本的初级特征带 入网络得到高级特征。这种方法虽然利用了样本结构特性，使得对样本高 级特征相对准确，且计算复杂度较低，但是，仍然存在的不足之处是：该 方法只用于了自然语言，其他领域(尤其是高光谱领域)并没有加以利用。

发明内容

本发明的目的是克服现有高光谱特征提取技术中类标比较难找的问 题。

为此，本发明提供了一种基于递归自动编码的高光谱特征学习方法， 包括以下步骤：

(1)输入高光谱遥感图像数据，每个像素即样本用光谱特征向量表示， 样本的特征维数为d，将样本集合归一化到0～1之间，所有归一化的样本 构成样本集，其中x_i为第i个样本，N为样本总个数，表示实数域；

(2)从样本集X中选出一定比例的样本作为训练集，剩余的样本作为测 试集；

(3)构建每个样本的邻域窗口块：在经过归一化样本集上，以各个样本 为中心，取其m×m邻域窗口内的所有样本，得到包含m²个样本的集合Y， 其中m为奇数；

(4)利用训练集及其邻域样本，训练递归自动编码网络：

4a)初始化网络参数θ＝(W₁,W₂,b₁,b₂)，W₁和b₁是输入层到隐层的权重和偏 置，W₂和b₂是隐层到重构层的权重和偏置，初始化b₁和b₂分别为全零d维和 2d维向量，W₁＝rand(2d,d)×2e-e，W₂＝rand(d,2d)×2e-e；其中，rand(2d,d)和rand(d,2d)分别是从0到1间随机取值而构成的2d×d矩阵和 d×2d矩阵；

4b)以一个训练样本及其邻域样本构成的集合Y作为输入，连接任意两 个样本的特征向量构成新特征向量z,根据自动编码方法计算z的d 维隐层特征向量z′,及其重构误差；

4c)合并重构误差最小的样本对，利用其隐层特征z′代替该样本对，此 时集合Y中的样本个数减1，重新计算Y中每对样本的隐层特征向量和重 构误差，重复上述过程，直至Y中仅包含一个样本为止，完成该训练样本 树结构递归网络拓扑的构建；

4d)按4b)和4c)中的步骤对每个训练样本构建树结构递归网络拓扑；

4e)以重构误差最小化作为目标函数，运用后向传播算法对训练集构建 的树结构递归网络参数θ求梯度，再用L-BFGS方法对树结构递归网络参数 θ进行训练；

(5)利用上述训练好的树结构递归网络提取训练样本和测试样本的新特 征向量，将每个样本及其邻域样本构成的集合Y输入到已训练好的树结构 递归网络，获得所有隐层特征向量，构成集合Q，计算Y和Q中所有特征 向量的平均值，作为该样本的新特征向量，从而得到新的训练集和测试集；

(6)分类：对新训练集和测试集，利用支撑矢量机进行分类；

(7)输出分类结果。

上述步骤4b)中自动编码方法的具体步骤为：

第1步，计算输入样本的隐层特征向量：z′＝f(W₁z+b₁)，其中， W₁,表示权值，b₁,表示偏置,f(·)表示tanh函数， $\tanh \left(a\right)=\frac{(e^a-e^a)}{(e^a+e^a)},$其中a＝W₁z+b₁；

第2步，利用隐层特征向量重构输入特征向量，计算公式如下： 其中，表示重构后的特征向量,W₂为d×2d维的权值矩 阵，b₂为一个2d维偏置向量；

第3步，计算重构误差

上述步骤4e)L-BFGS方法的目标函数为迭代公式为 θ^(k+1)＝θ^(k)+α^(k)d^(k)，其中，k为迭代次数，α^(k)为步长，d^(k)为搜索方向，具 体过程为：

第1步，给定初始化的网络参数θ＝(W₁,W₂,b₁,b₂)，转换成列向量 θ′＝{W₁(:)；W₂(:)；b₁；b₂}，其中W₁(:)和W₂(:)分别表示矩阵W₁和W₂按列扫描 成的3d维列向量，n＝2d+d+d+2d+d+2d＝9d，初始化对称正定矩阵为单位方阵，给定需要存储的迭代向量数r＝70，误差ε>0为一个很小的数， 令k＝0；

第2步，若||g^(k)||≤ε，则算法终止，否则，令d^(k)＝-B^(k)g^(k)，其中 $g^{\left(k\right)}=▿E\left({\theta }^{\left(k\right)}\right),$$▿E\left({\theta }^{\left(k\right)}\right)$是E对θ^(k)求梯度；

第3步，根据如下Wolfe-Powell条件确定步长α^(k)>0： E(θ^(k)+α^(k)d^(k))≤E(θ^(k))+δα^(k)(g^(k))^Td^(k),，(g^(k+1))^Td^(k)≥σ(g^(k))^Td^(k)；

其中θ^(k+1)＝θ^(k)+α^(k)d^(k)，σ∈(δ,1)，T表示矩正转置；

第4步，令r′＝min{k+1,r}，取B^(k)＝((s^(k))^Ty^(k)/||y^(k)||²)I，其中I为单位方 阵，T表示矩正转置，由以下公式确定B^(k+1)：

$\left(\begin{array}{c}{B}^{(k+1)}=({\left(V^{\left(k\right)}\right)}^T···{\left(V^{(k-r^{\prime })}\right)}^T)B^{\left(k\right)}(\left(V^{(k-r^{\prime })}\right)···\left(V^{\left(k\right)}\right))+\underset{j=0}{\overset{r^{\prime }}{\Sigma }}{\rho }^{(k-r^{\prime }+j)}\left(\underset{l=0}{\overset{r^{\prime }-j-1}{\Pi }}{\left(V^{(k-l)}\right)}^T\right)\\ \times S^{(k-r^{\prime }+j)}{\left(S^{(k-r^{\prime }+j)}\right)}^T\left(\underset{l=r^{\prime }-j-1}{\overset{0}{\Pi }}{V}^{(k-l)}\right)\end{array}\right)$

其中s^(k)＝θ^(k+1)-θ^(k)，ρ^(k)＝1/(s^(k))^Ty^(k),y^(k)＝g^(k+1)-g^(k),V^(k)＝(I-ρ^(k)y^(k)(s^(k))^T)；

第5步，令k＝k+1，转第2步。

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

第一，本发明把深度学习算法应用于高光谱领域，更重要的是，整个 特征提取的过程是无监督的，这对于高光谱图像来说有重大的意义。因为 高光谱数据中无类标数据很容易得到，但有类标数据集则非常困难得到， 并且造成巨大的时间消耗和经济消费，而本发明对高光谱数据集，克服了 现有技术中需要类标信息来进行特征提取的困难，使得本发明具有更经济 和容易实现的优点。

第二，本发明在特征提取阶段基本上是无参数的，除了在训练网络权 值求梯度时，需要一个控制权值正则的参数，这克服了现有技术中参数调 节而带来的复杂度高的不足，提高了特征提取的计算效率，使得本发明具 有体征提取速度更快的优点。

第三，本发明利用了高光谱图像中空谱结构信息，再利用深度学习的 方法，克服了深度学习方法没有应用于高光谱领域的缺点，使得本发明具 有对高光谱判别性高和识别精度高的优点。

以下将结合附图对本发明做进一步详细说明。

附图说明

图1为本发明方法的流程图示意图。

图2(a)为本发明仿真实验中Indian Pines高光谱数据图像。

图2(b)为本发明仿真实验中Pavia University高光谱数据图像。

图3为本发明方法自动编码网络图。

图4为本发明方法递归自动编码的结构图。

具体实施方式

实施例1：

为了克服现有高光谱特征提取技术中类标比较难找的问题，本实施例 提供了一种如图1，图3和图4所示的基于递归自动编码的高光谱特征学习 方法，包括以下步骤：

(1)输入高光谱遥感图像数据，每个像素即样本用光谱特征向量表示， 样本的特征维数为d，将样本集合归一化到0～1之间，所有归一化的样本 构成样本集，其中x_i为第i个样本，N为样本总个数，表示实数域；

(2)从样本集X中选出一定比例的样本作为训练集，剩余的样本作为测 试集；

(3)构建每个样本的邻域窗口块：在经过归一化样本集上，以各个样本 为中心，取其m×m邻域窗口内的所有样本，得到包含m²个样本的集合Y， 其中m为奇数；

(4)利用训练集及其邻域样本，训练递归自动编码网络：

4a)初始化网络参数θ＝(W₁,W₂,b₁,b₂)，W₁和b₁是输入层到隐层的权重和偏 置，W₂和b₂是隐层到重构层的权重和偏置，初始化b₁和b₂分别为全零d维和 2d维向量，W₁＝rand(2d,d)×2e-e，W₂＝rand(d,2d)×2e-e。其中，rand(2d,d)和rand(d,2d)分别是从0到1间随机取值而构成的2d×d矩阵和 d×2d矩阵；

4b)以一个训练样本及其邻域样本构成的集合Y作为输入，连接任意两 个样本的特征向量构成新特征向量z,根据自动编码方法计算z的d 维隐层特征向量z′,及其重构误差；

4c)合并重构误差最小的样本对，利用其隐层特征z′代替该样本对，此 时集合Y中的样本个数减1，重新计算Y中每对样本的隐层特征向量和重 构误差，重复上述过程，直至Y中仅包含一个样本为止，完成该训练样本 树结构递归网络拓扑的构建；

4d)按4b)和4c)中的步骤对每个训练样本构建树结构递归网络拓扑；

4e)以重构误差最小化作为目标函数，运用后向传播算法对训练集构建 的树结构递归网络参数θ求梯度，再用L-BFGS方法对树结构递归网络参数 θ进行训练；

(5)利用上述训练好的树结构递归网络提取训练样本和测试样本的新 特征向量：将每个样本及其邻域样本构成的集合Y输入到已训练好的树结 构递归网络，获得所有隐层特征向量，构成集合Q，计算Y和Q中所有特 征向量的平均值，作为该样本的新特征向量，从而得到新的训练集和测试 集；

(6)分类：对新训练集和测试集，利用支撑矢量机进行分类；

(7)输出分类结果。

根据上述基于递归自动编码的高光谱特征学习方法，结合图3，步骤(4) 中所述自动编码方法的具体步骤为：

第1步，计算输入样本的隐层特征向量：z′＝f(W₁z+b₁)，其中， W₁,表示权值，b₁,表示偏置,f(·)表示tanh函数， $\tanh \left(a\right)=\frac{(e^a-e^a)}{(e^a+e^a)},$其中a＝W₁z+b₁。

第2步，利用隐层特征向量重构输入特征向量，计算公式如下： 其中，表示重构后的特征向量,W₂为d×2d维的权值矩 阵，b₂为一个2d维偏置向量。

第3步，计算重构误差：

另外，根据所述的基于递归自动编码的高光谱特征学习方法步骤，结 合图4，步骤(4)中所述递归自动编码方法步骤为：

第1步，一个样本的邻域窗口块中的所有样本形成集合Y输入，假设窗 口大小为3×3，即集合Y为[x₁；x₃；x₂；…x₉]，其中所有样本对利用自动编码方 法计算出这对样本的隐层特征(隐层特征向量维数和样本的维数相同)，同 时得到重构误差。

第2步，把重构误差最小的样本对用隐层特征向量表示，并且此隐层 特征向量完全代替此样本对。假设样本对{x₁；x₂}是重建错误最小的，其得到 的隐层特征向量y₁，则这集合Y变成[y₁；x₃；x₄；…x₉]，重复这样的过 程，直到集合Y只剩下一个样本，这样递归地构建了一个树结构递归网络。

第3步，根据后向传播算法对树结构递归网络参数θ＝(W₁,W₂,b₁,b₂)求梯 度，再用L-BFGS算法训练这个树结构递归网络的网络参数θ＝(W₁,W₂,b₁,b₂)。

第4步，一个样本Y集合输入到这个训练好的树结构递归网络的网络， 然后用所有树结构结点的平均值来表示这个样本，从而提取到该样本的新 特征。

根据上述所述的基于递归自动编码的高光谱特征学习方法步骤，结合 图3，步骤(4)中所述反向传播算法对网络参数求梯度，其具体过程为：

第1步，重复步骤4a)，初始化网络参数θ＝(W₁,W₂,b₁,b₂)。

第2步，计算隐层特征向量z′和重构的特征向量。

第3步，计算其中λ为正则项，具体参数的梯度为：

$\frac{\partial E}{{\partial W}_2\left({\mathrm{jp}}^{\prime }\right)}=({\hat{z}}_{p^{\prime }}-z_{p^{\prime }})f^{\prime }\left(a_{p^{\prime }}^{\prime }\right)z_j^{\prime },\frac{\partial E}{{\partial W}_1\left(\mathrm{pj}\right)}=\underset{p^{\prime }=1}{\overset{2d}{\Sigma }}\left(({\hat{z}}_{p^{\prime }}-z_{p^{\prime }})f^{\prime }\left(a_{p^{\prime }}^{\prime }\right)z_j^{\prime }{W}_2\left({\mathrm{jp}}^{\prime }\right)\right)f^{\prime }\left(a_j\right)z_p,$

$\frac{\partial E}{{\partial b}_2\left(p^{\prime }\right)}=({\hat{z}}_{p^{\prime }}-z_{p^{\prime }})f^{\prime }\left(a_{p^{\prime }}^{\prime }\right),\frac{\partial E}{{\partial b}_1\left(j\right)}=\underset{p^{\prime }=1}{\overset{2d}{\Sigma }}\left(({\hat{z}}_{p^{\prime }}-z_{p^{\prime }})f^{\prime }\left(a_{p^{\prime }}^{\prime }\right)z_j^{\prime }{W}_2\left({\mathrm{jp}}^{\prime }\right)\right)f^{\prime }\left(a_j\right)$

其中f′(·)代表tanh函数的导数，j＝1,2,…,d，p＝1,2,…,2d，p′＝1,2,…,2d， W₂(jp′)为矩阵W₂中第j行p′列的值，W₁(pj)为矩阵W₁中第p行j列的值，为 向量z′第j维的值，z_p为向量z第p维的值，z_p′为向量z第p′维的值，为 向量第p′维的值，b₁(j)为向量b₁第j维的值，b₂(p′)为向量b₂第p′维的值， $a_{p^{\prime }}^{\prime }=\underset{j=1}{\overset{d}{\Sigma }}{W}_2\left({\mathrm{jp}}^{\prime }\right)z_j^{\prime },$$a_j=\underset{p=1}{\overset{2d}{\Sigma }}{W}_1\left(\mathrm{pj}\right)z_p.$

根据上面所述的基于递归自动编码的高光谱特征学习方法步骤，步骤(4) 中所述L-BFGS方法的目标函数为：迭代公式为 θ^(k+1)＝θ^(k)+α^(k)d^(k)，其中，k为迭代次数，α^(k)为步长，d^(k)为搜索方向；具 体过程为：

第1步，给定初始化的网络参数θ＝(W₁,W₂,b₁,b₂)，转换成列向量 θ′＝{W₁(:)；W₂(:)；b₁；b₂}，其中W₁(:)和W₂(:)分别表示矩阵W₁和W₂按列扫描 成的3d维列向量，n＝2d+d+d+2d+d+2d＝9d，初始化对称正定矩阵为单位方阵，给定需要存储的迭代向量数r＝70，误差ε>0为一个很小的数， 令k＝0。

第2步，若||g^(k)||≤ε，则算法终止，否则，令d^(k)＝-B^(k)g^(k)，其中 $g^{\left(k\right)}=▿E\left({\theta }^{\left(k\right)}\right),$$▿E\left({\theta }^{\left(k\right)}\right)$是E对θ^(k)求梯度。

第3步，根据如下Wolfe-Powell条件确定步长α^(k)>0： E(θ^(k)+α^(k)d^(k))≤E(θ^(k))+δα^(k)(g^(k))^Td^(k),，(g^(k+1))^Td^(k)≥σ(g^(k))^Td^(k)；

其中θ^(k+1)＝θ^(k)+α^(k)d^(k)，σ∈(δ,1)，T表示矩正转置。

第4步，令r′＝min{k+1,r}，取B^(k)＝((s^(k))^Ty^(k)/||y^(k)||²)I，其中I为单位方 阵，T表示矩正转置，由以下公式确定B^(k+1)：

$\left(\begin{array}{c}{B}^{(k+1)}=({\left(V^{\left(k\right)}\right)}^T···{\left(V^{(k-r^{\prime })}\right)}^T)B^{\left(k\right)}(\left(V^{(k-r^{\prime })}\right)···\left(V^{\left(k\right)}\right))+\underset{j=0}{\overset{r^{\prime }}{\Sigma }}{\rho }^{(k-r^{\prime }+j)}\left(\underset{l=0}{\overset{r^{\prime }-j-1}{\Pi }}{\left(V^{(k-l)}\right)}^T\right)\\ \times S^{(k-r^{\prime }+j)}{\left(S^{(k-r^{\prime }+j)}\right)}^T\left(\underset{l=r^{\prime }-j-1}{\overset{0}{\Pi }}{V}^{(k-l)}\right)\end{array}\right)$

其中s^(k)＝θ^(k+1)-θ^(k)，ρ^(k)＝1/(s^(k))^Ty^(k),y^(k)＝g^(k+1)-g^(k),V^(k)＝(I-ρ^(k)y^(k)(s^(k))^T)。

第5步，令k＝k+1，转第2步。

实施例2：

下面结合图2(a)和图2(b)对本发明的效果做进一步描述。

本实施例的仿真实验是在Intel Core(TM)2Duo CPU、主频2.33GHz，内 存2G，Windows7平台上的MATLAB7.14上实现的。

本实施例的仿真是在两个代表性的高光谱数据Indian Pines和Pavia  University上进行的，Indian Pines图像包含16类地物：苜蓿，玉米-未耕犁， 玉米-灌溉，玉米，牧草，树木，切割的牧草，干草料堆，荞麦，大豆-未耕 犁，大豆-灌溉，黄豆，小麦，树林，建筑物-草-树，石头-钢筋；Pavia University 图像包含9类地物：苜蓿，草甸，碎石，树木，着色的金属薄片，赤裸的土 壤，沥青，砖块，阴影。

图2(a)表示Indian Pines高光谱数据图像，图2(b)表示Pavia University 高光谱数据图像。一般情况下，两幅图用于分类，其中的边缘部分更易被 错分，一部分的匀质区域也会被错分；图2(a)和图2(b)为本实验的数据集。

本发明把未经过特征提取直接在原图上进行支持矢量机(SVM)分类 作为基准，与稀疏自编码单隐层网络(SAE)，PCA提取方法进行对比，都 是应用SVM分类得到分类精度，其中SAE的稀疏度参数为0.5。

在Indian Pines数据集上隐层单元数和PCA降维的维数都是100，本发 明方法窗口m大小为7，训练集为每类取10％，得到的分类精度如表1所示。

表1不同方法在Indian Pines的分类精度

方法 分类精度 SVM 84.13％ SAE+SVM 79.29％ PCA+SVM 82.78％ 本发明方法 95.61％

从表1可以看出，在Indian Pines上的仿真结果，本发明的分类精度为 95.61％，分类精度更高。

本发明在Pavia University数据集上窗口大小m为5，SAE隐层单元数 和PCA降维的维数都为50，训练集为每类取9％，得到的平均分类精度如 表2所示。

表2不同方法在Pavia University的分类精度

方法 分类精度 SVM 89.48％ SAE+SVM 90.29％ PCA+SVM 90.36％ 本发明方法 98.66％

从表2可以看出，在Pavia University上的仿真结果，本发明的分类精度 为98.66％，分类精度更高。

综上所述，通过在Indian Pines和Pavia University的仿真实验，可以看 出本发明方法能够获得较高的分类精度，因此，本发明是一种有效的高光 谱特征提取方法。

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