基于非负矩阵分解和核FCM的多时相遥感影像变化检测方法

摘要

本发明公开了基于非负矩阵分解和核FCM的多时相遥感影像变化检测方法，该方法首先基于非负矩阵分解算法融合多时相遥感影像的变化矢量幅值和多时相的光谱角映射图，然后将融合结果作为核FCM的输入，再基于核FCM结合空间邻域信息的方法，获取最终的变化检测结果。本发明在基于核FCM的变化检测中，联合变化矢量幅值和光谱角映射图作为输入，利用这两个特征，优于仅仅使用变化矢量幅值的FCM方法。在变化检测中，采用伪训练样本，并采用变化指数准则来自适应选择核FCM结合邻域空间信息的检测方法中的相关参数，可以使得变化检测的结果更加可靠，也更加具有稳健性。

说明书

技术领域

本发明涉及多时相高分辨率光学遥感影像变化检测方法，具体涉及基于非负矩阵分解和核FCM的多时相遥感影像变化检测方法，属于遥感影像处理技术领域。

背景技术

随着多时相高分辨率遥感数据的不断积累以及空间数据库的相继建立，如何从这些遥感数据中提取和检测变化信息已成为遥感科学和地理信息科学的重要研究课题。根据同一区域不同时相的遥感影像，可以提取城市、环境等动态变化的信息，为资源管理与规划、环境保护等部门提供科学决策的依据。我国“十二五”将加大拓展实施“十一五”已启动实施的高分辨率对地观测工程，关注包括高分辨率遥感目标与空间环境特征分析及高可靠性自动解译等基础理论与关键技术研究，正在成为解决国家安全和社会经济发展重大需求的研究焦点。

遥感影像的变化检测就是从不同时期的遥感数据中，定量地分析和确定地表变化的特征与过程。各国学者从不同的角度和应用研究提出了许多有效的检测算法，主要包括非监督变化检测和监督变化检测两大类算法。其中非监督变化检测算法比较直观，不需要先验信息，由于不需要训练样本，检测算法成本低。其中代表性的非监督变化检测算法如变化矢量分析法(Change Vector Analysis，CVA)、基于Fuzzy C-means(FCM)的聚类方法、基于水平集的方法等。其中，传统的基于FCM的多时相光学遥感变化检测，多先进行CVA变换，然后对变化矢量的幅值进行FCM聚类，进而得到变化检测结果。该类技术中，使用FCM的不足是仅适用于球状或椭球状聚类，且对噪声及其野值(Outlier)极为敏感。另外，仅仅使用变化矢量的幅值，使得原始多光谱信息没有得到充分的挖掘。

针对上述问题，许多学者试图通过在FCM目标函数中加上不同的空间邻域的约束来解决，但是高分辨率影像检测环境的复杂化以及目标先验信息匮乏等，导致这些算法都存在着一定的局限性。为此，有必要研究新的高分辨率可见光遥感图像变化检测技术来有效克服上述难点。

发明内容

本发明所要解决的技术问题是：提供基于非负矩阵分解和核FCM的多时相遥感影像变化检测方法，该方法解决了多时相高分辨率多光谱遥感影像变化检测精度不高的问题。

本发明为解决上述技术问题采用以下技术方案：

基于非负矩阵分解和核FCM的多时相遥感影像变化检测方法，包括如下步骤：

步骤1，输入两时相的高分辨率光学遥感影像，分别记为：X₁和X₂；

步骤2，对X₁和X₂进行影像配准；

步骤3，利用多元变化检测方法分别对影像配准后的X₁和X₂进行辐射归一化校正；

步骤4，对校正后的两时相影像分别计算变化矢量幅值和光谱角映射图；

步骤5，利用非负矩阵分解算法融合变化矢量幅值和光谱角映射图，得到一幅新的差异影像；

步骤6，利用核FCM结合空间邻域信息的方法对差异影像进行变化检测，得到最终的变化检测结果。

作为本发明的一种优选方案，步骤2所述影像配准包括几何粗校正和几何精校正两个步骤，其中，几何粗校正采用双线性差值法，几何精校正采用自动匹配与三角剖分法。

作为本发明的一种优选方案，步骤4所述变化矢量幅值的计算公式为：

$X_M(i,j)=\sqrt{\underset{b=1}{\overset{B}{\Sigma }}{\left(X_{1b}\right(i,j)-X_{2b}(i,j\left)\right)}^2}$

其中，B表示每一个时相遥感影像的波段数目，X_1b(i,j)、X_2b(i,j)分别表示X₁、X₂时相遥感影像b波段像素坐标(i,j)的像素值，X_M(i,j)表示像素坐标(i,j)的变化矢量幅值。

作为本发明的一种优选方案，步骤4所述光谱角的计算公式为：

$X_{SA}(i,j)=\arccos \left(\underset{b=1}{\overset{B}{\Sigma }}\right(X_{1b}\left(i,j\right)-X_{2b}\left(i,j\right))/\sqrt{\underset{b=1}{\overset{B}{\Sigma }}{X}_{1b}^2(i,j)}\sqrt{\underset{b=1}{\overset{B}{\Sigma }}{X}_{2b}^2(i,j)})$

其中，B表示每一个时相遥感影像的波段数目，X_1b(i,j)、X_2b(i,j)分别表示X₁、X₂时相遥感影像b波段像素坐标(i,j)的像素值，X_SA(i,j)表示像素坐标(i,j)的光谱角。

作为本发明的一种优选方案，步骤6所述核FCM结合空间邻域信息的模型为：

$J_m=\underset{t=1}{\overset{C}{\Sigma }}\underset{k=1}{\overset{N}{\Sigma }}{u}_{tk}^m(1-K(X_F\left(k\right),v_t\left)\right)+\alpha \underset{t=1}{\overset{C}{\Sigma }}\underset{k=1}{\overset{N}{\Sigma }}{u}_{tk}^m(1-K({\bar{X}}_F\left(k\right),v_t\left)\right)$

其中，C表示聚类数目，N表示样本的总数，表示第k样本X_F(k)对于第t类聚类中心v_t的模糊隶属度，m为隶属度的加权指数，α表示控制惩罚效果的参数，为差异影像X_F的局部均值影像，x,y表示任意变量，g表示核宽度。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：

1、本发明在基于核FCM的变化检测中，融合变化矢量幅值和光谱角映射图信息，优于仅仅使用变化矢量幅值的FCM方法。

2、本发明变化检测中，采用伪训练样本，并采用变化指数准则来自适应选择核FCM结合邻域空间信息的检测方法中的相关参数，可以使得变化检测的结果更加可靠，也更加具有稳健性。

附图说明

图1是本发明基于非负矩阵分解和核FCM的多时相遥感影像变化检测方法的实现流程示意图。

图2(a)是本发明所采用的2006年的SPOT5高分辨率图像三个波段合成的示意图。图2(b)是本发明所采用的2009年的SPOT5高分辨率图像三个波段合成的示意图。图2(c)是变化检测的参考图像。

图3(a)是CVA-EM算法检测结果图像。图3(b)是FCM-S算法的检测结果图像。图3(c)是本发明算法检测结果图像。

具体实施方式

下面详细描述本发明的实施方式，所述实施方式的示例在附图中示出。下面通过参考附图描述的实施方式是示例性的，仅用于解释本发明，而不能解释为对本发明的限制。

针对高空间分辨率遥感影像背景信息复杂、噪声干扰严重，使得变化检测所面临的问题采用常规的变化检测方法难以解决。本发明首先基于非负矩阵分解(Non-Negative Factorization，NMF)算法融合多时相遥感影像的变化矢量幅值(Magnitudes of Change Vectors，MCV)和多时相的光谱角映射图(Spectral Angle Mapper，SAM)，然后将融合结果作为核FCM的输入，再基于核FCM结合空间邻域信息的方法，获取最终的变化检测结果。

如图1所示，本发明的实现步骤如下：

步骤1，输入同一区域、不同时相的两幅高分辨率光学遥感影像，分别记为：X₁和X₂。

步骤2，利用ENVI遥感软件对X₁和X₂进行影像配准，分为几何粗校正和几何精校正两个步骤：

对于几何粗校正，利用ENVI4.8软件中的相关功能实现，具体操作步骤为：(1)选择X₁和X₂其中一个作为基准影像，另一个作为待校正影像；(2)采集地面控制点(GCPs)；GCPs应均匀分布在整幅图像内，GCPs的数目至少大于等于9；(3)计算误差；(4)选择多项式模型；(5)采用双线性插值进行重采样输出。

双线性差值法，若求未知函数f在点P＝(x,y)的值，假设我们已知函数f在Q₁₁＝(x₁,y₁)、Q₁₂＝(x₁,y₂)、Q₂₁＝(x₂,y₁)、及Q₂₂＝(x₂,y₂)四个点的值。如果选择一个坐标系统使得这四个点的坐标分别为(0,0)、(0,1)、(1,0)和(1,1)，那么双线性插值公式就可以表示为：

f(x,y)≈f(0,0)(1-x)(1-y)+f(1,0)x(1-y)+f(0,1)(1-x)y+f(1,1)xy (1)

对于几何精校正，将经过几何粗校正的多光谱遥感影像数据，利用自动匹配与三角剖分法进行几何精校正。

三角剖分法为，采用逐点插入法构建Delaunay三角网，对每一个三角形，利用其三个顶点的行列号与其对应的基准影像同名点的地理坐标来确定该三角形内部的仿射变换模型参数，对待校正影像进行纠正，得到校正后的遥感影。

步骤3，利用多元变化检测(Multivariate Alteration Detection，MAD)方法对X₁和X₂进行辐射归一化校正。该方法首先找到两期影像各波段亮度值的一个线性组合，得到变化信息增强的差异影像，通过阈值确定变化和未变化区域，然后通过未变化区域对应的两时相像元建立映射方程，完成相对辐射校正。

步骤4，对输入的多时相高分辨率影像分别进行变化矢量幅值X_M和光谱角信息X_SA的计算。

1)X_M的计算如下：

$X_M(i,j)=\sqrt{\underset{b=1}{\overset{B}{\Sigma }}{\left(X_{1b}\right(i,j)-X_{2b}(i,j\left)\right)}^2}---\left(2\right)$

其中，B表示每一个时相遥感影像的波段数目，(i,j)是影像的坐标。

2)X_SA的计算如下：

$X_{SA}(i,j)=\arccos \left(\underset{b=1}{\overset{B}{\Sigma }}\right(X_{1b}\left(i,j\right)-X_{2b}\left(i,j\right))/\sqrt{\underset{b=1}{\overset{B}{\Sigma }}{X}_{1b}^2(i,j)}\sqrt{\underset{b=1}{\overset{B}{\Sigma }}{X}_{2b}^2(i,j)})---\left(3\right)$

步骤5，利用NMF融合X_M和X_SA的信息，得到一幅新的差异影像X_F，具体过程如下：

1)将X_M和X_SA分别按列的顺序排成一列，并组成一个新的矩阵X，其维数为N×2，其中，N表示差异影像X_F的高和宽的乘积。

2)建立NMF的目标函数为：

$O=\frac{1}{2}\left|\right|X-AY|{|}_F^2---\left(4\right)$

式中，F表示Frobenius范数，A和Y是非负矩阵，A的维数是N×1，Y的维数是1×2。

3)采用交替最小二乘法求解式(4)的A和Y，迭代规则为：

式中，分别表示A和Y的广义逆，[x]₊＝max{ω,x}，ω表示任意大于0的常数，x表示任意变量。

4)将矩阵A重新排列为和原始差异影像X_M相同的维数，即为新的差异影像X_F。

步骤6，利用核FCM结合空间邻域信息的方法对X_F进行分割，具体过程如下：

1)首先构建核FCM结合空间邻域信息的模型如下：

$J_m=\underset{t=1}{\overset{C}{\Sigma }}\underset{k=1}{\overset{N}{\Sigma }}{u}_{tk}^m(1-K(X_F\left(k\right),v_t\left)\right)+\alpha \underset{t=1}{\overset{C}{\Sigma }}\underset{k=1}{\overset{N}{\Sigma }}{u}_{tk}^m(1-K({\bar{X}}_F\left(k\right),v_t\left)\right)---\left(7\right)$

式中，C是聚类数目，N是样本的总数，表示第k样本对于第t类聚类中心的模糊隶属度，m为隶属度的加权指数，u_tk∈[0,1]且参数α控制惩罚效果，为X_F的局部均值影像，x,y表示任意变量，g表示核宽度。

2)式(7)的目标函数最小化可以用下述的公式交替进行：

$u_{tk}=\frac{{\left(\right(1-K\left(X_F\left(k\right),v_t\right))+\alpha (1-K\left({\bar{X}}_F\left(k\right),v_t\right)\left)\right)}^{-1/(m-1)}}{\underset{t=1}{\overset{C}{\Sigma }}\left(\right(1-K\left(X_F\left(k\right),v_t\right)\left)\right)+\alpha \left(1-K\left({\bar{X}}_F\left(k\right),v_t\right)\right){)}^{-1/(m-1)}}---\left(8\right)$

$v_t=\frac{\underset{k=1}{\overset{N}{\Sigma }}{u}_{tk}^m\left(K\right(X_F\left(k\right),v_t\left)X_F\right(k)+\alpha K({\bar{X}}_F\left(k\right),v_t\left){\bar{X}}_F\right(k\left)\right)}{\underset{k=1}{\overset{N}{\Sigma }}{u}_{tk}^m\left(K\right(X_F\left(k\right),v_t)+\alpha K({\bar{X}}_F\left(k\right),v_t\left)\right)}---\left(9\right)$

步骤7，利用式(8)和(9)的交替迭代，获取最终的变化和未变化类的模糊隶属度，并根据隶属度的大小确定影像的变化区域和非变化区域。具体实现步骤如下：

1)设定C＝2，初始的未变化类和变化类的中心，选择和变化矢量幅值最小值和最大值相对应的矢量；设m＝2，ε为大于0的非常小的常量；

2)计算X_F的局部窗口均值，窗口大小设置为3×3；

3)采用式(8)更新模糊划分矩阵；

4)采用式(9)更新聚类中心；

5)重复3)和4)直到相邻两次迭代的聚类中心距离小于ε；

6)根据模糊划分矩阵u_tk确定最终的变化检测图。

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

如图2(a)和图2(b)所示，本发明的实验数据为法国Littoral地区的多时相SPOT高分辨影像数据，图像大小为400×400，使用B1、B2和B3三个波段。图2(c)是变化检测的参考图像。为了验证本发明的有效性，将本发明变化检测方法与下述变化检测方法进行比对：

如图3(a)所示，基于CVA的EM方法(CVA-EM)[意大利的Bruzzone L.等在文章“Automatic analysis of difference image for unsupervised change detection”(IEEE Transactions on Geoscience and Remote Sensing,2000,38(3):1171-1182)中所提的检测方法]。

如图3(b)所示，FCM结合空间邻域信息的分类方法(FCM-S)[Chen songchan等在文章“Robust Image Segmentation Using FCM With Spatial Constraints Based on New Kernel-Induced Distance Measure”(IEEE Transactions on Systems,Man,and Cybernetics-Part B:Cybernetics,2004,34(4):1907-1916)中所提的方法]。

如图3(c)所示，本发明方法的检测结果。

检测性能用错检数FP、漏检数FN、总错误数OE和Kappa系数四个指标来衡量。FP、FN和OE越接近于0、Kappa系数越接近于1，表明变化检测方法的性能越好。检测结果如表1所示。由表1可见，本发明所提的检测方法Kappa系数比其余两种算法要高，总体误差OE的数值也是最低的，这表明本发明技术性能优于其他两种检测方法，所提的变化检测方法是有效的。

表1 Littoral地区的多时相SPOT5影像变化检测结果比较

方法FPFNOEKappaCVA-EM79193882118010.705FCM-S1822692887500.737本发明方法2951436173120.797理想0001

以上实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。