CN103955926A - 基于Semi-NMF的遥感图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于Semi-NMF的遥感图像变化检测方法，处理对象同时包含了光学遥感和合成孔径雷达图像，主要解决现有遥感图像变化检测方法在获得较强变化区域的同时无法检测到微弱、较小的变化区域，以及不能有效保留更多细节和边缘信息的问题。其实现过程包括：(1)根据遥感图像类型产生差异图像。(2)通过PCA获得差异图像每个像素所对应的基于邻域信息的特征向量，并利用其构建特征矩阵X。(3)对X执行Semi-NMF算法，经过迭代运算将其分解为基矩阵F和系数矩阵G。(4)根据系数矩阵判定变化类ω_c和不变化类ω_u，实现软聚类的功能，得到二值变化检测结果。本发明减少边缘信息的丢失，同时检测较强和微弱、较小的变化区域，减少总错误率，保留了更多的细节信息，有效地获取变化结果。

Description

基于Semi-NMF的遥感图像变化检测方法

技术领域

本发明属于数字图像处理技术领域，主要涉及多时相遥感图像变化检测研究方向，具体地说是基于Semi-NMF(Semi-Nonnegative Matrix Factorization，半非负矩阵分解)的遥感图像变化检测方法，处理对象同时包含了多时相光学遥感图像和合成孔径雷达(Synthetic Aperture Radar，SAR)图像。该方法可应用于对地观测的诸多实际问题。

背景技术

近年来,随着航天技术、传感器技术、计算机技术及相关学科的迅猛发展,遥感技术亦得以不断进步，各种搭载光学和SAR传感器的遥感平台(如飞机或卫星等)相继运行，可大面积、迅速、动态地进行对地观测，获取宽幅、多时相的遥感图像，为地物感知和检测提供有效数据源。受制于人类活动加剧以及自然灾害频发，利用不同时间所获取的同一地区遥感图像进行地表变化检测一直是遥感应用研究的热点,并且已经广泛地应用于国民经济和国防建设的诸多方面，譬如土地利用、环境监测、森林资源调查、城市规划、灾情评估等。

多时相遥感图像变化检测是指通过对同一区域、不同时期获取的两幅或多幅遥感图像进行比较分析，进而根据图像之间差异来获取所感兴趣地物、场景或目标的变化信息，主要包括三个基本步骤：1)图像预处理；2)变化信息检测提取；3)后处理及性能评估。本发明主要是针对第2步变化信息检测提取的创新工作。

截至目前，国内外学者提出了很多有效的变化检测方法，概括起来可分为有监督变化检测和无监督变化检测。其中，有监督变化检测需要事先获得关于地面真实类别的样本进行训练，但在大多数情况下，关于地面的真实情况是难以或无法获取的，在很大程度上制约了其在实际中的应用。而无监督变化检测则是针对不同时期的图像直接进行比较分析，不需要其它信息，相对简单、直观，因此得到了广泛研究和应用。其大致分为：(1)基于分布模型差异的多时相遥感图像变化检测；(2)基于差异图像分析的多时相遥感图像变化检测；以及(3)基于马尔科夫融合的多时相遥感图像变化检测。特别是对第二类算法的研究最为普遍，核心思想是将变化检测问题视为图像的二元分类/分割问题，又可细分为聚类分析、智能优化、阈值分割、有限混合模型、马尔科夫随机场、主动轮廓和水平集等策略,涉及图像处理、模式识别与机器视觉领域的许多基本理论和方法。其中，聚类分析由其简单、有效而受到普遍认可。

T.Celik在文献[T.Celik,“Unsupervised change detection in satellite images usingprincipal component analysis and K-means clustering,”IEEE Geoscience and Remote SensingLetters,vol.6,no.4,pp.772-776,2009.]中最早提出聚类分析的遥感图像变化检测方法。在此基础上，申请人通过低阶次分数阶Fourier变换产生噪声抑制的差异图像，以降低噪声对检测性能的影响。两者均使用最简单的硬聚类算法(即K-means聚类)实现差异图像特征聚类，产生对应的变化类和不变化类，但限于K-means聚类算法的局限性，难以获得比较好的检测结果。后续，A.Ghosh等[A.Ghosh,N.S.Mishra,and S.Ghosh,“Fuzzy clustering algorithms for unsupervised change detectionin remote sensing images”,Information Sciences,vol.181,no.4,pp.699-715,2011.]以及M.Volpi等[M.Volpi,D.Tuia,G.Camps-Valls,and M.Kanevski,“Unsupervised change detectionwith kernels,”IEEE Geoscience and Remote Sensing Letters,vol.9,no.6,pp.1026-1030,2012]等分别采用模糊C均值聚类(Fuzzy C-means)、Gustafson-Kesse模糊聚类和核K-means聚类算法进行多时相遥感图像变化检测，在一定程度上提高了检测性能，但仍具有较高的错误检测概率，特别是针对高分辨率的遥感图像。因此，有必要发明一种简单且能获得更好检测效果的多时相遥感图像变化检测方法，以提高遥感图像变化检测在实际中的应用价值。

发明内容

本发明提出一种基于Semi-NMF的遥感图像变化检测新方法。该方法借助主成分分析(PrincipalComponent Analysis,PCA)提取差异图像每个像素对应的紧致特征矢量，以构建特征矩阵；然后执行Semi-NMF算法将其分解为基矩阵F和系数矩阵G。其中，系数矩阵为差异图像像素的隶属度矩阵，对其进行判决实现像素集的软聚类，得到最终变化检测结果。具体实现过程包括如下步骤：

(1)利用在同一地区不同时间所获得两幅大小相同、相互配准的遥感图像Y₀和Y₁以产生差异图像Y_D，对于光学遥感图像和SAR图像计算方式不同。

(2)构建特征矩阵X。这一步主要是利用PCA获得差异图像中每一个像素的特征向量，利用这些特征向量构建特征矩阵X。具体的执行过程可以分为如下四步。

(2a)将差异图像Y_D分成h×h非重叠子块，并且h≥2。按照一定的顺序将差异图像块排列为向量形式，表示为y_D(x，y)。

(2b)对非重叠子块所对应的向量集执行PCA以生成本征向量空间。首先计算上述列向量集的平均向量μ；然后将每一个向量减去平均向量μ得到差值向量集Δ_a，并构建协方差矩阵C；最后经奇异值分解得到协方差矩阵C的特征向量{e_s}和特征值{λ_s}，并将特征向量依据与它们相对应的特征值大小降序排列。

(2c)对于差异图像的每一个像素坐标(i，j)，映射y_D(x，y)到本征向量空间生成特征向量，

V(i，j)＝[v₁ v₂...v_s]^T

其中1≤S≤h²，并且参数S决定在空间位置(i，j)的特征向量V(i，j)的维数，y_D(i，j)是在差异图像上获得的重叠块，按照与(2a)中相同的规则排列成向量形式。

(2d)将差异图像所有像素的特征向量V(i，j)作为列向量构建大小为S×HW的特征矩阵X。

(3)执行Semi-NMF算法。该步主要是对步骤(2)中获得特征矩阵X执行Semi-NMF算法，将其分解为基矩阵F和系数矩阵G。

为了找到基矩阵F和系数矩阵G，我们必须求解代价函数(即最小化重构误差)使其收敛。鉴于代价函数对F和G而言不是联合凸的，可使用F和G交替更新的迭代算法解决上述最小化问题，包含下面两个步骤，迭代前首先对Semi-NMF算法进行随机初始化，

(3a)保持系数矩阵G不变，更新基矩阵F，

(3b)保持基矩阵F不变，更新系数矩阵G，

然后验证是否收敛。若达到收敛条件，则Semi-NMF算法运算结束。若没有达到收敛条件，则返回(3a)处继续迭代。

(4)获得变化检测结果。Semi-NMF算法能够实现软聚类功能，优于K-means与核K-means聚类算法，可以减弱噪声和异质像素影响，完整地从差异图像中提取变化信息。特征矩阵X经过Semi-NMF算法处理获得大小为2×HW的矩阵G^T，该矩阵就是差异图像像素归属于变化类ω_c和不变化类ω_u的隶属度指示矩阵。根据矩阵G^T每一列中隶属度大小是判定所代表像素属于变化类ω_c或者不变化类ω_u，

$c(i,j)=\left\{\begin{array}{cc}0,& G^T(1,(i-1)*H+j)\&GreaterEqual;G^T(2,(i-1)*H+j)\\ 1,& G^T(1,(i-1)*H+j)<G^T(2,(i-1)*H+j)\end{array}\right.$

然后分别计算两类所包含差异图像像素大小的平均值，对于均值较大的像素类设置为变化类ω_c，另外一个像素类设置为不变化类ω_u，获得二值的变化检测图。

对于不同的图像来源，差异图像Y_D有不同的获取方式：对于光学遥感图像，差异图像Y_D为输入图像Y₀和Y₁直接相减的绝对值，Y_D＝|Y₁-Y₀|；对于合成孔径雷达(SAR)图像，其差异图像Y_D取输入图像Y₀和Y₁作对数比的绝对值， $Y_D=\left|\log \frac{Y_1}{Y_0}\right|=|\log \left(Y_1\right)-\log \left(Y_0\right)|.$

有益效果

1、本发明结合每个像素的邻域信息，利用PCA提取与差异图像中每个像素对应的紧致特征向量，并用其构建具有优异判决能力的特征矩阵，可有效减弱噪声和异质像素影响。

2、本发明首次将Semi-NMF算法应用到遥感图像变化检测中，主要用于对特征矩阵分解，获得差异图像像素的隶属度矩阵，实现软聚类的功能。该算法能够有效区分变化类和不变化类像素，准确实现变化检测。

3、该方法能够在完整获得较强变化区域同时，检测到微弱、较小的变化区域，减少边缘信息的丢失。可兼顾错误检测率和丢失检测率，使得总错误率降低，保留了更多细节信息，有效地获取变化结果。

附图说明

图1基于Semi-NMF的遥感图像变化检测方法整体框架；

图2仿真实验使用的第一组遥感图像数据集Bangladesh，(a)2007年4月12日获得的ESA EnvisatASAR图像，(b)2007年7月26日获得的ESA Envisat ASAR图像，(c)变化检测参考图；

图3仿真实验使用的第二组遥感图像数据集lake，(a)1986年8月5日获得的光学遥感图像，(b)1992年8月5日获得的光学遥感图像，(c)变化检测参考图，(d)色基(Color Key)模型图；

图4二组遥感图像数据集分别对应的差异图像：(a)Bangladesh数据集差异图像，(b)lake数据集差异图像；

图5本发明与对比方法仿真SAR图像数据Bangladesh的变化检测结果，(a)对比方法1的变化检测结果，(b)对比方法2的变化检测结果，(c)本发明的变化检测结果，(d)变化检测参考图；

图6本发明与对比方法仿真光学遥感图像数据lake的变化检测结果，(a)对比方法1的变化检测结果，(b)对比方法2的变化检测结果，(c)本发明的变化检测结果，(d)变化检测参考图；

具体实施方式

下面参照附图描述本发明的实施例。应当理解，在此描述的本发明各个实施例仅仅是为了更好地解释本发明的原理和概念，而不是要限制本发明。在阅读这样的描述之后，本领域技术人员很容易构造其他修改或替换，这样的修改或替换应被理解为落入本发明的范围中。

图1给出了本发明实施例的流程框架，具体实施包含以下步骤：

(1)利用多时相遥感图像产生差异图像。两幅大小为H×W且相互配准的遥感图像Y₀＝{y₀(i，j)|1≤i≤H，1≤j≤W}和Y₁＝{y₁(i，j)|1≤i≤H，1≤j≤W}，它们是在同一地区不同时间t₀和t₁分别获得的遥感数据。根据遥感图像的类型生成差异图像Y_D。

对于光学图像(Optical Images)，输入图像Y₀和Y₁直接相减的绝对值就是差异图像Y_D，

Y_D＝|Y₁-Y₀|

而对于SAR(Synthetic Aperture Radar)图像，考虑乘性相干斑噪声，差异图像Y_D是输入图像Y₁和Y₀的比值图像，同时为了提高低幅度的像素值，我们对比值图像执行一个对数运算，

$Y_D=\left|\log \frac{Y_1}{Y_0}\right|=|\log \left(Y_1\right)-\log \left(Y_0\right)|$

其中log是自然对数算子，能够使差异图像中变化类ω_c和不变化类ω_u分布更加均匀。同时使Y₀和Y₁内在残留的乘性(×)斑点干扰转换成加性(+)干扰。

(2)构建特征矩阵X。该步主要是利用PCA(主成分分析)获得差异图像中每一个像素对应的特征向量，利用这些特征向量构建特征矩阵X。具体的执行过程可以分为如下四步。

(2a)将差异图像Y_D分成h×h非重叠子块，并且h≥2。以(x，y)为参考坐标的h×h差异图像块可以表示为当h是奇数时，(x，y)正好处于图像块的中心，其中是数学上的求上限操作符，得到比某一数大的最小整数，例如按照一定的顺序将差异图像块Y_D(x，y)排列为向量形式表示为，

(2b)对非重叠子块对应的向量集执行PCA以生成本征向量空间。为了简化数学表示用表示向量Y_D(x，y)，其中a代表一个数字索引， 是数学上的求下限操作符。

首先计算上述向量集的平均向量，即

$\mathrm{\&mu;}=\frac{1}{M}\underset{a=1}{\overset{M}{\mathrm{\&Sigma;}}}{y}_D^a$

然后将每一个向量减去平均向量μ：此处1≤a≤M，并构建协方差矩阵

$C=\frac{1}{M}\underset{a=1}{\overset{M}{\mathrm{\&Sigma;}}}{{\mathrm{\&Delta;}}_a\mathrm{\&Delta;}}_a^T$

其中T为转置符号，矩阵C大小为h²×h²。为了更好的描述这些数据，使用奇异值分解得到协方差矩阵C的h²个特征向量{e_s}和特征值{λ_s}。此处假定矩阵C得到的特征向量依据与它们相对应的特征值降序排列，例如λ_s≥λ_s+1。

(2c)对于差异图像的每一个像素坐标(i，j)，映射y_D(i，j)到本征向量空间生成特征向量，

V(i，j)＝[v₁ v₂...v_s]^T

其中1≤S≤h²，并且参数S决定在空间位置(i，j)的特征向量V(i，j)的维数，y_D(i，j)是在差异图像上获得的重叠块，按照与(2a)中相同的规则排列成向量形式。

(2d)将差异图像所有像素的特征向量V(i，j)作为列向量构建大小为s×HW的特征矩阵X。

(3)执行Semi-NMF算法。对步骤(2)中获得特征矩阵X执行Semi-NMF算法，将其分解为基矩阵F和系数矩阵G，如下所示，

X≈FG^T

根据Semi-NMF算法可知，特征矩阵X和基矩阵F没有非负的约束，允许存在负值，但是系数矩阵G必须满足非负的条件。

为了找到基矩阵F和系数矩阵G，我们必须求解如下代价函数(即最小化重构误差)，

arg min_F，G||X-FG^T||_F s.t.G≥0

其中，||·||_F表示Frobenius范数，G≥0表示矩阵G中的每一个元素都不小于0。

鉴于代价函数对F和G而言不是联合凸的，下面使用F和G交替更新的迭代算法解决上述最小化问题，包含下面两个步骤，迭代前首先对Semi-NMF算法进行随机初始化，

(3a)保持系数矩阵G不变，更新基矩阵F，

F＝XG(G^TG)^-1

(3b)保持基矩阵F不变，更新系数矩阵G，

$G_{\mathrm{ik}}\&LeftArrow;G_{\mathrm{ik}}\sqrt{\frac{{\left(X^{T}F\right)}_{\mathrm{ik}}^{+}+{\left[G{\left(F^{T}F\right)}^{-}\right]}_{\mathrm{ik}}}{{\left(X^{T}F\right)}_{\mathrm{ik}}^{-}+{\left[{G\left(F^{T}F\right)}^{+}\right]}_{\mathrm{ik}}}}$

然后验证是否收敛。若达到收敛条件，则Semi-NMF算法运算结束。若没有达到收敛条件，则返回(3a)处继续迭代。注意，矩阵A的正和负部分分别表示为，

$A_{\mathrm{ik}}^{+}=\left(\right|A_{\mathrm{ik}}|+A_{\mathrm{ik}})/2,A_{\mathrm{ik}}^{-}=\left(\right|A_{\mathrm{ik}}|-A_{\mathrm{ik}})/2$

(4)获得变化检测结果。Semi-NMF算法能够实现软聚类的功能，优于K-means与核K-means算法，可减弱噪声和异质像素影响，完整地从差异图像中提取变化信息。特征矩阵X经过Semi-NMF算法处理获得大小为2×HW的矩阵G^T，该矩阵就是差异图像像素归属于变化类ω_c和不变化类ω_u的隶属度指示矩阵。根据矩阵G^T每一列中隶属度大小是判定所代表像素属于变化类ω_c或者不变化类ω_u，

$c(i,j)=\left\{\begin{array}{cc}0,& G^T(1,(i-1)*H+j)\&GreaterEqual;G^T(2,(i-1)*H+j)\\ 1,& G^T(1,(i-1)*H+j)<G^T(2,(i-1)*H+j)\end{array}\right.$

然后分别计算两类所包含差异图像像素大小的平均值，对于均值较大的像素类设置为变化类ω_c，另外一个像素类设置为不变化类ω_u，获得二值的变化检测图。

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

1、仿真实验条件

本发明提出的基于Semi-NMF的遥感图像变化检测方法及对比方法都是在相同的实验仿真平台上实现的。实验计算机配置为AMD Athlon(tm)7750双核处理器，主频2.70GHz，RAM为4.00GB(3.25GB可用)，32位操作***，软件平台为MATLAB7.13。

2、仿真实验数据

本发明的仿真实验使用的是真实的多时相遥感数据，包括光学遥感图像和SAR图像。

第一组遥感图像数据集Bangladesh，图2中是两幅SAR图像数据和一幅通过人工分析产生的地面真实变化参考图。其中(a)和(b)是两幅大小相同的ASAR(Advanced Synthetic Aperture Radar)图像，由ESAEnvisat卫星在Bangladesh地区(印度)分别于2007年4月12日和2007年7月26日获得的。第一幅图像是2007年4月12日洪水发生前该区域的地表信息，第二幅图像是2007年7月26日洪水发生后同一区域的地表信息，两幅图像的大小为300×300，灰度级为256。第三幅图像是变化参考图，它是通过人工分析已知的图像数据结合当地真实的地面变化信息而得到的一幅二值图像，图中的黑色像素区域表示地面的变化情况。

第二组遥感图像数据集lake，图3中(a)和(b)是两幅光学遥感图像，(c)是变化参考图，(d)是一幅色基(Color Key)模型图。其中图像(a)、(b)和(c)的大小均为200×200，而且前两幅图像分别于1986年8月5日和1992年8月5日在Lake Tahoe，Reno，Nevada获得。这一区域受到了干旱的影响，导致湖面变化，同时还有土地使用和森林砍伐导致地表区域发生变化，第三幅图像是该区域的变化参考图。该区域的变化参考图是利用Landsat多光谱扫描器(Multi Spectral Scaner，MSS)和专题制图仪(ThematicMapper，TM)数据获得，用于识别和描述地表变化信息。

图3中(d)是对颜色进行编码的Color Key模型图，用于解释变化检测参考图。在Color Key图中，图像中绿色的增加和减少依次用绿色和品红(Magenta)来表示，而亮度的减弱显示为不同深浅的蓝色，亮度的增强显示为橙色/红色色调，表1给出了这些对应关系。任何亮度和绿色变化的组合都可以用ColorKey图中对应的颜色来表征。在图3(c)的变化检测参考图中，灰白色表示没有变化或者不重要的变化；湖的边界区域显示为橘红色，意味着由于湖面干涸而造成亮度的显著增强；橙色的小块区域可能是由于土地利用变化或木材砍伐等森林地区的活动导致的变化。

表1图像变化与Color Key模型图之间的对应关系

图像的变化	绿色增加	绿色减少	亮度增强	亮度减弱
					Color Key模型图	绿色	品红	橙色/红色	蓝色

3、结果评价指标

一般在实际变化已知的情况下，衡量变化检测算法性能有如下三个指标。

(1)P_FA：错误检测率，表示错误检测像素数(False Alarms，FA)占图像像素总数(N)的百分比，即P_FA＝FA/N×100。错误检测像素数表示实际未变化而检测成变化的像素总数。

(2)P_MA：丢失检测率，表示丢失检测像素数(Missed Alarms，MA)占图像像素总数(N)的百分比，即P_MA＝MA/N×100。丢失检测像素数表示实际变化而检测成未变化的像素总数。

(3)P_TE：总错误率，表示总错误像素数(Total Errors，TE)占图像像素总数(N)的百分比，即P_TE＝（FA+MA）/N×100。总错误像素数是上面的错误检测像素数(FA)和丢失检测像素数(MA)的总和。这是衡量变化检测算法最重要的指标，是证明结果优劣的主要判断依据。

4、对比仿真实验

为了说明本发明的有效性，我们将本发明的仿真效果与如下两个方法进行比较：

对比方法1，文献[T.Celik,“Unsupervised change detection in satellite images using principalcomponent analysis and K-means clustering,”IEEE Geoscience and Remote Sensing Letters,vol.6,no.4,pp.772-776,2009.]中利用主成分分析(PCA，Principal Component Analysis)和K-means聚类算法实现遥感图像的变化检测。该方法首先利用主成分分析提取每个像素对应的特征向量，再使用K-means算法将特征向量分成两类(即变化类和不变化类)，得到变化检测图。

对比方法2，文献[Y.Q.Cheng,H.C.Li,T.Celik,and F.Zhang,“FRFT-based improved algorithmof unsupervised change detection in SAR images via PCA and K-means clustering”,in Proceedingof IEEE International Geoscience and Remote Sensing Symposium,Melbourne,Australia,July2013:1952-1955.]是利用低阶次的分数阶Fourier变换对上面文献提出的方法进行改进。

5、仿真实验内容

实验1计算差异图像

本发明的仿真实验采用两组遥感图像数据集，包括SAR图像和光学遥感图像。根据遥感图像的类型依据不同的公式，生成差异图像Y_D，分别如图4(a)和(b)所示。

实验2SAR图像数据Bangladesh

本次仿真实验使用的是第一组遥感图像数据集Bangladesh。图5是本发明与对比方法仿真SAR图像数据Bangladesh的变化检测结果，其中图5(a)是对比方法1的变化检测结果，图5(b)是对比方法2的变化检测结果，图5(c)是本发明的变化检测结果，图5(d)是变化检测参考图。

除了依据视觉观察进行定性对比外，还进行定量比较。表2列出了仿真实验输出图像的精度评估数据。

表2本发明与对比方法仿真Bangladesh数据变化检测结果精度评估数据

Algorithm	PFA	PMA	PTE
				对比方法1	0％	4.84％	4.84％
对比方法2	0％	4.76％	4.76％
				本方法	1.39％	0.57％	1.96％

实验3光学遥感图像数据lake

本次仿真实验使用的是第一组遥感图像数据集lake。图6是本发明与对比方法仿真光学遥感图像数据lake的变化检测结果，其中图6(a)是对比方法1的变化检测结果，图6(b)是对比方法2的变化检测结果，图6(c)是本发明的变化检测结果，图6(d)是变化检测参考图。

6、仿真结果分析

通过分析本发明与对比方法仿真SAR图像数据Bangladesh的变化检测结果及精度评估数据可知，本发明仿真结果的总错误率为1.96％，相比于其它对比方法得到了较好的结果。而且本发明的丢失检测率为0.57％，相比于对比方法的4.84％和4.76％明显较低，检测性能更好。对比方法的错误检测率为0％，主要是因为对比算法不能完整检测变化区域导致的。

通过分析本发明与对比方法仿真光学遥感图像数据lake的变化检测结果及图像数据集变化参考图信息可知，该图像数据主要显示了湖面干涸、土地使用和森林砍伐导致的变化。在变化检测参考图中，灰白色表示没有变化或者不重要的变化；湖的边界区域显示为橘红色，意味着由于湖面干涸而造成亮度的显著增强；橙色的小块区域可能是由于土地利用变化或木材砍伐等森林地区的活动导致的变化。对比方法能够较好的检测湖面干涸导致的变化，但是无法有效检测土地使用和森林砍伐导致的变化，丢失大量微弱变化区域。本发明仿真结果相比于其它对比方法不仅能够较完整地检测湖面干涸导致的变化，同时保留了土地使用和森林砍伐导致的细微变化。

综上所述，本发明同时适用于检测多时相光学遥感图像和SAR图像，借助PCA得到紧致特征向量，并用其构建具有优异判决能力的特征矩阵。然后执行Semi-NMF算法，将特征矩阵分解为基矩阵F和系数矩阵G。系数矩阵就是差异图像像素的隶属度矩阵，对像素集实现软聚类的功能，减少噪声和异质像素的影响，得到较优的变化检测结果。通过定性和定量的比较，可以证明本发明所提出的基于Semi-NMF的遥感图像变化检测方法，能够在完整获得较强变化区域的同时检测到微弱、较小的变化区域，减少边缘信息的丢失。并且兼顾了丢失检测率和错误检测率，使得总错误率降低，保留了更多细节信息，有效地获取变化结果。

本领域的普通技术人员显然清楚并且理解，本发明方法所举的以上实施例仅用于说明本发明方法，而并不用于限制本发明方法。虽然通过实施例有效描述了本发明，本发明存在许多变化而不脱离本发明的精神。在不背离本发明方法的精神及其实质的情况下，本领域技术人员当可根据本发明方法做出各种相应的改变或变形，但这些相应的改变或变形均属于本发明方法要求的保护范围。

Claims (5)

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

步骤1，利用在同一地区不同时间所获得两幅大小相同、相互配准的遥感图像Y₀和Y₁以产生差异图像Y_D；

步骤2，构建特征矩阵X，该步主要是通过PCA获得差异图像中每一个像素对应的特征向量，利用这些特征向量构建特征矩阵X；

步骤3，对步骤2中获得特征矩阵X执行Semi-NMF算法，将其分解为基矩阵F和系数矩阵G；

步骤4，根据系数矩阵G判定变化类和不变化类，变化类像素位置置0，不变化类像素位置置1，得到最终的二值变化检测结果。

2.根据权利要求1所述的遥感图像变化检测方法，其特征在于步骤2所述的构建特征矩阵X，具体实施包括：

(a)将差异图像Y_D分成h×h非重叠子块，并且h≥2；按照一定的顺序将差异图像块排列为向量形式，表示为y_D(x，y)；

(b)对非重叠子块所对应的向量集执行PCA以生成本征向量空间；首先计算上述向量集的平均向量μ；然后将每一个向量减去平均向量μ得到差值向量集Δ_a，并构建协方差矩阵C；最后经奇异值分解得到协方差矩阵C的特征向量{e_s}和特征值{λ_s}，并将特征向量依据与它们相对应的特征值大小降序排列；

(c)对于差异图像的每一个像素坐标(i，j)，映射y_D(i，j)到本征向量空间生成特征向量，

V(i，j)＝[v₁ v₂...v_s]^T

其中1≤S≤h²，并且参数S决定在空间位置(i，j)的特征向量V(i，j)的维数，y_D(x，y)是在差异图像上获得的重叠块，按照与(a)中相同的规则排列成向量形式；

(d)将差异图像所有像素的特征向量V(i，j)作为列向量构建大小为S×HW的特征矩阵X。

3.根据权利要求1所述的遥感图像变化检测方法，其特征在于步骤3所述的对特征矩阵X执行Semi-NMF算法，将其分解为基矩阵F和系数矩阵G，为找到基矩阵F和系数矩阵G，使用F和G交替更新的迭代策略解决最小化重构误差，包含下面两个步骤，迭代前首先对Semi-NMF算法进行随机初始化，

(a)保持系数矩阵G不变，更新基矩阵F，

(b)保持基矩阵F不变，更新系数矩阵G，

然后验证是否收敛；若达到收敛条件，则Semi-NMF算法运算结束；若没有达到收敛条件，则返回(a)处继续迭代。

4.根据权利要求1所述的遥感图像变化检测方法，其特征在于步骤4所述的获得变化检测结果，具体实施为：特征矩阵X经过Semi-NMF算法处理获得大小为2×HW的矩阵G^T，该矩阵就是差异图像像素归属于变化类ω_c和不变化类ω_u的隶属度指示矩阵；根据矩阵G^T每一列中隶属度大小是判定所代表像素属于变化类ω_c或者不变化类ω_u，

$c(i,j)=\left\{\begin{array}{cc}0,& G^T(1,(i-1)*H+j)\&GreaterEqual;G^T(2,(i-1)*H+j)\\ 1,& G^T(1,(i-1)*H+j)<G^T(2,(i-1)*H+j)\end{array}\right.$

然后分别计算两类所包含差异图像像素大小的平均值，对于均值较大的像素类设置为变化类ω_c，另外一个像素类设置为不变化类ω_u，获得二值的变化检测图。

5.根据权利要求1所述的遥感图像变化检测方法，其特征在于所述的差异图像Y_D：对于光学遥感图像，差异图像Y_D为输入图像Y₀和Y₁直接相减的绝对值，Y_D＝|Y₁-Y₀|；对于合成孔径雷达(SAR)图像，其差异图像Y_D取输入图像Y₀和Y₁作对数比的绝对值，