CN101650439B - 基于差异边缘和联合概率一致性的遥感图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于差异边缘和联合概率一致性的遥感影像变化检测方法，主要解决传统方法在提高抗噪声能力时，无法较好地保持变化区完整性的问题。具体实现步骤为：对输入的两时相影像采用灰度一致性校正，并创建联合概率一致性图；求取两时相原始图像的canny边缘图，并获得差异边缘图；确定引起配准噪声的差异边缘的位置；结合联合概率一致性构建新差异图；对新差异图取OTSU阈值，获得最终的变化检测结果。本分明在消除噪声的同时实现了变化区块空间一致性的保持，提高了的变化检测结果的精度，可用于多时相遥感影像的火灾、水灾等灾害的评估。

Description

基于差异边缘和联合概率一致性的遥感图像变化检测方法

技术领域

本发明属于图像处理领域，涉及遥感图像的变化检测方法，具体地说是一种基于差异边缘和联合概率一致性的遥感图像变化检测方法。

背景技术

遥感图像变化检测技术旨在通过分析同一地区不同时相的遥感图像，提取遥感图像中地物的变化信息，可应用于森林资源调查、土地利用/覆盖动态监测、环境灾害评估、城市规划检测和国防军情监控等对地观测中。

遥感图像变化检测技术本质上是将遥感图像数据分为变化类和未变化类的两分类问题。遥感图像含有如光照变化、大气环境变化、几何形变等多种因素引起的噪声，这些噪声影响了变化类与非变化类的可分性，降低了变化检测的精度，成为无监督变化检测方法研究的难点。

针对无监督变化检测方法中低数据可分性的问题，研究人员主要尝试从变化阈值的选取、空间关系的有效利用方面进行相关工作。变化阈值的选取方面的研究工作一般是先提出概率模型分布的假设，然后在此前提条件下结合一些阈值方法进行变化检测。例如Bruzzone等提出的基于EM算法的贝叶斯阈值自动选取方法，在一定的分布模型假设前提条件下，通过EM算法来估计变化类和非变化类的分布模型参数，然后应用统计模式识别中的贝叶斯决策理论加以解决，但EM算法在进行变化分类的参数估计时，存在收敛速度过慢和估计出现较大偏差等问题，影响了变化检测的性能。且该类方法都存在一个前提假设，认为数据符合高斯分布、混合高斯分布、瑞利高斯分布、Gamma分布等等概率分布模型，但实际获取的遥感数据由于地物等各种因素千差万别，简单的概率模型分布很难对多种变化和噪声均获得良好的检测结果，且目标区域的一致性无法很好的保持。变化检测中空间关系的利用方面，目前的研究工作一般是从空间上下文关系，如隐马尔科夫链、马尔科夫随机场(Markov random field，MRF)等方法着手，如2000年Bruzzone和Prieto提出利用贝叶斯最小错误率阈值进行初始分类并结合Markov随机场的变化检测方法，在考虑概率统计的同时充分考虑了局部邻域的一致性，取得了较好的效果。但是能量函数的选择、区域边界的各向异性、迭代算法寻找最优解时的问题等等使其在达到较好的检测结果方面应用不便。

目前也有一些基于边缘信息的变化检测方法，但都是基于图像中的目标边缘检测分割的，检测的结果很大程度上依赖于目标边缘检测的精度，要求目标(如建筑物、道路、汽车等)边缘比较突出并对检测出目标区边缘的连续性要求较高。

本发明的目的在于克服变化检测中输入图像数据难以有效区分的不足，针对以上变化检测方法中的不足之处，提出了一种基于差异边缘和联合概率一致性的遥感图像变化检测方法，以减少配准噪声的影响，同时充分利用了联合概率一致性对变化区块边缘细节的保持，取得了较好的变化检测结果。

实现本发明目的的技术方案是：首先，准确定位由配准噪声引起的差异边缘，即利用差异边缘并结合联合概率一致性确定未变化区域差异边缘的位置；其次，对两时相图像中的未变化区差异边缘，做加权均值滤波以减少差异图中的噪声干扰；最后，利用处理后得到的差异图，并结合联合概率一致性图，进行阈值分割，获得精确的变化检测结果。具体步骤包括如下：

(1)对输入的两时相图像采用灰度一致性校正；

(2)取步骤(1)校正后图像I₁和I₂中位于(i，j)点处的灰度值v₁和v₂，组成灰度点对(v₁，v₂)，并统计图像I₁和I₂中点对(v₁，v₂)出现的概率P(v₁，v₂)，然后构建一幅与图像I₁和I₂相同大小的图像I_P，对应I₁和I₂出现灰度点对(v₁，v₂)的位置处的I_P的像素值赋为概率P(v₁，v₂)，得到一幅联合概率一致性图I_P；

(3)求取步骤(1)校正后图像I₁和I₂的canny边缘图，并对canny边缘图相减后取绝对值，获得差异边缘图I_D；

(4)从差异边缘图中粗提取可能噪声边缘；

(5)对图像I₁和I₂中1~2像素宽噪声边缘上的点和3像素宽噪声边缘上的点分别采用5*5的加权均值滤波和7*7的加权均值滤波，最后将滤波后的两时相图像相减后取绝对值，得到一幅滤波后差异图I_F；

(6)采用Otsu阈值法将滤波后差异图分成变化区和非变化区，依据联合概率空间一致图I_P的空间一致性对其中的变化区进行区域修正，得到初始分割结果图；

(7)对步骤(1)校正后图像I₁和I₂中处于初始分割结果图非变化区的可能噪声边缘上的点采用k×k加权模板均值滤波，并对滤波后的两时相图像相减后取绝对值，得到一幅修正后差异图I_M；

(8)对步骤(7)修正后差异图I_M处于初始分割结果图变化区的点采用5×5像素大小的加权模板均值滤波，对非变化区的点采用5×5像素大小的八方向加权均值滤波，将加权均值滤波后差异图中灰度值大于滤波前灰度值的点，用加权滤波前的灰度值替代加权滤波后的灰度值，得到一幅加权方向滤波差异图；

(9)根据Otsu阈值将步骤(8)得到的差异图分为变化类与非变化类，获得最终的变化检测结果；

本发明由于将差异边缘信息与联合概率的空间一致性相结合，创建出了新差异图，与现有技术相比较具有如下优点：

a、本发明对配准误差和雷达入射角度不同引起的投影干扰具有一定的鲁棒性。

b、本发明不受变化区域形状和边缘联通性的限制。

c、本发明对于变化区域内部的一致性保持较好。

d、本发明对变化区域边缘的细节信息保持较好。

e、提高了变化检测的精度。

附图说明

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

图2是本发明联合概率一致性图创建的示意图；

图3是本发明的8个方向加权均值滤波权值模板的示意图；

图4是本发明的第一组实验数据和结果图；

图5是本发明的第二组实验数据和结果图；

具体实施方式

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

步骤1，对输入的两时相图像采用灰度一致性校正，得到图像I₁和I₂。

输入图像为同一地理位置处已配准的、大小为M×N的两个不同时相的原始遥感图像，分别记为和

本发明的灰度一致性校正是对时相1图像

保持不变，并以此图像为基础对时相2图像进行灰度校正，具体公式如下：

$I_2=\left\{\begin{array}{cc}\frac{{\mathrm{\σ}}_1}{{\mathrm{\σ}}_2}({\tilde{I}}_2-{\widetilde{\mathrm{\μ}}}_2)+{\widetilde{\mathrm{\μ}}}_1& \mathrm{if}({\widetilde{\mathrm{\μ}}}_1-{\widetilde{\mathrm{\μ}}}_2)\·({\widetilde{\mathrm{\σ}}}_1-{\widetilde{\mathrm{\σ}}}_2)>0\\ {\tilde{I}}_2& \mathrm{others}\end{array}---\left(1\right)\right.$

其中

和

表示两个时相图像的标准方差，

和

表示两时相图像的均值，I₂为校正后的时相2图图像。

步骤2，对经过灰度一致性校正的两时相图像I₁和I₂，创建联合概率一致性图I_P。

如图2所示，取已经校正图像I₁和I₂中位于(i，j)点(其中i＝1，2…M，j＝1，2…N)处的灰度值v₁和v₂，组成灰度点对(v₁，v₂)，并统计图像I₁和I₂中点对(v₁，v₂)出现的概率P(v₁，v₂)；然后构建一幅与图像I₁和I₂相同大小的图像I_P，对应I₁和I₂出现灰度点对(v₁，v₂)的位置处的I_P的像素值赋为概率P(v₁，v₂)，得到一幅联合概率一致性图I_P。

步骤3，求取图像I₁和I₂的canny边缘图，将两时相canny边缘图相减后取绝对值，获得一幅差异边缘图I_D。

对经过步骤(1)处理后的两时相图像I₁和I₂利用canny边缘检测算子检测得到两时相图的canny边缘图，其中边缘点的值为1，非边缘点值为0；然后将两时相的canny边缘图相减，对其差值图再取绝对值，得到一幅差异边缘图I_D。

步骤4，从差异边缘图中提取可能噪声边缘。

利用5×5像素大小的滑窗从差异边缘图I_D中提取符合以下条件的边缘，本发明称之为可能噪声边缘。所要求的边缘条件符合以下三个之一即可：①2个素宽的边缘；②由相隔1个像素的错位边缘及其中间区域所组成的3个像素宽的边缘；③与2个像素或者与3个像素宽的边缘存在四连通的1个像素宽的边缘；

步骤5，对图像I₁和I₂中1~2像素宽噪声边缘上的点和3像素宽噪声边缘上的点分别采用5*5的加权均值滤波和7*7的加权均值滤波，最后将滤波后的两时相图像相减后取绝对值，得到一幅滤波后差异图I_F。

5a)对于两时相图像I₁和I₂中可能噪声边缘上的点，采用k×k像素大小的加权模板。这里k×k加权均值滤波模板是指：对两时相图像I₁和I₂所取的当前k×k像素大小的窗口计算其联合概率一致性图I_P对应位置处窗口内的均值μ_w和方差σ_w，将联合概率一致性图I_P的当前窗口内像素值不在(μ_w-2×σ_w，μ_w+2×σ_w)范围内的点置0，当前窗口内剩余像素值作归一化处理，得到k×k窗口的加权模板；

5b)对两时相图像I₁和I₂中可能噪声边缘上的点，依据k×k窗口的加权模板进行均值滤波。当噪声边缘为1～2个像素宽时滤波窗口大小k取值为5；当噪声边缘为3个像素宽时滤波窗口大小k取值为7。

5c)对滤波后的两时相图像相减的差值取绝对值，得到一幅滤波后差异图I_F；

步骤6，对滤波后差异图采用Otsu取阈法得到变化区和非变化区，对其中的变化区进行基于联合概率空间一致图I_P的区域修正，得到初始分割结果图。

对滤波后差异图采用Otsu取阈法得到变化区和非变化区，统计变化区域内各变化区块在联合概率一致性图I_P中对应位置区块的均值μ_r和方差σ_r，对其中属于(μ_r-2×σ_r，μ_r+2×σ_r)灰度范围内且与当前变化区块存在四连通的像素点标记为变化点，得到修正变化区后的图像。

步骤7，然后依据此修正的变化区对差异边缘图I_D中的边缘分类为变化区边缘和非变化区边缘，对非变化区边缘上的点再次用k×k加权模板均值滤波，并对滤波后的两时相图像相减后取绝对值，得到一幅修正后差异图I_M。

7a)将步骤(6)所得图像与差异边缘图ID做与运算，将差异边缘图中的边缘分为变化区边缘和非变化区边缘两类。

7b)将两时相图像I₁和I₂的非变化区边缘上的点，再次计算其k×k窗口的加权模板并进行均值滤波。当可能噪声边缘为1～2个像素宽时滤波窗口大小k取值为5；当可能噪声边缘为3个像素宽时滤波窗口大小k取值为7。

7c)将步骤(7b)处理后的两时相图像相减的差值取绝对值，得到修正后差异图I_M。

步骤8，对步骤(7)修正后差异图I_M处于初始分割结果图变化区的点采用5×5像素大小的加权模板均值滤波，对非变化区的点采用5×5像素大小的八方向加权均值滤波，将加权均值滤波后差异图中灰度值大于滤波前灰度值的点，用加权滤波前的灰度值替代加权滤波后的灰度值，得到一幅加权方向滤波差异图。这里5×5像素大小的八方向加权均值滤波的方向加权模板如图3所示，各个方向模板的中心位置处灰色点为当前像素点，黑色点为计算方向加权的点，其权值为联合概率一致性图中当前方向模板中所有黑色点对应值的归一化，模板内其他点的权值为零。

步骤9，根据Otsu阈值将步骤(6)得到的差异图分为变化类与非变化类，获得最终的变化检测结果。

本发明的效果通过以下实验具体说明：

本发明使用了两组真实遥感图像。第一组真实实验数据集是通过Landsat 7ETM+4波段获得的墨西哥郊外的两幅大小为512×512的真实图像，第一时相图像是在2000年4月获得，如图4(a)所示；第二时相SAR图像是在2002年5月获得，如图4(b)所示，变化区域主要为大火破坏了大面积的当地植被所致。图4(c)为图4(a)和图4(b)变化检测结果参考图，图4(d)为本发明关于图4(a)和图4(b)的变化检测结果。第二组真实实验数据集是通过Landsat 5 TM 4波段获得的意大利Elba岛的两幅大小为414×326的真实图像，第一时相图像是在1992年8月获得，如图5(a)所示；第二时相SAR图像是在1994年8月获得，如图5(b)所示，变化区域主要是野火破坏的植被。图5(c)为图5(a)和图5(b)的变化检测结果参考图，图5(d)为本发明关于图5(a)和图5(b)的变化检测结果。

本发明设计了两个实验，第一个实验比较了直接差值法手动阈值的最优结果和本发明的结果，用于证明去噪声的必要性及发明中去噪声的有效性；第二个实验比较了引入空间一致性和未引入一致性的结果，用于说明本发明提出的空间一致性的必要性和有效性。

表1给出了实验一中第一组实验数据结果技术指标，表2给出了实验一中第二组实验数据结果技术指标。从表1和表2可以看出，对两组真实的遥感图像，本发明的总错误数比直接差值法的总错误数分别少了25％和93％，明显的减少了噪声的影响。表3给出了实验二中第一组实验数据结果技术指标，表4给出了实验二中第二组实验数据结果技术指标。从表3和表4可以看出，本发明引入概率空间一致性后，在保持误检测点数基本不变的情况下，总错误数减分别减少了7.9％和42％，说明变化区的完整性得到了保持。总的来说，本发明实现了噪声的消除，同时较好的保持了变化区的完整性。

表1实验一的第一组实验数据结果

表2实验一的第二组实验数据结果

表3实验二的第一组实验数据结果

表4实验二的第二组实验数据结果

Claims (6)

1.一种基于差异边缘与联合概率一致性的遥感图像变化检测方法，包括如下步骤：

1)对输入的两时相图像采用灰度一致性校正；

2)取步骤1)校正后图像I₁和I₂中位于(i，j)点处的灰度值v₁和v₂，组成灰度点对(v₁，v₂)，并统计图像I₁和I₂中点对(v₁，v₂)出现的概率P(v₁，v₂)，然后构建一幅与图像I₁和I₂相同大小的图像I_P，对应I₁和I₂出现灰度点对(v₁，v₂)的位置处的I_P的像素值赋为概率P(v₁，v₂)，得到一幅联合概率空间一致性图I_P；

3)求取步骤1)校正后图像I₁和I₂的canny边缘图，并对canny边缘图相减后取绝对值，获得差异边缘图I_D；

4)从差异边缘图中粗提取可能的噪声边缘；

5)对图像I₁和I₂中1-2像素宽噪声边缘上的点和3像素宽噪声边缘上的点分别采用5*5的加权均值滤波和7*7的加权均值滤波，最后将滤波后的两时相图像相减后取绝对值，得到一幅滤波后差异图I_F；

6)采用Otsu阈值法将滤波后差异图分成变化区和非变化区，依据联合概率空间一致图I_P的空间一致性对其中的变化区进行区域修正，得到初始分割结果图；

7)对步骤1)校正后图像I₁和I₂中处于初始分割结果图非变化区的噪声边缘上的点采用k×k加权模板均值滤波，并对滤波后的两时相图像相减后取绝对值，得到一幅修正后差异图I_M；

8)对步骤7)修正后差异图I_M处于初始分割结果图变化区的点采用5×5像素大小的加权模板均值滤波，对非变化区的点采用5×5像素大小的八方向加权均值滤波，将加权均值滤波后差异图中灰度值大于滤波前灰度值的点，用加权滤波前的灰度值替代加权滤波后的灰度值，得到一幅加权方向滤波差异图；

9)根据Otsu阈值将步骤8)得到的差异图分为变化类与非变化类，获得最终的变化检测结果。

2.根据权利要求1所述遥感图像变化检测方法，其中步骤1)所述的对输入的两时相图像采用灰度一致性校正，是将同一地理位置处已配准的大小为M×N的两个时相的原始遥感图像，分别记为 

和 对时相1图像 

保持不变，并以此图像为基础，依据如下公式对第二时相 进行校正： 

其中 

和 

表示两个时相图像的标准方差， 

和 

表示两时相图的均值，I₂表示校正后的第二时相图。

3.根据权利要求1所述的遥感图像变化检测方法，其中步骤4)所述的从差异边缘图中提取可能噪声边缘，是利用5×5像素大小的滑窗从差异边缘图I_D中选取出三种边缘：第一种是2个像素宽的边缘；第二种是由相隔1个像素的错位边缘及其中间区域所组成的3个像素宽的边缘；第三种是与2个像素宽边缘和3个像素宽的边缘存在四连通的1个像素宽的边缘。

4.根据权利要求1所述的遥感图像变化检测方法，其中步骤7)所述k×k加权模板均值滤波，是将联合概率空间一致性图中当前k*k窗口内(μ-2*σ，μ+2*σ)范围内的点作归一化处理，并将其他点置零得到权值模板，计算加权均值得到当前点的加权均值滤波结果，其中μ和σ分别为窗口内点的均值和方差。

5.根据权利要求1所述的遥感图像变化检测方法，其中步骤6)所述采用Otsu阈值法将滤波后差异图分成变化区和非变化区，并依据联合概率空间一致图I_P的空间一致性对其中的变化区进行区域修正，是依次统计变化区内各变化区块在联合概率空间一致性图中对应区块的均值μ_r和方差σ_r，对于联合概率空间一致性图中属于(μ_r-2*σ_r，μ_r+2*σ_r)灰度范围内且与当前变化区块四联通的点标记为变化点，得到初始分割结果图。

6.根据权利要求1所述的遥感图像变化检测方法，其中步骤8)所述5×5像素大小的八方向加权均值滤波，是以联合概率空间一致性图中当前方向所有黑色点对应值的归一化值为权值，计算各方向黑色点的加权均值，选取与中心点值最近的加权均值为八方向加权均值滤波结果。 

Images