CN101694719A - 基于非参数密度估计的遥感图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于非参数密度估计的遥感图像变化检测方法，主要解决现有技术对差异图像中与变化类和非变化类相关的统计项的估计存在偏差的问题。其实现过程是：输入两幅不同时相的遥感图像，并对每幅图像的每个通道去噪，得到两时相的去噪后图像，并采用变化矢量分析法构造差异影像；应用K-均值聚类算法将差异图像聚成变化类和非变化类，得到初始分类结果，并采用非参数密度估计的方法估计差异影像中与变化类和非变化类相关的统计项；结合变权马尔科夫随机场模型进行自适应的空间约束，得到最终的变化检测结果。实验表明本发明能够有效地保持图像的结构信息，并去除孤立噪声，提高变化检测处理效率，可用于灾情监测、土地利用、农业调查领域。

Description

基于非参数密度估计的遥感图像变化检测方法

技术领域

本发明属于数字图像处理技术领域，涉及多时相遥感图像的变化检测，具体的说是一种基于非参数密度估计的遥感图像变化检测。

背景技术

变化检测技术是指通过分析在同一地区但不同时间获得的两幅图像来辨识其变化信息。随着遥感图像获取技术和手段的日益先进以及遥感图像数据的海量积累，变化检测技术在环境监测、土地利用/覆盖、森林/植被变化分析、灾情监测、农业调查、城市变化分析、军事侦察和打击效果评估等方面的应用越来越广泛。

在已发表的文献中，基于非监督的变化检测技术主要基于以下3个步骤：1)图像的预处理，包括辐射校正、几何配准、图像去噪等；2)差异图像的构建，具体指的是将两幅图像进行逐个像素的比较；3)变化区域的提取，主要包括阈值法和分类法，其中基于MRF(Markov Random Fields)模型的分类方法，由于顾及了上下文关系，有较强的抗噪性，得到了一些学者的关注。

Bruzzone和Prieto(2000)在文章“Automatic analysis of the difference image forunsupervised change detection”中提出了基于Bayes理论和MRF模型的非监督变化检测方法，假设差异图像中与变化类和非变化类相关的统计项符合高斯混合模型(GMM，Gaussian Mixture Models)，并采用期望最大化(EM，ExpectationMaximum)算法来估计其模型参数，最后分别采用贝叶斯最小错误率阈值和MRF对差异图像进行分类。作为该方法的进一步改进，2002年Bruzzone和Prieto又在文章“An adaptive semiparametric and context-based approach to unsupervisedchange detection in multitemporal remote-sensing images”采用了简化的Parzen估计和EM算法来估计差异图像中与变化类和非变化类像素灰度级相关的统计项，但由于统计项的非监督估计和MRF空间正则化的过程是分离的，所以变化检测处理效率低。

江利明，廖明生(2006)等在文章“顾及空间邻域关系的多时相SAR图像变化检测”提出基于EM-MPM模型的变化检测方法，并和双阈值EM算法进行了比较，有效地提高了变化区域提取的可靠性和准确性。

孙强(2007)在其博士论文“基于统计模型的SAR图像处理与解译”中提出了一种基于广义高斯混合模型的SAR图像变化检测方法。在GGM(General GaussMixture，GGM)的先验下，通过基于Gibbs采样估计方法的模型统计推断，对两幅相关SAR图像的对数比图像进行最大似然分类，并在此基础上基于MRF进行自适应的空间约束，完成检测结果的更新。

宋妍，袁修孝(2009)等在文章“基于混合高斯密度模型和空间上下文信息的遥感图像变化检测方法及扩展”中提出了一种运用遗传K均值算法与EM算法联合解算高斯混合密度模型参数的方法，该方法可以自动地解算出模型的统计参数；然后，比较概率松弛迭代法以及MRF模型法的图像变化检测效果；最后，对传统的基于模拟退火法的MRF方法进行改进，提出了一种变权MRF方法，检测结果能更好地保持图像的结构性，并有效地去除了孤立噪声。

以上方法假设了差异图像中与变化类和非变化类相关的统计项符合具体的模型如高斯混合模型、广义高斯混合模型等，需要进行复杂的参数估计过程，且参数估计的精确程度会影响变化检测的结果，而实际中差异图像的统计项不一定符合这些具体的模型，使得这些方法对差异图像中与变化类和非变化类相关的统计项的估计存在偏差，进而影响变化检测精度。

发明内容

本发明的目的在于克服上述已有的遥感图像变化检测技术的不足，提出一种基于非参数密度估计的遥感图像变化检测，以减小差异图像中与变化类和非变化类相关统计项的估计偏差，提高变化检测精度。

实现本发明目的的技术方案是采用非参数密度估计方法估计差异图像与变化类和非变化类相关的统计项，并结合宋妍，袁修孝(2009)等人提出的变权马尔科夫随机场(Markov Random Fields)进行自适应的空间约束，对遥感图像的变化进行检测，其实现步骤包括如下：

(1)输入两幅不同时相的遥感图像，并对每幅图像的每个通道分别进行窗口大小为3×3像素的中值滤波，得到两时相的去噪后图像；

(2)将去噪后的两幅图像应用变化矢量分析，得到一幅差异图像，并根据该差异图像计算变权马尔科夫随机场的权值因子W；

(3)应用K-均值聚类算法将差异图像聚成变化类和非变化类，得到初始分类结果；

(4)利用初始分类结果，采用非参数密度估计方法估计差异图像中变化类和非变化类的类条件概率密度，再对该类条件概率密度取负自然对数得到变化类和非变化类的似然能量；

(5)对初始分类结果利用马尔科夫随机场计算变化类和非变化类的先验能量；

(6)利用权值因子W、变化类和非变化类的似然能量及变化类和非变化类的先验能量计算变化类的总能量和非变化类的总能量，将总能量较小的那一类作为当前类别，得到类别更新后的结果；

(7)对类别更新后的结果，采用非参数密度估计方法重新估计差异图像中变化类和非变化类的类条件概率密度，再对该类条件概率密度取负自然对数得到变化类和非变化类的似然能量，并利用马尔科夫随机场重新计算变化类和非变化类的先验能量；

(8)重复步骤(6)及步骤(7)直至迭代终止，并存储每次类别更新后的结果，得到每个像素点的类别更新集合，该迭代终止条件为迭代次数不超过50次及两次迭代之间相异的像元数目比例小于给定的阈值；

(9)利用每个像素点的类别更新集合估计变化类的后验概率和非变化类的后验概率，将后验概率较大的那一类作为该像素点的最终变化检测结果。

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

(1)本发明由于采用非参数密度估计方法估计差异图像的类条件概率密度，克服了现有技术采用高斯混合模型和广义高斯混合模型假设的缺陷，不需要事先对遥感影像的类条件概率密度做出假设，能够得到精确的估计结果。

(2)本发明由于结合了变权马尔科夫随机场进行自适应的空间约束来迭代更新变化检测结果，使得检测结果能更好地保持图像的结构信息，并有效地去除孤立噪声。

(3)本发明由于将类别统计项的估计和自适应的空间约束融为一体，提高了变化检测处理效率。

附图说明

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

图2是本发明第一组实验的变化检测结果图；

图3是本发明第二组实验的变化检测结果图；

图4是本发明第三组实验的变化检测结果图；

具体实施方式

参照图1，本发明的实施如下：

步骤1，输入两幅不同时相的遥感图像，并对每幅图像的每个通道分别进行窗口大小为3×3像素的中值滤波，得到两时相的去噪后图像X₁和X₂。

步骤2，将去噪后的两幅图像X₁和X₂应用变化矢量分析，得到一幅差异图像X_d，并根据该差异图像计算变权马尔科夫随机场的权值因子W，具体步骤如下：

(2a)利用变化矢量分析法计算差异图像X_d，即

$X_d=\sqrt{\frac{{|X_{11}-X_{21}|}^2+{|X_{11}-X_{21}|}^2+{|X_{11}-X_{21}|}^2}{3}};---\left(1\right)$

其中，X₁₁、X₁₂及X₁₃为图像X₁的三个通道图像；X₂₁、X₂₂及X₂₃为图像X₂的三个通道图像。

(2b)计算权值因子W：

首先，计算每个像素点的特征值，即

$t(i,j)=\underset{m=1}{\overset{p}{\mathrm{\Σ}}}\underset{n=1}{\overset{p}{\mathrm{\Σ}}}|x(m,n)-u(i,j)|;---\left(2\right)$

其中p为像素点局部窗口的大小，x(m，n)为局部窗口内每个像素点的灰度值，u(i，j)为局部窗口像素的均值；

然后，利用整幅图像中像素特征值t(i，j)的最大值及最小值将t(i，j)映射到[V_min，V_max]区间上，得到每个像素点的权值因子W(i，j)，V_min＝0.5，V_max＝8。

步骤3，应用K-均值聚类算法将差异图像聚成两类，将均值较大的一类作为变化类，均值较小的一类作为非变化类，得到初始分类结果。

步骤4，利用初始分类结果，采用非参数密度估计方法估计差异图像中变化类和非变化类的类条件概率密度，再对该类条件概率密度取负自然对数，得到变化类和非变化类的似然能量，具体步骤如下：

(4a)采用非参数密度估计方法估计差异图像中变化类的类条件概率密度

和非变化类的类条件概率密度

即

$\hat{P}\left(X_{\mathrm{ij}}\right|{\mathrm{\ω}}_n)=\frac{1}{H_n{N}_n}\underset{k=1}{\overset{N_n}{\mathrm{\Σ}}}K\left(\frac{X_{\mathrm{ij}}-S_{n_k}}{H_n}\right)---\left(3\right)$

$\hat{P}\left(X_{\mathrm{ij}}{|\mathrm{\ω}}_c\right)=\frac{1}{H_c{N}_c}\underset{k=1}{\overset{N_c}{\mathrm{\Σ}}}K\left(\frac{X_{\mathrm{ij}}-S_{c_k}}{H_c}\right);---\left(4\right)$

其中，S_n和S_c分别表示非变化类和变化类的像素集合，N_n和N_c分别表示非变化类和变化类的像素数目，K(□)为高斯核函数，H_n和H_c分别表示非变化类和变化类的自适应窗宽平滑参数，与像素数目及像素点X_ij的频数f(X_ij)有关，通过如下公式计算：

H_n＝H₀(α/N_n-f(X_ij)·β)；                    (5)

H_c＝h₀(α/N_c-f(X_ij)·β)；                    (6)

其中，H₀，α和β均为经验常数，H₀＝1，α＝40000，β＝10；

(4b)对非变化类和变化类的类条件概率密度取负自然对数，得到的非变化类的似然能量LEu(i，j)及变化类的似然能量LEc(i，j)，即

$\mathrm{LEu}(i,j)=-\ln \left(\hat{p}\left(X_{\mathrm{ij}}\right|{\mathrm{\ω}}_n)\right);---\left(7\right)$

$\mathrm{LEc}(i,j)=-\ln \left(\hat{p}\left(X_{\mathrm{ij}}\right|{\mathrm{\ω}}_c)\right).---\left(8\right)$

步骤5，对初始分类结果利用马尔科夫随机场计算变化类和非变化类的先验能量，并对马尔科夫随机场采用各向同性的二阶马尔科夫随机场邻域，则变化类的先验能量PEc(i，j)及非变化类PEu(i，j)的先验能量为：

$\mathrm{PEc}(i,j)=-\underset{(p,q)\∈S}{\mathrm{\Σ}}V(C(i,j),C(p,q));---\left(9\right)$

PEu(i，j)＝-8-PEc(i，j)。                            (10)

其中，C(i，j)为像素点(i，j)处的类别，S为C(i，j)的二阶马尔科夫随机场邻域，C(p，q)为S中的类别，V(C(i，j)，C(p，q))为邻域势函数，通过狄拉克函数计算：

$V(C(i,j),C(p,q))=\left\{\begin{array}{cc}-1,& C(i,j)=C(p,q)\\ 0,& \mathrm{esle}\end{array}\right.;---\left(11\right)$

步骤6，利用权值因子W、变化类和非变化类的似然能量及变化类和非变化类的先验能量计算非变化类的总能量TEu(i，j)和变化类的总能量TEc(i，j)：

TEu(i，j)＝LEu(i，j)+W(i，j)×PEu(i，j)；                (12)

TEc(i，j)＝LEc(i，j)+W(i，j)×PEc(i，j)，                (13)

若TEu(i，j)＜TEc(i，j)，则将像素点(i，j)处的类别更新为非变化类，否则为变化类，得到类别更新后的结果。

步骤7，对类别更新后的结果，采用非参数密度估计方法重新估计差异图像中变化类和非变化类的类条件概率密度，再对该类条件概率密度取负自然对数得到变化类和非变化类的似然能量，并利用马尔科夫随机场重新计算变化类和非变化类的先验能量。

步骤8，重复步骤(6)及步骤(7)直至迭代终止，并存储每次类别更新后的结果，得到每个像素点的类别更新集合，该迭代终止条件有两种：一种是迭代次数不超过50次，另一种是两次迭代之间相异的像元数目比例小于给定的阈值T，T＝5×10^-8。

步骤9，利用每个像素点的类别更新集合估计变化类的后验概率和非变化类的后验概率，将后验概率较大的那一类作为该像素点的最终变化检测结果。

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

本发明的对比实验为宋妍和袁修孝(2009)等在文章“基于混合高斯密度模型和空间上下文信息的遥感图像变化检测方法及扩展”中提出的变化检测方法，变化检测结果的性能采用虚警数、漏检数及总错误数三个指标进行评价。

本发明所设计的三组实验；

第一组为ATM(Airborne Thematic Mapper)3波段图像和模拟变化图像构成的模拟数据集，分别如图2(a)和图2(b)所示。其中ATM图像位于英国Feltwell村庄的一个农田区，模拟变化图像是通过模拟地球的天气变化和电磁波的辐射特性等因素影响并人工地嵌入一些变化区域得到，图像大小均为470×335，256灰度级，两幅图像的配准误差为1.5个像元左右。图2(c)为变化参考图。对图2(a)和图2(b)应用变化矢量分析法得到的差异图像，如图2(d)所示。图2(e)为采用对比实验方法得到的变化检测结果，图2(f)为采用本发明方法得到的变化检测结果。

第二组为2000年4月和2002年5月的墨西哥郊外的两幅Landsat 7ETM+4波段遥感图像，分别如图3(a)和图3(b)所示。图像大小均为512×512，256灰度级，图像配准误差为1.5个像元左右，变化区域主要为大火破坏了大面积的当地植被所致，变化参考图如图3(c)所示。对图3(a)和图3(b)应用变化矢量分析法得到的差异图像，如图3(d)所示。图3(e)为采用对比实验方法得到的变化检测结果，图3(f)为采用本发明方法得到的变化检测结果。

第三组为1995年9月和1996年7月Landsat-5卫星TM(Thematic Mapper)传感器接收的两幅多光谱图像，分别如图4(a)和图4(b)所示。图像大小均为300×412，256灰度级。试验区为意大利撒丁岛包含湖泊的一部分，变化前后湖中水位上升，变化参考图如图4(c)所示。对图4(a)和图4(b)应用变化矢量分析法得到的差异图像，如图4(d)所示。图4(e)为采用对比实验方法得到的变化检测结果，图4(f)为采用本发明方法得到的变化检测结果。

表1为第一组实验结果，从表中可以看出：与对比实验相比，本发明方法的虚警数减少了1280个像元，漏检数增加了779个像元，但总的错误数减少了501个像元。从图2(e)和图2(f)可以看出：与对比实验方法结果相比，本发明方法减少了孤立噪声，有效地保持变化区域的结构信息，总体上说本发明方法是有效的。

表2为第二组实验结果，从表中可以看出：与对比实验相比，本发明方法的虚警数减少了623个像元，漏检数增加了500个像元，但总的错误数减少了123个像元。从图3(e)和图3(f)可以看出：与对比实验方法结果相比，本发明方法减少了孤立噪声，有效地保持变化区域的结构信息，总体上说本发明方法是有效的。

表3为第三组实验结果，从表中可以看出：与对比实验相比，本发明方法的虚警数减少了571个像元，漏检数增加了153个像元，但总的错误数减少了418个像元。从图4(e)和图4(f)可以看出：与对比实验方法结果相比，本发明方法减少了孤立噪声，有效地保持变化区域的结构信息，总体上说本发明方法是有效的。

表1第一组实验结果

表2第二组实验结果

表3第三组实验结果

Claims (3)

1.一种基于非参数密度估计的遥感图像变化检测方法，包括如下步骤：

(1)输入两幅不同时相的遥感图像，并对每幅图像的每个通道分别进行窗口大小为3×3像素的中值滤波，得到两时相的去噪后图像；

(2)将去噪后的两幅图像应用变化矢量分析，得到一幅差异图像，并根据该差异图像计算变权马尔科夫随机场的权值因子W；

(3)应用K-均值聚类算法将差异图像聚成变化类和非变化类，得到初始分类结果；

(4)利用初始分类结果，采用非参数密度估计方法估计差异图像中变化类和非变化类的类条件概率密度，再对该类条件概率密度取负自然对数得到变化类和非变化类的似然能量；

(5)对初始分类结果利用马尔科夫随机场计算变化类和非变化类的先验能量；

(6)利用权值因子W、变化类和非变化类的似然能量及变化类和非变化类的先验能量计算变化类的总能量和非变化类的总能量，将总能量较小的那一类作为当前类别，得到类别更新后的结果；

(7)对类别更新后的结果，采用非参数密度估计方法重新估计差异图像中变化类和非变化类的类条件概率密度，再对该类条件概率密度取负自然对数得到变化类和非变化类的似然能量，并利用马尔科夫随机场重新计算变化类和非变化类的先验能量；

(8)重复步骤(6)及步骤(7)直至迭代终止，并存储每次类别更新后的结果，得到每个像素点的类别更新集合，该迭代终止条件为迭代次数不超过50次及两次迭代之间相异的像元数目比例小于给定的阈值；

(9)利用每个像素点的类别更新集合估计变化类的后验概率和非变化类的后验概率，将后验概率较大的那一类作为该像素点的最终变化检测结果。

2.根据权利要求1所述的遥感图像变化检测方法，其中步骤(2)所述计算变权马尔科夫随机场的权值因子，按如下步骤计算：

首先，计算每个像素点的特征值，即

$t(i,j)=\underset{m=1}{\overset{p}{\mathrm{\Σ}}}\underset{n=1}{\overset{p}{\mathrm{\Σ}}}|x(m,n)-u(i,j)|;$

其中p为像素点局部窗口的大小，x(m，n)为局部窗口内每个像素点的灰度值，u(i，j)为局部窗口像素的均值；

然后，利用整幅图像中像素特征值t(i，j)的最大值及最小值将t(i，j)映射到[V_min，V_max]区间上，得到每个像素点的权值因子W(i，j)，V_min＝0.5，V_max＝8。

3.根据权利要求1所述的遥感图像变化检测方法，其中步骤(4)所述采用非参数密度估计方法估计差异图像中变化类和非变化类的类条件概率密度，通过如下公式进行：

$\hat{P}\left(X_{\mathrm{ij}}\right|{\mathrm{\ω}}_n)=\frac{1}{H_n{N}_n}\underset{k=1}{\overset{N_n}{\mathrm{\Σ}}}K\left(\frac{X_{\mathrm{ij}}-S_{n_k}}{H_n}\right);$

$\hat{P}\left(X_{\mathrm{ij}}\right|{\mathrm{\ω}}_c)=\frac{1}{H_c{N}_c}\underset{k=1}{\overset{N_c}{\mathrm{\Σ}}}K\left(\frac{X_{\mathrm{ij}}-S_{c_k}}{H_c}\right);$

其中，

和

分别表示非变化类和变化类的类条件概率密度，S_n和S_c分别表示非变化类和变化类的像素集合，N_n和N_c分别表示非变化类和变化类的像素数目，K(□)为高斯核函数，H_n和H_c分别表示非变化类和变化类的自适应窗宽平滑参数：

H_n＝H₀(α/H_n-f(X_ij)·β)；

H_c＝H₀(α/N_c-f(X_ij)·β)；

其中，H₀，α和β均为经验常数，H₀＝1，α＝40000，β＝10。

Images