CN103226825B - 基于低秩稀疏模型的遥感图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于低秩稀疏模型的遥感图像的变化检测方法，主要解决现有技术中构造差异图信息丢失严重的问题。其实现过程是：首先，用对数比值法从变化前后的遥感图像中提取出初始变化区域；其次，用该初始变化区域与变化前后图像，构造模拟遥感图像序列；然后，用低秩稀疏分解算法，将模拟遥感图像序列分解为低秩矩阵，稀疏矩阵和噪声矩阵三个部分；最后，用模糊C均值的方法对稀疏矩阵最后一列的列向量进行分类，得到最终的变化检测结果图。本发明具有检测精度高，信息丢失和累积误差少的优点，可用于目标检测与识别，图像分割以及机器学习领域。

Description

基于低秩稀疏模型的遥感图像变化检测方法

技术领域

本发明属于图像处理技术领域，特别涉及遥感图像的处理，可应用于生态和环境的监控以及自然灾害的评测与预防。

背景技术

遥感图像变化检测是指选用合适的检测方法，从不同时期获取的多幅遥感图像中提取出变化信息并加以分析，生成变化分布图和其他检测结果的技术。目前遥感图像变化检测技术变化检测已成为遥感图像处理研究的一个研究重点，被广泛的应用到社会经济的各个领域，比如灾害监测和评估、土地使用情况的分析、水资源质量和地理分布情况调查、城市的规划及布局、气候的变化监测、战场态势的评估等。

目前遥感图像变化检测还处于初步阶段，大致有两个路线，一条路线是分类后比较方法，该方法先对多时相遥感图像单独进行分类标记，然后根据相应像素类别的差异来提取发生变化的区域/像素，最后得到变化检测结果；另外一条路线是图像直接比较方法，该方法先对多时相遥感图像进行逐像素的比较，比如逐像素点的差值、比值、CVA等，再对比较得出的差异图像进行进一步地多种变换、概率分布等处理以达到两分类，最后得到变化检测结果。图像直接比较方法能够减少由于数据获取平台和环境的不同引起的伪变化信息，不需要数据的辐射校正等复杂的预处理，目前较多的研究是沿着图像直接比较方法这个路线进行的。

目前变化检测方法的缺陷主要有：对于图像直接比较方法而言，首先，差异图可分性的高低直接影响着变化检测结果，而高可分性差异图的构建至今是一个难题；其次，简单的图像直接比较方法，既没有利用到图像的空间信息，也没有考虑到噪声或光照变化造成的影响，很容易受噪声因素的影响，使检测不精确；最后，在构造差异图的过程中不可避免的造成了图像信息的丢失，图像直接比较方法把变化信息压缩到一维的差异图像上，损失了遥感图像的波段信息，从而使检测不准确。

发明内容

本发明的目的在于针对上述已有技术的不足，提出了一种基于低秩稀疏模型的遥感图像变化检测方法，以降低构造差异图降维造成的较高漏检率，并考虑图像像素邻域之间的信息，降低了遥感图像变化检测的漏检率。

为实现上述目的，本发明基于低秩稀疏模型的遥感图像变化检测方法，包括如下步骤：

（1）输入变化前的遥感图像P₁和变化后的遥感图像P₂，其中，变化前后图像大小均为m×n；

（2）利用变化前的遥感图像P₁和变化后的遥感图像P₂，构造k幅模拟变化后遥感图像X_i,i=1,2,…,k，且模拟变化后遥感图像X_i的大小为m×n；

（3）由变化前的遥感图像P₁，变化后的遥感图像P₂以及k幅模拟变化后遥感图像X_i,i=1,2,…,k，构造模拟遥感图像序列X=[p₁,…x_i,…p₂]，其中，p₁，x_i，p₂分别为将变化前图像P₁，模拟变化后图像X_i以及变化后图像P₂，拉成列向量后的列向量，则p₁，x_i，p₂∈R^mn×1，X∈R^mn×(k+2)，R^mn×(k+2)表示大小为mn×(k+2)的实数空间；

（4）用GODEC低秩稀疏分解算法，对模拟遥感图像序列X进行低秩稀疏分解：

X=S+L+G，

其中，S为稀疏矩阵，L为低秩矩阵，G为噪声矩阵，S,L,G∈R^mn×(k+2)，且S＝{s₁,s₂,…s_l,…,s_k+2}，l=1,2,…,(k+2)，s_l∈R^mn×1为稀疏矩阵S的列向量；

（5）用模糊C均值方法对稀疏矩阵S中的最后一列的列向量s_k+2∈R^mn×1进行聚类，得到每个像素对应的类别标签，最终得到变化检测结果图。

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

1.本发明针对多时相遥感图像之间的稀疏性，利用低秩稀疏模型提取出多时相遥感图像中的稀疏信息，较好保存了遥感图像之间的波段信息，提高了遥感图像的变化检测准确率；

2.本发明降低了构造差异图时，将多时相遥感图像压缩到一维差异图中所造成的信息丢失和累积误差问题，提高了变化检测的准确率；

3.本发明中的双边随机投影方法，既考虑了单时相遥感图像内的邻域信息，又考虑了多时相遥感图像之间的差异信息，降低了变化检测的漏检率；

附图说明

图1是本发明的流程图；

图2是本发明仿真所使用的模拟遥感图像数据集；

图3是通过人工标记得到图2的标准变化检测结果图；

图4是用本发明和传统变化检测方法对图2的变化检测结果图；

图5是本发明仿真所使用的墨西哥遥感图像数据集；

图6是通过人工标记得到图5的标准变化检测结果图；

图7是用本发明和传统变化检测方法对图5的变化检测结果图。

具体实施方式

参照图1，本发明基于低秩稀疏模型的遥感图像变化检测方法，包括如下步骤：

步骤1:输入变化前遥感图像P₁和变化后遥感图像P₂，其中，变化前后的遥感图像大小均为m×n。

步骤2:构造k幅模拟变化后遥感图像。

2a）用对数比值法，从变化前的遥感图像P₁和变化后的遥感图像P₂中，提取出变化后遥感图像P₂的初始变化区域P₀；

2b）利用变化前的遥感图像P₁，变化后的遥感图像P₂以及初始变化区域P₀，构造k幅模拟变化后遥感图像X_i：

$X_i=\left\{\begin{array}{cc}{F}_i& F_i\stackrel{F}{\←}\frac{i}{k}{P}_0\\ U_i& U_i\stackrel{U}{\←}{P}_1-\frac{i}{k}{P}_0\end{array}\right.,$

其中，模拟变化后遥感图像X_i,i=1,2,…,k，10≤k≤100，k为自然数，且模拟变化后遥感图像X_i的大小为m×n；

F_i表示模拟变化后遥感图像X_i中发生变化的区域；

U_i表示模拟变化后遥感图像X_i中没有发生变化的区域；

式中，表示将初始变化区域P₀划分为k等份，并取其中的i份，表示将变化后图像P₂中的部分变化区域的灰度值，迁移到模拟变化后遥感图像X_i中的对应位置上；

式中，表示变化前遥感图像P₁中除与部分变化区域位置相同的区域之外的所有区域，表示将变化前遥感图像P₁中除与部分变化区域位置相同的区域之外的所有区域的灰度值，迁移到模拟变化后遥感图像X_i中的对应位置上。

步骤3:由变化前的遥感图像P₁，变化后的遥感图像P₂以及k幅模拟变化后遥感图像X_i,i=1,2,…,k，构造模拟遥感图像序列X=[p₁,…x_i,…p₂]，其中，p₁，x_i，p₂分别为将变化前图像P₁，模拟变化后图像X_i以及变化后图像P₂，拉成列向量后的列向量，且p₁，x_i，p₂∈R^mn×1，X∈R^mn×(k+2)，R^mn×(k+2)表示大小为mn×(k+2)的实数空间。

步骤4:用GODEC低秩稀疏分解算法，对模拟遥感图像序列X进行低秩稀疏分解。

4a）用双边随机投影的方法，对模拟遥感图像序列X进行低秩逼近，得到迭代低秩矩阵L_t为：

$L_t=Q_1{\left[R_1{\left({A_{2t}}^T{Y}_{1t}\right)}^{-1}{R_2}^T\right]}^{1/(2q+1)}{Q_2}^T,$

其中，t≥1为迭代次数，A_2t＝XA_1(t-1)，当t＝1时，A_1(t-1)是初始化为秩为r，2≤r≤30的高斯矩阵，当t≥2时，A_1(t-1)＝X^TA_2(t-1)，这里，q≥0为收敛系数，且Q₁,R₁和Q₂,R₂分别为对左投影矩阵Y_1t和右投影矩阵Y_2t进行奇异值分解得到：Y_1t＝XA_1t＝Q₁R₁，Y_2t＝X^TA_2t＝Q₂R₂；

4b）用模拟遥感图像序列X与迭代低秩矩阵L_t去逼近迭代稀疏矩阵S_t：

S_t＝P_Ω(X-L_t)，

其中，Ω表示矩阵|X-L_t|的前K个最大非零值子集空间，K≥200，P_Ω表示将非零数据从Ω空间投影到迭代稀疏矩阵S_t上；

4c）重复步骤4a到4b，每重复一次，迭代次数t加1，直至迭代稀疏矩阵S_t和迭代低秩矩阵L_t逼近模拟遥感图像序列X为止，该终止条件为：

||X-L_t-S_t||² _F/||X||² _F＜ε，

其中，|| ||² _F表示2范数，ε为设定的阈值，且0＜ε1；

分别记满足该终止条件的迭代稀疏矩阵S_t和迭代低秩矩阵L_t，并将其作为稀疏矩阵S和低秩矩阵L；

4d）由模拟遥感图像序列X，稀疏矩阵S和低秩矩阵L，得到噪声矩阵G：

G=X-S-L，

其中，稀疏矩阵S＝{s₁,s₂,…s_l,…,s_k+2}，l=1,2,…,(k+2)，s_l∈R^mn×1为稀疏矩阵S的列向量，且S,L,G∈R^mn×(k+2)。

步骤5:用模糊C均值方法对稀疏矩阵S中的最后一列的列向量s_k+2∈R^mn×1进行聚类，得到每个像素对应的类别标签，最终得到变化检测结果图。

本发明的效果可以通过以下仿真实验进一步说明：

1.实验条件

实验环境为：windows XP，SPI，CPU Pentium（R）4，基本频率为2.4GHZ，软件平台为Matlab7.0.4。

仿真选用的第一个数据集为英格兰Feltwell村庄模拟遥感图像数据集,如图2所示，其中图2(a)为英格兰Feltwell村庄变化前的遥感图像，图2(b)为英格兰Feltwell村庄变化后的遥感图像。该图像位于英国Feltwell村庄的一个农田区，模拟变化图像通过模拟地球的天气变化和电磁波的辐射特性等因素影响并人工地嵌入一些变化区域得到，图像大小均为470×335像素，灰度级为256，包括153214个非变化像素和4236个变化像素。检测的标准结果图采用如图3所示的对Feltwell遥感图像数据集进行人工标记的标准变化检测结果图。

第二个数据集为墨西哥遥感图像数据集，如图5所示，其中图5(a)为墨西哥变化前的遥感图像，图5(b)为墨西哥变化后的遥感图像。该组真实遥感数据集由2000年4月和2002年5月的墨西哥郊外的两幅Landsat7ETM+第4波段遥感图像构成。图像大小均为512×512，256灰度级。包含了25599个变化像素和236545个非变化像素。检测的标准结果图采用如图6所示的对墨西哥遥感图像数据集进行人工标记的标准变化检测结果图。

2.实验内容和实验结果

实验一：用本发明方法和两种传统的变化检测方法：ID差值法和KI阈值分割法，对图2进行变化检测。实验结果如图4所示，其中4(a)为现有ID差值法对图2(b)进行变化检测的结果图，4(b)为现有KI阈值分割法对图2(b)进行变化检测的结果图，4(c)为本发明方法对图2(b)进行变化检测的结果图。

实验二：用本发明方法和两种变化检测传统方法：ID差值法和LR对数比值法，对图5进行变化检测，实验结果如图7所示。其中7(a)为现有ID差值法对图5(b)进行变化检测的结果图，7(b)为现有LR对数比值法对图5(b)进行变化检测的结果图，7(c)为本发明方法对图5(b)进行变化检测的结果图。

从图4(a)和图7(a)中可以看出，ID差值法对图像进行变化检测的实验结果中有很多的杂点；从图4(b)和图7(b)中可以看出，LR对数比值法对图像进行变化检测的实验结果产生了很多的细节丢失，边缘很模糊，图4(b)丢失了大部分的变化区域；从图4(c)和图7(c)中可以看出，本发明方法能够在较好保持边缘细节的情况下，减少杂点的个数。

本发明方法和所述两种传统变化检测方法，对图2和图5进行变化检测的结果数据，如表1所示。，表1中列出了四种评价指标：分别为虚警数FA，漏检数MA，错误数OE和正确率PCC，其中，虚警数FA为实际没有发生变化但被当作变化检测出来的像素，漏检数MA为没有检测出来的实际发生了变化的像素，错误数OE=FA+MA，正确率PCC=错误数/图像像素个数。

表1试验结果数据

从表1可以看出，本发明与两种传统变化检测方法相比：首先，从表1的虚检指标中可以看出，本发明方法与两种传统变化检测方法相比有相对较高的虚检数，这与本发明中的低秩稀疏分解好坏有一定的关系；其次，值得注意的是，从表1的漏检指标中可以看出，本发明方法在模拟数据集以及墨西哥遥感数据集上与两种传统变化检测方法相比，均获得较少的漏检数，这在变化检测问题中是比较好的，因此减少漏检数一直是许多算法的目标，而本发明方法在模拟数据集上与差值法相比，漏检数降低了44个像素点，以及在墨西哥数据集上降低了619个像素点。

综上所知，本发明方法与两种传统变化检测方法相比，均有较少的漏检数，虽然虚检率偏高，却在可观的范围内，从而使本发明方法无论是在模拟数据集以及墨西哥遥感数据集中，均获得较少的漏检数，错误数，以及较高的准确率。这说明本发明方法对提高变化检测精度是有所改善的。

Claims (3)

1.一种基于低秩稀疏模型的遥感图像的变化检测方法，包括如下步骤：

(1)输入变化前的遥感图像P₁和变化后的遥感图像P₂，其中，变化前后图像大小均为m×n；

(2)利用变化前的遥感图像P₁和变化后的遥感图像P₂，构造k幅模拟变化后遥感图像X_i,i＝1,2,…,k，且模拟变化后遥感图像X_i的大小为m×n，k为自然数，10≤k≤100；

(3)利用变化前的遥感图像P₁，变化后的遥感图像P₂以及k幅模拟变化后遥感图像X_i,i＝1,2,…,k，构造模拟遥感图像序列X＝[p₁,…x_i,…p₂]，其中，p₁，x_i，p₂分别为将变化前图像P₁，模拟变化后图像X_i以及变化后图像P₂，拉成列向量后的列向量，则p₁，x_i，p₂∈R^mn×1，X∈R^mn×(k+2)，R^mn×(k+2)表示大小为mn×(k+2)的实数空间；

(4)用GODEC低秩稀疏分解算法，对模拟遥感图像序列X进行低秩稀疏分解：

X＝S+L+G，

其中，S为稀疏矩阵，L为低秩矩阵，G为噪声矩阵，S,L,G∈R^mn×(k+2)，R^mn×(k+2)表示大小为mn×(k+2)的实数空间，且S＝{s₁,s₂,…s_l,…,s_k+2}，l＝1,2,…,(k+2)，s_l∈R^mn×1为稀疏矩阵S的列向量；

(5)用模糊C均值方法，对稀疏矩阵S中最后一列的列向量s_k+2∈R^mn×1进行分类，得到每个像素对应的类别标签，最终得到变化检测结果图。

2.根据权利要求1所述的方法，其中步骤(2)所述的利用变化前的遥感图像P₁和变化后的遥感图像P₂，构造k幅模拟变化后遥感图像X_i,i＝1,2,…,k，按如下步骤进行：

2a)用对数比值法，从变化前的遥感图像P₁和变化后的遥感图像P₂中，提取出变化后遥感图像P₂的初始变化区域P₀；

2b)利用变化前的遥感图像P₁，变化后的遥感图像P₂以及初始变化区域P₀，构造k幅模拟变化后遥感图像X_i,i＝1,2,…,k：

$X_i=\left\{\begin{array}{cc}{F}_i& F_i\stackrel{F}{\←}\frac{1}{k}{P}_0\\ U_i& U_i\stackrel{U}{\←}{P}_1-\frac{i}{k}{P}_0\end{array}\right.,$

其中，F_i表示模拟变化后遥感图像X_i中发生变化的区域，U_i表示模拟变化后遥感图像X_i中没有发生变化的区域,式中，表示将初始变化区域P₀划分为k等份，并取其中的i份，表示将变化后图像P₂中的部分变化区域的灰度值，迁移到模拟变化后遥感图像X_i中的对应位置上,式中，表示变化前遥感图像P₁中除与部分变化区域位置相同的区域之外的所有区域，表示将变化前遥感图像P₁中除与部分变化区域位置相同的区域之外的所有区域的灰度值，迁移到模拟变化后遥感图像X_i中的对应位置上。

3.根据权利要求1所述的方法，其中步骤(4)所述的对模拟遥感图像序列X，用GODEC低秩稀疏分解算法进行低秩稀疏分解，按如下步骤进行：

3a)用双边随机投影的方法，对模拟遥感图像序列X进行低秩逼近，得到迭代低秩矩阵L_t为：

$L_t=Q_1{\left[R_1{\left({A_{2t}}^T{Y}_{1t}\right)}^{-1}{R_2}^T\right]}^{1/(2q+1)}{Q_2}^T,$

其中，t≥1为迭代次数，A_2t＝XA_1(t-1)，当t＝1时，A_1(t-1)是初始化为秩为r的高斯矩阵，2≤r≤30，当t≥2时，A_1(t-1)＝X^TA_2(t-1)，这里，q为收敛系数，且Q₁,R₁和Q₂,R₂分别为对左投影矩阵Y_1t和右投影矩阵Y_2t进行奇异值分解得到：

Y_1t＝XA_1t＝Q₁R₁，Y_2t＝X^TA_2t＝Q₂R₂；

3b)用模拟遥感图像序列X与迭代低秩矩阵L_t去逼近迭代稀疏矩阵S_t：

S_t＝P_Ω(X-L_t)，

其中，Ω表示矩阵|X-L_t|的前K个最大非零值子集空间，K≥200，P_Ω表示将非零数据从Ω空间投影到迭代稀疏矩阵S_t上；

3c)重复步骤3a到3b，每重复一次，迭代次数t加1，直至迭代稀疏矩阵S_t和迭代低秩矩阵L_t逼近模拟遥感图像序列X为止，该终止条件为：

||X-L_t-S_t||² _F/||X||² _F<ε，

其中，||||² _F表示2范数，ε为设定的阈值；

分别记满足该终止条件的迭代稀疏矩阵S_t和迭代低秩矩阵L_t，并将其作为稀疏矩阵S和低秩矩阵L；

3d)由模拟遥感图像序列X，稀疏矩阵S和低秩矩阵L，得到噪声矩阵G：

G＝X-S-L。