基于空间关联条件概率融合的多时相SAR图像变化检测方法

摘要

本发明公开了一种基于空间关联条件概率融合的多时相SAR图像变化检测方法，主要解决传统变化检测方法在保持变化区域的完整性时，伪变化区域较多的问题，具体实现步骤为：(1)对增强Lee滤波后的两时相SAR图像进行分割；(2)创建条件概率差异图并求得用于融合的空间关联条件概率；(3)利用空间关联条件概率融合两时相分割结果，得到融合差异图；(4)对融合差异图取阈值得到变化检测结果。本发明在较好地保证变化区域完整性的同时，减少了伪变化区的出现，提高了变化检测的精度，可用于多时相SAR图像的火灾、水灾等灾害的评估。

说明书

技术领域

本发明属于图像处理领域，涉及多时相SAR图像的变化检测，具体地说是一种基于空间关联条件概率融合的多时相SAR图像变化检测方法。

背景技术

随着卫星遥感技术的发展，遥感数据飞速增加。其中的SAR数据由于不受大气环境和和云层遮挡的影响成为一种重要的遥感数据。如何准确快速自动地从不同时相的SAR图像中找出显著变化的区域则具有十分重要的意义。

在SAR图像中寻找“非变化”和“变化”两分类的研究目前还处于初步阶段，大致有两个路线，一条路线是分类后比较方法，也称后分类比较法，该方法先对两个时相的图像进行独立分类，再对两幅分类图像进行逐象素的比较，最后得到变化检测图；另外一条路线是差异图分类方法，该方法先对两个时相的图像进行逐象素的比较，比如逐象素点的差值、比值、CVA等，再对比较得出的差异图像进行进一步地多种变换、概率分布等处理以达到两分类，最后得到变化检测图。后分类比较法能够减少由于数据获取平台和环境的不同引起的伪变化信息，不需要数据的辐射校正等复杂的预处理，但是该方法是基于不同时相的分类图的，对分割结果准确性要求较高，象素值的微小偏差可能会被判为不同的类别，经差异比较后会产生非变化与变化间的巨大差异，直接影响最终的变化检测结果。目前较多的研究是沿着差异分类模型这个路线进行的。差异图分类方法简单直观，得到的变化细节较为显著。构造差异影像得到的结果与实际的变化和非变化趋势大体一致，但是它们之间反映出的地物的辐射能量的变化却并不总是相同的，因此反映出的地物变化程度及属性也不完全相同，往往含有较多的伪变化信息。同时该方法将两时相图像直接进行逐象素的差值或者比值运算并没有充分利用两时相图像的内在关联特性，无法较好的保持变化区块的完整性。

发明内容

本发明的目的在于针对上述两类变化检测方法的不足，提出一种基于空间关联条件概率融合的SAR图像变化检测方法，以较好的保持变化区块的完整性，并减少伪变化区块的出现，取得了较好的变化检测结果。

实现本发明目的的技术方案是：首先分别计算两时相SAR图像的边缘匹配最优分割图和空间关联条件概率图，然后逐点计算条件概率图与两时相SAR图像中对应3×3邻域的空间关联条件概率，据此条件概率逐点对两时相的最优分割图进行融合，得到融合差异图，最后对融合差异图分割得到变化检测结果。具体步骤如下：

(1)对输入的两时相图像分别采用两次增强Lee滤波，得到图像I₁和I₂；

(2)用Canny边缘算子分别计算两时相图像I₁和I₂的边缘，得到两幅Canny边缘图；

(3)计算两幅Canny边缘图中每个像素点与其周围每个边缘点的欧氏距离，选择每个像素与其周围边缘点的最小的距离作为该像素点的边缘距离值，得到两幅Canny边缘图中各点的边缘距离值；

(4)取图像I₂的每一个灰度级作为一个阈值对图像I₂分别二值化，得到与各灰度级对应的多幅二值图；计算每个二值图的边缘与图像I₂的Canny边缘的匹配度m₁(g₂)/m₂(g₂)，其中，m₁(g₂)为LE₂(g₂)/BE₂(g₂)，BE₂(g₂)是与灰度级g₂对应的二值图的边缘上的像素点数，LE₂(g₂)是二值图边缘与时相2图像I₂的Canny边缘存在重叠的像素点数，m₂(g₂)是灰度级g₂对应二值图的边缘上点在第二幅Canny边缘图中对应点的边缘距离值的和；选取匹配度最大的二值图作为图像I₂的最优分割图IE₂；

(5)取图像I₁的每一个灰度级作为一个阈值对图像I₁分别二值化，得到与各灰度级对应的多幅二值图；计算每个二值图的边缘与图像I₁的Canny边缘的匹配度m₃(g₁)·m₁(g₁)/m₂(g₁)，其中，m₁(g₁)为LE₁(g₁)/BE₁(g₁)，BE₁(g₁)是与灰度级g₁对应的二值图的边缘上的像素点数，LE₁(g₁)是二值图边缘与时相1图像I₁的Canny边缘图中的边缘存在重叠的像素点数，m₂(g₁)是灰度级g₁对应二值图边缘上点在第一幅Canny边缘图中对应点的边缘距离值的和，m₃(g₁)为第一时相图像I₁灰度级g₁对应的二值图的边缘与第二时相图像I₂的最优分割图IE₂边缘两者存在空间位置重叠的像素点数除以第二时相图像I₂的最优分割图IE₂边缘的像素点数的值；选取出匹配度最大的二值图作为图像I₁的最优分割图IE₁；

(6)取图像I₁和I₂中位于(i，j)处像素的灰度值v₁和v₂，组成灰度点对(v₁，v₂)，统计图像I₁和I₂中灰度点对(v₁，v₂)出现的概率P(v₁，v₂)及灰度值v₁出现的概率P(v₁)；

(7)构建一幅与图像I₁和I₂相同大小的M×N图像I_P，对应I₁和I₂出现灰度点对(v₁，v₂)的位置处的I_P的像素值赋为概率P(v₁，v₂)/P(v₁)，得到空间关联条件概率图I_P；

(8)统计空间关联条件概率图I_P中每个像素点与图像I₁和I₂中相同位置像素的3×3邻域内各点的条件概率其中，v_P(i，j)为I_P中(i，j)处点的值，v(s)为图像I₁和I₂中(i，j)处点3×3邻域内18个点的值，N_h，l(v_P(i，j)，v(s))为(v_P(i，j)，v(s))在由I_P中(h，l)处像素值与图像I₁和I₂中(h，l)处像素的3×3邻域内18个像素值组成的18组灰度对中出现的频数，P(v(s))为图像I₁和I₂中值v(s)出现的概率；

(9)依据步骤(8)图像I₁和I₂中(i，j)处3×3邻域内各点与I_P中(i，j)处的点的条件概率，对应地将步骤(4)和步骤(5)中两时相图像的最优分割结果图IE₁和IE₂中(i，j)处3×3邻域内各点融合，得到融合后差异图(i，j)处的值，得到一幅融合差异图I_F；

(10)利用Otsu阈值将步骤(9)所得融合差异图I_F分为变化类与非变化类，得到变化检测结果。

本发明由于采用基于边缘匹配的分割方法分割两时相图，并利用空间关联条件概率将两时相最优分割图融合，创建出了融合差异图，与现有技术相比具有如下优点：

(1)本发明实现了边缘硬匹配与软匹配的结合，减少了其他分割方法中过分割的现象，取得了较好的分割结果。

(2)本发明既充分考虑了两时相原图像的区域一致性，又保持了空间关联条件概率图相对于原图像的概率关系，达到了在尽量保证区域一致性的同时，减少检测出的伪变化区域的目的。

附图说明

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

图2是本发明条件概率差异图的创建流程示意图；

图3是本发明空间关联条件概率构建示意图；

图4是本发明基于空间关联条件概率的最优分割图四方向融合示意图；

图5是模拟SAR实验数据的Canny边缘匹配曲线图；

图6是模拟SAR实验数据及结果图；

图7是真实SAR实验数据及结果图。

具体实施方式：

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

步骤1，输入两时相图像，如图6(a)和图6(b)，对输入的两时相图像均采用两次增强Lee滤波，得到图像I₁和I₂，其中，增强Lee滤波的窗口大小为5×5像素。

步骤2，用Canny边缘算子分别计算两时相图像I₁和I₂的边缘，得到两幅Canny边缘图。

步骤3，依次取图像I₂的每一个灰度级作为阈值对图像I₂二值化，得到与各灰度级对应的多幅二值图，计算每个二值图的边缘与图像I₂的Canny边缘的匹配度m₁(g₂)/m₂(g₂)，选取匹配度最大的二值图作为图像I₂的最优分割图IE₂，其中，

m₁(g₂)定义为LE₂(g₂)/BE₂(g₂)，该BE₂(g₂)是与灰度级g₂对应的二值图的边缘上的像素点数，该LE₂(g₂)是二值图边缘与时相2图像I₂的Canny边缘存在重叠的像素点数；

计算图像I₂的Canny边缘图中每个像素点与其周围每个边缘点的欧氏距离，选择其中最小的距离作为该点的像素值，得到图像I₂的边缘距离图；

m₂(g₂)定义为灰度级g₂对应二值图的边缘上点在图像I₂的边缘距离图中对应点的边缘距离值的和。

步骤4，分别取图像I₁的每一个灰度级g₁作为阈值对图像I₁二值化，得到与各灰度级对应的多幅二值图，计算每个二值图的边缘与图像I₁的Canny边缘的匹配度m₃(g₁)·m₁(g₁)/m₂(g₁)，选取出匹配度最大的二值图作为图像I₁的最优分割图IE₁，其中，

m₁(g₁)定义为LE₁(g₁)/BE₁(g₁)，BE₁(g₁)是与灰度级g₁对应的二值图的边缘上的像素点数，LE₁(g₁)是二值图边缘与图像I₁的Canny边缘图中的边缘存在重叠的像素点数；

计算图像I₁的Canny边缘图中每个像素点与其周围每个边缘点的欧氏距离，选择其中最小的距离作为该点的像素值，得到图像I₁的边缘距离图；

m₂(g₁)定义为灰度级g₁对应二值图边缘上点在图像I₁的边缘距离图中对应点的边缘距离值的和；

m₃(g₁)定义为第一时相图像I₁灰度级g₁对应的二值图的边缘与第二时相图像I₂的最优分割图IE₂边缘两者存在空间位置重叠的像素点数除以第二时相图像I₂的最优分割图IE₂边缘的像素点数的值。

步骤6，空间关联条件概率图I_P。

参照图2，本步骤的具体实现如下：

首先，取步骤(1)滤波后图像I₁和I₂中位于(i，j)点处的灰度值v₁和v₂，组成灰度点对(v₁，v₂)，并统计图像I₁和I₂中点对(v₁，v₂)出现的概率P(v₁，v₂)及灰度值v₁出现的概率P(v₁)；

然后，构建一幅与图像I₁和I₂相同大小的M×N图像I_P，对应I₁和I₂出现灰度点对(v₁，v₂)的位置处的I_P的像素值赋为概率P(v₁，v₂)/P(v₁)，得到空间关联条件概率图I_P。

步骤7，统计空间关联条件概率图I_P中每个像素点对应图像I₁和I₂中相同像素位置的3×3邻域内各点的条件概率，如图3所示，具体步骤如下：

7a)记图像I₁和I₂中(i，j)处点的3×3邻域内的18个点的像素值构成集合其中图像I₁的像素值集合为图像I₂的像素值集合为

7b)空间关联条件概率图I_P中(i，j)处点的像素值为v_P，与该点对应的图像I₁和I₂的以(i，j)为中心的3×3像素大小范围内的共18个点的像素值集合为v^S，统计其中与v_P组成的像素值对在像素对集合(v_P，v^S)中出现的频数并在空间关联条件概率图I_P中统计该像素值对出现的总频数得到出现的概率

7c)统计灰度值在I₁和I₂两幅图中出现的总概率则可以得到I_P中(i，j)处灰度值v_P相对于I₁和I₂中(i，j)处3×3邻域内各灰度值的条件概率：$P(v_P/v^{S_l})=P(v_P,v^{S_l})/P\left(v^{S_l}\right);$

步骤8，依据步骤(6)图像I₁和I₂中(i，j)处3×3邻域内各点与I_P中(i，j)处点的条件概率，对应地将步骤(2)中两时相图像的最优分割图IE₁和IE₂中(i，j)处3×3邻域内各点融合，得到融合后差异图(i，j)处的值，得到一幅融合差异图I_F，具体步骤如下：

8a)分别选出两时相图像的最优分割图IE₁和IE₂在(i，j)点处的3×3邻域内0°、45°、90°、135°方向上大于各自最优阈值的像素点，计算四个方向上大于最优阈值的像素点在各自方向所占比例r₁(ori)和r₂(ori)及其在各自方向上的对应条件概率的均值u₁(ori)和u₂(ori)，其中ori分别为45°，90°，135°；

8b)分别计算两时相图的最优分割图IE₁和IE₂在(i，j)点处四个方向上的比例与对应方向条件概率的均值的乘积；

8c)将两时相的四个方向的乘积取平均求和得到(i，j)点处的融合像素值，得到一幅融合差异图I_F。

步骤9，利用Otsu阈值将步骤(8)所得融合差异图I_F中值大于OTSU阈值的点分成变化类，其他点分为非变化类，得到变化检测结果。

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

本发明分别使用了模拟SAR和真实SAR数据两组实验数据。

第一组模拟SAR数据是按照Gibbs随机场的分布形式获得具有两类场景的参考图像，然后按照SAR图像的成像机理获得带有相干斑的SAR图像，最后将仿真的变化区域嵌入到参考图像中得到的。图6给出了模拟SAR数据及其实验结果，其中图6(a)为模拟SAR数据在第一时相的原图，图6(b)为模拟SAR数据在第二时相的原图，图6(c)为第一时相图的最优分割图，图6(d)为第二时相图的最优分割图，图6(e)为两时相的变化检测参考图，图6(f)为本发明的变化检测结果图。

第二组的真实SAR数据是1999年4月和5月由ERS 2所携带的SAR在瑞士Bern城区附近获得。该组数据变化部分是在Thun市、Bern市和Bern机场段Aare河的洪水引起的。实验中所用图像大小为301×301像素，256灰度级，其中变化像元数为1155，未变化像素元数为89446。图7给出了真实SAR数据及实验结果，图7(a)为真实SAR第一时相的原图，图7(b)为真实SAR第二时相的原图，图7(c)为第一时相图的最优分割图，图7(d)为第二时相图的最优分割图，图7(e)为两时相的变化检测参考图，图7(f)为本发明的变化检测结果图。

本发明通过与对数比值法比较，利用误检率、漏检率和正检率三个变化检测领域的经典指标来说明本发明较好地保持了变化区块的完整性，并减少了伪变化区块的出现，取得了较好的变化检测结果。

表1分别给出了本发明和对数比值法的实验结果。从表1可以看出，相对于对数比值法，本发明对模拟SAR和真实SAR数据，误检率和漏检率都有了明显的降低，正检率明显提高，总的来说，变化检测精度有了很大提高。对于模拟SAR数据，本发明的误检率相对于对数比值法从25.24％降低到6.93％，说明伪变化区块出现率减少了17个百分点，漏检率从11.304％降低到6.93％，说明变化区域的完整性得到较好的保持，同时正检率从88.69％提高到95.063％，变化检测精度有了很大提高。对于真实SAR数据，本发明在保持误检率变化不大的情况下，漏检率大幅降低，正检率明显提高，变化检测精度有了很大提高。

表1SAR数据实验结果