CN102831598B - 多分辨率NMF和Treelet融合的遥感图像变化检测方法 - Google Patents

Abstract

本发明公开一种多分辨率NMF和Treelet融合的遥感图像变化检测方法，针对从单一分辨率考虑图像细节和平滑区域时常无法较好权衡的问题，使得图像变化检测既保留图像的细节信息又保留平滑区域信息。实现过程是：输入两时相图像，直接差值法构造差异图像和进行中值滤波；之后运用NMF算法提取不同分辨率图像；对滤波后差异图像和不同分辨率图像分别取阈值；用Treelet算法融合上述阈值后图像；用区域生长法对融合后图像分割，得到最终变化检测结果。本发明解决了图像的邻域结构易受到孤立噪声点影响的问题，能够保持图像的细节信息与平滑区域信息，且能够去除孤立噪声，提高了变化检测精度，可用于灾情监测、土地利用、农业调查等领域。

Description

多分辨率NMF和Treelet融合的遥感图像变化检测方法

技术领域

本发明属于图像处理技术领域，涉及遥感图像变化检测，主要涉及多分辨率非负矩阵分解(Non-negative Matrix Factorization，NMF)和Treelet融合的遥感图像变化检测，具体是一种多分辨率NMF和Treelet融合的遥感图像变化检测方法。用于对遥感图像变化区域的检测。

背景技术

遥感图像变化检测通过分析同一地区不同时相的遥感图像间光谱特征差异或空间结构特征的差异，从而得到地物的动态发展变化信息。遥感图像的变化检测已经得到了广泛的应用，如土地、森林、草场等的资源监测，农作物估产及病虫害监测，地图数据的校正更新，土地利用和覆盖变化的监测，海洋、湖泊水资源变化的监测，海岸线变化、湿地变化、城区变化等的监测以及地震、海啸、火灾等突发事件的评估；另外，军事目标动态侦查、战场动态情报获取、军事部署情况监测等。随着变化检测应用的发展，变化检测方法也成为重要的研究内容。

不要求地面真实数据和额外分类信息的无监督的变化检测方法是目前变化检测方法研究的热点。一般的无监督变化检测方法是对输入图像所构造的差异图像进行分析，以判断出二分类的变化和非变化区域。然而传统的基于马尔科夫随机场(MarkovRandom Fields，MRF)的方法在考虑图像的邻域结构或纹理信息时易受到孤立噪声点的影响，从单一分辨率角度考虑图像的细节和平滑区域的权衡的方法常无法较好的既保留图像的细节信息又保留平滑区域信息，导致最终的变化检测结果不够准确。

为了提高变化检测的精度，研究人员提出了许多方法。如：Celik等在文章“Unsupervised Change Detection in Satellite Images Using Principal ComponentAnalysis and K-Means Clustering，IEEE Geoscience and Remote Sensing Letters，2009，6(4)：772-776”中提出了基于主成分分析(Principle Component Analysis，PCA)和K-means聚类的变换检测方法，该方法提取数据特征时结合了像素的邻域信息，操作简单方便，但是PCA方法对于非线性相关的数据分类结果却不是很理想，且采用PCA方法提取特征数据时，对数据的邻域信息考虑有限，只在固定的邻域块内进行结构纹理分析。没有从多分辨率数据融合的角度对数据进行分析，没有考虑好图像的细节信息和平滑区域的信息和噪声的影响给变化检测带来错误的结果，从而影响了变化检测的准确性。

黄世奇等在文章“基于小波变换的多时相SAR图像变化检测技术，测绘学报，2010，39(2)：180-186”中提出了一种基于图像多分辨率分解的可靠尺度融合的遥感图像变化检测方法。该方法对差异图像采用2维离散平稳小波分解，通过判定局部变化系数和全局变化系数的关系确定可靠尺度，并对可靠尺度进行加入权重的特征级融合，对融合后的图像采用EM双阈值得到变化检测结果。该方法的优点在于选择可靠尺度达到了去除噪声和保留细节的均衡，同时能够获得发生变化和未发生变化的像素区域，还可以区分发生变化的类型，如变化区域增强类和变化减弱类。但是由于平稳小波对图像具有平滑作用，分解的层数越高，平滑效应越大，造成一定程度上扩散了变化区域的面积，使得检测结果中虚警较高。同时特征级融合时只采用图像的低频信息，造成变化区域的边缘保持不理想。此外，对特征融合后的图像采用EM双阈值进行分割，造成伪变化信息过多，虚警较高。

Li等在文章“Multitemporal Image Change Detection Using a Detail-EnhancingApproach With Nonsubsampled Contourlet Transform，IEEE Geoscience and RemoteSensing Letters，2012，9(5)：836-840”中提出了一种基于非下采样contourlet变换的变化检测方法。该方法对差异图采用非下采样contourlet变换分解得到低频子带和高频的方向子带，对高频方向子带采用尺度内和尺度间的方法进行融合，以此增强图像的细节信息，同时选出较优的低频子带，将该低频子带加上一定权重值的融合后的方向信息得到细节信息增强的差异图像，接着对差异图像采用PCA构造特征矢量空间，对该特征矢量空间采用基于PCA指导的K-means进行聚类得到变化检测结果。该方法采用尺度内和尺度间融合的方法增强了图像的细节信息，对于变化区域较为明显的图像能够检测出较好的结果，但对于变化区域不明显的图像，变化检测的精度大大降低。

综上所述，上述现有的方法均不能既处理好非线性相关的数据，又同时保持图像的细节信息和平滑区域的信息，且易受噪声影响，造成漏检或者虚警过多，降低变化检测性能。

本发明的发明人就本主题对国内外的专利文献和公开发表的期刊论文搜索，尚未发现与本发明完全相同的文献及报道。

发明内容

本发明的目的在于克服上述已有技术的不足，提出了一种多分辨率NMF和Treelet融合的遥感图像变化检测方法。本方法利用了NMF是非线性维数约减方法，能够处理非线性相关的数据，且本方法中多分辨分析既能保留图像的细节信息又能保留图像平滑区域信息的特点，受噪声影响较小，可提高后续的变化检测的性能。

本发明是一种多分辨率NMF和Treelet融合的遥感图像变化检测方法，其特征在于：包括有如下步骤：

(1)取同一地区在不同时间获取的已配准的两幅大小均为P×Q的遥感图像，将该两图像对应空间位置像素的灰度值相减取绝对值，得到一幅差异图像，对此差异图像进行大小为m×m像素的中值滤波，其中，m的取值范围为3、5、7、9，在该步骤中包含了两部分内容，首先是对输入的两幅图像构造差异图，然后对差异图进行中值滤波。

(2)对滤波后的差异图像运用NMF算法提取五幅不同分辨率图像F_r(r＝1，2....，5)；

2a)将滤波后的差异图像X_D分成大小均为h×h且不重叠的正方形图像块E，将每个小块E转成h²×1的列向量C_h，所有块的列向量合并构成矩阵V_h，其中，图像块尺寸h取2～～10中的所有偶数值；

2b)对矩阵V_h采用NMF进行分解，得到NMF基矩阵W_h和系数矩阵H_h；

2c)对滤波后的差异图像X_D边界进行扩展得到边界扩展图像D，即将差异图像X_D的首列向左扩展w列，末列向右扩展w列，对列扩展完成后的图像再进行首行向上扩展w行，末行向下扩展w行，即可得到边界扩展后的差异图像D。其中， 为向下取整符号；

2d)逐个选择图像D中的非边界扩展像素——即D中对应的差异图像X_D的每个像素点，作为中心像素点，以此为中心，也取大小为h×h的邻域块，将每个h×h邻域块转化成h²×1的列向量，这些列向量合并构成了矩阵V_vh，将矩阵V_vh在NMF的相应的基W_h上进行投影，得到特征数据集F_dh，将特征数据集F_dh转化成P×Q的图像大小，可得特征图像F_h；

2e)重复步骤(2a)至步骤(2d)直到图像块尺寸h依次取遍2，4，6，8，10，可得到五幅不同分辨率的图像F_r(r＝1，2....，5)，其中r是图像的标号。

(3)估计滤波后的差异图像X_D的噪声标准差对X_D取阈值得到图像Y，估计五幅不同分辨率图像F_r(r＝1，2，...，5)的噪声标准差对F_r(r＝1，2,...，5)取阈值得到图像Y_r(r＝1，2，...，5)，其中K为常数，K的取值为2，在该步骤中包含了两部分内容，一部分内容是对滤波后的差异图像取阈值，另一部分内容是对五幅不同分辨率图像取阈值，该步骤的目的是为了抑制图像的背景噪声的影响。

(4)用Treelet算法对阈值后的差异图像Y和阈值后的五幅不同分辨率图像Y_r进行融合，得到一幅融合后的图像A，该步骤包含三部分内容，首先将阈值后的差异图像Y和阈值后的五幅不同分辨率图像Y_r构成初始样本，然后对初始样本进行逐层聚类，直至最高层，得到最终的狄拉克基矩阵B，最后将初始样本与在狄拉克基矩阵B转置的方向进行投影，将投影得到的数据集转化为P×Q的图像大小，得到融合后的图像A。

(5)采用区域生长算法对融合后的图像A进行分割，得到最终的变化检测结果。

本发明将直接差值法构造的差异图像运用NMF算法在不同分辨率下提取特征图像，将差异图像和不同分辨率的特征图像取阈值，以此抑制背景噪声影响，然后运用Treelet算法将阈值后的差异图像和不同分辨率的特征图像进行融合，采用区域生长算法对融合后的图像进行分割，得到的二值图像为最终结果，即变化检测结果。

本发明采用的NMF方法能够较好地处理非线性相关的数据，同时沿用多分辨分析的方法既能够保留图像的细节信息又能保留图像的平滑区域的信息的特点，且采用的Treelet融合方法能够既保留图像的强变化区域的信息又能保留图像的弱变化区域的信息，综合利用了这些方法的优势，扬长避短，提高了变化检测的准确性。

本发明的实现还在于：步骤2b)中对矩阵V_h采用NMF进行分解的过程是：

2b1)随机初始化W_h和H_h，设置最大迭代次数e和停止精度ε，最大迭代次数e的取值范围为50～1000，停止精度ε的取值范围为10^-4～10^-6；

2b2)在W_h和H_h均为非负矩阵的约束下，极小化目标函数||V_h-W_hH_h||²，对矩阵H_h进行更新，

$H_h=H_h\frac{W_h^T{V}_h}{W_h^T{W}_h{H}_h}---\left(1\right)$

对矩阵W_h进行更新，

$W_h=W_h\frac{V_h{H}_h^T}{W_h{H}_h{H}_h^T}---\left(2\right)$

直到满足或者达到最大迭代次数e为止，可得到NMF基矩阵W_h和系数矩阵H_h，其中，b为当前迭代次数，迭代次数b的取值范围为50～1000。

本发明采用了NMF方法的优势在于，能够较好地逼近原矩阵，不存在负分量，能方便、合理的进行数据解释，同时具有稀疏性，能够抑制外界变化带来的不利影响，能够突出变化区域，有利于后续的变化检测工作。

本发明的实现还在于：步骤(4)中用Treelet算法对阈值后的差异图像Y和阈值后的五幅不同分辨率图像Y_r进行融合的过程是：

4a)将阈值后图像Y和Y_r合并表示为Y_t(f＝1，2,...，6)，t为图像标号，将Y_t转为大小均为P×Q的列向量Y₁、Y₂、Y₃、Y₄、Y₅、Y₆，所有列向量构成初始样本矩阵X＝[Y₁，Y₂，Y₃，Y₄，Y₅，Y₆]，阈值后图像Y为一幅图像，Y_r为阈值后的五幅不同分辨率图像，需要进行融合的图像为六幅。

4b)初始化Treelet变换的逐层聚类层数ι＝0，1，...，L-1，L为矩阵X的列向量的个数，即L＝6，在第0层，每个变量采用初始样本X的列向量表示，初始化和变量的下标集δ＝{1，2，...，L}，初始化狄拉克基矩阵B⁰为L×L的单位阵，计算矩阵X的协方差矩阵和相关系数矩阵计算公式如下：

${\widehat{\mathrm{\Σ}}}_{\mathrm{ij}}=E\left[(X_u-E{X}_v)(X_u-E{X}_v)\right]---\left(3\right)$

${\hat{M}}_{\mathrm{ij}}=\frac{{\widehat{\mathrm{\Σ}}}_{\mathrm{ij}}}{\sqrt{{\widehat{\mathrm{\Σ}}}_{\mathrm{ii}}{\widehat{\mathrm{\Σ}}}_{\mathrm{jj}}}}---\left(4\right)$

其中，表示初始协方差矩阵第i行第j列的值，i＝1，2,...，P×Q，j＝1，2，...，P×Q，X_u和X_v表示样本矩阵X中两个不同的列向量，u＝1，2,...，6，v＝1，2，...，6；表示初始化相关系数矩阵第i行第j列的值；

4c)当l≠0时，由公式(4)计算相关系数矩阵找出中最大的两个值，将最大值和次大值的对应位置序号分别记为α和β：

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

这里i＜j，分别表示相关系数矩阵中任意值的行和列，且只在变量下标集合δ内进行；

4d)对图像协方差矩阵进行局部主成分分析变换，得到第一主成分的和变量s₁和第二主成分的差变量d₁，且使得图像协方差矩阵中α行β列的值和β行α列的值都为零，即得到旋转角度θ₁，并由下式得到雅克比旋转矩阵J：

其中，|θ₁|≤π/4；

4e)根据雅克比旋转矩阵J计算当前聚类层级l的狄拉克基矩阵B_l＝B_l-1J，狄拉克基矩阵B_l的第α和β列分别是尺度函数Φ_l和细节函数Ψ_l，前l层级的尺度向量{Φ_l}是尺度函数Φ_l和上一层的尺度向量{Φ_l-1}的合集，同时更新相关系数矩阵 ${\hat{M}}^{\left(l\right)}=J^T{\hat{M}}^{(l-1)}J$ 和协方差矩阵 ${\widehat{\mathrm{\Σ}}}^{\left(l\right)}=J^T{\widehat{\mathrm{\Σ}}}^{(l-1)}J;$

4f)将差变量的序号β的下标从和变量下标集δ中去除，即δ＝δ\{β}；

4g)重复步骤(4c)至步骤(4f)直至分解l＝L-1层，得到最终Treelet分解的基矩阵：

B＝[Φ_L-1，Ψ₁，...，Ψ_L-1]^T                                      (7)

其中，Φ_L-1∈{Φ_L-1}，Ψ₁为第一层Treelet变换得到的细节函数，Ψ_L-1是最高层Treelet变换得到的细节函数；

4h)将矩阵X与在狄拉克基矩阵B转置的方向进行投影，即R＝X×B^T，得到融合后的数据集R，将数据集R转化为P×Q的图像大小，得到融合后的图像A。

本发明采用Treelet融合阈值后的差异图像和阈值后的多分辨率图像的优势在于，能够较好地保留图像的强变化区域的信息，同时也能够保留图像的弱变化区域的信息，提高了变化检测的准确性。

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

1、NMF是非线性维数约减方法，能够处理非线性相关的数据，采用NMF提取数据特征，由于不存在负分量，能方便、合理的进行数据解释，解决了PCA线性维数约减方法不能很好地处理非线性相关的数据，且存在负分量的问题。同时由于NMF具有一定的稀疏性，在一定程度上抑制了外界变化(如光照变化等)给特征提取带来的不利影响，从而提高了变化检测的性能。

2、本方法中的多分辨分析考虑了像素点的邻域信息，同时既保留了图像的细节信息又保留了图像的平滑区域信息。克服了只考虑固定大小的邻域块容易丢失细节信息或者平滑区域信息，并且具有抗噪性，一定程度上提高了变化检测的准确性和可靠性。

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

附图说明

图1是本发明遥感图像变化检测流程框图；

图2是用于实验的第一组真实遥感图像数据集原始图像及变化参考图像；

图3是用于实验的第二组真实遥感图像数据集原始图像及变化参考图像；

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

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

具体实施方式

结合附图对本发明详细说明：

随着遥感图像获取技术和手段的日益先进以及遥感图像数据的海量积累，遥感图像变化检测技术已成为遥感图像处理与分析领域的研究热点，被广泛应用于土地、森林资源的监测，农作物估产及病虫害监测，土地利用和覆盖变化的监测，城区变化等的监测以及地震、海啸、火灾等突发事件的评估等。

实施例1

本发明是一种多分辨率NMF和Treelet融合的遥感图像变化检测方法，实验仿真环境为：MATLAB R2010a，CPU Inter Core2 1.80GHz，内存2G，Windows XPProfessional，参照图1，遥感图像变化检测方法包括有如下步骤：

(1)取同一地区在不同时间获取的已配准的两幅大小均为P×Q的遥感图像，时相1图像为X₁＝{X₁(p，q)，1≤p≤P，1≤q ≤Q}，参见图2(a)和图3(a)，时相2图像为X₂＝{X₂(p，q)，1≤p≤P，1≤q≤Q}，参见图2(b)和图3(b)，其中X₁(p,q)和X₂(p,q)分别为时相1图像和时相2图像在空间位置(p,q)处的像素点灰度值。将输入的两时相遥感图像对应空间位置像素的灰度值相减取绝对值，得到一幅差异图像X’_D(X’_D＝{X'_D(p,q)＝|X₂(p,q)-X₁(p,q)|})。对此差异图像X’_D进行大小为m×m像素的中值滤波，得到滤波后差异图像X_D。其中，m的取值范围为3、5、7、9，本例中m＝5。

(2)对滤波后的差异图像运用NMF算法提取五幅不同分辨率图像F_r(r＝1，2，...，5)。

2a)将滤波后的差异图像X_D分成大小均为h×h且不重叠的正方形图像块E，将每个小块E转成h²×1的列向量C_h，所有块的列向量合并构成矩阵V_h，其中，图像块尺寸h取2～10中的所有偶数值；

2b)对矩阵V_h采用NMF进行分解，得到NMF基矩阵W_h和系数矩阵H_h，对矩阵V_h采用NMF进行分解的过程是：

2b1)随机初始化W_h和H_h，设置最大迭代次数e和停止精度ε，对于不同的h的值，采用相同的最大迭代次数e和停止精度ε，分解得到的W_hH_h都能够较好地逼近矩阵V_h，因此本发明中采用相同的最大迭代次数e和停止精度ε。最大迭代次数e的取值范围为50～1000，停止精度ε的取值范围为10^-4～10^-6，本例中e的取值为100，ε的取值为10^-5。

2b2)在W_h和H_h均为非负矩阵的约束下，极小化目标函数||V_h-W_hH_h||²，对矩阵H_h进行更新，

$H_h=H_h\frac{W_h^T{V}_h}{W_h^T{W}_h{H}_h}---\left(1\right)$

对矩阵W_h进行更新，

$W_h=W_h\frac{V_h{H}_h^T}{W_h{H}_h{H}_h^T}---\left(2\right)$

直到满足或者达到最大迭代次数e为止，可得到NMF基矩阵W_h和系数矩阵H_h，其中，b为当前迭代次数，迭代次数b的取值范围为50～1000。对于不同的h的值，采用相同的当前迭代次数b，分解得到的W_hH_h都能够较好地逼近矩阵V_h，因此本发明中采用相同的当前迭代次数b。本例中b的取值为100。

2c)对滤波后的差异图像X_D边界进行扩展得到边界扩展图像D，将滤波后的差异图像X_D的首列向左扩展w列，末列向右扩展w列，对列扩展完成后的图像再进行首行向上扩展w行，末行向下扩展w行，即可得到边界扩展后的差异图像D，其中， 为向下取整符号。

2d)逐个选择图像D中的非边界扩展像素——即D中对应的差异图像X_D的每个像素点，作为中心像素点，以此为中心，也取大小为h×h的邻域块，将每个h×h邻域块转化成h²×1的列向量，这些列向量合并构成了矩阵V_vh，将矩阵V_vh在NMF的相应的基Wh上进行投影，得到特征数据集F_dh F_dh的行数为1，列数为P×Q。取特征数据集F_dh的前P个数据作为第一列，取P+1～2P个数据作为第二列，依此类推，取完特征数据集的数据，即可得到特征图像F_h。

2e)重复步骤(2a)至步骤(2d)直到图像块尺寸h依次取遍2，4，6，8，10，可得到五幅不同分辨率的图像F_r(r＝1，2,...，5)，其中r是图像的标号。

(3)估计滤波后差异图像X_D的噪声标准差对X_D取阈值得到图像Y；估计五幅不同分辨率图像F_r(r＝1，2,...，5)的噪声标准差对F_r(r＝1，2,...，5)取阈值得到图像Y_r(r＝1，2，...，5)，其中K为常数，K的取值为2。

(4)用Treelet算法对阈值后的差异图像Y和阈值后的五幅不同分辨率图像Y_r进行融合，得到一幅融合后的图像A。

4a)将阈值后图像Y和Y_r合并表示为Y_t(t＝1，2，...，6)，t为图像标号，将Y_t转为大小均为P×Q的列向量Y₁、Y₂、Y₃、Y₄、Y₅、Y₆，所有列向量构成初始样本矩阵X＝[Y₁，Y₂，Y₃，Y₄，Y₅，Y₆]；

4b)初始化Treelet变换的逐层聚类层数l＝0，1，...，L-1，L为矩阵X的列向量的个数，即L＝6，在第0层，每个变量采用初始样本X的列向量表示，初始化和变量的下标集δ＝{1，2，...，L}，初始化狄拉克基矩阵B⁰为L×L的单位阵，计算矩阵X的协方差矩阵和相关系数矩阵计算公式如下：

${\widehat{\mathrm{\Σ}}}_{\mathrm{ij}}=E\left[(X_u-E{X}_v)(X_u-E{X}_v)\right]---\left(3\right)$

${\hat{M}}_{\mathrm{ij}}=\frac{{\widehat{\mathrm{\Σ}}}_{\mathrm{ij}}}{\sqrt{{\widehat{\mathrm{\Σ}}}_{\mathrm{ii}}{\widehat{\mathrm{\Σ}}}_{\mathrm{jj}}}}---\left(4\right)$

其中，表示初始协方差矩阵第i行第j列的值，i＝1，2,...，P×Q，j＝1，2,...，P×Q，X_u和X_v表示样本矩阵X中两个不同的列向量，u＝1，2,...，6，v＝1，2,...，6；表示初始化相关系数矩阵第i行第j列的值；

4c)当l≠0时，由公式(4)计算相关系数矩阵找出最大的两个值，将最大值和次大值的对应位置序号分别记为α和β：

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

这里i＜j，分别表示相关系数矩阵中任意值的行和列，且只在变量下标集合δ内进行；

4d)对图像协方差矩阵进行局部主成分分析变换，得到第一主成分的和变量S1和第二主成分的差变量d₁，且使得图像协方差矩阵中α行β列的值和β行α列的值都为零，即得到旋转角度θ₁，并由下式得到雅克比旋转矩阵J：

其中，|θ₁|≤π/4；

4e)根据雅克比旋转矩阵J计算当前聚类层级l的狄拉克基矩阵B_l＝B_l-1J，狄拉克基矩阵B_l的第α和β列分别是尺度函数Φ_l和细节函数Ψ_l，前l层级的尺度向量{Φ_l}是尺度函数Φ_l和上一层的尺度向量{Φ_l-1}的合集，同时更新相关系数矩阵 ${\hat{M}}^{\left(l\right)}=J^T{\hat{M}}^{(l-1)}J$ 和协方差矩阵 ${\widehat{\mathrm{\Σ}}}^{\left(l\right)}=J^T{\widehat{\mathrm{\Σ}}}^{(l-1)}J;$

4f)将差变量的序号β的下标从和变量下标集δ中去除，即δ＝δ\{β}；

4g)重复步骤(4c)至步骤(4f)直至分解l＝L-1层，得到最终Treelet分解的基矩阵：

B＝[Φ_L-1，Ψ₁，...，Ψ_L-1]^T                                       (7)

其中，Φ_L-1∈{Φ_L-1}，Ψ₁为第一层Treelet变换得到的细节函数，Ψ_L-1是最高层Treelet变换得到的细节函数；

4h)将矩阵X与在狄拉克基矩阵B转置的方向进行投影，即R＝X×B^T，得到融合后的数据集R，将数据集R转化为P×Q的图像大小，得到融合后的图像A。

(5)采用区域生长算法对融合后的图像A进行分割，得到最终的变化检测结果。

本发明采用多分辨率NMF方法，不但能够利用NMF方法的优势，能够较好地逼近原矩阵，同时能够突出变化区域，而且多分辨率分析方法考虑了像素点的邻域信息，既保留了图像的细节信息又保留了图像的平滑区域信息，并且具有抗噪性，一定程度上提高了变化检测的准确性和可靠性。采用Treelet融合方法不但能够保留图像的强变化区域的信息还保留了图像的弱变化区域的信息，提高了变化检测的性能。

实施例2

多分辨率NMF和Treelet融合的遥感图像变化检测方法同实施例1，其中步骤(3)中对于阈值的计算是：

首先，估计滤波后的差异图像X_D的噪声标准差

${\mathrm{\σ}}_{X_D}=\mathrm{Median}\left(\right|X_D\left|\right)/0.6745---\left(8\right)$

其中，Median(·)表示取中值操作。

选取阈值若图像的灰度值小于阈值T_x，则将该值置为0，若图像的灰度值大于阈值T_x，则保留该值。得到一幅背景像素灰度值为0，变化区域保留的图像Y。其中K的取值为2。

然后，估计不同分辨率的特征图像F_r(r＝1，2,...，5)的噪声标准差

${\mathrm{\σ}}_{F_r}=\mathrm{Median}\left(\right|F_r\left|\right)/0.6745,r=\mathrm{1,2},...,5---\left(9\right)$

选取阈值对五幅图像进行阈值，若图像的灰度值小于阈值T_F，则将该值置为0，若图像的灰度值大于阈值T_F，则保留该值。得到五幅背景像素灰度值为0，变化区域保留的图像Y_r(r＝1，2，...，5)。其中K的取值为2。

步骤(5)对融合后的图像A采用区域生长算法进行分割的具体实现过程是：

首先，将阈值后的差异图Y和不同分辨率的特征图像Y_r(r＝1，2,...，5)进行对应点像素值相乘，选取灰度值大于0的像素点作为种子点；

然后，从图像A的第一个点开始，判断当前点是否为种子点，若当前点为种子点，则在当前点的8邻域内寻找与当前点灰度差小于阈值T_R的点，并把寻找到的点合并到种子点集合中，重复上述步骤，遍历全图。阈值T_R取值范围为10～30，本例中取为10；

最后，重复上一个步骤，直到找不到满足条件的点为止。所有的种子点集合就是变化区域。

实施例3

多分辨率NMF和Treelet融合的遥感图像变化检测方法同实施例1-2，其中步骤2b)中的最大迭代次数e的取值为50，停止精度ε的取值为10^-4，当前迭代次数b的取值为50，这种情况下，运算速度较快，检测精度稍微降低一些。

另外，步骤(5)中的阈值T_R的取值为20，也可以得到较好的变化检测结果。

实施例4

多分辨率NMF和Treelet融合的遥感图像变化检测方法同实施例1-2，其中步骤2b)中的最大迭代次数e的取值为1000，停止精度ε的取值为10^-6，当前迭代次数b的取值为1000，这种情况下，效率有所降低，但是也可以取得相近的变化检测效果。

另外，步骤(5)中的阈值T_R的取值为30，也可以得到较好的变化检测结果。

实施例5

多分辨率NMF和Treelet融合的遥感图像变化检测方法同实施例1-2

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

1.实验数据

本发明的实验数据是2组真实的遥感图像数据，第一组真实的遥感图像数据集是由2000年4月和2002年5月的墨西哥郊外的两幅Landsat 7 ETM+第4波段光谱图像组成，分别如图2(a)和图2(b)所示，图像大小均为512×512，256灰度级，配准误差为1.5个像素左右。图2(c)为参考变化图，其包含25599个变化像素和236545个非变化像素，白色区域表示变化的区域。第二组真实的遥感图像数据集是由1994年8月和1994年9月意大利Elba岛西部地区的两时相Landsat-5TM第4波段光谱图像组成，分别如图3(a)和图3(b)所示。图像大小均为326×414，灰度级为256，其参考变化图如图3(c)所示，包含2415个变化像素和132549个非变化像素，白色区域表示变化的区域。

2.实验方法

为了验证多分辨率NMF和Treelet融合的遥感图像变化检测方法的实验效果，我们将本发明方法与以下方法进行对比。

方法1：Celik等学者在文章“Unsupervised Change Detection in Satellite ImagesUsing Principal Component Analysis and K-Means Clustering，IEEE Geoscience andRemote Sensing Letters，2009，6(4)：772-776”中提出的方法。本发明与方法1进行对比是为了验证NMF对数据的非线性分类效果优于PCA对数据的线性分类效果。

方法2：黄世奇等学者在文章“基于小波变换的多时相SAR图像变化检测技术，测绘学报，2010，39(2)：180-186”中提出的方法。本发明与方法2进行对比是为了验证方法2中采用的平稳小波分解方法对图像具有平滑作用，造成一定程度上扩散了变化区域的面积，使得检测结果中虚警较高。本发明方法采用空域处理的多分辨率方法，达到了保留图像细节信息和平滑区域信息的均衡，减少了虚警点。

方法3：Shutao Li等学者在文章“Multitemporal Image Change Detection Using aDetail-Enhancing Approach With Nonsubsampled Contourlet Transform，IEEE GeoscienceandRemote Sensing Letters，2012，9(5)：836-840”中提出的方法。本发明与方法3进行对比是为了验证方法3的尺度内和尺度间融合方法能够保留强变化区域，但是丢失了弱变化区域，本发明方法中的Treelet融合方法既能保留强变化区域，也能保留弱变化区域，从而提高了变化检测的可靠性和准确性。

3.实验评价指标

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

4.实验结果与分析

图4是对第一组真实遥感图像实验数据采用不同方法的变化检测结果比较图。图4(a)是现有方法1的变化检测结果，从图中可看出方法1的检测结果图中存在较多漏检。图4(b)是现有方法2的变化检测结果，明显看出方法2对噪声敏感，存在较多虚警点，伪变化信息较多。图4(c)是现有方法3的变化检测结果，从图中可看出检测结果的边缘没有得到较好得保持。图4(d)是本发明的变化检测结果，从图中可看出本发明方法的检测结果，含有很少的杂点，伪变化信息较少，边缘保持得较好，较接近实际变化参考图。

实施例6

多分辨率NMF和Treelet融合的遥感图像变化检测方法同实施例1-2

图5是对第二组真实遥感图像实验数据采用不同方法的变化检测结果比较图。图5(a)是方法1的变化检测结果，图5(b)是方法2的变化检测结果，图5(c)方法3的变化检测结果，图5(d)是本发明的变化检测结果。从图中可看出方法1和方法2的检测结果，都产生了较多的虚警点，存在较多的伪变化信息，降低了变化检测的可靠性。方法3的检测结果更是产生了大范围的虚警点，导致变化区域提取失败，大大降低了变化检测的可靠性。本发明方法含有的伪变化信息少，边缘保持也较好，能够准确地检测出有效的变化区域，更接近实际的变化参考图。

实施例7

多分辨率NMF和Treelet融合的遥感图像变化检测方法同实施例1-2

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

从表中可以看出，对于第一组真实遥感图像数据集，方法1和方法3都存在较多的漏检，方法2和方法3都存在较多的虚警，本发明方法相对于方法1、方法2和方法3，总错误数分别减少了1458点、3089点、889点，正确率分别提高了0.55％、1.17％、0.33％。

对于第二组真实遥感图像数据集，方法1、方法2、方法3的虚警数都较高，总错误数也较高，本发明方法相对于方法1、方法2和方法3，总错误数分别减少了2820点、4123点、54595点，正确率分别提高了2.42％、3.05％、40.45％。

从以上两组数据集的实验结果可以看出，本发明方法的变化检测效果都优于方法1、方法2和方法3的变化检测效果，尤其当两时相数据集的变化和非变化区域的界限不是十分明显时更能体现本发明的优势。

综上，本发明的多分辨率NMF和Treelet融合的遥感图像变化检测方法，解决了传统方法在考虑图像的邻域结构时易受到孤立噪声点的影响，从单一分辨率角度考虑图像的细节和平滑区域的权衡的方法常无法较好的既保留图像的细节信息又保留平滑区域信息的问题。其实现过程是：输入两时相图像，采用直接差值法构造差异图像，对此差异图像进行大小为m×m像素的中值滤波；对滤波后差异图像运用NMF算法提取不同分辨率图像；对滤波后差异图像和不同分辨率特征图像取阈值；用Treelet算法融合阈值后的差异图像和不同分辨率的特征图像；采用区域生长算法对融合后的图像进行分割，得到最终的变化检测结果。本发明能够保持图像的细节信息与平滑区域信息，且能够去除孤立噪声，提高了变化检测精度，可用于灾情监测、土地利用、农业调查等领域。

Claims (2)

1.一种多分辨率NMF和Treelet融合的遥感图像变化检测方法，其特征在于：包括有如下步骤：

(1)取同一地区在不同时间获取的已配准的两幅大小均为P×Q的遥感图像，将获取的已配准的两图像对应空间位置像素的灰度值相减取绝对值，得到一幅差异图像，对此差异图像进行大小为m×m像素的中值滤波，其中，m的取值为3、5、7、9中的任意一值；

(2)对滤波后的差异图像运用NMF算法提取五幅不同分辨率图像F_r，r＝1，2，...，5；

2a)将滤波后的差异图像X_D分成大小均为h×h且不重叠的正方形图像块E，将每个小块E转成h²×1的列向量C_h，所有块的列向量合并构成矩阵V_h，其中，图像块尺寸h取2～10中的所有偶数值；

2b)对矩阵V_h采用NMF进行分解，得到NMF基矩阵W_h和系数矩阵H_h；

2c)对滤波后的差异图像X_D边界进行扩展得到边界扩展图像D，将滤波后的差异图像X_D的首列向左扩展w列，末列向右扩展w列，对列扩展完成后的图像再进行首行向上扩展w行，末行向下扩展w行，即可得到边界扩展后的差异图像D，其中， 为向下取整符号；

2d)逐个选择图像D中的非边界扩展像素——即D中对应的滤波后的差异图像X_D的每个像素点，作为中心像素点，以此为中心，也取大小为h×h的邻域块，将每个h×h邻域块转化成h²×1的列向量，这些列向量合并构成了矩阵V_vh，将矩阵V_vh在NMF的相应的基W_h上进行投影，得到特征数据集F_dh，将特征数据集F_dh转化成P×Q的图像大小，可得特征图像F_h；

2e)重复步骤2a)至步骤2d)直到图像块尺寸h依次取遍2，4，6，8，10，可得到五幅不同分辨率的图像F_r，r＝1，2，...，5，其中r是图像的标号；

(3)估计滤波后的差异图像X_D的噪声标准差，对X_D取阈值得到取阈值后的差异图像Y，估计五幅不同分辨率图像F_r的噪声标准差对F_r取阈值r＝1，2，...，5，得到取阈值后的五幅不同分辨率图像Y_r，r＝1，2，...，5，其中K为常数，K的取值为2；

(4)用Treelet算法对取阈值后的差异图像Y和取阈值后的五幅不同分辨率图像Y_r进行融合，得到一幅融合后的图像A；用Treelet算法对取阈值后的差异图像Y和取阈值后的五幅不同分辨率图像Y_r进行融合的过程是：

4a)将取阈值后的差异图像Y和取阈值后的五幅不同分辨率图像Y_r合并表示为Y_t，t＝1，2，...，6，t为图像标号，将Y_t转为大小均为P×Q的列向量Y₁、Y₂、Y₃、Y₄、Y₅、Y₆，所有列向量构成初始样本矩阵X＝[Y₁，Y₂，Y₃，Y₄，Y₅，Y₆]；

4b)初始化Treelet变换的逐层聚类层数l＝0，1，...，L-1，L为矩阵X的列向量的个数，即L＝6，在第0层，每个变量采用初始样本X的列向量表示，初始化和变量的下标集δ＝{1，2，...，L}，初始化狄拉克基矩阵B⁰为L×L的单位阵，计算矩阵X的初始协方差矩阵和初始相关系数矩阵计算公式如下：

${\widehat{\mathrm{\Σ}}}_{\mathrm{ij}}=E\left[(X_u-{\mathrm{EX}}_v)(X_u-{\mathrm{EX}}_v)\right]$

${\hat{M}}_{\mathrm{ij}}=-\frac{{\widehat{\mathrm{\Σ}}}_{\mathrm{ij}}}{\sqrt{{\widehat{\mathrm{\Σ}}}_{\mathrm{ii}}{\widehat{\mathrm{\Σ}}}_{\mathrm{jj}}}}$

其中，表示初始协方差矩阵第i行第j列的值，i＝1，2，...，P×Q，j＝1，2，...，P×Q，X_u和X_v表示样本矩阵X中两个不同的列向量，u＝1，2，...，6，v＝1，2，...，6；表示初始化相关系数矩阵第i行第j列的值；

4c)当l≠0时，寻找相关系数矩阵中最大的两个值，将最大值和次大值的对应位置序号分别记为α和β：

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

这里i＜j，分别表示相关系数矩阵中任意值的行和列，且只在变量下标集合δ内进行；

4d)对图像协方差矩阵进行局部主成分分析变换，得到第一主成分的和变量s₁和第二主成分的差变量d₁，且使得图像协方差矩阵中α行β列的值和β行α列的值都为零，即得到旋转角度θ_l，并由下式得到雅克比旋转矩阵J：

其中，|θ_l|≤π/4；

4e)根据雅克比旋转矩阵J计算当前聚类层级l的狄拉克基矩阵B_l＝B_l-1J，狄拉克基矩阵B_l的第α和β列分别是尺度函数Φ_l和细节函数Ψ_l，前l层级的尺度向量{Φ_l}是尺度函数Φ_l和上一层的尺度向量{Φ_l-1}的合集，同时更新相关系数矩阵 ${\hat{M}}^l=J^T{\hat{M}}^{l-1}J$ 和图像协方差矩阵 ${\widehat{\mathrm{\Σ}}}^l=J^T{\widehat{\mathrm{\Σ}}}^{l-1}J;$

4f)将差变量的序号β的下标从和变量下标集δ中去除，即δ＝δ\{β}；

4g)重复步骤4c)至步骤4f)直至分解l＝L-1层，得到最终Treelet分解的基矩阵：

${B=[{\mathrm{\Φ}}_{L-1},{\mathrm{\Ψ}}_1,...,{\mathrm{\Ψ}}_{L-1}]}^T$

其中，Φ_L-1∈{Φ_L-1}，Ψ₁为第一层Treelet变换得到的细节函数，Ψ_L-1是最高层Treelet变换得到的细节函数；

4h)将矩阵X与在狄拉克基矩阵B转置的方向进行投影，即R＝X×B^T，得到融合后的数据集R，将数据集R转化为P×Q的图像大小，得到融合后的图像A；

(5)采用区域生长算法对融合后的图像A进行分割，得到最终的变化检测结果。

2.根据权利要求1所述的多分辨率NMF和Treelet融合的遥感图像变化检测方法，其特征在于：步骤2b)中对矩阵V_h采用NMF进行分解的过程是：

2b1)随机初始化W_h和H_h，设置最大迭代次数e和停止精度ε，最大迭代次数e的取值范围为50～1000，停止精度ε的取值范围为10^-4～10^-6；

2b2)在W_h和H_h均为非负矩阵的约束下，极小化目标函数||V_h-W_hH_h||²，对矩阵H_h进行更新，

$H_h=H_h\frac{W_h^T{V}_h}{W_h^T{W}_h{H}_h}$

对矩阵W_h进行更新，

$W_h=W_h\frac{V_h{H}_h^T}{W_h{H}_h{H}_h^T}$

直到满足或者达到最大迭代次数e为止，可得到NMF基矩阵W_h和系数矩阵H_h，其中，b为当前迭代次数，迭代次数b的取值范围为50～1000。