CN102867187A - Nsst域mrf与自适应阈值融合的遥感图像变化检测方法 - Google Patents

Abstract

本发明公开一种NSST域MRF与自适应阈值融合的遥感图像变化检测方法，解决了现有变化检测方法中不能在有效去除杂点的同时保持变化区域边缘信息的缺点。其实现过程是：输入两幅不同时相的遥感图像，用差值法构造差异图像；对差异图像进行非下采样Shearlet分解，将每一层方向子带合并为一个高频子带；对各层高频子带和低频子带分别自适应阈值分类，每层得到一幅高频自适应阈值分类图和一幅低频自适应阈值分类图；分别对各层的低频子带MRF分类，每层得到一幅MRF分类图；融合分类结果得到变化检测结果。本发明既有较强的抗噪性，又有边缘信息保持好的优点，检测结果误检少，准确率高。用于城区变化监测、森林和植被变化监测、军事目标监测等领域。

Description

NSST域MRF与自适应阈值融合的遥感图像变化检测方法

技术领域

本发明属于数字图像处理领域，涉及遥感图像变化检测，主要涉及一种非下采样Shearlet变换(Nonsubsampled Shearlet Transform，NSST)域马尔科夫随机场(MarkovRandom Field，MRF)与自适应阈值融合的遥感图像变化检测，具体是一种NSST域MRF与自适应阈值融合的遥感图像变化检测方法。该方法用于对遥感图像变化检测中的差异图像分类。

背景技术

遥感图像的变化检测主要依据同一地区不同时相的遥感图像存在电磁波谱特征差异或空间结构特征差异，通过分析和提取这些差异特征来实现识别地物类型的转变或内部条件和状态的变化。在国民经济和国防建设的诸多领域已得到广泛应用，如农作物的生长监测、森林和植被变化监测、城区变化监测、军事目标监测等。

遥感图像变化检测中变化/非变化信息检测的关键步骤是选择合适的阈值或二值分类方法。其中，无监督变化检测由于不需要昂贵且难以获得、甚至是常常无法获得的地面真实数据，因此在许多应用中是基本的普遍采用的变化检测方法。基于贝叶斯理论的阈值方法是无监督变化检测中的常用方法，MRF分类方法则由于考虑了像素点的邻域信息常用于提高分类准确性。

Bruzzone等学者(2000)在文章“Automatic Analysis of the Difference Image forUnsupervised Change Detection，IEEE Transactions on Geoscience and Remote Sensing，38(3)：1171-1182”中提出了两种基于贝叶斯理论的无监督变化检测方法。一种方法是在假设差值差异图的像素在空间上是独立的条件下，最小化全局变化检测误差来自动选取判决阈值的方法。该方法能保持较好的边缘信息，但是抗噪性能较差。另一种方法是考虑差值差异图中各个像素的邻域信息来分析差值图，该方法将差值差异图看作一个马尔科夫随机场，在基于贝叶斯准则下对差异图像分类。该方法有较强的抗噪性能，但是得到的变化检测结果并不能保持很好的边缘。

为了进一步提高变化检测的精度，一些学者提出了变化域的无监督变化检测方法，Bovolo等学者(2005)在文章“A Detail-Preserving Scale-Driven Approach to ChangeDetection in Multitemporal SAR Images，IEEE Transactions on Geoscience and RemoteSensing，2005，43(12)：2963-2972”中提出了一种尺度驱动融合的SAR图像变化检测方法，该方法首先对两时相图的对数比值图进行二维离散平稳小波分解，然后通过局部变化系数和全局变化系数的关系确定可靠尺度，并采用可靠尺度特征级融合策略进行变化检测。该方法虽然根据一定准则对变化域各尺度低频信息加以利用，提高了变化检测的准确性，但是用该方法提出的准则融合各尺度信息时，如果融合零尺度信息(即原对数比值图)，则变化检测结果边缘保持较好但是误检高，杂点多，如果不融合零尺度信息，则变化检测结果误检少但边缘保持较差，所以该方法并不能很好的权衡误检和边缘保持特性。

Celik(2010)在文章“A Bayesian approach to unsupervised multiscale changedetection in synthetic aperture radar images，Signal Processing，2010，90：1471-1485”中提出了一种变化域的自适应阈值变化检测方法，该方法先将差值差异图进行双树复小波变换，然后在贝叶斯最小错误率准则下对各子带系数进行分类，最后融合各子带系数分类结果得到变化检测结果图，该方法中双树复小波是下采样的，在进行尺度间融合时要对粗尺度结果进行插值，从而使得最终结果图的边缘呈块状，边缘保持很差。

Moser等学者(2011)在文章“Multiscale Unsupervised Change Detection on OpticalImages by Markov Random Fields and Wavelets，IEEE Geoscience and Remote SensingLetters，2011，8(4)：725-729”中提出了一种多尺度上下文的无监督变化检测方法，该方法首先对差值差异图进行离散小波分解，将各个尺度的低频子带分别重构回原图像域，提取多尺度特征，然后用马尔科夫链数据融合方法来融合这些特征信息和空间上下文信息，最后最小化能量函数得到变化检测结果。该方法用马尔科夫链将多个尺度信息和空间上下文信息融合起来，提高了变化检测结果的准确性，但其得到的最终的似然能量函数仅是各个尺度似然能量函数的简单叠加，并没考虑不同尺度特征之间的关系，并不能严格满足条件独立的假设，在一定程度上影响变化检测结果的精度。

上述的变化检测方法并不能很好的权衡降低误检和保持边缘信息这两方面特性，使得检测的结果或者边缘信息保持好，但误检多，杂点多；或者误检少，杂点少，但边缘信息保持不好，不能得到理想的变化检测结果。

发明内容

本发明的目的在于克服上述已有变换检测技术的不足，提出了一种NSST域MRF与自适应阈值融合的遥感图像变化检测方法。该方法有效地融合了基于贝叶斯最小错误率的自适应阈值方法与MRF分类方法，既有较强的抗噪性又使得检测的结果误检少，边缘信息保持好，检测准确率高。

本发明是一种NSST域MRF与自适应阈值融合的遥感图像变化检测方法，包括有如下步骤：

步骤1，输入同一地区不同时相的已配准的两幅遥感图像，对该两幅图对应像素点灰度值做差得到一幅差值差异图像。

步骤2，对差值差异图像进行非下采样Shearlet分解，分解的每层有一个低频子带和多个高频方向子带，将每层中所有高频方向子带合并为一个高频子带，也就是说非下采样Shearlet分解后，每层有一个低频子带和一个由该层所有高频方向子带合并的高频子带。

步骤3，应用K-均值聚类算法分别将各层的低频子带和高频子带聚类，在每一层分别得到一幅低频初始分类图和一幅高频初始分类图。

步骤4，利用每一层的低频初始分类图，对对应层中的低频子带在基于贝叶斯最小错误率下进行自适应阈值分类，在每一层得到一幅低频自适应阈值分类图。

步骤5，利用每一层的高频初始分类图，对对应层中的高频子带在基于贝叶斯最小错误率下进行自适应阈值分类，在每一层得到一幅高频自适应阈值分类图。

在差值差异图的各分解层中，每一层都有一幅低频自适应阈值分类图和一幅高频自适应阈值分类图。

步骤6，对非下采样Shearlet分解得到的每一层的低频子带进行MRF建模，求满足最小化能量函数的类别标记，得到该层的一幅MRF分类图。

这样各分解层的每一层都得到一幅杂点少，误检少的MRF分类图。

步骤7，忽略第1层的高频自适应阈值分类图，将第1层的低频自适应阈值分类图作为第1层的总自适应阈值分类图，在其余各分解层中，对每一层的高频自适应阈值分类图和低频自适应阈值分类图进行并集融合，得到该层的一幅总自适应阈值分类图。

第1层的高频自适应阈值分类图中主要分布的是杂点，包含的边缘信息很少，几乎没有，所以本发明忽略了第1层高频自适应阈值分类图。并集运算使得每层得到的总自适应阈值分类图的变化信息更完整，边缘信息保持更好。

步骤8，对所有分解层的总自适应阈值分类图进行交集融合，得到一幅最终自适应阈值分类图，对所有分解层的MRF分类图进行交集融合，得到一幅最终MRF分类图。

交集融合的作用主要是提取进行融合的图像的公共信息，在一定程度上减少杂点，通过交集融合使得最终自适应阈值分类图和最终MRF分类图中的杂点有所减少。

步骤9，将最终自适应阈值分类图和最终MRF分类图进行融合得到变化检测结果图。

为了提高变化检测的精度，本发明对差值差异图进行了非下采样Shearlet分解，分解的低频子带保留了差值差异图的大部分信息，高频方向子带保留了差值差异图各个方向的边缘信息，本发明按照步骤2中的方法合并高频方向子带为一个高频子带，使得该高频子带包含差值差异图的所有边缘信息，除此以外，采用非下采样Shearlet分解克服了融合时需对粗尺度进行插值的缺点，保持了变化检测结果的边缘，提高了变化检测结果的正确性。本发明考虑到变换域低频子带保留了差值差异图的大部分信息，同时去除了大部分噪声，且越高层低频子带越平滑，噪声越少，所以本发明分别对各层低频子带进行MRF建模分类，然后对所有层的MRF分类结果进行交集融合，得到一幅杂点很少，误检少的最终MRF分类图。本发明还考虑到基于贝叶斯最小错误率的自适应阈值方法使检测结果的边缘保持好的优点，将该方法在变换域的分类结果与最终MRF分类结果进行了融合，因此，本发明同时具有使变化检测结果误检少，漏检少，边缘保持好的优点。

本发明的实现还在于：步骤2所述的将所有高频方向子带合并为一个高频子带是通过如下公式进行：

$D_{h,s}=\frac{1}{N_s}\underset{b\∈N_s}{\mathrm{\Σ}}{\left|D_{b,s}\right|}^2$

其中，D_b，s表示第s层第b个高频方向子带，D_h，s表示第s层的高频子带，下标h表示为高频，分解层s＝1，2，3，4，N_s表示高频方向子带个数，s＝1和2时，N_s＝10，s＝3和4时，N_s＝6。

其他方法中合并高频方向子带是将同层所有高频方向子带相加，然后除以方向子带个数，再开方，这样抑制了变化区域的边缘信息，不利于提取边缘信息。本发明中对同层的高频方向子带的绝对值进行了平方和，然后除以方向子带个数，该合并方法对变化区域的边缘信息进行了增强，有利于提取边缘信息。

本发明的实现还在于：步骤6所述的得到MRF分类图的步骤包括：

(6a)将步骤3所得的某一层的低频初始分类图作为该层的初始标记图；如将第s层的低频初始分类图作为第s层的初始标记图，变化类像素

用1表示，非变化类像素

用0表示，设迭代次数K＝1，最终迭代次数为K，K的取值范围为K＞＝5；

(6b)第s层低频子带D_l,s中像素点(i，j)处变化类的似然能量

和非变化类的似然能量

由二维高斯概率密度的负自然对数计算得到，即

${\mathrm{Uc}}_{\mathrm{data}}^{l,s}(i,j)=\frac{1}{2}\ln \left|2\mathrm{\π}{\left({\mathrm{\σ}}_{l,c}^s\right)}^2\right|+\frac{1}{2}{(d_l^s(i,j)-{\mathrm{\μ}}_{l,c}^s)}^2{\left[{\left({\mathrm{\σ}}_{l,c}^s\right)}^2\right]}^{-1};$

${\mathrm{Uu}}_{\mathrm{data}}^{l,s}(i,j)=\frac{1}{2}\ln \left|2\mathrm{\π}{\left({\mathrm{\σ}}_{l,u}^s\right)}^2\right|+\frac{1}{2}{(d_l^s(i,j)-{\mathrm{\μ}}_{l,u}^s)}^2{\left[{\left({\mathrm{\σ}}_{l,u}^s\right)}^2\right]}^{-1}$

(6c)计算第s层低频子带内像素点(i，j)处变化类的先验能量

和非变化类的先验能量

${\mathrm{Uc}}_{\mathrm{context}}^{l,s}(i,j)=\underset{(p,q)\∈N(i,j)}{\mathrm{\Σ}}\mathrm{\δ}(C^{l,s}(i,j),C^{l,s}(p,q));$

${\mathrm{Uu}}_{\mathrm{context}}^{l,s}(i,j)=-8-{\mathrm{Uc}}_{\mathrm{context}}^{l,s}(i,j)$

其中，C^l，s(i，j)为第s层低频子带内像素点(i，j)处的类别标记，N(i，j)为像素点(i，j)的二阶邻域，N(i，j)＝{(i，j)+(v，ζ)，(v，ζ)∈N}，N＝{(±1，0)，(0，±1)，(1，±1)，(-1，±1)}，C^l,s(p，q)为二阶邻域N(i，j)中(p,q)处的类别标记，p和q分别为二阶邻域N(i，j)中的行序号和列序号，δ(C^l,s(i，j)，C^l，s(p，q))为狄拉克函数，其表达式为

(6d)用第s层低频子带内像素点(i，j)处变化类的似然能量

和先验能量

计算第s层低频子带内像素点(i，j)处变化类的总能量

用非变化类的似然能量和先验能量计算非变化类的总能量 $U_u^{l,s}(i,j):$

$u_c^{l,s}(i,j)={\mathrm{Uc}}_{\mathrm{data}}^{l,s}(i,j)+\mathrm{\β}{\mathrm{Uc}}_{\mathrm{context}}^{l,s}(i,j)$

$u_u^{l,s}(i,j)={\mathrm{Uu}}_{\mathrm{data}}^{l,s}(i,j)+\mathrm{\β}{\mathrm{Uu}}_{\mathrm{context}}^{l,s}(i,j)$

其中，参数β为常数，β＞0，用来调整先验能量和似然能量在变化检测过程中的重要性；

(6e)若

则将像素点(i，j)处的类别标记更新为变化类，否则，更新为非变化类，对图像所有像素点做同样处理，得到更新的类别标记图像，更新迭代次数k，令k＝k+1，当k＜＝K时，转到(6b)；否则，转到(6f)；

(6f)将更新的类别标记图像作为第s层的MRF分类图M_s＝{M_s(i，j)|1≤i≤I，1≤j≤J}。这样就得到了第s层的MRF分类图。

本发明对变换域中每层的低频子带都进行马尔科夫建模分类，有效的利用了变换域低频子带中像素之间的依赖关系，有利于去除分类结果中的杂点。

本发明的实现还在于：步骤9所述的将最终自适应阈值分类图B和最终MRF分类图M进行融合，融合步骤包括：

(9a)设融合循环次数z＝1，融合终止循环次数为Z，Z＞50；

(9b)检测最终MRF分类图M的像素值是否等于1，如果M(i，j)＝1，则将最终自适应阈值分类图B中(i，j)处及其二阶邻域N(i，j)中的值赋给融合图像CM的(i，j)处及其二阶邻域的对应位置，否则令CM(i，j)的值为0；

(9c)按照步骤(9b)从左到右从上到下检测M中的所有像素值，令融合图像CM中像素值大于0的像素值为1；

(9d)更新融合循环次数z，令z＝z+1，当z＞Z或两次迭代之间类别发生变化的像元总数小于预先设定的阈值T，则停止，融合图像CM为变化检测结果图；否则将融合图像CM赋给M，转到(9b)。

最终自适应阈值分类图B的边缘保持好，但误检多，杂点多，而最终MRF分类图M的误检少，杂点少，但边缘保持较差，该步骤检测最终MRF分类图M中变化区域的位置，将自适应阈值分类图B中的对应位置及其邻域像素值赋给融合图像CM，得到变化检测结果，这样能有效的去除杂点，使得变换检测结果具有误检少，漏检少，边缘保持好的特点。

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

1、本发明由于采用非下采样Shearlet变换分解差值差异图，克服了融合时需对粗尺度进行插值的缺点，使检测结果的边缘不再呈现块状，提高了变化检测结果的正确率。

2、本发明对变换域中每层的低频子带都进行马尔科夫建模分类，充分利用了变换域低频子带中像素之间的依赖关系，去除了变化检测结果中的杂点，减少了变化检测结果的误检数。

3、本发明由于有效地融合了基于贝叶斯最小错误率的自适应阈值分类结果和MRF分类结果，不仅减少了变化检测结果的误检数，还保持了变化区域的边缘。

附图说明

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

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

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

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

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

具体实施方式

下面结合附图对本发明详细说明

实施例1

本发明是一种NSST域MRF与自适应阈值融合的遥感图像变化检测方法，主要是对遥感图像进行处理，对图像的处理需要大容量的计算机作为硬件支持，使用matlab等软件来实现图像处理，参见图1，遥感图像变化检测包括有如下步骤：

步骤1，对输入的两幅大小均为I×J的同一地区不同时相已配准的遥感图像X₁和X₂，如图2(a)和图2(b)所示，将图像X₁和X₂空间对应位置(m,n)处的像素点灰度值X₁(m，n)和X₂(m，n)进行差值计算，得到差值X_d(m,n)＝|X₁(m,n)-X₂(m，n)|，其中，m和n分别为遥感图像的行和列序号，m＝1，2，...，I，n＝1，2，…，J，由此得到一幅差异图像X_d。

步骤2，对差异图像X_d进行非下采样Shearlet分解，分解层数取值范围为1～5层，为了在去除噪声和信息保留两方面取得较好的权衡，本例中进行了4层分解。分解后每层得到一个低频子带和多个高频方向子带，其中，第1层到第4层的高频方向子带个数分别为10、10、6、6。第s层的低频子带D_l，s表示为

其中，

为第s层的低频子带中(i，j)处的系数值，i和j分别为变换域图像的行和列序号，下标l表示为低频，分解层s＝1，2，3，4；第s层的一个高频方向子带表示为D_b，s＝{d_b，s(i，j)|1≤i ≤I，1≤j≤J}，其中，d_b，s(i，j)表示第s层的第b个高频方向子带中(i，j)处的系数值，下标b表示为方向子带序号，当s＝1和2时，b＝1，2，…，10，当s＝3和4时，b＝1，2，…，6。

将第s层的高频方向子带根据下式合并为一个高频子带D_h，s， $D_{h,s}=\left\{d_h^s(i,j)\right|1\≤i\≤I,1\≤j\≤J\}$

$D_{h,s}=\frac{1}{N_s}\underset{b\∈N_s}{\mathrm{\Σ}}{\left|D_{b,s}\right|}^2$

其中，

表示第s层的高频子带中(i，j)处的系数值，下标h表示为高频，N_s表示高频方向子带个数，当s＝1和2时，N_s＝10，当s＝3和4时，N_s＝6。

步骤3，应用K-均值聚类算法将第s层的低频子带D_l，s和高频子带D_h，s分别聚为两类，即变化类和非变化类，变化类用下标c表示，非变化类用下标u表示，得到第s层的一幅低频初始分类图和一幅高频初始分类图。

步骤4，依据第s层的低频初始分类图，对第s层的低频子带D_l,s在基于贝叶斯最小错误率下进行自适应阈值分类，得到一幅低频自适应阈值分类图B_l，s；

(4a)利用步骤3所得的第s层的低频初始分类图，计算该层低频子带D_l，s中变化类的先验概率

和非变化类的先验概率

$p\left({\mathrm{\ω}}_{l,c}^s\right)=\frac{N_{l,c}^s}{N_{l,c}^s+N_{l,u}^s};$

$p\left({\mathrm{\ω}}_{l,u}^s\right)=\frac{N_{l,u}^s}{N_{l,c}^s+N_{l,u}^s}$

其中，

和

分别表示第s层的低频初始分类图中的变化类标记和非变化类标记，

和

分别表示第s层中低频初始分类图中变化类和非变化类的像素数目；

(4b)计算第s层低频子带D_l,s中处于(i，j)点的变化类的类条件概率密度

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

即

$p(d_l^s(i,j)/{\mathrm{\ω}}_{l,c}^s)=\frac{1}{\sqrt{2\mathrm{\π}{\left({\mathrm{\σ}}_{l,c}^s\right)}^2}}\exp \{-\frac{{(d_l^s(i,j)-{\mathrm{\μ}}_{l,c}^s)}^2}{2{\left({\mathrm{\σ}}_{l,c}^s\right)}^2}\};$

$p(d_l^s(i,j)/{\mathrm{\ω}}_{l,u}^s)=\frac{1}{\sqrt{2\mathrm{\π}{\left({\mathrm{\σ}}_{l,u}^s\right)}^2}}\exp \{-\frac{{(d_l^s(i,j)-{\mathrm{\μ}}_{l,u}^s)}^2}{2{\left({\mathrm{\σ}}_{l,u}^s\right)}^2}\}$

其中，

为第s层低频子带D_l,s中(i，j)处的系数值，

和

分别表示第s层低频子带D_l,s中变化类和非变化类的系数的均值，

和

分别表示第s层低频子带D_l,s中变化类和非变化类的系数的标准差；

(4c)用第s层的低频子带D_l,s中变化类的先验概率

和(i，j)点的类条件概率

及非变化类的先验概率

和(i，j)点的类条件概率

计算第s层的低频子带D_l，s中(i，j)点的概率密度函数

$p_l^s(i,j)=p(d_l^s(i,j)/{\mathrm{\ω}}_{l,c}^s)p\left({\mathrm{\ω}}_{l,c}^s\right)+p(d_l^s(i,j)/{\mathrm{\ω}}_{l,u}^s)p\left({\mathrm{\ω}}_{l,u}^s\right)$

(4d)对第s层的低频子带的所有像素点重复步骤(4b)至(4c)，用期望最大化(expectation maximization，EM)算法估计第s层低频子带D_l,s中变化类的先验概率均值

和方差

及非变化类的先验概率

均值和方差

(4e)用步骤(4d)中估计的变化类的均值

和方差

更新第s层低频子带D_l，s中每一个像素点的变化类的类条件概率密度，用估计的非变化类的均值

和方差

更新第s层低频子带D_l，s中每一个像素点的非变化类的类条件概率密度；

(4f)基于贝叶斯最小错误率准则，对(4a)至(4e)计算得到的第s层的低频子带D_l,s中(i，j)点的变化类的类条件概率密度和非变化类的类条件概率密度

变化类的先验概率

和非变化类的先验概率

计算(i，j)点的阈值分类结果B_l,s(i，j)：

$B_{l,s}(i,j)=\left\{\begin{array}{cc}1,& p\left(d_l^s(i,j)\right|{\mathrm{\&omega;}}_{l,c}^s)p\left({\mathrm{\&omega;}}_{l,c}^s\right)>=p\left(d_l^s(i,j)\right|{\mathrm{\&omega;}}_{l,u}^s)p\left({\mathrm{\&omega;}}_{l,u}^s\right)\\ 0,& p\left(d_l^s(i,j)\right|{\mathrm{\&omega;}}_{l,c}^s)p\left({\mathrm{\&omega;}}_{l,c}^s\right)<p\left(d_l^s(i,j)\right|{\mathrm{\&omega;}}_{l,u}^s)p\left({\mathrm{\&omega;}}_{l,u}^s\right)\end{array}\right.$

其中，1表示变化类，0表示非变化类。

(4g)对第s层的低频子带的所有像素点重复步骤(4f)得到第s层的一幅低频自适应阈值分类图B_l,s＝{B_l,s(i，j)|1≤i≤I，1≤j≤J}。

步骤5，依据第s层的高频初始分类图，对第s层的高频子带D_h，s在基于贝叶斯最小错误率下进行自适应阈值分类，得到第s层的一幅高频自适应阈值分类图B_h,s

步骤6，对第s层的低频子带D_l,s进行MRF建模，求满足最小能量函数的类别标记，得到第s层的MRF分类图。

(6a)将步骤3所得的某一层的低频初始分类图作为该层的初始标记图；如将第s层的低频初始分类图作为第s层的初始标记图，变化类像素

用1表示，非变化类像素用0表示，设迭代次数k＝1，最终迭代次数为K，K的取值范围为K＞＝5，本例中K＝5；

(6b)第s层低频子带D_l,s中像素点(i，j)处变化类的似然能量

和非变化类的似然能量

由二维高斯概率密度的负自然对数计算得到，即

${\mathrm{Uc}}_{\mathrm{data}}^{l,s}(i,j)=\frac{1}{2}\ln \left|2\mathrm{\&pi;}{\left({\mathrm{\&sigma;}}_{l,c}^s\right)}^2\right|+\frac{1}{2}{(d_l^s(i,j)-{\mathrm{\&mu;}}_{l,c}^s)}^2{\left[{\left({\mathrm{\&sigma;}}_{l,c}^s\right)}^2\right]}^{-1};$

${\mathrm{Uu}}_{\mathrm{data}}^{l,s}(i,j)=\frac{1}{2}\ln \left|2\mathrm{\&pi;}{\left({\mathrm{\&sigma;}}_{l,u}^s\right)}^2\right|+\frac{1}{2}{(d_l^s(i,j)-{\mathrm{\&mu;}}_{l,u}^s)}^2{\left[{\left({\mathrm{\&sigma;}}_{l,u}^s\right)}^2\right]}^{-1}$

(6c)计算第s层低频子带内像素点(i，j)处变化类的先验能量和非变化类的先验能量

${\mathrm{Uc}}_{\mathrm{context}}^{l,s}(i,j)=\underset{(p,q)\&Element;N(i,j)}{\mathrm{\&Sigma;}}\mathrm{\&delta;}(C^{l,s}(i,j),C^{l,s}(p,q));$

${\mathrm{Uu}}_{\mathrm{context}}^{l,s}(i,j)=-8-{\mathrm{Uc}}_{\mathrm{context}}^{l,s}(i,j)$

其中，C^l,s(i，j)为第s层低频子带内像素点(i，j)处的类别标记，N(i，j)为像素点(i，j)的二阶邻域，N(i，j)＝{(i，j)+(v，ζ)，(v，ζ)∈N}，N＝{(±1，0)，(0，±1)，(1，±1)，(-1，±1)}，C^l，s(p，q)为二阶邻域N(i，j)中(p,q)处的类别标记，p和q分别为二阶邻域N(i，j)中的行序号和列序号，δ(C^l,s(i，j)，C^l,s(p，q))为狄拉克函数，其表达式为

(6d)用第s层低频子带内像素点(i，j)处变化类的似然能量

和先验能量

计算第s层低频子带内像素点(i，j)处变化类的总能量用非变化类的似然能量

和先验能量

计算非变化类的总能量 $U_u^{l,s}(i,j):$

$u_c^{l,s}(i,j)={\mathrm{Uc}}_{\mathrm{data}}^{l,s}(i,j)+\mathrm{\&beta;}{\mathrm{Uc}}_{\mathrm{context}}^{l,s}(i,j);$

$u_u^{l,s}(i,j)={\mathrm{Uu}}_{\mathrm{data}}^{l,s}(i,j)+\mathrm{\&beta;}{\mathrm{Uu}}_{\mathrm{context}}^{l,s}(i,j)$

其中，参数β为常数，B＞0，用来调整先验能量和似然能量在变化检测过程中的重要性，本例中β＝5；

(6e)若

则将像素点(i，j)处的类别标记更新为变化类，否则，更新为非变化类，对图像所有像素点做同样处理，得到更新的类别标记图像，更新迭代次数k，令k＝k+1，当k＜＝K时，转到(6b)；否则，转到(6f)；

(6f)将更新的类别标记图像作为第s层的MRF分类图M_s＝{M_s(i，j)|1≤i ≤I，1≤j≤J}。这样就得到了第s层的MRF分类图。

步骤7，忽略第1层的高频自适应阈值分类图，将第1层的低频自适应阈值分类图作为第1层的总自适应阈值分类图，在其余各分解层中，对每一层的高频自适应阈值分类图和低频自适应阈值分类图进行并集融合，得到该层的一幅总自适应阈值分类图。

步骤8，对所有分解层的总自适应阈值分类图进行交集融合，得到一幅最终自适应阈值分类图B，对所有分解层的MRF分类图进行交集融合，得到一幅最终MRF分类图M。

步骤9，对最终自适应阈值分类图B和最终MRF分类图M进行融合。

(9a)设融合循环次数z＝1，融合终止循环次数为Z，Z＞50，本例中Z＝100；

(9b)检测最终MRF分类图M的像素值是否等于1，如果M(i，j)＝1，则将最终自适应阈值分类图B中(i，j)处及其二阶邻域N(i，j)中的值赋给融合图像CM的(i，j)处及其二阶邻域的对应位置，否则令CM(i，j)的值为0；

(9c)按照步骤(9b)从左到右从上到下检测M中的所有像素值，令融合图像CM中像素值大于0的像素值为1；

(9d)更新融合循环次数z，令z＝z+1，当z＞Z或两次迭代之间类别发生变化的像元总数小于预先设定的阈值T，则停止，融合图像CM为变化检测结果图；否则将融合图像CM赋给M，转到(9b)。本例中T＝1，T还可以有不同的取值，但T＝1为最佳阈值。

这样就得到了变化检测结果图CM。本发明的融合方法不是对最终自适应阈值分类图B和最终MRF分类图简单叠加，而是最佳融合，完全融合了两幅图的优点。

实施例2

NSST域MRF与自适应阈值融合的遥感图像变化检测方法同实施例1，参照图1，对本发明的实现步骤进一步说明如下：

其中步骤5中得到第s层的一幅高频自适应阈值分类图B_h,s的具体操作是：

(5a)利用步骤3所得的第s层的高频初始分类图，计算该层高频子带D_h,s中变化类的先验概率和非变化类的先验概率

$p\left({\mathrm{\&omega;}}_{h,c}^s\right)=\frac{N_{h,c}^s}{N_{h,c}^s+N_{h,u}^s};$

$p\left({\mathrm{\&omega;}}_{h,u}^s\right)=\frac{N_{h,u}^s}{N_{h,c}^s+N_{h,u}^s}$

其中，

和

分别表示第s层的高频初始分类图中的变化类标记和非变化类标记，

和

分别表示第s层中高频初始分类图中变化类和非变化类的像素数目；

(5b)计算第s层高频子带D_h，s中(i，j)处变化类的类条件概率密度

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

即

$p(d_h^s(i,j)/{\mathrm{\&omega;}}_{h,c}^s)=\frac{1}{\sqrt{2\mathrm{\&pi;}{\left({\mathrm{\&sigma;}}_{h,c}^s\right)}^2}}\exp \{-\frac{{(d_h^s(i,j)-{\mathrm{\&mu;}}_{h,c}^s)}^2}{2{\left({\mathrm{\&sigma;}}_{h,c}^s\right)}^2}\};$

$p(d_h^s(i,j)/{\mathrm{\&omega;}}_{h,u}^s)=\frac{1}{\sqrt{2\mathrm{\&pi;}{\left({\mathrm{\&sigma;}}_{h,u}^s\right)}^2}}\exp \{-\frac{{(d_h^s(i,j)-{\mathrm{\&mu;}}_{h,u}^s)}^2}{2{\left({\mathrm{\&sigma;}}_{h,u}^s\right)}^2}\}$

其中为第s层高频子带D_h，s中(i，j)处的系数值，

和分别表示第s层高频子带D_h，s中变化类和非变化类的系数的均值，和分别表示第s层高频子带D_h,s中变化类和非变化类的系数的标准差；

(5c)用第s层的高频子带D_h，s中变化类的先验概率和(i，j)点的类条件概率及非变化类的先验概率

和(i，j)点的类条件概率

计算第s层的高频子带D_h,s中(i，j)点的概率密度函数

$p_h^s(i,j)=p(d_h^s(i,j)/{\mathrm{\&omega;}}_{h,c}^s)p\left({\mathrm{\&omega;}}_{h,c}^s\right)+p(d_h^s(i,j)/{\mathrm{\&omega;}}_{h,u}^s)p\left({\mathrm{\&omega;}}_{h,u}^s\right)$

(5d)对第s层的高频子带的所有像素点重复步骤(5b)至(5c)，用期望最大化算法估计第s层高频子带D_h，s中变化类的先验概率

均值

和方差及非变化类的先验概率

均值

和方差

(5e)用步骤(5d)中估计的变化类的均值

和方差

更新第s层高频子带D_h，s中每一个像素点的变化类的类条件概率密度，用估计的非变化类的均值

和方差

更新第s层高频子带D_h，s中每一个像素点的非变化类的类条件概率密度；

(5f)基于贝叶斯最小错误率准则，对(5a)至(5d)计算得到的第s层的高频子带D_h,s中(i，j)点的变化类的类条件概率密度

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

和非变化类的先验概率

计算(i，j)点的阈值分类结果B_h，s(i，j)：

$B_{h,s}(i,j)=\left\{\begin{array}{cc}1,& p\left(d_h^s(i,j)\right|{\mathrm{\&omega;}}_{h,c}^s)p\left({\mathrm{\&omega;}}_{h,c}^s\right)>=p\left(d_h^s(i,j)\right|{\mathrm{\&omega;}}_{h,u}^s)p\left({\mathrm{\&omega;}}_{h,u}^s\right)\\ 0,& p\left(d_h^s(i,j)\right|{\mathrm{\&omega;}}_{h,c}^s)p\left({\mathrm{\&omega;}}_{h,c}^s\right)<p\left(d_h^s(i,j)\right|{\mathrm{\&omega;}}_{h,u}^s)p\left({\mathrm{\&omega;}}_{h,u}^s\right)\end{array}\right.$

其中，1表示变化类，0表示非变化类。

(5g)对第s层的高频子带的所有像素点重复步骤(5f)，得到第s层的一幅高频自适应阈值分类图B_h，s＝{B_h，s(i，j)|1≤i ≤I，1≤j≤J}。

其中步骤7中忽略第1层的高频自适应阈值分类图B_h，1，将第1层低频自适应阈值分类图B_l，1作为第1层的总自适应阈值分类图B₁，对于第2、3、4层，对每一层的高频自适应阈值分类图B_h，s和低频自适应阈值分类图B_l，s进行并集融合，得到该层的一幅总自适应阈值分类图B_s，并集融合公式为：

B_s＝B_l，s∪B_h，s

其中，s＝2，3，4，符号∪表示并运算，B_l，s和B_h，s空间对应位置(i，j)处的像素值B_l，s(i，j)和B_h，s(i，j)有一个为1或两个同时为1时，则B_s的(i，j)处的像素值为1。

其中步骤8中对所有分解层的总自适应阈值分类图B_s进行交集融合得到一幅最终自适应阈值分类图B，对所有分解层的MRF分类图M_s进行交集融合得到一幅最终MRF分类图M，交集融合公式为：

B＝B₁∩B₂∩B₃∩B₄

M＝M₁∩M₂∩M₃∩M₄

其中符号∩表示交运算，当且仅当B₁、B₂、B₃、B₄空间对应位置(i，j)处的像素值B₁(i，j)、B₂(i，j)、B₃(i，j)、B₄(i，j)都为1时，B中(i，j)处的像素值为1，当且仅当M₁、M₂、M₃、M₄空间对应位置(i，j)处的像素值M₁(i，j)、M₂(i，j)、M₃(i，j)、M₄(i，j)都为1时，M中(i，j)处的像素值为1。

实施例3

NSST域MRF与自适应阈值融合的遥感图像变化检测方法同实施例1-2，其中步骤(6a)中最终迭代次数K的取值选定为15，这样得到的各层MRF分类图与K＝5时得到的各层MRF分类图同层相比，漏检减少，误检增加；K的取值选定为15，对所有分解层MRF分类图进行交集融合得到的最终MRF分类图与K＝5时的最终MRF分类图相比，漏检减少，误检增加；K的取值选定为15对最终的变化检测结果图影响不大，但其效率有所降低，也就是说K的取值增大时会对效率略有影响。

其中步骤(9a)中融合终止循环次数Z的取值选定为200，这样当最终MRF分类图的漏检很多时，也能很好的与最终自适应阈值分类图融合，使得融合图像保持较好的边缘信息。

实施例4

NSST域MRF与自适应阈值融合的遥感图像变化检测方法同实施例1-2，本发明的效果可通过以下实验结果与分析进一步说明：

1.实验数据

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

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

2.对比试验

本发明中的最终自适应阈值分类结果具有漏检少，边缘保持较好的优点，MRF分类结果具有误检少的优点。本发明将最终自适应阈值分类结果与MRF分类结果进行融合，能很好的在保持漏检少的情况下减少误检，并较好的保持了变化区域的边缘信息，具有较高的检测可靠性和准确性。

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

对比方法1，是Bruzzone等学者(2000)在文章“Automatic Analysis of the DifferenceImage for Unsupervised Change Detection”中提出的对差值差异图进行马尔科夫建模的分类方法得到的结果与本发明中最终自适应阈值分类结果进行融合的方法。本发明与对比方法1对比是为了验证将马尔科夫随机场模型引入到变换域，减少了误检数，提高了变化区域检测的可靠性和准确性。

对比方法2，是Celik(2010)在文章“A Bayesian approach to unsupervised multiscalechange detection in synthetic aperture radar images”中提出的一种变换域的无监督变化检测方法。本发明与对比方法2对比是为了验证非下采样shearlet变换能较好的保持边缘信息，提高变化区域检测精确性的优点。

对比方法3，本发明方法中不融合马尔科夫分类结果的方法。本发明与对比方法3对比是为了验证本发明将最终自适应阈值分类结果与MRF分类结果进行融合，能很好的在保持漏检少的同时减少误检的优点。

对比方法4，是Moser等学者在文章“Multiscale Unsupervised Change Detection onOptical Images by Markov Random Fields and Wavelets”中提出的方法。

3.实验结果与分析

图4是本发明与对比方法仿真第一组遥感图像得到的变化检测结果图。图4(a)是对比方法1的变化检测结果图。图4(b)是对比方法2的变化检测结果图。图4(c)是对比方法3的变化检测结果图。图4(d)是对比方法4的变化检测结果图。图4(e)是本发明的变化检测结果图。

图5是本发明与对比方法仿真第二组遥感图像得到的变化检测结果图。图5(a)是对比方法1的变化检测结果图。图5(b)是对比方法2的变化检测结果图。图5(c)是对比方法3的变化检测结果图。图5(d)是对比方法4的变化检测结果图。图5(e)是本发明的变化检测结果图。

从说明书附图4和附图5可以看出本发明比四种对比方法的变化区域的形状保持好，且本发明的变化检测结果图中杂点最少。本发明与对比方法1、对比方法3、对比方法4相比，杂点少；本发明与对比方法2相比，边缘保持好，没有块效应。图4中的(a)、(b)、(c)、(d)分别与图2(c)对比，明显可见图4(a)-图4(d)这四幅图都分布着很多杂点，而本发明的结果图，即图4(e)，几乎没有杂点。从图中分布的各个小图形可见，图4中(e)图的边缘明显优于图4(b)和图4(d)。

实施例5

NSST域MRF与自适应阈值融合的遥感图像变化检测方法同实施例1-2，表1列出了本发明与对比方法仿真第一组遥感图像和第二组遥感图像的定量结果。

从表1可以看出，对于第一组遥感图像，本发明的误检数最少，本发明的漏检数比对比方法2和对比方法4少，比对比方法1和对比方法3多。对于第二组遥感数图像，本发明的误检数最少，本发明的漏检数比对比方法2和对比方法4少，与对比方法1和对比方法3相等。本发明的总错误数最少，性能最优。综上所述，本发明能够获得比较准确的变化信息，具有较强的抗噪性，能有效的去除杂点，同时能保持较好的变化区域的边缘信息。

表1

综上，本发明的NSST域MRF与自适应阈值融合的遥感图像变化检测方法，主要解决了现有变化检测方法中不能在有效去除杂点的同时保持变化区域边缘信息的缺点。其实现过程是：输入两幅不同时相的遥感图像，用差值法构造差异图像；对差异图像进行非下采样Shearlet分解，在每一层将该层所有方向子带合并为一个高频子带；对各层的高频子带和低频子带分别进行自适应阈值分类，每层得到一幅高频自适应阈值分类图和一幅低频自适应阈值分类图；同时对各层的低频子带进行MRF建模，每层得到一幅MRF分类图；融合分类结果得到变化检测结果。实验表明本发明能够在保持变化区域边缘的同时有效的去除杂点，提高变化检测的精度，可用于城区变化监测、森林和植被变化监测、军事目标监测等领域。

Claims (4)

1.一种NSST域MRF与自适应阈值融合的遥感图像变化检测方法，其特征在于：包括如下步骤：

步骤1，输入同一地区不同时相的已配准的两幅遥感图像，对该两幅图对应像素点灰度值做差得到一幅差值差异图像；

步骤2，对差值差异图像进行非下采样Shearlet分解，分解的每层有一个低频子带和多个高频方向子带，将每层中所有高频方向子带合并为一个高频子带；

步骤3，应用K-均值聚类算法分别将各层的低频子带和高频子带聚类，在每一层分别得到一幅低频初始分类图和一幅高频初始分类图；

步骤4，利用每一层的低频初始分类图，对对应层中的低频子带在基于贝叶斯最小错误率下进行自适应阈值分类，在每一层得到一幅低频自适应阈值分类图；

步骤5，利用每一层的高频初始分类图，对对应层中的高频子带在基于贝叶斯最小错误率下进行自适应阈值分类，在每一层得到一幅高频自适应阈值分类图；

步骤6，对非下采样Shearlet分解得到的每一层的低频子带进行MRF建模，求满足最小化能量函数的类别标记，得到该层的一幅MRF分类图；

步骤7，忽略第1层的高频自适应阈值分类图，将第1层的低频自适应阈值分类图作为第1层的总自适应阈值分类图，在其余各分解层中，对每一层的高频自适应阈值分类图和低频自适应阈值分类图进行并集融合，得到该层的一幅总自适应阈值分类图；

步骤8，对所有分解层的总自适应阈值分类图进行交集融合，得到一幅最终自适应阈值分类图，对所有分解层的MRF分类图进行交集融合，得到一幅最终MRF分类图；

步骤9，将最终自适应阈值分类图B和最终MRF分类图M进行融合得到变化检测结果图。

2.根据权利要求1所述的遥感图像变化检测方法，其特征在于：步骤2所述的将所有高频方向子带合并为一个高频子带是通过如下公式进行：

$D_{h,s}=\frac{1}{N_s}\underset{b\&Element;N_s}{\mathrm{\&Sigma;}}{\left|D_{b,s}\right|}^2$

其中，D_h，s表示第s层的高频子带，下标h表示为高频，分解层s＝1，2，3，4，N_s表示高频方向子带个数，s＝1和2时，N_s＝10，s＝3和4时，N_s＝6。

3.根据权利要求1所述的遥感图像变化检测方法，其特征在于：步骤6所述的得到MRF分类图的步骤包括：

(6a)将步骤3所得的某一层的低频初始分类图作为该层的初始标记图；如将第s层的低频初始分类图作为第s层的初始标记图，变化类像素

用1表示，非变化类像素

用0表示，设迭代次数k＝1，最终迭代次数为K，K的取值范围为K＞＝5；

(6b)第s层低频子带D_l，s中像素点(i，j)处变化类的似然能量

和非变化类的似然能量由二维高斯概率密度的负自然对数计算得到，即

${\mathrm{Uc}}_{\mathrm{data}}^{l,s}(i,j)=\frac{1}{2}\ln \left|2\mathrm{\&pi;}{\left({\mathrm{\&sigma;}}_{l,c}^s\right)}^2\right|+\frac{1}{2}{(d_l^s(i,j)-{\mathrm{\&mu;}}_{l,c}^s)}^2{\left[{\left({\mathrm{\&sigma;}}_{l,c}^s\right)}^2\right]}^{-1};$

${\mathrm{Uu}}_{\mathrm{data}}^{l,s}(i,j)=\frac{1}{2}\ln \left|2\mathrm{\&pi;}{\left({\mathrm{\&sigma;}}_{l,u}^s\right)}^2\right|+\frac{1}{2}{(d_l^s(i,j)-{\mathrm{\&mu;}}_{l,u}^s)}^2{\left[{\left({\mathrm{\&sigma;}}_{l,u}^s\right)}^2\right]}^{-1}$

(6c)计算第s层低频子带内像素点(i，j)处变化类的先验能量

和非变化类的先验能量

${\mathrm{Uc}}_{\mathrm{context}}^{l,s}(i,j)=\underset{(p,q)\&Element;N(i,j)}{\mathrm{\&Sigma;}}\mathrm{\&delta;}(C^{l,s}(i,j),C^{l,s}(p,q));$

${\mathrm{Uu}}_{\mathrm{context}}^{l,s}(i,j)=-8-{\mathrm{Uc}}_{\mathrm{context}}^{l,s}(i,j)$

其中，C^l,s(i，j)为第s层低频子带内像素点(i，j)处的类别标记，N(i，j)为像素点(i，j)的二阶邻域，N(i，j)＝{(i，j)+(v，ζ)，(v，ζ)∈N}，N＝{(±1，0)，(0，±1)，(1，±1)，(-1，±1)}，C^l,s(p，q)为二阶邻域N(i，j)中(p,q)处的类别标记，p和q分别为二阶邻域N(i，j)中的行序号和列序号，δ(C^l，s(i，j)，C^l,s(p，q))为狄拉克函数，其表达式为

(6d)用第s层低频子带内像素点(i，j)处变化类的似然能量和先验能量

计算第s层低频子带内像素点(i，j)处变化类的总能量

用非变化类的似然能量

和先验能量

计算非变化类的总能量 $U_u^{l,s}(i,j):$

$u_c^{l,s}(i,j)={\mathrm{Uc}}_{\mathrm{data}}^{l,s}(i,j)+\mathrm{\&beta;}{\mathrm{Uc}}_{\mathrm{context}}^{l,s}(i,j)$

$u_u^{l,s}(i,j)={\mathrm{Uu}}_{\mathrm{data}}^{l,s}(i,j)+\mathrm{\&beta;}{\mathrm{Uu}}_{\mathrm{context}}^{l,s}(i,j)$

其中，参数β为常数，β＞0，用来调整先验能量和似然能量在变化检测过程中的重要性；

(6e)若

则将像素点(i，j)处的类别标记更新为变化类，否则，更新为非变化类，对图像所有像素点做同样处理，得到更新的类别标记图像，更新迭代次数k，令k＝k+1，当k＜＝K时，转到(6b)；否则，转到(6f)；

(6f)将更新的类别标记图像作为第s层的MRF分类图M_s＝{M_s(i，j)|1≤i≤I，1≤j≤J}。

4.根据权利要求1所述的遥感图像变化检测方法，其特征在于：其中步骤9所述的将最终自适应阈值分类图B和最终MRF分类图M进行融合，融合步骤包括：

(9a)设融合循环次数z＝1，融合终止循环次数为Z，Z＞50；

(9b)检测最终MRF分类图M的像素值是否等于1，如果M(i，j)＝1，则将最终自适应阈值分类图B中(i，j)处及其二阶邻域N(i，j)中的值赋给融合图像CM的(i，j)处及其二阶邻域的对应位置，否则令CM(i，j)的值为0；

(9c)按照步骤(9b)从左到右从上到下检测M中的所有像素值，令融合图像CM中像素值大于0的像素值为1；

(9d)更新融合循环次数z，令z＝z+1，当z＞Z或两次迭代之间类别发生变化的像元总数小于预先设定的阈值T，T＝1为最佳阈值，则停止，融合图像CM为变化检测结果图；否则将融合图像CM赋给M，转到(9b)。

Images