CN103198482A - 基于差异图模糊隶属度融合的遥感图像变化检测方法 - Google Patents

Abstract

本发明公开一种基于差异图模糊隶属度融合的遥感图像变化检测方法，主要解决现有变化检测方法不能既有效去除伪变化信息又保持边缘信息的问题。其实现过程是：输入两幅不同时相的遥感图像，计算其对应像素点的结构相似度系数，得到一幅相似度差异图；对两幅遥感图像做差得到一幅差值图像；对差值图的像素进行类别标记得到一幅类别标记图；根据类别标记图对差值图进行滤波处理得到一幅去噪差值图；对相似度差异图和去噪差异图进行模糊隶属度融合并分类得到变化检测结果。本发明具有较强的抗噪性，能有效去除伪变化信息，同时保留较好的边缘信息，检测结果准确率高，可用于城区扩展监测、森林和植被变化监测。

Description

基于差异图模糊隶属度融合的遥感图像变化检测方法

技术领域

本发明属于数字图像处理领域，主要涉及遥感图像变化检测，具体是一种基于差异图模糊隶属度融合的遥感图像变化检测方法，可用于遥感图像分析和处理。

背景技术

遥感图像的变化检测是通过分析和提取同一地区不同时相的遥感图像之间存在的电磁波谱特征差异或空间结构特征差异，来识别物体的状态变化或现象变化的过程。在国民经济和国防建设的诸多领域已得到广泛应用，如农业调查、森林和植被变化监测、城区扩展监测、军事目标监测等。

常见的遥感图像变化检测方法的一般步骤是先构造差异图，然后选取适当的阈值将差异图分为变化类和非变化类。其中差异图的构造和对差值图的处理是图像变化检测的重要步骤。比较简单的差异图像构造方法有差值法，这种方法易于实现但是构造出的差异图噪声较多，需要有效的方法对差异图中的噪声进行处理。将差异图像中化较剧烈的噪声去除，将灰度值不大的变化类像素进行增强，能有效提高差异图的质量，使检测结果更准确。

为了构造较好的差异图，一些学者通过度量两时相图像对应像素灰度值的相似度来构造差异图，Inglada和Mercier（2007）在文章“A New Statistical Similarity Measurefor Change Detection in Multitemporal SAR Images and Its Extension to MultiscaleChange Analysis,IEEE Transactions on Geoscience and Remote Sensing,2007,45(5):1432-1445”中提出了一种基于统计相似度的SAR图像变化检测方法，该方法用KL散度衡量两时相图像对应像素邻域的统计相似度来构造差异图，然后阈值分割差异图得到变化结果，该方法是对局部直方图建模，但是局部区域的像素很少，很难有效的对其建模，所以该方法的检测结果较差。He（2010）在文章“Application of EuclideanNorm in Multi-temporal Remote Sensing Image Change Detection,International Congresson Image and Signal Processing(CISP’2010),2010,5:2111-2115”中提出了一种基于欧式距离的变化检测方法，该方法通过计算多个波段的两时相遥感图像的欧式距离来构造差异图，然后对差异图进行阈值分割得到变化结果，该方法能够有效地减少部分噪声的影响但检测结果仍存在较多伪变化信息。

为了合并不同差异图的优点，一些学者对不同的差异图进行了融合，马国锐等学者（2006）在文章“基于融合和广义高斯模型的遥感影像变化检测,遥感学报,2006,10(6):847-853”中提出了用乘积融合策略来融合差值图和比值图的方法，该方法能够抑制背景，在一定程度上增强变化区域，但是该方法并不稳定，有时也会抑制变化区域。汪闽和张星月（2010）在文章“多特征证据融合的遥感图像变化检测.遥感学报,2010,14(2):1-7”中提出了用证据理论融合方法融合多种特征差异图的方法，该方法提高了单一特征检测方法的检测精度，但该方法采用的结构相似度度量方法不稳定，从而降低了检测结果的正确率。Du等（2012）在文章“Fusion of Difference Images forChange Detection over Urban Areas,IEEE Journal of Selected Topics in Applied EarthObserbations and Remote Sensing,2012,5(4):1076-1086”中提出了对多种差异图进行特征级和决策级融合的变化检测方法，该方法合并了多种差异图的优点，提高了变化检测的正确率，但该方法并没有根据差异图的优缺点有针对性的融合差异图，所以检测结果的正确率提高不多，有时可能会下降。

发明内容

本发明的目的在于克服上述变换检测技术的不足，提出了一种基于差异图模糊隶属度融合的遥感图像变化检测方法，以降低噪声对检测结果的影响，减少检测结果中的伪变化信息，提高检测结果的正确率。

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

（1）输入的两幅大小均为I×J的同一地区不同时相的已配准的遥感图像X₁和X₂，计算这两幅遥感图像X₁和X₂对应像素点的结构相似度系数SIM(m,n)，得到一个相似度系数矩阵SIM：

SIM＝{SIM(m,n)|1≤m≤I,1≤n≤J}

其中, $\mathrm{SIM}(m,n)=\mathrm{\λ}\frac{2{\mathrm{\μ}}_1(m,n)\·{\mathrm{\μ}}_2(m,n)+C}{{\mathrm{\μ}}_1^2(m,n)+{\mathrm{\μ}}_2^2(m,n)+C}+(1-\mathrm{\λ})\frac{2{\mathrm{\σ}}_1(m,n)\·{\mathrm{\σ}}_2(m,n)+C}{{\mathrm{\σ}}_1^2(m,n)+{\mathrm{\σ}}_2^2(m,n)+C},$

式中，m和n分别为图像的行序号和列序号，m=1,2,...,I，n=1,2,...,J，μ₁(m,n)、σ₁(m,n)和

(m,n)分别为遥感图像X₁中以像素点(m,n)为中心、窗口大小为w×w的局部区域像素值的均值、标准差和方差，μ₂(m,n)、σ₂(m,n)和

(m,n)分别为遥感图像X₂中以像素点(m,n)为中心、窗口大小为w×w的局部区域像素值的均值、标准差和方差，窗口w的取值范围为3～9，C是用来避免分母接近零时产生不稳定现象的常数，C的取值范围为C>0，λ为权值系数，结构相似度系数的取值范围为0≤SIM(m,n)≤1；

（2）将相似度系数矩阵SIM在位置(m,n)处的结构相似度系数SIM(m,n)线性映射到区间[0,255]，得到相似度差异图X_S在位置(m,n)处的像素灰度值X_S(m,n)：

X_S(m,n)＝(1-SIM(m,n))×255，

（3）将遥感图像X₁和X₂空间对应位置(m,n)处的像素点灰度值X₁(m,n)和X₂(m,n)进行差值计算，得到差值X_d(m,n)=|X₁(m,n)-X₂(m,n)|，按照从左到右、从上到下的顺序依次计算遥感图像X₁和X₂空间对应位置的像素点灰度值的差值，得到一幅差值差异图X_d：

X_d＝{X_d(m,n)|1≤m≤I,1≤n≤J}，

（4）对差值差异图X_d中的像素进行窗口大小为3×3的中值滤波，得到滤波后差值图X_f，对滤波后差值图X_f进行统计直方图阈值分割，得到初始分类图X_m；

（5）根据滤波后差值图X_f的灰度值范围和初始分类图X_m，对差值差异图X_d中的像素点进行类别标记，得到一幅类别标记图X_b；

（6）根据类别标记图X_b中位置(m,n)处的标记，对差值差异图X_d中位置(m,n)处的像素进行滤波，得到去噪差值图X_N；

（7）计算相似度差异图X_S中位置(m,n)处像素的变化类隶属度

(m,n)和非变化类隶属度

(m,n)，计算去噪差值图X_N在位置(m,n)处像素的变化类隶属度(m,n)和非变化类隶属度

(m,n)；

（8）计算相似度差异图X_S和去噪差值图X_N的对应像素点(m,n)的变化类隶属度(m,n)和(m,n)的融合隶属度值H_c(m,n)，再计算该像素点的非变化类隶属度

(m,n)和

(m,n)的融合隶属度H_u(m,n)；

（9）建立一幅与相似度差异图X_S相同大小的融合图像X_I，如果H_c(m,n)＞H_u(m,n)，则将该像素点处的融合图像X_I值的标记为1，否则记为0，由此得到一幅变化检测结果图像。

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

1、本发明提出的结构相似度度量方法能稳定有效地度量两时相图像像素之间的相似性，构造的相似度差异图能够抑制背景噪声，增强目标与背景的对比度。

2、本发明对差值图中的像素进行类别标记，能够在保持边缘信息的同时有效抑制差值图中较剧烈的噪声。

3、本发明融合相似度差异图和去噪差值图，有效地合并了两种差异图的优点，不仅降低了噪声对检测结果的影响，还保持了变化区域的边缘，提高了变化检测结果的正确率。

附图说明

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

图2是用于实验的第一组遥感图像和对应的变化参考图像；

图3是用于实验的第二组遥感图像和对应的变化参考图像；

图4是本发明与对比方法仿真第一组遥感图像得到的变化检测结果图；

图5是本发明与对比方法仿真第二组遥感图像得到的变化检测结果图。

具体实施方式

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

步骤1，输入的两幅大小均为I×J的同一地区不同时相的已配准的遥感图像X₁和X₂，如图2(a)和图2(b)所示，计算这两幅遥感图像X₁和X₂对应像素点的结构相似度系数SIM(m,n)，得到一个相似度系数矩阵SIM。

1.1）分别计算时相1和时相2的相遥感图像X₁和X₂中以像素点(m,n)为中心、窗口大小为w×w像素的局部区域的均值μ₁(m,n)和μ₂(m,n)、标准差σ₁(m,n)和σ₂(m,n)及方差

(m,n)和

(m,n)，m和n分别为图像的行序号和列序号，m=1,2,…,I，n=1,2,…,J，为了在定位精度和抗噪能力方面进行折中，窗口w的取值范围为3～9，本发明实例中w=5；

1.2）计算时相1和时相2的遥感图像X₁和X₂对应像素点(m,n)处的结构相似度系数SIM(m,n)：

$\mathrm{SIM}(m,n)=\mathrm{\λ}\frac{2{\mathrm{\μ}}_1(m,n)\·{\mathrm{\μ}}_2(m,n)+C}{{\mathrm{\μ}}_1^2(m,n)+{\mathrm{\μ}}_2^2(m,n)+C}+(1-\mathrm{\λ})\frac{2{\mathrm{\σ}}_1(m,n)\·{\mathrm{\σ}}_2(m,n)+C}{{\mathrm{\σ}}_1^2(m,n)+{\mathrm{\σ}}_2^2(m,n)+C}$

式中，C是用来避免分母接近零时产生不稳定现象的常数，C的取值范围为C>0，本发明实例中取为1，λ为权值系数，本发明实例中λ=0.95，结构相似度系数计算所得范围0≤SIM(m,n)≤1，SIM(m,n)越接近于0，表明遥感图像X₁和X₂在像素点(m,n)处越不相似，属于变化类的可能性越大；

1.3）按照步骤（1.1）和步骤（1.2）计算图像中所有像素的结构相似度系数，得到一个相似度系数矩阵SIM＝{SIM(m,n)|1≤m≤I,1≤n≤J}。

步骤2，将相似度系数矩阵SIM在位置(m,n)处的结构相似度系数SIM(m,n)线性映射到区间[0,255]，得到相似度差异图X_S在位置(m,n)处的像素灰度值X_S(m,n)：

X_S(m,n)＝(1-SIM(m,n))×255。

步骤3，将时相1和时相2的遥感图像X₁和X₂空间对应位置(m,n)处的像素点灰度值X₁(m,n)和X₂(m,n)进行差值计算，得到差值X_d(m,n)=|X₁(m,n)-X₂(m,n)|，按照从左到右、从上到下的顺序依次计算遥感图像X₁和X₂空间对应位置的像素点灰度值的差值，得到一幅差值差异图X_d：

X_d＝{X_d(m,n)|1≤m≤I,1≤n≤J}。

步骤4，对差值差异图X_d进行窗口大小为3×3的中值滤波，得到滤波后差值图X_f，对滤波后差值图X_f进行统计直方图阈值分割，得到初始分类图X_m。

4.1）对差值差异图X_d进行窗口大小为3×3的中值滤波，得到滤波后差值图X_f；

4.2）对滤波后差值图X_f，将满足条件

的灰度级中的最小灰度级设置为低阈值T_l，其中，L表示进行灰度级像素个数累加的最大像素灰度级，其取值范围为0到255，N_x表示灰度级为x的像素的总个数，λ为一个比例常数，取值范围为0.5～0.6，本发明实例中λ=0.5；

4.3）对滤波后差值图X_f，将大于低阈值T_l且满足条件N_x＜(1-λ)×I×J/(255-T_l)的灰度级中的最小灰度级设定为统计直方图阈值

4.4）对滤波后差值图X_f以统计直方图阈值进行分割，得到初始分类图X_m：

$X_m(m,n)=\left\{\begin{array}{cc}1,& X_f(m,n)>=\stackrel{\&OverBar;}{T}\\ 0,& X_f(m,n)<\stackrel{\&OverBar;}{T}\end{array}\right.,$

其中，X_m(m,n)为初始分类图X_m中(m,n)处的像素值，X_f(m,n)为滤波后差值图X_f中(m,n)处的灰度值，标记1表示为变化类，标记0表示非变化类。

步骤5，根据滤波后差值图X_f的灰度值范围和初始分类图X_m，对差值差异图X_d中的像素点进行类别标记，得到一幅类别标记图X_b。

5.1）计算初始分类图X_m中变化类像素形成的各个区域的面积，如果面积小于70，则将初始分类图X_m中该区域内的像素标记赋值为2，否则像素标记不变；

5.2）计算滤波后差值图X_f的最大灰度值N_max，如果N_max大于阈值T，初始分类图X_m则为类别标记图X_b；否则，转到步骤（5.3），其中灰度级阈值T为一个常数，取值范围为100～150，本发明实例中T=100；

5.3）对初始分类图X_m中标记为1的像素形成的区域，进行结构元素为3×3的方形窗的数学形态学膨胀，得到小扩展图像X_m1，并小扩展图像X_m1中点(m,n)的像素值记为X_m1(m,n)；

5.4）对初始分类图X_m中标记为1的像素形成的区域，进行结构元素为7×7的方形窗的数学形态学膨胀，得到大扩展图像X_m2，并大扩展图像X_m2中点(m,n)的像素值记为X_m2(m,n)；

5.5）将大扩展图像X_m2和小扩展图像X_m1空间对应位置(m,n)处的像素点灰度值进行差值计算，得到差值X_m3(m,n)=X_m2(m,n)-X_m1(m,n)，由此得到一幅扩展差值图像X_m3={X_m3(m,n)}，其中X_m3(m,n)为扩展差值图像X_m3中点(m,n)处的像素值；

5.6）建立一幅与初始分类图像相同大小的类别标记图X_b，并按以下五种情况对该图像中点(m,n)处的像素值X_b(m,n)进行赋值：

对满足条件X_m(m,n)=2的像素点(m,n)，将类别标记图X_b中该像素点的值X_b(m,n)标记为2；

对满足条件X_m(m,n)≠2且X_m3(m,n)=1的像素点(m,n)，将类别标记图X_b中该像素点的值X_b(m,n)标记为3；

对满足条件X_m(m,n)≠2且X_m1(m,n)=1的像素点(m,n)，将类别标记图X_b中该像素点的值X_b(m,n)标记为1；

对满足条件X_m(m,n)=0的像素点(m,n)，将类别标记图X_b中该像素点的值X_b(m,n)标记为0；

对不满足以上四种条件的像素点(m,n)，将类别标记图X_b中该像素点的值X_b(m,n)标记为0。

步骤6，根据类别标记图X_b中位置(m,n)处的像素值，对差值差异图X_d中位置(m,n)处的像素进行滤波，得到一幅去噪差值图X_N。

滤波操作按如下规则进行：

对满足条件X_b(m,n)=0的像素点，计算差值差异图X_d中该像素点的9×9大小的窗口内所有像素的中值，将该中值赋给去噪差值图X_N中点(m,n)处的像素值X_N(m,n)；

对满足条件X_b(m,n)=1的像素点，将差值差异图X_d中该像素点的值X_d(m,n)赋给去噪差值图X_N中点(m,n)处的像素值X_N(m,n)；

对满足条件X_b(m,n)=2的像素点，计算差值差异图X_d中该像素点的11×11大小的窗口内所有像素的中值，将该中值赋给去噪差值图X_N中点(m,n)处的像素值X_N(m,n)；

如果存在满足条件X_b(m,n)=3的像素点，则计算差值差异图X_d中该像素点的高斯尺度为3×3的高斯核函数的值，将该值赋给去噪差值图X_N中点(m,n)处的像素值X_N(m,n)；

所有像素点(m,n)的值X_N(m,n)构成去噪差值图X_N={X_N(m,n)}。

步骤7，计算相似度差异图X_S中位置(m,n)处像素点的变化类隶属度(m,n)和非变化类隶属度(m,n)，再计算去噪差值图X_N在位置(m,n)处像素点的变化类隶属度

(m,n)和非变化类隶属度

(m,n)。

7.1）计算相似度差异图X_S中位置(m,n)处像素点的变化类隶属度

(m,n)和非变化类隶属度(m,n)：

$H_{{c,X}_S}(m,n)=\left\{\begin{array}{cc}0,& X_S(m,n)<T_1\\ \frac{1}{2}\&times;{\left(\frac{X_S(m,n)-T_1}{T_2-T_1}\right)}^2,& T_1\&le;X_S(m,n)<T_2\\ 1-\frac{1}{2}\&times;{\left(\frac{X_S(m,n)-T_3}{T_2-T_3}\right)}^2,& T_2\&le;X_S(m,n)<T_3\\ 1& X_S(m,n)\&GreaterEqual;T_3\end{array}\right.,$

$H_{{u,X}_S}(m,n)=1-H_{{c,X}_S}(m,n),$

其中，T₁、T₂和T₃是三个不同的阈值，且T₁<T₂<T₃，

T₂=(T₁+T₃)/2，其中，

为最大熵阈值分割方法分割相似度差异图X_S的最佳阈值；

7.2）计算去噪差值图X_N在位置(m,n)处像素点的变化类隶属度

(m,n)和非变化类隶属度

(m,n)：

$H_{{c,X}_N}(m,n)=\left\{\begin{array}{cc}0,& X_N(m,n)<T_1\\ \frac{1}{2}\&times;{\left(\frac{X_N(m,n)-T_1}{T_2-T_1}\right)}^2,& T_1\&le;X_N(m,n)<T_2\\ 1-\frac{1}{2}\&times;{\left(\frac{X_N(m,n)-T_3}{T_2-T_3}\right)}^2,& T_2\&le;X_N(m,n)<T_3\\ 1& X_N(m,n)\&GreaterEqual;T_3\end{array}\right.,$

$H_{{u,X}_N}(m,n)=1-H_{{c,X}_N}(m,n),$

其中，X_N(m,n)为去噪差值图X_N中(m,n)处的像素值。

步骤8，计算相似度差异图X_S和去噪差值图X_N的对应像素点(m,n)的变化类隶属度

(m,n)和

(m,n)的融合隶属度值H_c(m,n)，再计算该像素点的非变化类隶属度

(m,n)和(m,n)的融合隶属度H_u(m,n)。

8.1）计算相似度差异图X_S和去噪差值图X_N的对应像素点(m,n)的变化类隶属度(m,n)和

(m,n)的变化类的融合隶属度值H_c(m,n)：

$H_c(m,n)=\mathrm{\&beta;}\&times;H_{{c,X}_S}(m,n)+(1-\mathrm{\&beta;})\&times;H_{{c,X}_N}(m,n),$

其中，β为相似度差异图的变化类隶属度

(m,n)和非变化类隶属度(m,n)的融合权值，本发明实例中β=0.4。

8.2）计算像素点(m,n)的相似度差异图X_S的非变化类隶属度(m,n)和去噪差值图X_N的非变化类隶属度

(m,n)的融合隶属度H_u(m,n)：

$H_u(m,n)=\mathrm{\&beta;}\&times;H_{{u,X}_S}(m,n)+(1-\mathrm{\&beta;})\&times;H_{{u,X}_N}(m,n).$

步骤9，建立一幅与相似度差异图X_S相同大小的融合图像X_I，如果H_c(m,n)>H_u(m,n)，则将该处像素点(m,n)的融合图像X_I的值X_I(m,n)赋为1，否则赋为0；由此得到一幅变化检测结果图像X_I={X_I(m,n)}。

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

1.实验数据

图2(a)和图2(b)是第一组遥感数据，该组真实遥感数据是由2000年4月和2002年5月的墨西哥郊外的两幅Landsat7ETM+（Enhanced Thematic Mapper Plus）第4波段遥感图像组成。图2(a)和图2(b)大小均为512×512，256个灰度级，变化区域主要是大火破坏了大面积的当地植被所致，图2(c)为对应的变化参考图，图2(c)中白色区域表示变化区域，其中，变化的像元数为25599，非变化像元数为236545。

图3(a)和图3(b)是第二组遥感数据，该组真实遥感数据由1994年8月和1994年9月意大利Elba岛西部地区的两时相Landsat-5卫星TM第4波段光谱图像组成。图3(a)和图3(b)大小为326×414，256个灰度级，变化区域是由森林火灾引起的，图3(c)为对应的变化参考图，图3(c)中白色区域表示变化区域，其中，变化像素数为2415，非变化像素数为99985。

2.对比试验

为了说明本发明的有效性，本发明与如下三个对比方法进行了对比。

对比方法1，是差值差异图和相似度差异图模糊隶属度融合的方法。本发明与对比方法1对比是为了验证本发明根据标记图对差值图进行处理的有效性。

对比方法2，是将He（2010）在文章“Application of Euclidean Norm inMulti-temporal Remote Sensing Image Change Detection”提出的计算两时相遥感图像多个波段对应像素的欧式距离的方法换成计算两时相图单波段像素邻域的欧式距离来构造差异图的方法。本发明与对比方法2对比是为了验证本发明中构造的相似度差异图的有效性。

对比方法3，是Du等（2012）在文章“Fusion of Difference Images for ChangeDetection over Urban Areas”中提出的对多种差异图进行特征级融合的变化检测方法，其中特征级融合方法也采用模糊隶属度融合。本发明与对比方法3对比是为了验证本发明融合相似度差异图和去噪差值图进行变化检测的有效性。

3.实验内容与分析

实验1，用本发明和对比方法1、对比方法2和对比方法3对第一组遥感图像进行变化检测，结果如图4所示，其中：

图4(a)为对比方法1的变化检测结果图；

图4(b)为对比方法2的变化检测结果图；

图4(c)为对比方法3的变化检测结果图；

图4(d)为本发明的变化检测结果图。

从图4可以看出，三种对比方法的检测结果图中除了变化区域，整幅图还分布着面积很小的杂点；本发明与三种对比方法相比，结果图中变化区域的形状保持好，保留了较完整的边缘信息，且杂点等伪变化信息很少。

实验2，用本发明和对比方法1、对比方法2和对比方法3对第一组遥感图像进行变化检测，结果如图5所示，其中：

图5(a)是对比方法1的变化检测结果图；

图5(b)是对比方法2的变化检测结果图；

图5(c)是对比方法3的变化检测结果图；

图5(d)是本发明的变化检测结果图。

从图5可以看出，对比方法1的检测结果图中只检测出变化区域的一部分，漏检较多；对比方法2的检测结果图中的变化区域的整体形状保持较好，但结果图中有一些小面积的伪变化信息；对比方法3检测出的变化区域是一片白色区域，几乎覆盖了整幅图像，检测结果不能保持变化区域的基本形状，检测结果非常差；本发明与三种对比方法相比，检测结果图中变化区域的形状保持较好，且没有杂点等伪变化信息。

表1列出了本发明与三种对比方法仿真第一组遥感图像和第二组遥感图像的定量结果。

表1

从表1可以看出，对于第一组遥感图像，本发明的漏检数比对比方法1少，误检数比其稍多，总错误数少于对比方法1；本发明的漏检数比对比方法2多，但误检数比其少很多，总错误数少于对比方法2；本发明的漏检数比对比方法3稍多，但误检数远远少于对比方法3，总错误数少于对比方法3。对于第二组遥感图像，本发明的总错误数很少，正确率最高。

所以整体来看，本发明能够获得比较准确的变化信息，具有较强的抗噪性，能有效的去除杂点，同时能保持较好的变化区域的边缘信息；无论是视觉效果还是性能指标，本方法都较优。

Claims (4)

1.一种基于差异图模糊隶属度融合的遥感图像变化检测方法，其特征在于：包括如下步骤：

（1）输入的两幅大小均为I×J的同一地区不同时相的已配准的遥感图像X₁和X₂，计算这两幅遥感图像X₁和X₂对应像素点的结构相似度系数SIM(m,n)，得到一个相似度系数矩阵SIM：

SIM＝{SIM(m,n)|1≤m≤I,1≤n≤J}

其中, $\mathrm{SIM}(m,n)=\mathrm{\&lambda;}\frac{2{\mathrm{\&mu;}}_1(m,n)\&CenterDot;{\mathrm{\&mu;}}_2(m,n)+C}{{\mathrm{\&mu;}}_1^2(m,n)+{\mathrm{\&mu;}}_2^2(m,n)+C}+(1-\mathrm{\&lambda;})\frac{2{\mathrm{\&sigma;}}_1(m,n)\&CenterDot;{\mathrm{\&sigma;}}_2(m,n)+C}{{\mathrm{\&sigma;}}_1^2(m,n)+{\mathrm{\&sigma;}}_2^2(m,n)+C},$

式中，m和n分别为图像的行序号和列序号，m=1,2,...,I，n=1,2,...,J，μ₁(m,n)、σ₁(m,n)和

(m,n)分别为遥感图像X₁中以像素点(m,n)为中心、窗口大小为w×w的局部区域像素值的均值、标准差和方差，μ₂(m,n)、σ₂(m,n)和

(m,n)分别为遥感图像X₂中以像素点(m,n)为中心、窗口大小为w×w的局部区域像素值的均值、标准差和方差，窗口w的取值范围为3～9，C是用来避免分母接近零时产生不稳定现象的常数，C的取值范围为C>0，λ为权值系数，结构相似度系数的取值范围为0≤SIM(m,n)≤1；

（2）将相似度系数矩阵SIM在位置(m,n)处的结构相似度系数SIM(m,n)线性映射到区间[0,255]，得到相似度差异图X_S在位置(m,n)处的像素灰度值X_S(m,n)：

X_S(m,n)＝(1-SIM(m,n))×255

（3）将遥感图像X₁和X₂空间对应位置(m,n)处的像素点灰度值X₁(m,n)和X₂(m,n)进行差值计算，得到差值X_d(m,n)=|X₁(m,n)-X₂(m,n)|，按照从左到右、从上到下的顺序依次计算遥感图像X₁和X₂空间对应位置的像素点灰度值的差值，得到一幅差值差异图X_d：

X_d＝{X_d(m,n)|1≤m≤I,1≤n≤J}

（4）对差值差异图X_d中的像素进行窗口大小为3×3的中值滤波，得到滤波后差值图X_f，对滤波后差值图X_f进行统计直方图阈值分割，得到初始分类图X_m；

（5）根据滤波后差值图X_f的灰度值范围和初始分类图X_m，对差值差异图X_d中的像素点进行类别标记，得到一幅类别标记图X_b；

（6）根据类别标记图X_b中位置(m,n)处的标记，对差值差异图X_d中位置(m,n)处的像素进行滤波，得到去噪差值图X_N；

（7）计算相似度差异图X_S中位置(m,n)处像素的变化类隶属度

(m,n)和非变化类隶属度

(m,n)，计算去噪差值图X_N在位置(m,n)处像素的变化类隶属度

(m,n)和非变化类隶属度

(m,n)；

（8）计算相似度差异图X_S和去噪差值图X_N的对应像素点(m,n)的变化类隶属度(m,n)和

(m,n)的融合隶属度值H_c(m,n)，再计算该像素点的非变化类隶属度

(m,n)和

(m,n)的融合隶属度H_u(m,n)；

（9）建立一幅与相似度差异图XS相同大小的融合图像XI，如果H_c(m,n)＞H_u(m,n)，则将该像素点处的融合图像X_I值的标记为1，否则记为0，由此得到一幅变化检测结果图像。

2.根据权利要求1所述的遥感图像变化检测方法，其特征在于：步骤（4）所述的对滤波后差值图X_f进行统计直方图阈值分割，得到初始分类图X_m，按如下步骤进行：

（4a）对差值差异图X_d进行窗口大小为3×3的中值滤波，得到滤波后差值图X_f；

（4b）对滤波后差值图X_f，将满足条件

的灰度级中的最小灰度级设置为低阈值T_l，其中，L表示进行灰度级像素个数累加的最大素灰度级，其取值范围为0到255，N_x表示灰度级为x的像素的总个数，λ为一个比例常数，取值范围为0.5～0.6，I和J分别表示图像的行数和列数；

（4c）对滤波后差值图X_f，将大于低阈值T_l且满足条件N_x＜(1-λ)×I×J(25-5T_l)的灰度级中的最小灰度级设定为统计直方图阈值

；

（4d）对滤波后差值图X_f以统计直方图阈值

进行分割，得到初始分类图X_m：

$X_m(m,n)=\left\{\begin{array}{cc}1,& X_f(m,n)>=\stackrel{\&OverBar;}{T}\\ 0,& X_f(m,n)<\stackrel{\&OverBar;}{T}\end{array}\right.$

其中，X_m(m,n)为初始分类图X_m中(m,n)处的像素值，X_f(m,n)为滤波后差值图X_f中(m,n)处的灰度值，标记1表示为变化类，标记0表示非变化类。

3.根据权利要求1所述的遥感图像变化检测方法，其特征在于：步骤（5）所述的根据滤波后差值图X_f的灰度值范围和初始分类图X_m，对差值差异图X_d中的像素点进行类别标记，得到一幅类别标记图X_b，按如下步骤进行：

（5a）计算初始分类图X_m中变化类像素形成的各个区域的面积，如果面积小于70，则将初始分类图X_m中该区域内的像素标记赋值为2，否则像素标记不变；

（5b）计算滤波后差值图X_f的最大灰度值N_max，如果N_max大于灰度级阈值T，初始分类图X_m则为类别标记图X_b，结束；否则，转到（5c），其中灰度级阈值T为一个常数，取值范围为100～150；

（5c）对初始分类图X_m中标记为1的像素形成的区域，进行结构元素为3×3的方形窗的数学形态学膨胀，得到小扩展图像X_m1，并将小扩展图像X_m1中点(m,n)的像素值记为X_m1(m,n)；

（5d）对初始分类图X_m中标记为1的像素形成的区域，进行结构元素为7×7的方形窗的数学形态学膨胀，得到大扩展图像X_m2，并将大扩展图像X_m2中点(m,n)的像素值记为X_m2(m,n)；

（5e）将大扩展图像X_m2和小扩展图像X_m1空间对应位置(m,n)处的像素点灰度值进行差值计算，得到差值X_m3(m,n)=X_m2(m,n)-X_m1(m,n)，由此得到一幅扩展差值图像X_m3={X_m3(m,n)}，其中X_m3(m,n)为扩展差值图像X_m3中点(m,n)处的像素值；

（5f）建立一幅与初始分类图像相同大小的类别标记图X_b，并按以下五种情况对该图像中点(m,n)处的像素值X_b(m,n)进行赋值：

对满足条件X_m(m,n)=2的像素点(m,n)，将类别标记图X_b中该像素点的值X_b(m,n)标记为2；

对满足条件X_m(m,n)≠2且X_m3(m,n)=1的像素点(m,n)，将类别标记图X_b中该像素点的值X_b(m,n)标记为3；

对满足条件X_m(m,n)≠2且X_m1(m,n)=1的像素点(m,n)，将类别标记图X_b中该像素点的值X_b(m,n)标记为1；

对满足条件X_m(m,n)=0的像素点(m,n)，将类别标记图X_b中该像素点的值X_b(m,n)标记为0；

对不满足以上四种条件的像素点(m,n)，将类别标记图X_b中该像素点的值X_b(m,n)标记为0。

4.根据权利要求1所述的遥感图像变化检测方法，其特征在于：步骤（6）所述的根据类别标记图X_b中位置(m,n)处的像素值，对差值差异图X_d中位置(m,n)处的像素进行滤波，得到去噪差值图X_N，按如下规则进行：

对满足条件X_b(m,n)=0的像素点，计算差值差异图X_d中该像素点的9×9大小的窗口内所有像素的中值，将该中值赋给去噪差值图X_N中点(m,n)处的像素值X_N(m,n)；

对满足条件X_b(m,n)=1的像素点，将差值差异图X_d中该像素点的值X_d(m,n)赋给去噪差值图X_N中点(m,n)处的像素值X_N(m,n)；

对满足条件X_b(m,n)=2的像素点，计算差值差异图X_d中该像素点的11×11大小的窗口内所有像素的中值，将该中值赋给去噪差值图X_N中点(m,n)处的像素值X_N(m,n)；

如果存在满足条件X_b(m,n)=3的像素点，则计算差值差异图X_d中该像素点的高斯尺度为3×3的高斯核函数的值，将该值赋给去噪差值图X_N中点(m,n)处的像素值X_N(m,n)；

所有像素点(m,n)的值X_N(m,n)构成去噪差值图X_N={X_N(m,n)}。

Images