CN103258324A - 基于可控核回归和超像素分割的遥感图像变化检测方法 - Google Patents

Abstract

本发明公开一种基于可控核回归和超像素分割的遥感图像变化检测方法，主要解决构造差异图只考虑图像的灰度信息，对其他特征信息利用不足以及直接对差异图采用k-means聚类容易造成弱变化区域无法检测的问题。其实现过程是：将输入的两时相图像采用可控核回归分别提取结构特征矩阵；将邻域特征矩阵分别与结构特征矩阵相结合，得到局部结构特征矩阵，利用非负矩阵分解算法分解局部结构特征矩阵，将得到的系数矩阵构造差异图；最后用超像素分割方法分割差异图得到过分割图像；对过分割图像采用K-means聚类，获得变化检测结果。本发明能够保持图像的边缘信息，抗噪性能较好，提高了变化检测精度，可用于灾情监测、土地利用、农业调查等领域。

Description

基于可控核回归和超像素分割的遥感图像变化检测方法

技术领域

本发明属于图像处理技术领域，涉及遥感图像变化检测，具体地说是一种基于可控核回归SKR和超像素分割的遥感图像变化检测方法，可用于对遥感图像变化区域的检测。

背景技术

遥感图像变化检测是对同一地区，不同时相获取的遥感影像进行对比、分析，得到两时相的变化区域。遥感图像变化检测在诸多行业中已经具有广泛应用，例如：水资源、土地、森林、草场等资源变化监测，海洋、湖泊、湿地、城区等的变化监测，海啸、地震、火灾、农作物病虫害等自然灾害评估，农作物生长状况评估，地球数据校正更新等方面；除此以外，遥感图像变化检测在军事方面具有重要应用，例如战场情报获取、军事目标侦查、军力部署情况调查等。随着遥感图像变化检测应用的日趋广泛，遥感图像变化检测方法的研究也成为国内外学者研究的重要方面。

首先对同一地区不同时相两幅遥感图像进行对比构造差异图，再对构造的差异图进行分类从而找到变化区域的先比较后分类方法是遥感图像变化检测的一类重要方法。一般的先分类后比较方法，对两幅图像进行直接相减构造差异图，再对构造的差异图利用分割、聚类方法进行分类。然而此类方法一般只使用图像灰度信息，没有考虑图像的其他信息，从而造成检测结果对噪声敏感，对弱变化区域漏检高等问题，同时，此类方法在对差异图进行聚类的时候，直接采用k-means聚类，容易造成弱变化区域难以检测的问题。

许多学者对先比较后分类的方法，进行了深入的研究。Wang等在文章“Unsupervised Change Detection for Remote Rensing Images Using multiscaleDecomposition and Treelet Fusion：A Level Set Approach，Proceedings of2011IEEECIE International Conference on Radar(RADAR2011),2011:1558-1561.”中提出了基于Treelet融合和水平集分割的变化检测方法。该方法对两个时相的图像分别进行小波分解，对分解后的对应子代系数进行做差，再对具有边缘信息的垂直和水平子代系数差进行soble算子边缘增强，然后进行重构和treelet融合，该方法能够减少图像中的伪变化信息，一定程度上保留了图像的边缘信息，但该方法没有考虑图像的局部信息，对于弱变化区域的边缘保持并不理想，同时treelet融合的结果在保留图像信息的同时也引入了冗余信息，导致最后结果虚警较高。

Zhang等人在文章“A Novel SAR Image Change Detection Based on Graph-cut andGeneralized Gaussian Model，IEEE Geoscience and Remote Sensing Letters，2012”中提出了基于graph-cut和广义高斯模型的变化检测算法。该方法首先对差异图进行graph-cut进行过分割，然后用模糊C均值聚类(Fuzzy C-means,FCM)算法进行过分割图像的初分割，在FCM初分割的基础上利用广义高斯模型进行差异图像分割，得到最后结果。该方法利用graph-cut和FCM进行初分割，考虑了图像的局部特性，初分割结果较好的反映了图像的实际情况，提高了最终分割结果的准确率，但该方法直接对差异图进行分割，使该方法对于噪声敏感，对于弱变化区域出现漏检，导致最终检测结果准确率不高。

可控核回归是一种自适应的核回归方法，综合考虑了图像的灰度信息和结构信息。该方法首先对图像局部结构进行初始估计，得到图像的梯度估计，然后采用该估计进行自适应地控制局部核，因此可以得到一个沿着图像的局部边缘结构进行扩散的狭长的具有椭圆形状的核函数，即自适应可控核核函数。对于图像的不同结构，该自适应可控核函数能够产生不同的核的形状，能够较好地对图像结构进行表示。可控核回归目前主要应用在图像去噪、图像插值和图像恢复等方面。

超像素分割在目标识别、图像分割等方面具有重要的应用，超像素表示的是一个空间上相邻且具有相似特征的集合，即表示图像中的一个区域块。相对于采用像素点来表示图像，采用超像素表示图像，图像的基元个数会大大减少，因此能够减少计算复杂度，运算速度很快。Liu等人在文章“Entropy Rate Superpixel Segmentation.2011IEEE Conference on Computer Vision and Pattern Recognition(CVPR),2097-2104”中提出了一种基于熵率的超像素分割方法。首先对图像构造无向图，利用无向图的随机游走的熵率来划分无向图，从而得到图像的分割结果，熵率能够使得图像的聚类具有紧密且匀质的聚类结果，同时引入了一个平衡因子项，能够使得图像的每个超像素具有相似的大小。该方法分割效果准确，边界保持效果好，且计算速度快。

发明内容

本发明的目的在于针对上述现有遥感图像变化检测方法的不足，提出了一种基于可控核回归和超像素分割的遥感图像变化检测方法，以保持变化区域的边缘，抑制图像的背景噪声，提高变化检测的准确性。

实现本发明目的技术思路是，将可控核回归和超像素分割方法用于变化检测中，利用可控核回归能较好地对图像结构进行表示的特点，将可控核回归用于提取图像结构特征，较好地保持变化区域的边缘，同时，利用超像素分割方法计算速度快，边界保持效果好的特点，保持弱化变化区域的边缘，减少背景噪声的影响。其具体步骤包括如下：

（1）输入同一地区不同时间获取的已配准的两幅遥感图像X₁和X₂，对图像X₁采用可控核回归提取结构特征矩阵W_S1，对图像X₂采用可控核回归提取结构特征矩阵W_S2；

（2）提取图像X₁的邻域特征矩阵N₁和图像X₂的邻域特征矩阵N₂，将邻域特征矩阵N₁与结构特征矩阵W_S1相结合，得到图像X₁的局部结构特征矩阵F₁，将邻域特征矩阵N₂与可控核提取的结构特征矩阵W_S2相结合，得到图像X₂的局部结构特征矩阵F₂；

（3）对图像X₁的局部结构特征矩阵F₁，采用非负矩阵分解方法进行分解，得到基矩阵B₁和系数矩阵H₁，对图像X₂的局部结构特征矩阵F₂，采用非负矩阵分解方法进行分解，得到基矩阵B₂和系数矩阵H₂，利用图像X₁的基矩阵B₁、系数矩阵H₁和图像X₂的基矩阵B₂、系数矩阵H₂构造差异图X_D；

（4）利用基于熵率的超像素分割方法分割差异图X_D得到分割后的图像X_D′；

（5）对分割后的图像X_D′采用k-means聚类，得到变化检测结果。

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

1、本发明使用可控核回归提取图像的邻域结构特征，考虑了图像的邻域信息，能够自适应的进行图像结构特征提取，较好的保持了图像的边缘，抑制了图像的噪声，保留了变化区域。

2、本发明采用超像素进行特征图像的过分割，然后再聚类，考虑了图像的局部性，有利于局部信息的保留，同时能够抑制杂点，避免了直接对特征图像聚类造成弱变化区域错分割的情况。

3、实验结果仿真表明，本发明方法与现有方法相比，减少了总错误数，提高了正确率。

附图说明

图1是本发明的实现流程图；

图2是本发明实验使用的一组模拟遥感图像数据集原始图像及变化参考图像；

图3是本发明实验使用的一组真实遥感图像数据集原始图像及变化参考图像；

图4是本发明对第一组实验数据采用不同方法进行变化检测的结果对比图；

图5是本发明对第二组实验数据采用不同方法进行变化检测的结果对比图。

具体实施方式

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

步骤1，输入同一地区不同时间获取的已配准的两幅遥感图像X₁和X₂。

输入同一地区获取的已配准的两幅大小均为I×J的两时相遥感图像：第一时相图像为X₁＝{X₁(i,j)|1≤i≤I,1≤j≤J}，第二时相图像为X₂＝{X₂(i,j)|1≤i≤I,1≤j≤J}，其中X₁(i,j)和X₂(i,j)分别为第一时相图像和第二时相图像在空间位置(i,j)处的像素点灰度值，i和j分别为图像的行序号和列序号，I为图像的总行数，J为图像的总列数。

步骤2，计算第一时相图像X₁和第二时相图像X₂的结构特征矩阵W_S1和W_S2。

计算结构特征矩阵的方法有非线性自适应插值方法，各向异性的方向扩散PDE方法，可控核回归方法等，本实施例给出的可控核回归方法，具体步骤如下：

2a）将第一时相图像X₁首列向左扩展w列，末列向右扩展w列，将列扩展后的图像再将其首行向上扩展w行，末行向下扩展w行，即可得到第一时相图像X₁的扩展后的图像Y₁，其中，

为向下取整符号，k是邻域滤波的图像块的大小，k的取值范围是11～17，本发明实例中k的取值为15。

2b）计算扩展后的第一时相图像Y₁的水平梯度图像Zx₁和垂直梯度图像Zy₁：

2b1）计算扩展后的第一时相图像Y₁的像素点(m,n)沿水平方向的导数，得到水平梯度Zx₁(m,n)，遍历扩展后的第一时相图像Y₁的所有的像素点，得到水平梯度图像Zx₁，其中m和n分别为扩展后的第一时相图像Y₁的行序号和列序号，m和n的取值范围分别为1≤m≤I+2w和1≤n≤J+2w，扩展后的第一时相图像Y₁的第m=1+w行为第一时相图像X₁的第i=1行，扩展后的第一时相图像Y₁的第n=1+w列为第一时相图像X₁的第j=1列；

2b2）计算扩展后的第一时相图像Y₁的每个像素点(m,n)沿垂直方向的导数，得到垂直梯度Zy₁(m,n)，遍历扩展后的第一时相图像Y₁的所有的像素点，得到水平梯度图像Zy₁；

2c）以水平梯度图像Zx₁中的像素点(i,j)为中心，选取大小为k×k的邻域块，采用半径为w的圆形均值滤波器对k×k的邻域块滤波，得到大小为k×k的水平梯度图像Zx₁的滤波后的矩阵

2d）对垂直梯度图像Zy₁采用与2c）步骤相同的操作得到垂直梯度图像Zy₁的滤波后的矩阵

2e）分别将水平梯度图像Zx₁的滤波后的矩阵

和垂直梯度图像Zy₁的滤波后的矩阵

拉成列向量，合并成大小为k²×2的局部梯度矩阵

2f）对局部梯度矩阵进行奇异值分解，得到k²×k²的对角矩阵

k²×2的特征值矩阵

及2×2特征向量矩阵

2g）计算第一时相图像X₁的协方差矩阵

$C_1^{\mathrm{ij}}={\mathrm{\&gamma;}}_1{\mathrm{\&alpha;}}_1{\mathrm{ra}}_1^{\mathrm{ij}}{\left({\mathrm{ra}}_1^{\mathrm{ij}}\right)}^T+{\mathrm{\&gamma;}}_1{\mathrm{\&beta;}}_1{\mathrm{rb}}_1^{\mathrm{ij}}{\left({\mathrm{rb}}_1^{\mathrm{ij}}\right)}^T---<1>$

其中，

为特征向量矩阵

的第一列，

为特征向量矩阵

的第二列，上标T为转置操作，γ₁为尺度因子，α₁为缩放因子，β₁为另一个缩放因子，

为能量矩阵，该矩阵是由特征值矩阵

的前2行构成，

为能量矩阵

的第一行第一列的元素，

为能量矩阵

的第二行第二列的元素，λ为正则参数，λ＞0，本发明实例中取值为1，ρ为结构因子，0＜ρ＜0.5，本发明实例取值为0.1；

2h）将2×2的协方差矩阵的第一行第一列的元素

作为第一个局部协方差矩阵C₁₁中元素C₁₁(i,j)的值，协方差矩阵的第一行第二列的元素

作为第二个局部协方差矩阵C₂₁中元素C₂₁(i,j)的值，协方差矩阵

的第二行第二列的元素作为第三个局部协方差矩阵C₃₁中元素C₃₁(i,j)的值，计算协方差矩阵

的行列式的平方根并将该值作为第四个局部协方差矩阵C₄₁中元素C₄₁(i,j)的值；

2i）重复步骤2c）至步骤2h）直到计算完水平梯度图像Zx₁和垂直梯度图像Zy₁的所有像素点(i,j)，得到大小均为I×J的第一个局部协方差矩阵C₁₁、第二个局部协方差矩阵C₂₁、第三个局部协方差矩阵C₃₁和第四个局部协方差矩阵C₄₁；

2j）将四个局部协方差矩阵C₁₁、C₂₁、C₃₁和C₄₁分别进行首列向左扩展o列，末列向右扩展o列，将列扩展完成后的矩阵再将其首行向上扩展o行，末行向下扩展o行，得到四个扩展后的局部协方差矩阵C′₁₁、C′₂₁、C′₃₁和C′₄₁；

2k）分别以扩展后的局部协方差矩阵C′₁₁、C′₂₁、C′₃₁和C′₄₁的像素点(i,j)为中心，选取大小为z×z的邻域块，分别得到四个邻域块协方差矩阵

和

其中，

z是可控核的大小，z的取值范围是3～7，本发明实例中取值为5；

2l）计算大小为z×z的可控核矩阵

$W_1^{\mathrm{ij}}=M_{41}^{\mathrm{ij}}\&CenterDot;\exp (-(0.5/h^2)\&times;(b_r^{\mathrm{ij}}\&CenterDot;M_{11}^{\mathrm{ij}}\&CenterDot;b_r^{\mathrm{ij}}+M_{21}^{\mathrm{ij}}\&CenterDot;b_c^{\mathrm{ij}}+b_c^{\mathrm{ij}}\&CenterDot;M_{21}^{\mathrm{ij}}\&CenterDot;b_r^{\mathrm{ij}}+M_{31}^{\mathrm{ij}}\&CenterDot;b_c^{\mathrm{ij}}))---<2>$

其中，h为平滑参数，0＜h＜10，本发明实例中取值为3，exp(·)表示对目标函数求指数，

为z×z大小的行距离矩阵，其元素值的邻域块为z×z中像素点到其中心像素点(i,j)的行距离，

$b_r^{\mathrm{ij}}=\left[\begin{array}{ccccc}-2& -2& -2& -2& -2\\ -1& -1& -1& -1& -1\\ 0& 0& 0& 0& 0\\ 1& 1& 1& 1& 1\\ 2& 2& 2& 2& 2\end{array}\right]---<3>$

为z×z大小的列距离矩阵，其元素值为z×z的邻域块中像素点到其中心像素点(i,j)的列距离，

$b_c^{\mathrm{ij}}=\left[\begin{array}{ccccc}-2& -1& 0& 1& 2\\ -2& -1& 0& 1& 2\\ -2& -1& 0& 1& 2\\ -2& -1& 0& 1& 2\\ -2& -1& 0& 1& 2\end{array}\right]---<4>$

2m）重复步骤2k）至步骤2l），逐像素(i,j)进行计算，得到第一时相图像的像素点(i,j)的可控核矩阵将所有像素点的可控核矩阵

均转化为z²×1的列向量，合并所有I×J个列向量，构成大小为z²×(I×J)的结构特征矩阵W_S1；

2n）对于第二时相图像，重复步骤2a）至步骤2m），得到第二时相图像X₂的所有像素点的可控核矩阵

合并所有I×J个列向量，构成大小为z²×(I×J)的结构特征矩阵W_S2。

步骤3，利用第一时相图像X₁的邻域特征矩阵N₁和结构特征矩阵W_S1构造第一时相图像X₁的局部结构特征矩阵F₁，对第二时相图像X₂采用同样地操作得到第二时相图像X₂的局部结构特征矩阵F₂。

3a)将第一时相图像X₁首列向左扩展o列，末列向右扩展o列，对列扩展完成后的图像再将其首行向上扩展o行，末行向下扩展o行，得到第一时相图像X₁扩展后的图像X_s1；

3b）逐个以扩展后的第一时相图像X_s1中的非边界扩展像素(i,j)为中心，选取大小为z×z的邻域块

并将

转化成z²×1的列向量

将所有非边界扩展像素(i,j)的邻域块的列向量合并构成大小为z²×(I×J)的邻域特征矩阵 $N_1=\left\{D_1^{\mathrm{ij}}\right|1\&le;i\&le;I,1\&le;j\&le;J\};$

3c）生成一个大小为z²×(I×J)的局部结构特征矩阵F₁，该矩阵中的元素值赋为邻域特征矩阵N₁与结构特征矩阵W_S1中对应元素位置的矩阵元素值的乘积；

3d）对于第二时相图像X₂，重复步骤3a）至步骤3c），得到第二时相图像的局部结构特征矩阵F₂。

步骤4，分解局部结构特征矩阵F₁和局部结构特征矩阵F₂。

分解局部结构特征矩阵的方法有主成分分析方法、二维主成分分析方法、独立成分分析方法、非负矩阵分解方法等，本实施例采用参考文献“Daniel D L.,H SebastianS..Algorithms for non-negative matrix factorization.Advances in Neural InformationProcessing Systems,2001,13:556-562.”中的非负矩阵分解方法，利用该方法分解局部结构特征矩阵F₁和局部结构特征矩阵F₂，得到第一时相图像X₁的基矩阵B₁、系数矩阵H₁和第二时相图像X₂的基矩阵B₂、系数矩阵H₂。

步骤5，利用第一时相图像X₁的基矩阵B₁、系数矩阵H₁和第二时相图像X₂的基矩阵B₂、系数矩阵H₂构造差异图X_D。

5a）令第二时相图像X₂的基矩阵B₂和第一时相图像X₁的基矩阵B₁之间的变换矩阵为P，即B₂＝B₁·P；

5b）将局部结构特征矩阵F₁和F₂的差在基矩阵B₁上投影，得到系数差异列向量：H＝H₁-PH₂＝H₁-B₁ ^-1B₂H₂，B₁ ^-1为基矩阵B₁的逆矩阵；

5c）将列向量H转化为大小为I×J的图像，得到差异图像X_D。

步骤6，利用基于熵率的超像素分割差异图像X_D得到分割后的图像X′_D。

6a）构造初始无向图G₀＝(V,E)，其中V表示初始无向图中所有节点的集合，E={e_mn}表示初始无向图中所有边的集合，e_mn为m节点和n节点之间的连线形成的边；

6b）对初始无向图G₀进行分割，得到分割后的无向图G＝(V,A)，使得差异图像X_D划分为K个超像素区域，其中A为子集，该子集中的元素是从集合E中选择边构成的，A∈E，K为超像素区域的个数，该分割的具体步骤如下：

6b1）首先初始化超像素区域的个数为K＝I×J/20；

6b2）计算差异图像X_D的两个像素u和v间的相似度权重w(u,v)：

$w(u,v)=\exp (-\frac{{(X_{\mathrm{Du}}-X_{\mathrm{Dv}})}^{2}d(X_{\mathrm{Du}},X_{\mathrm{Dv}})}{{2\mathrm{\&sigma;}}^2})---<5>$

其中X_Du和X_Dv分别为像素点u和v的灰度值，d(X_Du,X_Dv)为像素点u和v的灰度值之间的欧氏距离，u＝1,...,I×J，v＝1,...,I×J，σ为宽度参数，0＜σ＜5，本发明实例中取值为3；

6b3）计算分割后的无向图G＝(V,A)的随机游走的熵率ψ(A)：

$\mathrm{\&psi;}\left(A\right)=-\underset{u}{\mathrm{\&Sigma;}}{\mathrm{\&mu;}}_u\underset{v}{\mathrm{\&Sigma;}}{p}_{\mathrm{uv}}\left(A\right)\log \left(p_{\mathrm{uv}}\left(A\right)\right)---<6>$

其中，p_uv(A)为转移概率，其计算公式如下，

其中，w_u表示的是所有与点u相连接的点的权重和，即e_uv为子集A中的像素点u和像素点v之间的边，l为与像素点u相连接的像素点，e_ul为子集A中的像素点u和像素点l之间的边，w(u,l)为像素点u和l的相似度权重值；

μ_u为平稳分布因子，

上标T为转置操作，

是归一化常数，表示的是对所有节点的w_u进行求和，|V|表示节点的个数；

6b4）将使下式的值最大权值w(u,v)对应的两个像素点u和v合并，即将像素点u和像素点v之间的边e_uv加入到子集A中，使像素点u和像素点v有相同的标记，

max ψ(A)    <8>

s.t.A∈E and NA≥K

其中，s.t.是约束条件的英文缩写，NA是每次迭代得到的差异图像X_D的超像素区域的个数；

6b5）如果像素点u和像素点v之间的边e_uv被选择加入到子集A中，则像素点u和像素点v的自身的权重值保持不变，u和v的边的权重值也保持不变；如果像素点u和像素点v之间的边e_uv没有被选择，则调整像素点u和像素点v的自身的权重值 $\left\{\begin{array}{c}w(u,u)=w(u,u)+w(u,v)\\ w(v,v)=w(v,v)+w(u,v)\end{array}\right.,$ 并调整像素点u和v的边的权重值w(u,v)=0；

6b6）重复迭代步骤6b2）至步骤6b5），直到NA＝K停止迭代，使差异图像X_D分割为K个超像素区域；

6c）计算每个超像素区域中的所有像素点的灰度均值，在差异图像X_D中将每个超像素区域的像素点的灰度值用该超像素的灰度均值代替得到分割后的图像X′_D。

步骤7，对分割后的图像X′_D进行聚类，得到变化检测结果。

聚类的方法有FCM方法、KFCM方法、k-means方法等，本实施例采用参考文献“MacQueen,J.B..Some Methods for classification and Analysis of MultivariateObservations.Proceedings of5th Berkeley Symposium on Mathematical Statistics andProbability.Berkeley,CA,1967,1:281–297.”中的k-means聚类方法。

本发明的效果可以通过以下实验结果和分析进一步说明：

1.实验数据

本发明的实验数据是一组模拟实验数据和一组真实的遥感图像数据，其中模拟数据集的原始图像为ATM(Airborne Thematic Mapper)3波段位于英国Feltwell村庄的一个农田区的图像，如图2(a)所示，其模拟变化图像是通过模拟地球的天气变化和电磁波的辐射特性等因素影响并人工地嵌入一些变化区域得到的，如图2(b)所示。两幅图像大小均为470×335像素，灰度级为256，配准误差为1.5个像素左右，变化参考图如图2(c)所示，图2(c)中白色像素区表示变化的区域，变化的像素数为4236，非变化像素数为153214。真实的遥感图像数据集是由1995年9月和1996年7月在意大利撒丁岛Mulargia湖泊区域的Landsat5TM第5波段光谱图像组成，分别如图3(a)和3(b)所示。两幅图像大小均为300×412像素，灰度级为256，它们之间的变化为湖水水位上涨引起的，变化区域的参考图如图3(c)所示，图3(c)中白色像素区表示变化的区域，非变化像素为115974，变化的像素数为7626。

2.实验方法

为了验证基于可控核回归和超像素分割的遥感图像变化检测方法的实验效果，我们将本发明方法与以下现有方法进行对比。

方法1：Wang等在文章“Unsupervised Change Detection for Remote Rensing ImagesUsing multiscale Decomposition and Treelet Fusion：A Level Set Approach，Proceedings of2011IEEE CIE International Conference on Radar(RADAR2011),2011:1558-1561.”中提出的基于Treelet融合和水平集分割的变化检测方法。

方法2：Wang等在文章“Change Detection in SAR Image Based oN MultiscaleProduct and PCA,The Second International Asia-Pacific Conference on Synthetic ApertureRadar(APSAR2009),Xian,China.2009:872-875.”中提出的基于多尺度积和主成分分析的变化检测方法。

方法3：Zhang等人在文章“A Novel SAR Image Change Detection Based onGraph-cut and Generalized Gaussian Model，IEEE Geoscience and Remote SensingLetters，2012”中提出的基于graph-cut和广义高斯模型的变化检测算法。

3.实验评价指标

实验使用的评价指标是虚警数、漏检数和总错误数。虚警数是实际没有发生变化但被当作变化检测出来的像素的总和，漏检数是没有检测出来的实际发生了变化的像素的总和，总错误数是虚警数和漏检数之和。

4.实验内容与结果

仿真实验一：采用现有方法1、现有方法2、现有方法3和本发明方法对模拟实验数据集进行变化检测，得到的变化检测结果如图4所示，其中：

图4(a)是现有方法1的变化检测结果，图4(b)是现有方法2的变化检测结果，图4(c)是现有方法3的变化检测结果，图4(d)是本发明的变化检测结果。

从图4(a)可以看出，现有方法1的检测结果图中存在较多漏检，从图4(b)看出，现有方法2的检测结果对噪声敏感，存在较多虚警点，从图4(c)看出，现有方法3的检测结果存在较多漏检，图像的边缘没有得到较好地保持，从图4(d)看出，本发明方法的检测结果漏检少，边缘保持得较好，较接近实际变化参考图。

仿真实验二：采用现有方法1、现有方法2、现有方法3和本发明方法对真实遥感图像实验数据集进行变化检测，得到的变化检测结果如图5所示，其中：

图5(a)是现有方法1的变化检测结果，图5(b)是现有方法2的变化检测结果，图5(c)是现有方法3的变化检测结果，图5(d)是本发明的变化检测结果。

从图5(a)可看出，现有方法1的检测结果图中存在较多的漏检，从图5(b)看出，现有方法2的检测结果存在较多的杂点，虚警较多，从图5(c)看出，现有方法3的检测结果也存在较多的杂点，变化区域边缘较粗，虚警较多。从图5(d)看出，本发明方法的检测结果的边缘保持较好，准确地检测出变化区域，较接近实际的变化参考图。

下面从漏检数、虚警数、总错误数和正确率四个方面客观评价本发明方法，两组实验数据集的结果如表1所示。

表1两组遥感图像数据集采用不同方法变化检测结果的性能评价

从表1可以看出，对于模拟实验数据集，本发明方法的漏检最少，总错误数最少，正确率最高，整体性能最好，本发明方法的漏检数比现有方法1的漏检数少了223个像素点，比现有方法2的漏检数少了52个像素点，比现有方法3的漏检数少了234个像素点；本发明方法的总错误数比现有方法1的总错误数少了75个像素点，比现有方法2的总错误数少了582个像素点，比现有方法3的总错误数少了420个像素点。

对于真实遥感图像实验数据集，本发明方法虚警最少，总错误数最少，正确率最高，整体性能最优，本发明方法的虚警数比现有方法1的虚警数少了45个像素点，比现有方法2的虚警数少了1211个像素点，比现有方法3的虚警数少了2548个像素点；总错误数比现有方法1的总错误数少了305个像素点，比现有方法2的总错误数少了914个像素点，比现有方法3的总错误数少了2276个像素点。

综上，本发明方法能够较全面地检测出变化区域信息，较好地保持变化区域边缘，具有较高的检测精度。

Claims (4)

1.一种基于可控核回归和超像素分割的遥感图像变化检测方法，其特征在于：包括有如下步骤：

（1）输入同一地区不同时间获取的已配准的两幅遥感图像X₁和X₂，对图像X₁采用可控核回归提取结构特征矩阵W_S1，对图像X₂采用可控核回归提取结构特征矩阵W_S2；

（2）提取图像X₁的邻域特征矩阵N₁和图像X₂的邻域特征矩阵N₂，将邻域特征矩阵N₁与结构特征矩阵W_S1相结合，得到图像X₁的局部结构特征矩阵F₁，将邻域特征矩阵N₂与可控核提取的结构特征矩阵W_S2相结合，得到图像X₂的局部结构特征矩阵F₂；

（3）对图像X₁的局部结构特征矩阵F₁，采用非负矩阵分解方法进行分解，得到基矩阵B₁和系数矩阵H₁，对图像X₂的局部结构特征矩阵F₂，采用非负矩阵分解方法进行分解，得到基矩阵B₂和系数矩阵H₂，利用图像X₁的基矩阵B₁、系数矩阵H₁和图像X₂的基矩阵B₂、系数矩阵H₂构造差异图像X_D；

（4）利用基于熵率的超像素分割方法分割差异图像X_D得到分割后的图像X_D′；

（5）对分割后的图像X_D′采用k-means聚类，得到变化检测结果。

2.根据权利要求1所述的基于可控核回归和超像素分割的遥感图像变化检测方法，其特征在于：所述步骤（1）中对输入的两时相图像X₁、X₂采用可控核回归分别提取结构特征矩阵，按如下步骤进行：

1a）对输入的图像X₁的首列向左扩展w列，末列向右扩展w列；对列扩展后的图像再将其首行向上扩展w行，末行向下扩展w行，即可得到图像X₁的扩展后的图像Y₁，其中，

k是邻域滤波的图像块的大小，

为向下取整符号；

1b）计算图像X₁的扩展后的图像Y₁的水平梯度图像Zx₁和垂直梯度图像Zy₁；

1c）以水平梯度图像Zx₁的像素点(i,j)为中心，选取大小为k×k的邻域块，采用半径为w的圆形均值滤波器对k×k的邻域块进行滤波，得到大小为k×k的所述水平梯度图像Zx₁的滤波后矩阵

对垂直梯度图像Zy₁采用相同的操作得到图像Zy₁的滤波后的矩阵

分别将图像Zx₁的滤波后的矩阵

和图像Zy₁的滤波后的矩阵

拉成列向量，合并成大小为k²×2的局部梯度矩阵

1d）对局部梯度矩阵

进行奇异值分解，得到k²×k²的对角矩阵

k²×2的特征值矩阵

及2×2特征向量矩阵取特征值矩阵

的前2行构成2×2的能量矩阵

计算协方差矩阵

$C_1^{\mathrm{ij}}={\mathrm{\&gamma;}}_1{\mathrm{\&alpha;}}_1{\mathrm{ra}}_1^{\mathrm{ij}}{\left({\mathrm{ra}}_1^{\mathrm{ij}}\right)}^T+{\mathrm{\&gamma;}}_1{\mathrm{\&beta;}}_1{\mathrm{rb}}_1^{\mathrm{ij}}{\left({\mathrm{rb}}_1^{\mathrm{ij}}\right)}^T,$

其中，

为特征向量矩阵

的第一列，

为特征向量矩阵

的第二列，上标T为转置操作，γ₁为尺度因子， ${\mathrm{\&gamma;}}_1={\left(\frac{s_1^{\mathrm{ij}}\left(\mathrm{1,1}\right)\&times;s_1^{\mathrm{ij}}\left(\mathrm{2,2}\right)+\mathrm{\&lambda;}}{k^2}\right)}^{\mathrm{\&rho;}},$ α₁为缩放因子， ${\mathrm{\&alpha;}}_1=\frac{s_1^{\mathrm{ij}}\left(\mathrm{1,1}\right)+\mathrm{\&lambda;}}{s_1^{\mathrm{ij}}\left(\mathrm{2,2}\right)+\mathrm{\&lambda;}},$ β₁是另一个缩放因子，

为能量矩阵

的第一行第一列的元素，

为能量矩阵

的第二行第二列的元素，λ为正则参数，λ＞0，ρ为结构因子，0＜ρ＜0.5；

1e）取2×2的协方差矩阵

的第一行第一列的元素

作为第一个局部协方差矩阵C₁₁中元素C₁₁(i,j)的值，协方差矩阵

的第一行第二列的元素

作为第二个局部协方差矩阵C₂₁中元素C₂₁(i,j)的值，协方差矩阵

的第二行第二列的元素

作为第三个局部协方差矩阵C₃₁中元素C₃₁(i,j)的值，计算协方差矩阵

的行列式的平方根并将该值作为第四个局部协方差矩阵C₄₁中元素C₄₁(i,j)的值；

1f）重复步骤1c）至步骤1e）直到计算完水平梯度图像Zx₁和垂直梯度图像Zy₁的所有像素1≤i≤I且1≤j≤J的邻域块的协方差矩阵

得到大小均为I×J的第一个局部协方差矩阵C₁₁、第二个局部协方差矩阵C₂₁、第三个局部协方差矩阵C₃₁和第四个局部协方差矩阵C₄₁，其中I为图像的总行数，J为图像的总列数；

1g）对四个局部协方差矩阵C₁₁、C₂₁、C₃₁和C₄₁分别进行首列向左扩展o列，末列向右扩展o列，对列扩展完成后的矩阵再进行首行向上扩展o行，末行向下扩展o行，得到局部协方差矩阵C₁₁、C₂₁、C₃₁和C₄₁的扩展后局部协方差矩阵C′₁₁、C′₂₁、C′₃₁和C′₄₁；分别对局部协方差矩阵C₁₁、C₂₁、C₃₁和C₄₁的扩展后的局部协方差矩阵C′₁₁、C′₂₁、C′₃₁和C′₄₁以像素点(i,j)为中心，选取大小为z×z的邻域块，分别得到四个邻域块协方差矩阵

和

其中，

z是可控核的大小；

1h）计算大小为z×z的可控核矩阵

$W_1^{\mathrm{ij}}=M_{41}^{\mathrm{ij}}\&CenterDot;\exp (-(0.5/h^2)\&times;(b_r^{\mathrm{ij}}\&CenterDot;M_{11}^{\mathrm{ij}}\&CenterDot;b_r^{\mathrm{ij}}+M_{21}^{\mathrm{ij}}\&CenterDot;b_c^{\mathrm{ij}}+b_c^{\mathrm{ij}}\&CenterDot;M_{21}^{\mathrm{ij}}\&CenterDot;b_r^{\mathrm{ij}}+M_{31}^{\mathrm{ij}}\&CenterDot;b_c^{\mathrm{ij}})),$

其中，h为平滑参数，exp(·)表示求目标函数的指数，为z×z大小的行距离矩阵，其元素值为z×z的邻域块中的像素点到其中心像素点(i,j)的行距离，

$b_r^{\mathrm{ij}}=\left[\begin{array}{ccccc}-2& -2& -2& -2& -2\\ -1& -1& -1& -1& -1\\ 0& 0& 0& 0& 0\\ 1& 1& 1& 1& 1\\ 2& 2& 2& 2& 2\end{array}\right],$

为z×z大小的列距离矩阵，其元素值为z×z的邻域块中的像素点到其中心像素点(i,j)的列距离，

$b_c^{\mathrm{ij}}=\left[\begin{array}{ccccc}-2& -1& 0& 1& 2\\ -2& -1& 0& 1& 2\\ -2& -1& 0& 1& 2\\ -2& -1& 0& 1& 2\\ -2& -1& 0& 1& 2\end{array}\right],$

1i）重复步骤（1g）至步骤（1h），逐像素(i,j)进行计算，得到第一时相图像X₁的像素点(i,j)的可控核矩阵将所有像素点的可控核矩阵均转化为z²×1的列向量，合并所有I×J个列向量构成大小为z²×(I×J)的结构特征矩阵W_S1；

1j）对于第二时相图像X₂，重复步骤（1a）至步骤（1i），得到第二时相图像X₂的所有像素点的可控核矩阵

合并所有I×J个由可控核矩阵

转化成的列向量，构成大小为z²×(I×J)的结构特征矩阵W_S2。

3.根据权利要求1所述的基于可控核回归和超像素分割的遥感图像变化检测方法，其特征在于：步骤（3）所述的利用图像X₁的基矩阵B₁、系数矩阵H₁和图像X₂的基矩阵B₂、系数矩阵H₂构造差异图像X_D，按如下步骤进行：

3a）令图像X₂的基矩阵B₂和图像X₁的基矩阵B₁之间的变换矩阵为P，即B₂＝B₁·P；

3b）将局部结构特征矩阵F₁和F₂的差在图像X₁的基矩阵B₁上投影，得到系数差异列向量：H＝H₁-PH₂＝H₁-B₁ ^-1B₂H₂，B₁ ^-1为基矩阵B₁的逆矩阵；

3c）将列向量H转化为大小为I×J的图像，得到差异图像X_D。

4.根据权利要求1所述的基于可控核回归和超像素分割的遥感图像变化检测方法，其特征在于：步骤（4）所述的利用基于熵率的超像素分割方法分割差异图像X_D得到分割后图像X_D′，按如下步骤进行：

4a）构造初始无向图G₀＝(V,E)，其中V表示初始无向图G₀中所有节点的集合，E={e_mn}表示初始无向图G₀中所有边的集合，这里e_mn为m节点和n节点之间的连线形成的边，将图的分割问题看作是子集选择的问题，从集合E中选择边构成子集A，A∈E，得到分割无向图G＝(V,A)，使得差异图像X_D分割为K个超像素区域，初始化超像素区域的个数为K＝I×J/50；

4b）计算差异图像X_D的两个像素u和v间的相似度权重w(u,v)：

$w(u,v)=\exp (-\frac{{(X_{\mathrm{Du}}-X_{\mathrm{Dv}})}^{2}d(X_{\mathrm{Du}},X_{\mathrm{Dv}})}{{2\mathrm{\&sigma;}}^2}),$

其中X_Du为像素点u的灰度值，X_Dv为像素点v的灰度值，d(X_Du,X_Dv)为像素点u和v的灰度值之间的欧氏距离，u＝1,...,I×J，v＝1,...,I×J，σ为宽度参数，1＜σ＜5；

4c）计算分割无向图G＝(V,A)的随机游走的熵率ψ(A)：

$\mathrm{\&psi;}\left(A\right)=-\underset{u}{\mathrm{\&Sigma;}}{\mathrm{\&mu;}}_u\underset{v}{\mathrm{\&Sigma;}}{p}_{\mathrm{uv}}\left(A\right)\log \left(p_{\mathrm{uv}}\left(A\right)\right),$

其中，p_uv(A)为转移概率，

其计算公式如下：

其中，w_u表示的是与像素点u相连接的所有像素点的权重和，

e_uv为子集A中的像素点u和像素点v之间的边，l为与像素点u相连接的像素点，e_ul为子集A中的像素点u和像素点l之间的边，w(u,l)为像素点u和l的相似度权重值；

μ_u为平稳分布因子，

其中，上标T为转置操作，

是归一化常数，表示的是对所有节点的w_u进行求和，|V|为节点的个数；

4d）将使下式的值最大的权值w(u,v)对应的两个像素点u和v合并，即将像素点u和像素点v之间的边e_uv加入到子集A中，使像素点u和像素点v有相同的标记，

max ψ(A)

s.t.A∈E and NA≥K，

其中，s.t.是约束条件的英文缩写，NA是每次迭代得到的差异图像X_D的超像素区域的个数；

4e）如果像素点u和像素点v之间的边e_uv被选择加入到子集A中，则像素点u和像素点v的自身的权重值保持不变，u和v的边的权重值也保持不变；如果像素点u和像素点v之间的边e_uv没有被选择，则调整像素点u和像素点v自身的权重值， $\left\{\begin{array}{c}w(u,u)=w(u,u)+w(u,v)\\ w(v,v)=w(v,v)+w(u,v)\end{array}\right.,$ 并调整像素点u和v的边的权重值w(u,v)=0；

4f）重复迭代步骤（4b）至步骤（4e），直到NA＝K停止迭代，使差异图像X_D分割为K个超像素区域；

4g）计算每个超像素区域中的所有像素点的灰度均值，在差异图像X_D中将每个超像素区域的像素点的灰度值用该超像素的灰度均值代替，得到分割后图像X_D′。

Images