CN106384352A - 基于融合策略和fcm的多时相遥感影像变化检测方法 - Google Patents

Abstract

本发明公开了基于融合策略和FCM的多时相遥感影像变化检测方法，该方法首先联合多时相遥感影像的变化矢量幅值和多时相的光谱角映射图将检测范围分为确定区域和非确定区域。在非确定区域，基于模糊划分矩阵相加融合变化矢量幅值和SAM的信息。最后结合确定区域和非确定区域的结果，获取最终的变化检测结果。其中，非确定区域的FCM目标函数中的模糊指数，通过基于MCV和SAM在非确定区域的冲突指数来选择，以获得更加稳健和精度较高的变化检测结果。本发明将变化检测区域分为确定和非确定区域两个部分，采用融合的策略分别获取两个区域的检测结果，可以使得变化检测的结果更加可靠，也更加具有稳健性。

Description

基于融合策略和FCM的多时相遥感影像变化检测方法

技术领域

本发明涉及多时相光学遥感影像变化检测方法，具体涉及基于融合策略和FCM的多时相遥感影像变化检测方法，属于遥感影像处理技术领域。

背景技术

随着多时相遥感数据的不断积累以及空间数据库的相继建立，如何从这些遥感数据中提取和检测变化信息已成为遥感科学和地理信息科学的重要研究课题。根据同一区域不同时相的遥感影像，可以提取城市、环境等动态变化的信息，为资源管理与规划、环境保护等部门提供科学决策的依据。

遥感影像的变化检测就是从不同时期的遥感数据中，定量地分析和确定地表变化的特征与过程。各国学者从不同的角度和应用研究提出了许多有效的检测算法，如变化矢量分析法(Change Vector Analysis，CVA)、基于Fuzzy C-means(FCM)的聚类方法等。其中，传统的基于FCM的多时相光学遥感变化检测，多先进行CVA变换，然后对变化矢量的幅值进行FCM聚类，进而得到变化检测结果。

该类技术中，使用FCM的不足是仅仅使用变化矢量的幅值，使得原始多波段信息没有得到充分的挖掘。针对上述问题，许多学者试图通过在FCM目标函数中加上不同的空间邻域的约束来解决，但是遥感影像检测环境的复杂化以及目标先验信息匮乏等，导致这些算法都存在着一定的局限性。为此，有必要研究新的多时相可见光遥感图像变化检测技术来有效克服上述难点。

发明内容

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

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

基于融合策略和FCM的多时相遥感影像变化检测方法，包括如下步骤：

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

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

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

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

步骤5，针对变化矢量幅值，利用Bayes原理，并基于最大期望算法获取最优阈值T_M；针对光谱角映射图，采用Otsu阈值法获取阈值T_S；根据T_M和T_S构建确定区域和非确定区域；

步骤6，在非确定区域，利用FCM算法对变化矢量幅值和光谱角映射图分别进行变化检测，并对检测结果进行融合，获取非确定区域的检测结果；

步骤7，利用步骤6的检测结果，并根据确定区域的标记，确定最终的变化检测结果。

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

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

$X_M(i,j)=\sqrt{\underset{b=1}{\overset{B}{\Σ}}{\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_S(i,j)=\arccos \left(\underset{b=1}{\overset{B}{\Σ}}\left(X_{1b}\left(i,j\right)X_{2b}\left(i,j\right)\right)/\sqrt{\underset{b=1}{\overset{B}{\Σ}}{X}_{1b}^2(i,j)}\sqrt{\underset{b=1}{\overset{B}{\Σ}}{X}_{2b}^2(i,j)}\right)$

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

作为本发明的一种优选方案，步骤5所述最优阈值T_M的计算公式为：

$({\mathrm{\σ}}_n^2-{\mathrm{\σ}}_c^2){T_M}^2+2(m_n{\mathrm{\σ}}_c^2-m_c{\mathrm{\σ}}_n^2)T_M+m_c^2{\mathrm{\σ}}_n^2-m_n^2{\mathrm{\σ}}_c^2-2{\mathrm{\σ}}_n^2{\mathrm{\σ}}_c^{2}ln\[\frac{{\mathrm{\σ}}_{c}p\left({\mathrm{\ω}}_n\right)}{{\mathrm{\σ}}_{n}p\left({\mathrm{\ω}}_c\right)}\]=0$

其中，m_n和σ_n分别表示变化矢量幅值影像上未变化类ω_n的均值和方差，m_c和σ_c分别表示变化矢量幅值影像上变化类ω_c的均值和方差，p表示先验概率。

作为本发明的一种优选方案，步骤6所述FCM算法的模型为：

$J_m=\underset{c=1}{\overset{C}{\Σ}}\underset{k=1}{\overset{N}{\Σ}}{u}_{ck}^m\left|\right|X_M\left(k\right)-v_c|{|}^2=\underset{c=1}{\overset{C}{\Σ}}\underset{k=1}{\overset{N}{\Σ}}{u}_{ck}^m{d}_{kc}^2$

其中，C表示聚类数目，N表示样本的总数，表示第k样本X_M(k)对于第c类聚类中心v_c的模糊隶属度，m表示隶属度的加权指数，d_kc表示样本X_M(k)和v_c的距离。

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

1、本发明在基于核FCM的变化检测中，联合变化矢量幅值和光谱角映射图作为输入，利用这两个特征，优于仅仅使用变化矢量幅值的FCM方法。

2、本发明将变化检测区域分为确定和非确定区域两个部分，采用融合的策略分别获取两个区域的检测结果，可以使得变化检测的结果更加可靠，也更加具有稳健性。

附图说明

图1是本发明基于融合策略和FCM的多时相遥感影像变化检测方法的实现流程示意图。

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

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

具体实施方式

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

针对多时相多光谱遥感影像背景信息复杂、噪声干扰严重，使得变化检测所面临的问题采用常规的变化检测方法难以解决。本发明首先联合多时相遥感影像的变化矢量幅值(Magnitudes of Change Vectors，MCV)和多时相的光谱角映射图(Spectral AngleMapper，SAM)将检测范围分为确定区域和非确定区域。在非确定区域，基于模糊划分矩阵相加融合变化矢量幅值和SAM的信息。最后结合确定区域和非确定区域的结果，获取最终的变化检测结果。其中，非确定区域的FCM目标函数中的模糊指数，通过基于MCV和SAM在非确定区域的冲突指数来选择，以获得更加稳健和精度较高的变化检测结果。

如图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_S的计算，将X_M和X_S组合作为核FCM的输入。

1)X_M的计算如下：

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

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

2)X_S的计算如下：

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

步骤5，针对X_M，利用Bayes原理，并基于最大期望(Expectation-Maximization，EM)算法获取最优分割阈值T_M。最优分割阈值T_M的计算过程如下：

1)假设X_M影像上未变化类ω_n和变化类ω_c服从如下的高斯分布，即：

$p\left(X_M\right|{\mathrm{\ω}}_l)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_l}\exp \left\{-\frac{{(X_M-m_l)}^2}{2{\mathrm{\σ}}_l^2}\right\},{\mathrm{\ω}}_l\∈\{{\mathrm{\ω}}_n,{\mathrm{\ω}}_c\}---\left(4\right)$

其中，m_n和σ_n分别表示未变化类的均值和方差，m_c和σ_c分别表示变化类的均值和方差。

2)采用EM算法可估计m_n、σ_n、m_c和σ_c这四个参数，下面仅以未变化类的参数估计为例进行说明，变化类参数估计类似。

$p^{t+1}\left({\mathrm{\ω}}_n\right)=\frac{\underset{X_M(i,j)\∈X_M}{\Σ}\frac{p^t\left({\mathrm{\ω}}_n\right)p^t\left(X_M\right(i,j\left|{\mathrm{\ω}}_n\right))}{p^t\left(X_M\right(i,j))}}{IJ}---\left(5\right)$

$m_n^{t+1}=\frac{\underset{X_M(i,j)\∈X_M}{\Σ}\frac{p^t\left({\mathrm{\ω}}_n\right)p^t\left(X_M\right(i,j\left|{\mathrm{\ω}}_n\right))}{p^t\left(X_M\right(i,j))}{X}_M(i,j)}{\underset{X_M(i,j)\∈X_M}{\Σ}\frac{p^t\left({\mathrm{\ω}}_n\right)p^t\left(X_M\right(i,j\left|{\mathrm{\ω}}_n\right))}{p^t\left(X_M\right(i,j))}}---\left(6\right)$

${\left({\mathrm{\σ}}_n^2\right)}^{t+1}=\frac{\underset{X_M(i,j)\∈X_M}{\Σ}\frac{p^t\left({\mathrm{\ω}}_n\right)p^t\left(X_M\right(i,j\left|{\mathrm{\ω}}_n\right))}{p^t\left(X_M\right(i,j))}\[X_M(i,j)-m_n^t\]}{\underset{X_M(i,j)\∈X_M}{\Σ}\frac{p^t\left({\mathrm{\ω}}_n\right)p^t\left(X_M\right(i,j\left|{\mathrm{\ω}}_n\right))}{p^t\left(X_M\right(i,j))}}---\left(7\right)$

其中，I和J分别表示影像的行数和列数，t表示迭代次数。

3)依据Bayes最小误差准则，求解变化矢量幅值图像X_M的分割阈值T_M，在高斯分布的情况下，等同于求解下式：

$({\mathrm{\σ}}_n^2-{\mathrm{\σ}}_c^2){T_M}^2+2(m_n{\mathrm{\σ}}_c^2-m_c{\mathrm{\σ}}_n^2)T_M+m_c^2{\mathrm{\σ}}_n^2-m_n^2{\mathrm{\σ}}_c^2-2{\mathrm{\σ}}_n^2{\mathrm{\σ}}_c^{2}ln\[\frac{{\mathrm{\σ}}_{c}p\left({\mathrm{\ω}}_n\right)}{{\mathrm{\σ}}_{n}p\left({\mathrm{\ω}}_c\right)}\]=0---\left(8\right)$

式中，p(ω_n)是未变化类的先验概率，p(ω_c)是变化类的先验概率。

4)针对X_S，利用Otsu阈值法确定阈值T_S，计算过程如下：

4a)将X_S转化为8比特灰度图像；

4b)对每一灰度级，将其视为阈值，将影像分为变化和非变化区域，并计算此阈值下的w₀、u₀、w₁、u₁，其中，w₀表示变化区域的像素点数占图像的比例，u₀表示变化区域的平均灰度，w₁表示非变化区域的像素点数占图像的比例，u₁表示非变化区域的平均灰度；

4c)计算图像的总平均灰度：u＝w₀×u₀+w₁×u₁；

4d)计算变化和非变化区域的方差如下：

g＝w₀×(u₀-u)²+w₁×(u₁-u)² (9)

4e)遍历影像的每一个灰度级，重复步骤4b)-步骤4d)；

4f)找出和最大g对应的灰度级，即为X_S的Otsu阈值T_S。

5)根据EM算法所估计的阈值T_M和Otsu算法估计的阈值T_S，构建确定和非确定两个区域如下：

其中，确定区域包括两部分非变化区域和变化区域

非确定区域为式(10)和式(11)的补集。δ与X_M的动态范围相关，动态范围是指最大和最小值之间的差。一般取δ＝15％X_M的动态范围。

步骤6，在非确定区域，利用FCM算法分别对X_M、X_S进行分割(将变化检测视为分割)，以对X_M的分割为例，此时式(12)-(14)中的m＝m₁(X_S分割类似，此时式(12)-(14)中的m＝m₂)，具体过程如下：

1)首先构建FCM的模型如下：

$J_m=\underset{c=1}{\overset{C}{\Σ}}\underset{k=1}{\overset{N}{\Σ}}{u}_{ck}^m\left|\right|X_M\left(k\right)-v_c|{|}^2=\underset{c=1}{\overset{C}{\Σ}}\underset{k=1}{\overset{N}{\Σ}}{u}_{ck}^m{d}_{kc}^2---\left(12\right)$

式中，C是聚类数目，N是样本的总数，表示第k样本对于第c类聚类中心v_c的模糊隶属度，m为隶属度的加权指数，u_ck∈[0,1]且

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

$u_{ck}=\frac{1}{\underset{r=1}{\overset{C}{\Σ}}{\left(\frac{d_{kc}}{d_{kr}}\right)}^{2/(m-1)}}---\left(13\right)$

$v_c=\frac{\underset{k=1}{\overset{N}{\Σ}}{u}_{ck}^m{X}_M\left(k\right)}{\underset{k=1}{\overset{N}{\Σ}}{u}_{ck}^m}---\left(14\right)$

3)设定模糊指数m的取值范围，并以如下的准则搜索最优的参数m。

首先，定义非确定区域的冲突指数为：

$CI=\frac{n_1+n_2}{N_u}---\left(15\right)$

式中，N_u是非确定区域的像素总数，n₁和n₂的定义如下：

n₁为依据FCM算法对X_M判别为变化而对X_S判别为非变化的像素点总数，n₂为依据FCM算法对X_M判别为非变化而对X_S判别为变化的像素点总数。

然后，搜索模糊指数的组合m₁(对应于X_M)和m₂(对应于X_S)，最优的组合为使得式(15)取值最小。

步骤7，利用步骤6获取的最优参数，对X_M和X_S，分别利用式(13)和(14)的交替迭代，获取最终的变化和未变化类的模糊划分矩阵，并根据两类信息的模糊划分矩阵相加得到新的模糊划分矩阵，并根据隶属度u_ck的大小确定影像的变化区域和非变化区域。

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

如图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结合Markov随机场模型的分类方法(FCM-MRF)[Maoguo Gong等在文章“Fuzzy clustering with a modified mrf energy function for changedetection in synthetic aperture radar images”(IEEE Transactions on FuzzySystems,2014,22(1):98-109)中所提的方法]。

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

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

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

方法	FP	FN	OE	Kappa
					CVA-EM	7919	3882	11801	0.705
FCM-MRF	3419	6057	9476	0.730
					本发明方法	2255	6558	8813	0.739
理想	0	0	0	1

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

Claims (6)

1.基于融合策略和FCM的多时相遥感影像变化检测方法，其特征在于，包括如下步骤：

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

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

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

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

步骤5，针对变化矢量幅值，利用Bayes原理，并基于最大期望算法获取最优阈值T_M；针对光谱角映射图，采用Otsu阈值法获取阈值T_S；根据T_M和T_S构建确定区域和非确定区域；

步骤6，在非确定区域，利用FCM算法对变化矢量幅值和光谱角映射图分别进行变化检测，并对检测结果进行融合，获取非确定区域的检测结果；

步骤7，利用步骤6的检测结果，并根据确定区域的标记，确定最终的变化检测结果。

2.根据权利要求1所述基于融合策略和FCM的多时相遥感影像变化检测方法，其特征在于，步骤2所述影像配准包括几何粗校正和几何精校正两个步骤，其中，几何粗校正采用双线性差值法，几何精校正采用自动匹配与三角剖分法。

3.根据权利要求1所述基于融合策略和FCM的多时相遥感影像变化检测方法，其特征在于，步骤4所述变化矢量幅值的计算公式为：

$X_M(i,j)=\sqrt{\underset{b=1}{\overset{B}{\Σ}}{\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.根据权利要求1所述基于融合策略和FCM的多时相遥感影像变化检测方法，其特征在于，步骤4所述光谱角的计算公式为：

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

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

5.根据权利要求1所述基于融合策略和FCM的多时相遥感影像变化检测方法，其特征在于，步骤5所述最优阈值T_M的计算公式为：

$({\mathrm{\σ}}_n^2-{\mathrm{\σ}}_c^2){T_M}^2+2(m_n{\mathrm{\σ}}_c^2-m_c{\mathrm{\σ}}_n^2)T_M+m_c^2{\mathrm{\σ}}_n^2-m_n^2{\mathrm{\σ}}_c^2-2{\mathrm{\σ}}_n^2{\mathrm{\σ}}_c^{2}ln\[\frac{{\mathrm{\σ}}_{c}p\left({\mathrm{\ω}}_n\right)}{{\mathrm{\σ}}_{n}p\left({\mathrm{\ω}}_c\right)}\]=0$

其中，m_n和σ_n分别表示变化矢量幅值影像上未变化类ω_n的均值和方差，m_c和σ_c分别表示变化矢量幅值影像上变化类ω_c的均值和方差，p表示先验概率。

6.根据权利要求1所述基于融合策略和FCM的多时相遥感影像变化检测方法，其特征在于，步骤6所述FCM算法的模型为：

$J_m=\underset{c=1}{\overset{C}{\Σ}}\underset{k=1}{\overset{N}{\Σ}}{u}_{ck}^m\left|\right|X_M\left(k\right)-v_c|{|}^2=\underset{c=1}{\overset{C}{\Σ}}\underset{k=1}{\overset{N}{\Σ}}{u}_{ck}^m{d}_{kc}^2$

其中，C表示聚类数目，N表示样本的总数，表示第k样本X_M(k)对于第c类聚类中心v_c的模糊隶属度，m表示隶属度的加权指数，d_kc表示样本X_M(k)和v_c的距离。