CN102063720A - 基于Treelets的遥感图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于Treelets的遥感图像变化检测方法，它属于遥感图像分析与处理领域，主要解决现有遥感图像变化检测方法存在很多伪变化信息的问题。其实现步骤是：(1)对输入的时相1遥感图像进行Treelets滤波；(2)对输入的时相2遥感图像进行Treelets交叉滤波；(3)对滤波后的两时相图像对应像素点灰度值进行差值计算，得到一幅差异图；(4)对差异图像再进行Treelets滤波得到新的差异图像；(5)对新差异图采用K-means聚类分为变化类和非变化类，得到最后的变化检测结果图。本发明能有效减小辐射校正不理想和光照不均对检测效果的影响，提高变化检测的精度，可用于对灾情监测、土地利用及农业调查。

Description

基于Treelets的遥感图像变化检测方法

技术领域

本发明属于图像处理技术领域，具体地说是一种基于Treelets的遥感图像变化检测的方法，适用于遥感图像分析和处理。

背景技术

遥感图像变化检测的研究始于20世纪70年代，随着图像处理技术的不断发展，遥感图像变化检测的研究逐渐成为了一个热点，并广泛应用于国民经济和国防建设的诸多领域，如森林资源的动态监测，土地覆盖、利用的变化监测，农业资源调查，城市规划布局，环境监测分析，自然灾害评估，地理数据更新以及军事侦察中道路、桥梁、机场等战略目标的动态监视等。

早期的遥感图像变化检测方法由于受到当时技术条件的限制，采用的方法较为简单。其中很多是通过人工目测来进行解释的，但人工目测受到个人主观因素的影响较大，容易产生错误，且效率低下，因此在实际应用方面受到很大的限制。一般而言，遥感图像变化检测的流程主要包括：遥感图像的预处理、变化区域检测、后处理和变化检测性能评价几个步骤。随着图像处理技术的不断发展，遥感图像变化检测的研究逐渐成为了一个热点，近几年经过国内外遥感图像变化检测研究的蓬勃发展，涌现出了很多新的有效的方法。但宏观上讲主要分为两条研究路线：第一条是先比较后分类的方法，即首先构造两幅图像的差异图像，然后再对差异图像进行分类确定变化区域和非变化区域；第二条是先分类后比较的方法，即首先对两幅图像进行分类，然后再对分类后的两幅图像进行比较确定变化区域和非变化区域。

针对基于先比较后分类的遥感图像变化检测方法，它主要涉及差异图像的构建和分类两个方面。其中对于差异图像的构建方法有很多，其中主要有：差值法(灰度差值法、纹理特征差值法、图像回归法)、比值法、相关系数法、变化向量法和主分量分析法等。基于差异图像分析的方法由于简单直接，不同于其他复杂的方法，不会由于方法本身而对检测前的变化信息加以改变，忠于原始数据，保证了变化检测结果的可靠性。然而不同时相的遥感图像之间由于不同季节和情况下的光照、辐射等因素造成不同时相间的图像灰度值存在整体或者部分的差异，因此简单地对灰度值运算得到的差异图进行阈值分割，所得到的变化检测结果中存在很多伪变化信息，使得变化检测的精度较低。

发明内容

本发明的目的在于克服上述已有技术的不足，提出一种基于Treelets的遥感图像变化检测方法，以减少伪变化信息，提高遥感图像变化检测精度。

为实现上述目的，本发明的检测方法包括如下步骤：

(1)输入时相1遥感图像，并对该时相1遥感图像进行Treelets滤波；

(2)输入时相2遥感图像，并对时相2图像每个像素点取5×5滑动窗拉成列向量投影至步骤1中得到的尺度向量Φ中，用该结果取代时相2的像素点灰度值，得到交叉滤波后的时相2图像；

(3)将步骤1和步骤2得到的滤波后的两个时相图像对应点进行差值计算，得到一幅差异图；

(4)对步骤3得到的差异图像按步骤1再进行一次Treelets滤波得到一幅新的差异图像；

(5)对步骤4得到的新差异图采用K-means聚类分为变化类和非变化类，得到最后的变化检测结果图。

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

a、本发明采用了Treelets滤波和Treelets交叉滤波，避免了两时相图像因噪声、光照、辐射等因素对变化检测结果的影响。

b、本发明对于变化区域的一致性保持较好。

c、本发明减少了伪变化信息，提高了变化检测的精确性。

附图说明

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

图2是本发明的第一组实验数据；

图3是本发明的第二组实验数据；

图4是本发明的第三组实验数据；

图5是本发明的第四组实验数据；

图6是本发明的第一组实验数据的差异图像和变化检测结果图；

图7是本发明的第二组实验数据的差异图像和变化检测结果图；

图8是本发明的第三组实验数据的差异图和变化检测结果图像；

图9是本发明的第四组实验数据的差异图和变化检测结果图像。

具体实施方式

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

步骤1，输入时相1遥感图像，并对该时相1遥感图像进行Treelets滤波。

la)对时相1图像取5×5像素的滑动窗，计算该图像的初始协方差矩阵

${\widehat{\mathrm{\Σ}}}^{\left(0\right)}=\left[\begin{array}{cccc}{\mathrm{\σ}}_{11}& {\mathrm{\σ}}_{12}& L& {\mathrm{\σ}}_{1v}\\ {\mathrm{\σ}}_{21}& {\mathrm{\σ}}_{22}& L& {\mathrm{\σ}}_{2v}\\ M& M& L& M\\ {\mathrm{\σ}}_{u1}& {\mathrm{\σ}}_{u2}& L& {\mathrm{\σ}}_{\mathrm{uv}}\end{array}\right]$

其中，X为5×5像素的滑动窗向量，u，v分别为初始协方差矩阵

的位置索引，u＝1，2，L，25，v＝1，2，L，25；

1b)由初始协方差矩阵

计算图像的初始相似度矩阵

${\hat{M}}^{\left(0\right)}=\left|\frac{{{\widehat{\mathrm{\Σ}}}_{\mathrm{st}}}^{\left(0\right)}}{\sqrt{{{\widehat{\mathrm{\Σ}}}_{\mathrm{ss}}}^{\left(0\right)}{{\widehat{\mathrm{\Σ}}}_{\mathrm{tt}}}^{\left(0\right)}}}\right|+\left|{{\widehat{\mathrm{\Σ}}}_{\mathrm{st}}}^{\left(0\right)}\right|$

其中，s，t分别为初始协方差矩阵

的位置索引；

1c)初始化基矩阵B₀为25×25的单位矩阵：

$B_0=\left[\begin{array}{cccc}1& 0& L& 0\\ 0& 1& L& 0\\ M& M& O& M\\ 0& 0& L& 1\end{array}\right]$

初始化Treelets层级：l＝{0，1，2，K 24}以及和变量的下标集：δ＝{1，2，L，25}；1d)由l层的图像相似度矩阵

找出最相似的两个变量α和β：

$(\mathrm{\α},\mathrm{\β})=\arg \underset{i,j\∈\mathrm{\δ}}{\max }{{\hat{M}}_{\mathrm{ij}}}^{\left(l\right)}$

其中，α和β分别代表图像的协方差矩阵中的行变量和列变量，i和j是图像相似度矩阵

的位置索引，且必须属于和变量下标集δ，并且i＜j；

1e)对图像的协方差矩阵

进行局部主成分分析变换，得到第一主成分的和变量s_l和第二主成分的差变量d_l，且使得图像协方差矩阵

中α行β列的值和β行α列的值都为零，即

得到旋转角度为θ_l，并由下式得到雅克比旋转矩阵J：

$J=\left[\begin{array}{ccccccc}1& L& 0& L& 0& L& 0\\ M& O& M& & M& & M\\ 0& L& c& L& -s& L& 0\\ M& & M& O& M& & M\\ 0& L& s& L& c& L& 0\\ M& & M& & M& O& M\\ 0& L& 0& L& 0& L& 1\end{array}\right]$

其中，c和s为雅克比旋转矩阵J中的两个不同变量，c＝cos(θ_l)，s＝sin(θ_l)，|θ_l|≤π/4；

1f)根据雅克比旋转矩阵J计算当前层级的基矩阵：B_l＝B_l-1J，尺度函数φ_l和细节函数ψ_l分别为基矩阵B_l的第α和β列，当前层级的尺度基Φ是尺度函数φ_l和上一层的尺度向量集{φ_l-1，jt}_jt≠α，β的合集，计算相似度矩阵：

和协方差矩阵：

1g)将差变量的下标β从和变量的下标集δ中去除，即δ＝δ\{β}；

1h)重复步骤1d)至步骤1g)直至到l＝24层，得到基矩阵为：

B＝[Φ ψ₁ L ψ_l-1]^T

其中，Φ和ψ分别是基矩阵B的尺度向量和细节向量；

1i)对时相1图像每个像素点取5×5滑动窗拉成列向量，然后投影至尺度向量Φ中，用得到的结果取代时相1的像素点灰度值，得到滤波后的时相1图像；

步骤2，输入时相2遥感图像，并对时相2图像每个像素点取5×5滑动窗拉成列向量投影至步骤1h)中得到的尺度向量Φ中，用得到的结果取代时相2的像素点灰度值，得到交叉滤波后的时相2图像；

步骤3，将步骤1和步骤2得到的滤波后的两个时相图像对应点进行差值计算，得到一幅差异图；

步骤4，对步骤3得到的差异图像按步骤1a)至步骤1i)再进行一次Treelets滤波得到一幅新的差异图像；

步骤5，对步骤4得到的新差异图采用k-means聚类算法将差异图像分为两类{ω_c，ω_n}，其中ω_c代表变化类，ω_n代表非变化类，然后对差异图像的像素依据分类的结果进行二值化，得到最后的变化检测结果图。

本发明的效果可以通过以下内容进行说明：

1.实验数据

(1)模拟数据集

实验所用的一组模拟数据集原始图像及参考变化图如图2所示。其中模拟数据集的原始图像是ATM(Airborne Thematic Mapper)3波段，位于英国Feltwell村庄的一个农田区的图像，如图2(a)所示，其模拟变化图像是通过模拟地球的天气变化和电磁波的辐射特性等因素影响并人工地嵌入一些变化区域得到的，如图2(b)所示。两幅图像大小均为470×335像素，灰度级为256，配准误差为1.5个像素左右。参考变化图如图2(c)所示，其中包含变化的像素数为4236，非变化像素数为153214，白色像素区域表示变化的区域。

(2)真实遥感图像数据集

真实遥感图像数据集共三组。第一组真实遥感图像数据集原始图像及参考变化图如图3所示。该组真实遥感数据集是分别于1994年8月和1994年9月在意大利Elba岛西部地区的两时相Landsat-5TM第4波段多光谱图像组成，分别如图3(a)和图3(b)所示，两幅图像大小均为326×414，灰度级为256，它们之间发生的变化是由于森林火灾破坏了大量植被所致，如图3(b)中左上角较暗的区域，其参考变化图如图3(c)所示，图中包含2415个变化像素和132549个非变化像素，白色像素区域表示变化的区域。第二组真实遥感图像数据集原始图像及其参考变化图如图4所示。该组真实遥感图像数据集是由2000年4月和2002年5月的墨西哥郊外的两幅Landsat7ETM第4波段光谱图像组成，分别如图4(a)和(b)所示。两幅图像大小均为512×512像素，灰度级为256，配准误差为1.5个像素左右，它们之间的变化区域主要由于火灾破坏了大量的植被所致，如图4(b)中较暗的区域。其参考变化图如图4(c)所示错误！未找到引用源。，共包含25599个变化像素和236545个非变化像素，白色像素区域表示变化的区域。第三组真实遥感图像数据集原始图像及参考变化图如图5所示。该组真实遥感数据集由1995年9月和1996年7月在意大利撒丁岛Mulargia湖泊区域的Landsat5TM第5波段光谱图像组成，分别如图5(a)和(b)所示。两幅图像大小均为300×412像素，灰度级为256，它们之间的变化为湖水水位上涨引起的，发生变化区域的参考图如图5(c)所示，其中包括115974个非变化像素和7626个变化像素，图中白色像素区表示变化的区域。

2.实验评价指标

评价方法是检测变化检测方法好坏的依据，本发明中将采用变化检测方法研究中常用的三个评价指标：虚警数、漏检数和总错误数。由于变化检测类似于图像分类，所不同的是变化检测的结果只有两类，即变化类和非变化类。因此对于变化检测效果的评价类似

于对遥感图像分类精度的评价，即选取一定数量的样本构造一个2×2的变换检测混淆矩阵，如表1所示。

表1变化检测混淆矩阵

表1中的C₁为变化检测方法结果中的非变化像素数，C₀为变化像素数，变化参考图中非变化像素数为R₁，变化像素数为R₀，则真正变化的像素数D₀＝{R₀∩C₀}，真正非变化的像素数D₁＝{R₁∩C₁}；虚警数C_fa为变化参考图中像素为非变化类而变化检测方法所得结果中将该像素为变化类的像素的个数，即

C_fa＝{C|C∈C₀∩C∈R₁}

漏检像素数C_mil为变化参考图中像素为变化类而变化检测方法所得结果中将该像素为非变化类的像素的个数，即

C_mil＝{C|C∈R₀∩C∈C₁}

错误像素数C_e为漏检像素数与虚警数之和，即C_e＝C_fa+C_mil。

3.实验内容与结果

(1)差异图像的实验结果

为了验证基于基于Treelets构建的差异图像的实验效果，我们与基于像素差值的差异图像，即直接对两时相图差运算得到的差异图像进行了实验对比。其中：

模拟遥感图像数据集的差值法差异图和本发明差异图分别如图6(a)和6(b)所示；第一组真实遥感图像数据集的差值法差异图和本发明差异图分别如图7(a)和7(b)所示；第二组真实遥感图像数据集的差值法差异图和本发明差异图分别如8(a)和8(b)所示；第三组真实遥感图像数据集的差值法差异图和本发明差异图分别如图9(a)和9(b)所示。

从四组实验数据集采用不同方法构建的差异图像可以看出，本发明差异图像能够很好的将两时相遥感图像差异图中的变化区域和非变化区域的像素灰度值进行拉伸，更加凸显了变化的区域，而差值法差异图像中非变化区域中存在较多像素的灰度值与变化区域的灰度值相近，这将会导致后续的变化检测结果中含有较多伪变化信息。

(2)变化检测的实验结果与分析

采用k-means聚类方法分别计算本发明方法和差值法得到的差异图像进行分割，得到变化检测结果图。

模拟数据集采用不同方法得到的变化检测结果如图6所示，其中差值法方法的变化检测结果如图6(c)所示，本发明方法的变化检测结果如图6(d)所示。由于该模拟数据集变化信息很明显，所以采用像素差值法就可以得到较好的变化检测图，本发明方法也得到了同样的效果。第一组真实遥感图像数据集采用不同方法得到的变化检测结果如图7所示，差值法方法的变化检测结果如图7(c)所示，本发明方法的变化检测结果如7(d)所示。

从图7(c)中可以看出非变化区域中存在大量伪变化的区域，这是由于图3(a)、(b)两时相图像中辐射定标校正效果不理想，像素差值法容易受噪声、光照等因素的影响，导致变化检测结果中存在很多的伪变化信息，变化检测的效果较差。

第二组真实遥感图像数据集的变化检测结果如图8所示，其中，差值法方法的变化检测结果如图8(c)所示，本发明方法的变化检测结果如图8(d)所示。从如图8(c)可以看出，差值法方法得到的变化检测结果中含有较多的虚警点，这些孤立的像素点很多是由于噪声引起的，从而导致伪变化信息增加，而本发明方法得到的变化检测结果图中，含有较少的虚警点，更接近参考变化图。

第三组真实遥感图像数据集的变化检测结果如图9所示，差值法方法的变化检测结果如图9(c)所示，本发明方法的变化检测结果如图9(c)所示。从图9(c)可以看出，差值法方法得到的变化检测结果中含有很多的杂点，含有较多的伪变化信息，从图9(d)可以看出，本发明方法得到的变化检测结果图中，含有很少的杂点，伪变化信息较少，更好的反应了变化的区域。

从以上给出了四组实验数据集的变化检测结果，可以从主观上看出本发明方法的有效性。下面将从虚警数、漏检像素数和总错误像素数三个方面客观的评价本发明方法。

四组实验数据集的结果如表2所示。

表2四组实验数据的变化检测结果评价指标

从表2中可以看出，对模拟数据集本发明方法得到的变化检测结果与差值法方法得到的变化检测结果中总错误数增加了24个像素点，这是由于变化区域的边缘不精确所致，但并不存在其他多余的错误区域。

第一组真实遥感图像数据集的差值法方法得到的变化检测结果中虚警数为35031，总错误数为35039，而本发明法得到的变化检测结果中虚警数为41，比差值法方法减少了34990个像素点，总错误数为281，比差值法方法减少了34758个像素点。

第二组真实遥感图像数据集本发明法得到的变化检测结果比差值法方法虚警数减少了461个像素点，总错误数减少了1142个像素点。

第三组真实遥感图像数据集本发明方法得到的变化检测结果中比差值法方法虚警数减少了2908个像素点，总错误数减少了2604个像素点。

从以上四组数据集的实验结果可以客观的看到，除模拟实验数据集的本发明方法变化检测效果略低于差值法方法变化检测效果外，其它四组真实遥感图像实验数据集的本发明方法变化检测效果均好于差值法方法变化检测效果，尤其当两时相数据集辐射校正不理想，光照不均时更能体现本发明的优势。

Claims (1)

1.一种基于Treelets的遥感图像变化检测方法，包括如下步骤：

(1)输入时相1遥感图像，并对该时相1遥感图像进行如下Treelets滤波：

1a)对时相1图像取5×5像素的滑动窗，计算该图像的初始协方差矩阵

1b)由协方差矩阵

计算图像的初始相似度矩阵

1c)初始化基矩阵B₀为25×25的单位矩阵：

$B_0=\left[\begin{array}{cccc}1& 0& L& 0\\ 0& 1& L& 0\\ M& M& O& M\\ 0& 0& L& 1\end{array}\right]$

初始化Tree1ets层级：l＝{0，1，2，K 24}以及和变量的下标集：δ＝{1，2，L，25}；1d)由l层的图像相似度矩阵

找出最相似的两个变量α和β：

$(\mathrm{\α},\mathrm{\β})=\arg \underset{i,j\∈\mathrm{\δ}}{\max }{{\hat{M}}_{\mathrm{ij}}}^{\left(l\right)}$

其中，α和β分别代表图像的协方差矩阵

中的行变量和列变量，i和j是图像相似度矩阵

的位置索引，且必须属于和变量下标集δ，并且i＜j；

1e)对图像的协方差矩阵进行局部主成分分析变换，得到第一主成分的和变量s_l和第二主成分的差变量d_l，且使得图像协方差矩阵

中α行β列的值和β行α列的值都为零，即得到旋转角度为θ_l，并由下式得到雅克比旋转矩阵J：

$J=\left[\begin{array}{ccccccc}1& L& 0& L& 0& L& 0\\ M& O& M& & M& & M\\ 0& L& c& L& -s& L& 0\\ M& & M& O& M& & M\\ 0& L& s& L& c& L& 0\\ M& & M& & M& O& M\\ 0& L& 0& L& 0& L& 1\end{array}\right]$

其中，c和s为雅克比旋转矩阵J中的两个不同变量，c＝cos(θ_l)，s＝sin(θ_l)，|θ_l|≤π/4；

1f)根据雅克比旋转矩阵J计算当前层级的基矩阵：B_l＝B_l-1J，以及相似度矩阵：

和协方差矩阵：

1g)将差变量的下标β从和变量的下标集δ中去除，即δ＝δ\{β}；

1h)重复步骤1d)至步骤1g)直至到l＝24层，得到基矩阵为：

B＝[Φ ψ₁ L ψ_l-1]^T

其中，Φ和ψ分别是基矩阵B的尺度向量和细节向量；

1i)对时相1图像每个像素点取5×5滑动窗拉成列向量，然后投影至尺度向量Φ中，用得到的结果取代时相1的像素点灰度值，得到滤波后的时相1图像；

(2)输入时相2遥感图像，并对时相2图像每个像素点取5×5滑动窗拉成列向量投影至步骤1h)中得到的尺度向量Φ中，用该结果取代时相2的像素点灰度值，得到交叉滤波后的时相2图像；

(3)将步骤1和步骤2得到的滤波后的两个时相图像对应点进行差值计算，得到一幅差异图；

(4)对步骤3得到的差异图像按步骤1a)至步骤1i)再进行一次Treelets滤波得到一幅新的差异图像；

(5)对步骤4得到的新差异图采用K-means聚类分为变化类和非变化类，得到最后的变化检测结果图。

Images