CN102013017B - 一种高分辨率遥感图像场景粗分类方法 - Google Patents

Abstract

本发明公开了一种高分辨率遥感图像的场景粗分类方法，用于为遥感图像的目标识别提供上下文信息，其特征在于，该方法包括如下步骤：（1）图像初始分割；（2）上下文分类；（3）区域融合和（4）后处理。本发明将快速初分割、融入语义信息的上下文分类以及面向对象区域分类结合起来对遥感图像场景进行粗分类，可以用于星上多类目标识别的图像预处理。

Description

一种高分辨率遥感图像场景粗分类方法

技术领域

本发明属于遥感图像处理领域，具体涉及一种对高分辨率遥感图像进行场景分类的方法，用于为遥感图像的目标识别提供上下文信息。

背景技术

高分辨率遥感图像一般指1m～5m的可见光航空遥感图像，高分辨率遥感图像的数据量很大，如果直接对其进行特征表达、特征聚类及分割的方法进行目标识别，将耗费大量的计算资源，导致处理效率低下。因此，提出大尺度下对遥感图像进行场景粗分类的思想，将一幅图像分割为跟我们所感兴趣的目标相关联的有限区域，然后根据目标的性质在具有依赖关系的区域中提取目标。这样可以提高目标提取的效率和精度。这就是场景分类的目的所在。

遥感图像分类包含特征计算方法和分类方法两大主干。特征计算有多种方法。文献表明，当前应用较为广泛的有统计方法、随机场模型方法和信号处理方法。统计方法特征计算的计算量大、分割精度差，受到噪声影响较大；随机场模型对大尺寸、灰度级较多的图像分割计算量是极大的，不满足快速处理的需求。而信号处理中的小波变换方法具有多尺度特性，可以对图像进行不同分辨率上的处理，在遥感图像处理中应用成熟。

现有的分类技术大都应用于针对图像的单个像素进行分类，强调像素级分类的精度。这在小尺度目标分类中具有实际意义。但遥感图像在小尺度上不具备语义上的意义，因为地物场景通常都是大尺度区域，因此像素级分类无论是在效率上还是精度上都不符合实际需求。因此一种可以对大尺度对象进行分类的方法对遥感图像分类比较有效。

发明内容

本发明的目的在于提出一种高分辨率遥感图像的场景粗分类方法，为基于遥感图像的目标识别提供上下文信息。

为实现本发明的目的所采用的具体技术方案为：

一种高分辨率遥感图像的场景粗分类方法，包括如下步骤：

1.图像初始分割

图像初始分割包含特征提取、特征数据聚类以及图像标号。

特征提取采用基于小波分解的Laws纹理算子法，首先对图像进行二维离散小波分解，得到多组小波系数矩阵，再利用Laws算子对得到的小波系数矩阵进行平滑处理，得到表征图像每个像素点的特征数据。

接着，采用模糊C均值聚类方法对特征数据进行类别划分，计算出隶属度矩阵和聚类中心点并保存数据。

最后，采用最大化隶属度矩阵的准则对图像的每个像素点进行类别标记，由此得到初始分割结果。

采用小波分解进行特征描述是基于遥感图像的数据量较大而产生的计算效率低下的考虑，因为小波分解直接降低了图像的分辨率。同时由于遥感图像表现出来的纹理特性，采用Laws算子则正好体现了纹理信息在遥感图像特征提取中的应用。采用模糊C均值聚类方法有两个目的：一是对特征进行聚类；二是计算出隶属度矩阵可以作为后面的上下文分类处理的输入。

2.上下文分类

首先，将步骤1中利用模糊C均值聚类计算所得到的隶属度矩阵作为初始后验概率；再利用概率松弛迭代算法对每个像素8邻域内的中心点做概率更新，直至收敛，得到最终的后验概率；最后，根据最终的后验概率，采用最大化后验概率准则对图像的每个像素点进行类别标记，得到调整后的分割结果，即完成对图像的上下文分类。

其中，概率更新的公式定义如下：

$p_i^{(n+1)}(l_i=\mathrm{\α})=\frac{p_i^{\left(n\right)}(l_i=\mathrm{\α})\underset{(v_i,v_j)\∈E}{\mathrm{\Π}}{Q}_{\mathrm{ij}}^{\left(n\right)}\left(\mathrm{\α}\right)}{\underset{\mathrm{\λ}}{\mathrm{\Σ}}{p}_i^{\left(n\right)}(l_i=\mathrm{\λ})\underset{(v_i,v_j)\∈E}{\mathrm{\Π}}{Q}_{\mathrm{ij}}^{\left(n\right)}\left(\mathrm{\λ}\right)}---\left(1\right)$

其中，

表示像素点i属于类别α的概率；n表示迭代更新的次数，(v_i，v_j)表示像素点i及其8邻域中的任一像素点j组成的相邻像素对，且所有符合条件的像素对组成一个集合E；λ表示类别属性，Q为邻域函数，其作为邻域像素对中心点的类别属性的支撑度的度量，某个像素点的邻域函数Q是该像素点的8邻域内所有像素点的先验概率的加权和，邻域函数

表示第n步迭代时，邻域像素点j对中心像素点i的类别为α的支撑度，其数学表达为：

$Q_{\mathrm{ij}}^{\left(n\right)}\left(\mathrm{\α}\right)=\underset{\mathrm{\β}\∈\mathrm{\Λ}}{\mathrm{\Σ}}{p}_j^{\left(n\right)}\left(\mathrm{\β}\right)r_{\mathrm{ij}}(\mathrm{\α},\mathrm{\β})---\left(2\right)$

上式(2)中，Λ是类别属性的集合且α∈Λ，β∈Λ。

表示第n步迭代时类别β的先验概率分布。r_ij(α，β)为权重，量化地表示像素i属于类别α，且像素j属于类别β的可能性，称之为兼容性系数，其数学表达如下：

$r_{\mathrm{ij}}(\mathrm{\α},\mathrm{\β})=\frac{p(l_i=\mathrm{\α},l_j=\mathrm{\β})}{p(l_i=\mathrm{\α})p(l_j=\mathrm{\β})}---\left(3\right)$

在公式(3)中，p(l_i＝α，l_j＝β)表示像素点i的类别为α且像素点j的类别为β的联合概率。

将公式(1)作为更新准则，迭代至收敛，得到的结果为每个像素点的新的类别属性。然后根据后验概率最大化准则，得到分割图像。

3.区域融合

实验表明上面得到的分割图像包含了较多的零碎的区域块，这些区域块在尺度上尚不足以提供语义上下文信息，因此我们可以将它们合并至更大尺度的区域中去，从而在全图范围内形成少数只包含我们所关注的地物类别的并且相互之间轮廓明晰的区域。

以上下文分类所得到的分割图像作为输入，经过区域合并处理，得到一幅彩色图像，其中所包含RGB值的个数代表了合并之后的区域数总和。从实验可以看出，区域数目与处理前相比减少了许多，并且更能反映真实图像中地物分布。

接下来，对上述彩色图像做多灰度值标记，然后对标记图建立邻接矩阵，并进行如下处理：1)首先计算各区域面积，若其面积小于设定的阈值，则将其列入待合并区域序列，否则不作处理；2)对于合并序列中的孤立区域，将其合并至其邻接区域，迭代执行至不再存在孤立区域；对于在合并序列中非孤立区域的区域，分别计算该区域与其所有相邻的区域之间的公共边界的长度，将其合并至拥有最长公共边界的那个邻接区域中去。

经过上述处理即完成区域融合，得到融合后的分割图。

4.后处理

这一部分主要是采用面向对象的分类思想。首先对步骤3得到的分割图进行多值标记，得到若干个连通区域。然后把每个连通区域映射到原始图像，得到对应位置上的局部原始图像。

然后利用边缘密度、亮度均值、方差、能量等特征来表征每个连通区域的属性，分别计算每个区域对应的局部原始图像的上述特征值。

接着进行样本采集和样本特征计算处理。样本选取典型的，能够代表不同地物特点的图像，计算前面提到的特征。然后对采集的样本进行训练，得到区分不同地物的统计知识。

最后，根据训练得到的统计知识，利用不同地物类型的边缘密度的差异性，把区域划分为水域、陆地(山地、草地、人工建筑等)。再利用灰度差异(亮度均值、方差、能量等)将陆地上的不同地物(人工建筑、自然场景)区分开，最终达到场景分类的目的。

其中，所述利用特征表征每个区域的属性的具体过程为：首先获取每个区域相对应的原始图像数据，然后对获得的原始图像数据进行平滑滤波，边缘提取，以及边缘密度、亮度均值、方差、能量等特征计算的操作。

面向对象分类的具体流程如图5所示。

本方案的创新点为将快速初分割、融入语义信息的上下文分类以及面向对象区域分类结合起来对遥感图像场景进行粗分类。

(1)快速初分割

通常的遥感图像分类方法是对图像的每个像素点进行特征计算，然后通过样本学习及分类器训练，为每个像素分配类别标号，从而得到分类结果。考虑到遥感图像的大尺寸特性，这一过程需要消耗很大的计算量。如果能先把图像划分为若干区域，然后再对区域进行分类，那么计算量将大大降低。本发明正是基于这样的思路，首先对图像进行初始分割，得到若干同质区域。

如前所述，图像的特征提取部分占主要的计算量。因此在提取特征部分要实现快速计算。为此，本方案采用小波分解的方法，使用小波函数族及其相应的尺度函数将原始信号分解成不同的频带，其标准的分解过程是采用金字塔算法。这样就可以把高分辨率的图像降低分辨率为低分辨率图像，随着分辨率的降低，图像的纹理特性得到了增强，并且图像尺寸降低，带来了计算效率的提升。

待特征计算完毕之后，利用模糊C均值聚类算法，将图像分割为若干区域。在所有的处理完成之后，将图像恢复至原始的分辨率即可。

(2)融入语义信息的上下文分类

传统的像素分类仅根据像素的特征对其分类，而并没有考虑该像素周围的像素所属类别对它的影响。然而在任何实际的图像中，符合人的视觉认知的事实是相邻像素之间总存在一定的联系，这种相邻像素之间的联系可以为分类提供信息。在本发明中，我们利用邻域上下文约束对初始分类进行改进，得到语义上更为正确的结果。

(3)面向对象分类

经过区域合并之后得到若干形状不规则的区域。我们认为区域内部的像素点是同质的，相邻区域之间是互斥的。但是从地物类型的角度，不相邻的区域可能同属一种地物类型。对于判别区域的地物属性，传统的分类是基于像素的分类，不仅计算量大，而且分割结果语义性较差。面向对象的分类思想，是把每个连通区域作为一个对象，直接判定对象的属性。

附图说明

图1：本发明总体结构示意图

图2：图像初始分割流程图

图3：上下文分类流程图

图4：区域合并流程图

图5：面向对象分类流程图

具体实施方式

本实施例以一幅分辨率在3000*3000左右的可见光遥感图像为实例，对本发明的场景分类进行具体描述。

(1)图像初始分割

首先对原始图像进行小波分解，采用Daubechies4点小波，做一级分解，得到的是降低一倍分辨率的4组小波系数矩阵，分别是原始图像的近似系数以及3组细节系数。

接下来对这4组小波系数矩阵进行Laws纹理算子运算，得到描述原图的特征向量组。这一过程大致如下：已知了图像中包含的纹理类别数，然后从小波分解最大尺度的所有频带LL，HL，LH，HH开始，每个位置构造一个4维向量，再利用Laws算子定义大小为(2n₁+1)×(2n₁+1)的窗口上的标准差作为特征：

$s(i,j)=\frac{1}{{(2n_1+1)}^2}\underset{k=i-n_1}{\overset{i+n_1}{\mathrm{\Σ}}}\underset{l=j-n_1}{\overset{j+n_1}{\mathrm{\Σ}}}|g(k,l)-m(i,j)|---\left(4\right)$

式中，g(k，l)为小波分解子图，m(i，j)为窗口内g(k，l)的平均值。(k，l)及(i，j)表示像素点在图像中的坐标。

为了消除边缘效应，对特征图像作进一步平滑，计算大小为(2n₂+1)×(2n₂+1)的窗口内的均值：

$S(i,j)=\frac{1}{{(2n_2+1)}^2}\underset{k=i-n_2}{\overset{i+n_2}{\mathrm{\Σ}}}\underset{l=j-n_2}{\overset{j+n_2}{\mathrm{\Σ}}}s(k,l)---\left(5\right)$

其中，n₁和n₂分别是定义窗口的大小的值。

然后进入模糊C均值聚类处理，这个过程中的输入是前面得到的特征向量组以及参数值。参数值包括聚类数目以及迭代次数、误差阈值等，通过实验调整到较为优化的数值之后就可以保持参数值不变。输出结果为聚类中心点和隶属度矩阵，隶属度矩阵中每个元素的值表示每个像素点属于某一类别的隶属度大小。

最后根据最大化准则，将每个像素点的类别标记为具有最大隶属度值的那个类别，完成标号操作。同时，隶属度矩阵将作为上下文分类的输入，用于迭代求解引入邻域上下文信息的像素点的类别后验概率，最后的标记结果根据后验概率最大化准则得到。

(2)上下文分类

在图像中，定义每个像素点的8邻域。每个像素代表图中的一个节点，弧连接每个像素点与它的直接相邻点，并且节点标记表示可能属于的类别。令

为节点v_i的初始类别分布，即为节点v_i的类别的初始先验概率，它可以通过图像初始分类获取，即前面介绍的模糊C聚类分割。

根据公式(1)，分别求解类别先验概率，相邻像素点之间的联合概率以及兼容性系数。其中类别先验概率和联合概率在迭代过程中保持不变，它们仅与初始后验概率有关。具体计算方式为：

$p^{*}\left(\mathrm{\α}\right)=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{p}_i^{\left(0\right)}(l_i=\mathrm{\α})---\left(6\right)$

$p_k^{*}(\mathrm{\α},\mathrm{\β})=\frac{1}{\left|\right|R_k\left|\right|}\underset{(v_i,v_j)\∈R_k}{\mathrm{\Σ}}{p}_i^{\left(0\right)}(l_i=\mathrm{\α})p_j^{\left(0\right)}(l_j=\mathrm{\β})---\left(7\right)$

$r_k^{*}(\mathrm{\α},\mathrm{\β})=\frac{p_k^{*}(\mathrm{\α},\mathrm{\β})}{p^{*}\left(\mathrm{\α}\right)p^{*}\left(\mathrm{\β}\right)}---\left(8\right)$

其中，p^*(α)为类别α的先验概率，N为图像中总的像素点个数。为满足邻接关系为R_k的像素点i和它的邻接点的联合概率分布。R_k用于描述8邻域内中心像素点与其邻接点的邻接关系，||R_k||为满足关系R_k的像素点对的个数。

为满足关系R_k的兼容性系数。α与β表示类别属性。

然后用先验概率和兼容性系数的乘积累加和代替邻域函数Q，从而得到公式(1)的变换形式：

$p_i^{(n+1)}(l_i=\mathrm{\α})=\frac{1}{8}\underset{k=1}{\overset{8}{\mathrm{\Σ}}}\frac{\underset{\mathrm{\β}\∈\mathrm{\Λ}}{\mathrm{\Σ}}{p}_i^{\left(n\right)}(l_i=\mathrm{\α})p_j^{\left(n\right)}(l_j=\mathrm{\β})r_k(\mathrm{\α},\mathrm{\β})}{\underset{\mathrm{\λ}\∈\mathrm{\Λ}}{\mathrm{\Σ}}\underset{\mathrm{\β}\∈\mathrm{\Λ}}{\mathrm{\Σ}}{p}_i^{\left(n\right)}(l_i=\mathrm{\λ})p_j^{\left(n\right)}(l_j=\mathrm{\β})r_k(\mathrm{\λ},\mathrm{\β})}---\left(9\right)$

迭代的过程是用第k步迭代计算得到的后验概率作为第k+1步迭代的输入，计算新的后验概率，如此反复直至收敛。假设经过n+1步迭代完成，则得到表示像素点i的最终后验概率向量为

。再利用后验概率最大化准则：

${\mathrm{\θ}}^{*}=\underset{k}{\arg \max }\left(p_i^{\left(n\right)}(l_i=\mathrm{\θ})\right)---\left(10\right)$

得到每个像素点i的类别l_i＝θ^*。λ表示类别属性，Λ是类别属性的集合，且α∈Λ，β∈Λ。

(3)区域融合

在区域融合阶段，处理对象为一个包含多个标记值的矩阵。每个标记值代表一个连通的区域，需要将符合某些条件的区域合并至与其符合某种关系的区域中去。条件和关系均事先定义好，定义的依据是分割符合语义上的正确性以及保留我们所关注的大尺度场景。

区域合并的具体实施方式为：

首先以上下文分类得到的分割图像作为输入，进行区域合并处理，得到一幅彩色图像。其中所包含的RGB值的个数代表了合并之后的区域数总和。然后对这些区域分别做如下处理：

1)标记分割图像，得到若干连通区域。结果是一个矩阵，其中每个值代表一个连通区域；

2)建立这些区域的邻接矩阵，该矩阵仅包含0和1；

3)按区域的序号循环，分别计算连通区域R_i的面积S_i；

4)设定一个表示面积的阈值T，比较S_i和T的大小，若S_i＞T，则返回步骤3，反之则进入步骤5；

5)判断区域R_i是否为孤立区域。若是，则进入步骤6，反之进入步骤7；

6)将R_i合并至其邻接区域中，并更新邻接矩阵；找到更快速的更新邻接矩阵的方法；

7)计算R_i与其各个邻接区域的公共边界长度，找出具有最长边界的那个区域；

8)将R_i合并至该区域，并更新邻接矩阵；

9)判断是否遍历完所有的区域，若遍历完，则输出最终结果，若没有，返回步骤3；

(4)后处理，即进行基于面向对象的分类

这里的对象是指通过步骤3所做的区域融合得到的连通区域。首先把每个区域映射到原始图像中，得到每个区域对应于原始图像位置上的数据。定义边缘密度、亮度均值、方差以及能量等4种特征，计算每个区域对应的原始图像的特征值。

然后对图像数据进行采样，得到表示不同地物类型的样本集，用前面定义的相同的特征来计算这些样本的特征值。采用支撑向量机训练样本集，得到一个分类器模型。

最后，根据训练得到的分类器，对每个区域进行类别的判定。其中判定水域和陆地依据的是不同地物类型的边缘密度的差异。设定一个阈值，把边缘密度小于该阈值的区域划分为水域，反之则为陆地(山地、草地、人工建筑等)。再利用灰度差异(亮度均值、方差、能量等)，同样地设置一个阈值，将陆地上的不同地物划分为人工建筑和自然场景。从而最终得到一幅水陆分割图像，并且陆地又细分为人工建筑以及自然场景。

本发明是在MS VC 6.0编程环境下实现的。算法部分将利用C语言实现，便于移植到DSP上进行硬件实现。同时在实验阶段，借助于Matlab的Image ProcessingToolbox编写了主要的算法代码，并针对Matlab对于循环运算的效率低下的缺点，优化了部分算法以及代码，达到快速处理的目的。

Claims (4)

1.一种高分辨率遥感图像的场景粗分类方法，用于为遥感图像的目标识别提供上下文信息，其特征在于，该方法包括如下步骤：

(1)图像初始分割

首先对遥感图像进行二维离散小波分解，得到多组小波系数矩阵，再利用Laws算子对得到的小波系数矩阵进行平滑处理，得到表征图像每个像素点的特征数据；

接着，采用模糊C均值聚类方法对所述特征数据进行类别划分，计算出隶属度矩阵和聚类中心点并保存数据；

最后，采用最大化隶属度矩阵准则对遥感图像的每个像素点进行类别标记，由此得到初始分割结果；

(2)上下文分类

首先，将步骤(1)中得到的所述隶属度矩阵作为初始后验概率，再利用概率松弛迭代算法对每个像素8邻域内的中心点做概率更新，直至收敛，得到最终的后验概率；最后，根据最终的后验概率，采用最大化后验概率准则对图像的每个像素点进行类别标记，得到调整后的分割结果，即完成对图像的上下文分类；

(3)区域融合

以上下文分类所得到的分割图像作为输入，经过区域合并处理，得到一幅彩色图像，其中所包含RGB值的个数代表了合并之后的区域数总和；

对上述彩色图像进行多灰度值标记，然后对标记后的图像建立邻接矩阵，并进行如下处理：1)首先计算各区域面积，若其面积小于设定的阈值，则将其列入待合并区域序列，否则不作处理；2)对于待合并区域序列中的孤立区域，将其合并至其邻接区域，迭代执行至不再存在孤立区域；对于待合并区域序列中的非孤立区域，分别计算各孤立区域与其所有相邻的区域之间的公共边界的长度，将其合并至拥有最长公共边界的那个邻接区域中去；

经过上述处理即完成区域融合，得到融合后的分割图；

(4)后处理

对步骤(3)得到的分割图进行多灰度值标记，得到若干个连通区域，然后把每个连通区域映射到原始图像，得到对应位置上的局部原始图像，再利用特征来表征每个连通区域的属性，分别计算每个区域对应的局部原始图像的所述特征的特征值，其中上述特征包括边缘密度、亮度均值、方差和能量；

接着进行样本采集和样本特征计算处理，即选取能够代表不同地物特点的样本图像，计算上述各类特征的特征值，然后对采集的样本进行训练，得到区分不同地物的统计知识；

最后，根据训练得到的统计知识，利用不同地物类型的边缘密度的差异性，把区域划分为水域和陆地，再利用灰度差异将区域上的不同地物区分开，实现场景分类。

2.根据权利要求1所述的方法，其特征在于，所述的步骤(2)中的概率更新的公式如下：

$p_i^{(n+1)}(l_i=\mathrm{\α})=\frac{p_i^{\left(n\right)}(l_i=\mathrm{\α})\underset{(v_i,v_j)\∈E}{\mathrm{\Π}}{Q}_{\mathrm{ij}}^{\left(n\right)}\left(\mathrm{\α}\right)}{\underset{\mathrm{\λ}}{\mathrm{\Σ}}{p}_i^{\left(n\right)}(l_i=\mathrm{\λ})\underset{(v_i,v_j)\∈E}{\mathrm{\Π}}{Q}_{\mathrm{ij}}^{\left(n\right)}\left(\mathrm{\λ}\right)}---\left(1\right)$

其中，

表示像素点i属于类别α的概率；n表示迭代更新的次数，(v_i，v_j)表示像素点i及其8邻域中的任一像素点j组成的相邻像素对，且所有符合条件的像素对组成一个集合E；λ表示类别属性，Q为邻域函数，其作为邻域像素对中心点的类别属性的支撑度的度量，某个像素点的邻域函数Q是该像素点的8邻域内所有像素点的先验概率的加权和，邻域函数

表示第n步迭代时，邻域像素点j对中心像素点i的类别为α的支撑度。

3.根据权利要求2所述的方法，其特征在于，邻域函数

数学表达为：

$Q_{\mathrm{ij}}^{\left(n\right)}\left(\mathrm{\α}\right)=\underset{\mathrm{\β}\∈\mathrm{\Λ}}{\mathrm{\Σ}}{p}_j^{\left(n\right)}\left(\mathrm{\β}\right)r_{\mathrm{ij}}(\mathrm{\α},\mathrm{\β})---\left(2\right)$

其中，Λ是类别属性的集合且α∈Λ，β∈Λ，

表示第n步迭代时类别β的先验概率分布，r_ij(α，β)为权重，量化地表示像素i属于类别α，且像素j属于类别β的可能性，称之为兼容性系数。

4.根据权利要求3所述的方法，其特征在于，所述权重r_ij(α，β)的数学表达如下：

$r_{\mathrm{ij}}(\mathrm{\α},\mathrm{\β})=\frac{p(l_i=\mathrm{\α},l_j=\mathrm{\β})}{p(l_i=\mathrm{\α})p(l_j=\mathrm{\β})}---\left(3\right)$

其中，p(l_i＝α，l_j＝β)表示像素点i的类别为α且像素点j的类别为β的联合概率。

Images