基于非监督特征学习的高分辨率遥感影像场景分类方法

摘要

一种基于非监督特征学习的高分辨率遥感影像场景分类方法，包括对输入的原始高分辨率遥感图像进行划分得到场景，从每个场景中随机地提取若干个训练图像块，将训练图像块聚集起来做预处理操作；计算所有训练图像块的低维流形表示，聚类得到一组聚类中心；对每一幅场景密集采样得到局部图像块，对每个局部图像块做预处理操作后映射到相同的低维流形空间中，然后进行编码获得场景的所有局部特征；将所有场景的局部特征集合起来进行特征量化，统计每一幅场景的局部特征直方图，得到场景的全局特征表达；随机挑选若干幅场景作为训练样本，由分类器得到每一幅场景的预测类别标号，完成原始高分辨率遥感场景的标注任务。

说明书

技术领域

本发明属于遥感图像智能化分析领域，特别涉及用于高分辨率遥感图像的场景分类，是 一种基于非监督特征学习的遥感场景分类方法。

背景技术

随着摄影测量与遥感成像技术的快速发展，高分辨遥感影像的数量也迅猛增长，对遥感 图像解译的需求也变得越来越迫切。其中遥感图像的场景分类是智能化遥感信息处理领域的 一项重要的任务，因此对高分辨率遥感图像场景的快速理解和自动标注问题也受到了更多的 关注。

遥感场景通常是指遥感影像中具有语义信息的局部区域，例如，对于一幅城区的遥感影 像，它可能包含居民区、商业区、工业区、绿地、水域等多种不同的场景。由于遥感影像在 成像过程中受到环境因素的影响，使得同一类别的场景在光照、方向、尺度上呈现出较大的 区别。另一方面不同的场景类别可能包含相同的目标类别，例如商业区和居民区场景中都包 含道路、建筑、树木等一些具体的语义类别，只是由于这些语义目标在空间位置和疏密程度 上的不同导致了商业区和居民区的场景语义区别。因此，由于遥感场景复杂性和多变性，使 得智能化的场景分类成为了一项挑战性的工作。

基于词袋模型(Bag of Words)[1][2]的特征表示方法一直是有效解决场景分类问题的主 流。词袋模型将场景描述成若干无序视觉词汇(Visual Words)的集合，通过统计场景中视觉词 汇直方图来达到提取场景全局特征的目的。词袋模型中的视觉词汇是通过对场景中大量的局 部特征聚类量化得到，局部特征的选择的优劣会直接影响最后的场景分类精度。一般情况这 些局部特征都是人工设计得到的，然而设计一个稳定的多方面不变性（如尺度不变性、旋转 不变性、平移不变性等）的局部特征不仅需要十分专业的领域知识，更需要大量的反复试验， 因此人工设计特征显得耗时耗力。近些年来，很多研究致力于利用非监督特征学习 (Unsupervised Feature Learning)[3][4][5][6]的方法从海量无标号的数据中自动学习得到局部 低层和中层特征，代替这些人工设计的特征，并且在机器学习，计算机视觉和数据挖掘等领 域取得了重大的突破。但是这些算法都是直接利用数据的原始形式作为训练数据，并没有挖 掘这些训练数据在原始数据空间上的本征结构和信息，使其在训练和编码过程中的计算代价 太大。如何利用非监督特征学习的方法框架，训练编码得到能够有效描述遥感场景的局部特 征成为了提高场景分类精度的关键。

有关参考文献：

[1]J.Sivic and A.Zisserman,"Video Google:a text retrieval approach to object matching in  videos,"in Computer Vision,2003.Proceedings.Ninth IEEE International Conference o,2003,pp. 1470-1477vol.2.

[2]S.Lazebnik,C.Schmid,and J.Ponce,"Beyond Bags of Features:Spatial Pyramid Matching  for Recognizing Natural Scene Categories,"in Computer Vision and Pattern Recognition,2006 IEEE Computer Society Conference on,2006,pp.2169-2178.

[3]G.E.Hinton and R.R.Salakhutdinov,"Reducing the dimensionality of data with neural  networks,"Science,vol.313,pp.504-507,2006.

[4]A.Coates,A.Y.Ng,and H.Lee,"An analysis of single-layer networks in unsupervised feature  learning,"in International Conference on Artificial Intelligence and Statistics,2011,pp.215-223.

[5]Y.L.Boureau,F.Bach,Y.LeCun,and J.Ponce,"Learning mid-level features for recognition," in Computer Vision and Pattern Recognition(CVPR),2010IEEE Conference on,2010,pp. 2559-2566.

[6]A.M.Cheriyadat,"Unsupervised Feature Learning for Aerial Scene Classification," Geoscience and Remote Sensing,IEEE Transactions on,vol.52,pp.439-451,2014.

发明内容

针对现有技术不足，本发明提出一种基于非监督特征学习的遥感场景的分类方法，通过 研究遥感场景中局部区域（图像块，patch）的低维流形结构，使非监督特征学习方法在离线 训练和在线编码的过程都在图像块的低维流形中执行，从而达到提高非监督学习计算效率的 目的；然后利用编码得到的局部特征，结合词袋模型获取遥感场景的全局特征表达，从而实 现遥感场景的分类。

本发明技术方案提供一种基于非监督特征学习的高分辨率遥感影像场景分类方法，包括 以下步骤，

步骤1，对输入的原始高分辨率遥感图像进行均匀网格划分，每一个子网格代表一个场景；

步骤2，从每个场景中随机地提取若干个训练图像块，将训练图像块聚集起来做预处理操作；

步骤3，计算所有训练图像块的低维流形表示，并以训练图像块的低维矢量表示为基础，通 过聚类的方法得到一组聚类中心；

步骤4，对每一幅场景密集采样，得到若干个局部图像块，对每个局部图像块做预处理操作 后，映射到相同的低维流形空间中，并获得相应的低维流形表示，对局部图像块的低维流形 表示进行编码，获得描述该幅场景的所有局部特征；

步骤5，将所有场景的局部特征集合起来进行特征量化，统计每一幅场景的局部特征直方图， 得到场景的全局特征表达；

步骤6，随机挑选若干幅场景作为训练样本，并把每幅场景对应的全局特征表达和真实的类 别标号送入分类器中进行训练，得到分类器模型参数；将原始高分辨率遥感影像的所有场景 的全局特征表达送入分类器，由分类器得到每一幅场景的预测类别标号，完成原始高分辨率 遥感场景的标注任务。

而且，步骤2中，采用均值方差归一化和ZCA白化的方法对训练图像块做预处理；步骤 4中，采用均值方差归一化和ZCA白化的方法对局部图像块做预处理。

而且，步骤3中通过聚类的方法得到一组聚类中心时，所述聚类的方法采用K-means算 法，提取聚类中心的过程包括以下步骤，

1）采用K近邻算法以所有训练图像块为顶点构造一个无向邻接图G；

2）为无向邻接图G中每一条边赋予一个权重，来衡量由边连接的两个顶点v_i、v_j之间的 相互度W_ij，所得权重矩阵W中相似度的定义如下，

$W_{\mathrm{ij}}=\left(\begin{array}{cc}{e}^{(-{\left|\right|v_i-v_j\left|\right|}_2^2/t)}& \mathrm{if}{v}_i\in N\left(v_j\right)\\ 0& \mathrm{otherwise}\end{array}\right)$

其中，N(v_j)表示与顶点v_j最相邻的K个顶点的集合，t为权重衰减常数，e为数学常数；

3）计算Laplacian矩阵L,L＝A-W,其中A是一个对角矩阵，满足其中元素 A_jj＝Σ_iW_ij；

4）解决一个一般化的特征值问题如下，

$X_w{\mathrm{LX}}_w^{T}m=\lambda X_{w}A{X}_w^{T}m$

其中，λ和m表示待求的特征值和特征向量，将计算所得最小的d个特征值对应的特征向量 m₀,m₁,...,m_d-1按照特征值由小到大的顺序排列起来，得到一个线性映射矩阵M，然后利用M 将所有的预处理后的训练图像块相应向量映射到一个低维流形上，并得到对应的低维流形表 示如下，

p⁽ⁱ⁾＝M^Ty⁽ⁱ⁾,M＝(m₀,m₁,...,m_d-1)

其中，y⁽ⁱ⁾表示第i个训练图像块预处理后的向量，为第i个训练图像块的低维流形表 示，d<<n；

5）在4）所得训练图像块的低维流形表示的基础上，采用K-means聚类算法对所有训练 patch在低维流形上进行聚类，通过解如下优化问题得到一组聚类中心记为字典D，完成非监 督参数训练，

$\underset{D,c}{\min }\underset{i}{\Sigma }{\left|\right|{\mathrm{Dc}}^{\left(i\right)}-p^{\left(i\right)}\left|\right|}_2^2$

$s.t.{\left|\right|D^{\left(k\right)}\left|\right|}_2=1,\forall k$

$\mathrm{and}{\left|\right|c^{\left(i\right)}\left|\right|}_0\le 1,\forall i$

其中，c⁽ⁱ⁾表示第i个训练图像块的聚类中心分配向量，c表示任意训练图像块的聚类中心分 配向量，D^(k)表示字典D中第k个聚类中心。

而且，步骤4中，采用三角编码方法对场景中的所有局部图像块进行特征编码。

而且，步骤5中，通过采用K-means算法对局部图像块的局部特征进行聚类，对每一场 景类别都训练一个码本，将所有场景类别的码本串联起来构成一个联合码本，利用获得的联 合码本对每个局部图像块做量化操作。

而且，步骤6中，所述分类器采用直方图交叉核的SVM分类器。

本发明通过挖掘遥感场景中的局部图像块在在低维本征空间中的优良性质，利用非监督 特征学习的方法从大量维标号图像块中自动地学习得到场景局部特征；并进一步对模型参数 训练和图像块特征编码的过程进行优化。本发明在保证场景分类精度的前提下，大大降低了 非监督特征学习的计算代价，节约了系统资源。

附图说明

图1为本发明实施例的非监督特征学习整体框架流程图；

图2为本发明实施例的场景全局直方图表达流程示意图；

具体实施方式

本发明技术方案可采用计算机软件方式支持自动运行流程。以下结合附图和实施例详细 说明本发明技术方案。

参见图1，本发明实施例包括以下步骤：

步骤1：遥感影像的场景划分

要对遥感影像进行场景分类，首先要定义场景单元与场景类别个数，实施例在一副大尺 寸遥感影像中选择合适尺寸的矩形区域作为一个场景，最终的目标就是给这幅遥感影像中所 有的场景赋予一个场景类别标号，每一类场景用不同的颜色区分。具体实施过程中，采用均 匀网格对大幅遥感影像进行划分，每一个子网格就代表一个场景，相邻的场景之间没有重叠。 例如，大幅遥感影像的尺寸为6000×6000像素，每一个场景的尺寸为150×150像素，将影像 划分成了1600个子场景。划分得到测试遥感场景集。

步骤2：提取训练图像块并预处理

从每个场景中随机地提取若干个图像块作为训练图像块，将训练图像块聚集起来做预处 理操作。可对每一个训练图像块首先进行均值方差归一化操作，然后利用ZCA(Zero  Component Analysis)白化的方法除去训练图像块中每一维之间的相关性。

实施例的步骤2实现方式如下：

在步骤1获得大幅影像的所有场景后，需要计算每幅场景的局部特征描述，而局部特征 是从场景中的局部区域（小图像块，或者称为patch）提取得到的。因此在接下来的步骤中， 研究对象主要集中在以patch为基本单元的参数学习和编码上。参见图1，非监督特征学习算 法可分为两个阶段：模型参数学习和特征编码，分别负责非监督参数训练、场景局部特征编 码与全局特征表示。在参数学习阶段，首先要选取大量的无标号的patch。实施例中，根据步 骤1初始化分割所得未标号的遥感场景，对每一幅场景随机地抽取大量尺度较小的patch，将 所有patch在像素级空间上展开成一个向量，分别作均值方差的归一化，然后将所有归一化 后的patch向量组x⁽¹⁾,x⁽²⁾,...,x^(M)合成一个矩阵X＝[x⁽¹⁾,x⁽²⁾,...,x^(M)],表示n维实 空间，n＝r×r×l,r为patch的尺寸，l为patch的通道数（对于一般的彩色图像，通道数为 RGB三通道）；M表示训练patch的个数，i的取值为1,2,…,M。对矩阵X做ZCA白化操作， 得到白化后的矩阵X_w＝BX，其中矩阵V和D_Eig分别表示矩阵X的协方差矩 阵的特征向量和特征值构造的对角矩阵。

步骤3：低维流形空间聚类算法

经流形学习，计算所有训练图像块的低维流形表示，并以训练图像块的低维矢量表示为 基础，通过聚类的方法得到一组聚类中心（字典）。可采用LPP(Locality Preserving Projection, 局部保留投影)算法计算所有训练图像块的低维流形表达，并且得到相应的线性映射矩阵；采 用K-means算法对训练图像块的低维流形表达进行聚类，得到一组聚类中心。

在对所有训练patch进行预处理后，得到了训练patch新的表达X_w＝[y⁽¹⁾,y⁽²⁾,...,y^(M)],但是 对patch的预处理仅仅能够去除每一维之间的相关性、平衡每一维的尺度，每个patch向量的维 度并没有发生变化，也没有挖掘patch向量空间的本质结构。因此，本发明实施例提出了在低 维本征流形空间中做聚类操作的方法，得到模型参数D。具体而言，需要如下几步：

1）采用K近邻算法以所有训练patch为顶点构造一个无向邻接图G，K近邻算法为现有技 术；

2）为图G中每一条边赋予一个权重，来衡量由边连接的两个顶点v_i、v_j之间的相互关系 （相似度W_ij）。这个权重矩阵W（相似度矩阵）中相似度的定义为：

$W_{\mathrm{ij}}=\left(\begin{array}{cc}{e}^{(-{\left|\right|v_i-v_j\left|\right|}_2^2/t)}& \mathrm{if}{v}_i\in N\left(v_j\right)\\ 0& \mathrm{otherwise}\end{array}\right)$

其中，N(v_j)表示与顶点vj最相邻的K个顶点的集合，t为一个权重衰减常数，e为数学常数；

3）计算Laplacian矩阵L,L＝A-W,其中A是一个对角矩阵，满足其中元素A_jj＝Σ_iW_ij；

4）解决一个一般化的特征值问题，可以描述为：

$X_w{\mathrm{LX}}_w^{T}m=\lambda X_{w}A{X}_w^{T}m$

其中，λ和m表示待求的特征值和特征向量。将计算得到的特征向量m₀,m₁,...,m_d-1（只选择 最小的d个特征值对应的特征向量，d为预设数值）按照所对应的特征值由小到大的顺序排列 起来，得到一个线性映射矩阵M。然后利用M将所有的预处理后的patch向量映射到一个低 维流形上实现低维映射，并得到对应的低维流形表示：

p⁽ⁱ⁾＝M^Ty⁽ⁱ⁾,M＝(m₀,m₁,...,m_d-1)

其中，y⁽ⁱ⁾表示第i个训练图像块预处理后的向量，为第i个训练patch的低维 流形表示。

5）在获得训练patch的低维流形的基础上，采用K-means聚类算法对所有训练patch在低维 流形上进行聚类，通过解如下优化问题得到一组聚类中心D，完成非监督参数训练：

$\underset{D,c}{\min }\underset{i}{\Sigma }{\left|\right|{\mathrm{Dc}}^{\left(i\right)}-p^{\left(i\right)}\left|\right|}_2^2$

$s.t.{\left|\right|D^{\left(k\right)}\left|\right|}_2=1,\forall k$

$\mathrm{and}{\left|\right|c^{\left(i\right)}\left|\right|}_0\le 1,\forall i$

其中，c⁽ⁱ⁾表示第i个训练patch的聚类中心分配向量，c表示任意训练patch的聚类中心分配向 量；这里得到的聚类中心D既是由K-means算法得到的字典，又是非监督学习方法中的模型参 数，在接下来的特征编码过程中有着至关重要的作用。D^(k)表示字典D中第k个聚类中心。

基于从各场景中提取的所有训练图像块，按照以上方式进行训练后，即可得到一个线性 映射矩阵M。

步骤4、场景中patch特征编码

对步骤1所得每一幅场景密集采样，得到若干个局部图像块，构成场景中局部patch集合， 对每个局部图像块做预处理操作后，将其映射到相同的低维流形空间中，并获得相应的低维 流形表示，对局部图像块的低维流形表示进行编码，获得描述该幅场景的所有局部特征。

实施例的步骤4实现方式如下：

对每幅场景密集采样可采用等间隔方式，即每隔若干个像素点提取一个固定大小的局部 图像块，可提取与训练patch相同尺度的局部patch。每隔若干个像素点采样一次（间隔点数 可由本领域技术人员自行指定，实施例中选择的采样间隔为1个像素点），这样相邻的局部 patch之间会有部分重叠像素。这样的采样方式可以使得对场景的局部描述更加精细准确。

对每幅场景中的所有提取出的局部图像块做与步骤2中相同的预处理操作，即在采样得 到场景中所有局部patch后，仍然需要对这局部patch做均值归一化和ZCA白化预处理。

利用步骤3中由LPP(Locality Preserving Projections)得到的线性映射矩阵M，将场景中的 每一个局部patch都映射到同一个patch低维流形中，并同时获得相应的流形表示；采用三角 编码（Triangle encoding）的方法对局部patch的低维流形表示进行编码，得到每个局部patch 对应的特征表示。编码方式如下，这里需要用到步骤3中训练出的模型参数D：

φ_k＝max{0,mean(d)-d_k}

其中，中间变量d_k＝||p⁽ⁱ⁾-D^(k)||₂，中间变量N为字典D的大小（聚 类中心的个数），p⁽ⁱ⁾表示第i个待编码的局部patch的低维表示，D^(k)表示字典D中第k个聚 类中心，φ_k表示局部patch的第k维的编码特征，k的取值为1,2,…,N。由这个编码规则可以 发现，每个局部patch编码特征的维度是由K-means聚类得到的聚类中心的个数决定的。通 过特征编码，从而得到了每幅场景中所有局部patch的特征表达（特征矢量），也就是场景的 局部特征描述。

步骤5：场景全局特征表达

将所有场景的局部特征集合起来进行特征量化，统计每一幅场景的局部特征直方图，以 此作为场景的全局特征表达。可再次采用K-means算法对局部图像块的编码特征进行聚类， 利用形成的联合码本对每个局部图像块做量化操作；统计每幅场景中的所有局部图像块的量 化直方图，将直方图表示作为该幅场景的全局特征表示。

实施例的步骤5实现方式如下：

为了继续场景分类的工作，仅仅获得场景的局部特征描述是不够的，因为局部特征只能 反映场景局部的结构、纹理、方向等底层属性，要进行中高层的场景理解需要在底层特征的 基础上做全局性的分析，提取出能够全面有效地描述整幅场景的全局特征表达。基于词袋模 型的全局特征表达已经在场景分类中取得了很好的效果，在本发明中优化了传统的词袋模型， 并且进一步提高了场景分类的准确性。传统的词袋模型将所有场景局部特征集合起来聚类得 到一组聚类中心（码本），然后将每个局部特征量化成码本中的一个码字（通常计算特征到每 个码字的欧式距离，然后选择与之距离最小的那个码字），最后统计每幅场景中所有码字出现 的频率，将频率直方图作为场景全局特征。这样做忽略掉了场景类别之间的局部特征的区别， 使得训练出来的码字并不专属于某一场景类别。鉴于此，本发明中对每一场景类别都训练一 个码本，假设场景类别数为K，第i类场景训练得到的码本为i的取值为1,2,…,K， 码字个数为L,将所有场景类别的码本串联起来构成一个联合码本 对于第j幅场景，计算该场景中所有局部特征在联合码本C^J中 对应码字的频率直方图最后将这个KL维的直方图向量作为场景j的全局特征表 达。本发明中训练码本采用经典K-means聚类算法即可。本发明中的场景全局特征表达方法 流程如附图2所示，第1,2,…,K类场景分别为飞机场、高尔夫场…立交桥，相应码本 C⁽¹⁾,C⁽²⁾,...,C^(K)。

步骤6、利用SVM分类器实现场景分类

随机挑选若干幅场景作为训练样本，构成训练样本集，并把每幅场景对应的直方图特征 （即训练集特征）和真实的类别标号送入分类器中进行训练，得到分类器模型参数；将原大 幅遥感影像的所有场景作为测试样本集，将各场景的直方图特征（即测试集特征）送入分类 器，由分类器得到每一幅场景的预测类别标号，从而完成大幅遥感场景的标注任务。

本发明进一步提出，可采用SVM(Support Vector Machines)分类器对选取的训练样本进行 训练，从而得到分类器模型参数，其中训练SVM分类器时，采用直方图交叉核(Histogram  Intersection Kernel,HIK)。

实施例的步骤6实现方式如下：

通过以上步骤，已经获得了大幅遥感影像中所有场景的全局特征描述，在此基础上，本 发明采用SVM分类器来完成场景分类工作。场景分类过程分为两步：训练SVM分类器和测 试场景标号预测。首先对于预先定义好的每一类场景类别，选择若干幅场景，并通过步骤1-5 计算其全局特征，利用获得的全局特征和对应的标号训练SVM分类器，得到分类器的模型 参数；然后利用已训练的SVM分类器对大幅遥感影像中的所有场景进行标号预测，最终得 到每幅场景的场景类别，并将不同的场景类别用不同的颜色加以区分。

本发明中使用的SVM分类器是以最大分类间隔为准则，通过寻找最优分类超平面来对 样本进行空间划分，可参见现有SVM分类器技术。SVM分类器的训练过程可以最终转化成 解决如下的优化问题：

$\min \frac{1}{2}{\left|\right|w\left|\right|}^2+C\underset{i=1}{\overset{n}{\Sigma }}{\xi }_i$

s.t.,y_i(w^Tφ(x_i)+b)≥1-ξ_i,i＝1,...,n

ξ_i≥0,i＝1,...,n

其中w和b为定义分类超平面w^Tφ(x_i)+b＝0的参数向量，C是一个常量用于控制目标函数中 两项之间的权重，ξ_i表示松弛变量，φ(x_i)表示训练样本x_i的非线性映射，y_i表示第i个训 练样本x_i的类别标号,n为训练样本的个数。采用拉格朗日乘子法，将该带约束的优化问题转 化成无约束优化问题，其代价函数为：

$L(w,b,\xi ,\alpha ,r)=\frac{1}{2}{\left|\right|w\left|\right|}^2+C\underset{i=1}{\overset{n}{\Sigma }}{\xi }_i-\underset{i=1}{\overset{n}{\Sigma }}{\alpha }_i(y_i(w^T\phi \left(x_i\right)+b)-1+{\xi }_i)-\underset{i=1}{\overset{n}{\Sigma }}{r}_i{\xi }_i$

α_i,r_i均为需要优化的变量。在满足Karush-Kuhn-Tucker（KKT）条件的基础上，通过推到将 上面的无约束问题转化成如下优化问题：

$\max \underset{i=1}{\overset{n}{\Sigma }}{\alpha }_i-\frac{1}{2}\underset{i,j=1}{\overset{n}{\Sigma }}{\alpha }_i,{\alpha }_j,y_i,y_{j}K(x_i,x_j)$

s.t.,0≤α_i≤C,i＝1,...,n

$\underset{i=1}{\overset{n}{\Sigma }}{\alpha }_i{y}_i=0$

其中参数向量w和变量α之间满足K(x_i,x_j)为自定义的一种核函数，其目的 是将样本在特征空间中的内积操作转化成样本在原始空间中的核函数映射，y_j表示第j个训 练样本x_j的类别标号。本发明实施例采用直方图交叉核（HIK）作为核函数，定义为：

$K(x_i,x_j)=\underset{s}{\Sigma }\min (x_i\left(s\right),x_j\left(s\right))$

其中s表示样本向量的第s维。通过解上述优化问题得到变量{α_i}_i＝1,...,n的最优值，完成SVM 分类器模型参数的训练过程。

利用已训练好的SVM进行场景分类时，只需要将测试场景的特征向量带入SVM分类器 的分类函数中即可得到最终的预测类别标号，分类函数如下所示：

$f\left(H^{\left(t\right)}\right)=\underset{i=1}{\overset{n}{\Sigma }}{\alpha }_i{y}_{i}K(H_i,H^{\left(t\right)})+b$

其中H_i为第i个训练场景的全局特征，H^(t)为待分类场景的全局特征。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技 术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不 会偏离本发明的精神或者超越所附权利要求书所定义的范围。