CN108764119A - 基于迭代最大类间方差的sar图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于迭代最大类间方差的SAR图像变化检测方法，主要解决现有技术对于变化区域占整体区域的小部分情况下检测会产生过多虚警点的问题，其实现方案是：1.获取不同时相同一地区的SAR图像并进行预处理；2.生成对数比差异图；3.迭代使用最大类间方差法遍寻差异图查找变化元素与非变化元素相对均匀的区域；4.将这一区域的分割阈值应用到全图，得到分割结果图；5.对分割结果图进行中值滤波，得到最终变化检测结果图。本发明解决了直方图为单峰时最大类间方差法虚警过多的问题，提高了运行效率较高，可用于快速SAR图像变化检测，尤其是对于小面积变化区域的检测。

Description

基于迭代最大类间方差的SAR图像变化检测方法

技术领域

本发明属于图像处理技术领域，特别涉及SAR图像变化检测，可用于快速SAR图像变化检测，尤其是对于小面积变化区域的检测。

背景技术

合成孔径雷达SAR因其具有全天时、全天候的工作能力，且不受天气和光照的影响等优点，近年来得到了快速的发展。SAR图像变化检测就是获取不同时刻同一区域的变化情况，在环境监测、农业调查、城市研究、灾害监测等方面有着广泛的应用。尤其是自然灾害，如火灾，洪灾……发生时，高效的变化检测技术能够及时的避免或减少人身和财产的损失。

常见的变化检测方法最主要分为两类，一类是对图像分类后进行变化检测；另一类是对影像直接进行变化检测。为了避免因分类不精准而累积的误差，目前大部分变化检测工作都是直接针对影像进行研究。对影像直接进行变化检测包括3个部分：图像预处理、差异图生成与差异图分析。

差异图的分析方法常用的有四种：阈值分析、聚类分析、图切分析和水平集分析。最常用的差异图分析方法为阈值法或者聚类法。其中阈值分割算法的重点是阈值的确定。常用的阈值技术有全局阈值法，自适应阈值法以及最佳阈值发等等。阈值法由于其运算速度快，方法简明等优点得到了广泛应用。然而目前已有的全局阈值分割方法如最大类间方差法OSTU、峰谷法、迭代式阈值选择这些方法对于变化区域只占小部分的情况并不能进行有效分割。因为变化区域仅占整体区域的小部分，这表现到直方图上就是差异图像的直方图只出现一个明显的波峰，直方图将不具有明显的可分离性，最大类间方差法获取的自动阈值会偏向方差较小的一类也就是非变化类，进而导致分割结果虚警点过多，将大量非变化类元素划分到变化类。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于迭代最大类间方差的SAR图像变化检测方法，以解决在变化区域占整体区域比例低于10％时，使用最大类间方差法对差异图进行分割产生过多虚警点的问题，提高变化检测精度。

为实现上述技术目的，本发明的技术方案包括如下：

(1)通过合成孔径雷达获取在不同时相拍摄的同一地区的两幅SAR图像，对这两幅SAR图像依次进行辐射校正、几何校正、图像配准的预处理，获得预处理后的适合变化检测使用的两幅SAR图像I₁和I₂；

(2)计算第一幅SAR图像I₁在坐标(i,j)处的灰度值I₁(i,j)和第二幅SAR图像I₂在坐标(i,j)处的灰度值I₂(i,j)的比值，对比值运算的结果再进行一步对数运算，得到坐标(i,j)处的对数比值DI(i,j)：

DI(i,j)＝log(I₁(i,j)/I₂(i,j))，

遍历两幅SAR图像所有的像素点，得到对数比差异图：DI＝{DI(i,j)}；

(3)使用最大类间方差法对对数比差异图DI进行初分割，得到初步分割阈值T₁；

(4)用阈值T₁对对数比差异图DI进行二值分割得到初步分割图SI₁；

(5)设置大小为w*w的窗口对初步分割图SI₁进行滑窗，统计窗口内零元素与非零元素的个数，直至小窗口内零元素与非零元素个数差小于等于3时，停止滑窗，得到满足条件的窗口中心元素(p,q)，其中w*w窗口大小取值为待处理图像大小的0.5％-1.5％；

(6)取对数比差异图DI中中心像素为(p,q)，窗口大小为w*w的区域P，使用最大类间方差法获取区域P的分割阈值T，用阈值T对对数比差异图DI进行二值分割，得到新的分割图SI；

(7)计算区域P的分割阈值T与对数比差异图DI的预分割阈值T₁的差值，若差值小于10^-4则认为阈值收敛，执行步骤(8)，否则，将区域P的分割阈值T的值赋给对数比差异图预分割阈值T₁，将新的分割图SI的值赋给初步分割图SI₁，返回步骤(5)；

(8)对分割图SI进行3*3中值滤波，消除点状噪声对变化检测结果的影响，生成变化检测结果图。

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

1)本发明由于在图像处理的最后通过中值滤波方法去抑制图像中的杂噪点，无需对SAR图像进行斑噪滤波预处理，避免了因滤波程度设置不当而可能引起后续差异图生成、分析过程中效果降低的问题，同时提高了执行效率。

2)本发明由于应用最大类间方差法迭代获取最佳阈值，解决了当变化区域只占整体小比例，即差异图的直方图表征为单峰时，最大类间方差法分割结果产生过多虚警点的问题，提高变化检测精度。

3)本发明由于是在现有阈值方法上的进一步改进，保留了阈值分割方法的优势，计算公式简单，相对于其他差异图分析方法，例如：聚类法，水平集法具有很高的运行效率，能实现快速的变化检测。

附图说明

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

图2为用本发明对两时相Bern地区的SAR图像变化检测结果图；

图3为用本发明对两时相Ottawa地区的SAR图像变化检测结果图。

具体实施方式

下面结合附图对本发明做进一步的详细说明：

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

步骤1，获取适合变化检测使用的两幅SAR图像。

通过合成孔径雷达获取在不同时相拍摄的同一地区的两幅SAR图像，对这两幅SAR图像依次进行辐射校正、几何校正、图像配准的预处理，获得预处理后的适合变化检测使用的两幅SAR图像I₁和I₂；

所述预处理依次包括：辐射校正、几何校正、图像配准以及斑噪滤波；

由于斑噪滤波操作在有效的消除乘性噪声的影响的同时也会造成图像细节信息的丢失，进而影响定性性能和视觉效果，因而本发明方法将不涉及降斑处理，避免了因滤波程度设置不当引起后续差异图生成、分析过程中效果降低的可能性。

步骤2，计算生成对数比差异图。

对第一幅SAR图像I₁在坐标(i,j)处的灰度值I₁(i,j)和第二幅SAR图像I₂在坐标(i,j)处的灰度值I₂(i,j)进行比值运算，对比值运算的结果再进行一步对数运算，得到坐标(i,j)处的对数比值DI(i,j)：

DI(i,j)＝log(I₁(i,j)/I₂(i,j))，

遍历两幅SAR图像所有的像素点，得到对数比差异图：DI＝{DI(i,j)}。

步骤3，使用最大类间方差法对对数比差异图DI进行初分割，得到初步分割阈值T1。

(3a)设门限t，并根据该门限将包含S个灰度级的对数比差异图DI划分为C₁和C₂两类，C₁表示对数比差异图DI中灰度值小于等于t的像素点集合，C₁＝{1,2,K,t}，C₂表示对数比差异图DI中灰度值大于t的像素点集合，C₂＝{t+1,t+2,K,S-1}，0≤t≤S-1；

(3b)按照下式，计算C₁和C₂的类间方差σ²(t)：

σ²(t)＝p₁×(μ₁-μ_T)²+p₂×(μ₂-μ_T)²，

其中，p₁表示C₁出现的概率，n_λ表示灰度值为λ的像素点总数，W表示对数比差异图DI的像素点总数，p₂表示C₂出现的概率，μ₁表示C₁的灰度均值，μ₂表示C₂的灰度均值，μ(t)表示对数比差异图DI的灰度均值，

(3c)按照下式，计算类间方差σ²(t)取得最大值时的阈值T：

T＝Argmax{σ²(t)}，

其中，Argmax表示取最大值对应的变量操作；

(3d)将阈值T作为对数比差异图DI的初步分割阈值。

步骤4，用阈值T₁对对数比差异图DI进行二值分割，得到初步分割图SI₁。

(4a)判断对数比差异图DI中任意像素点的像素值是否大于阈值T₁，若是，则执行(4b)；否则，执行(4c)；

(4b)将对数比差异图DI中像素点的像素值大于阈值T₁的像素点的像素值置为0；

(4c)将对数比差异图DI中像素点的像素值小于等于阈值T₁的像素点的像素值置为1；

(4d)用对数比差异图DI中所有像素点的像素值置0置1后的像素点组成初步分割图SI₁。

步骤5，对初步分割结果图SI₁进行滑窗，获取所需的区域位置。

设置大小为w*w的窗口对初步分割图SI₁进行滑窗，统计窗口内零元素与非零元素个数直至小窗口内零元素与非零元素个数差小于等于3，得到满足条件的窗口中心元素记为(p,q)，其中w*w窗口大小取值为待处理图像大小的0.5％-1.5％，本实例的w取值为21。

步骤6，对对数比差异图DI上目标区域应用最大类间方差法进行分割。

取对数比差异图DI中心像素为(p,q)，窗口大小为w*w的区域P，使用最大类间法获取区域P的分割阈值T，用阈值T对对数比差异图DI进行二值分割得到新的分割图SI。

步骤7，判断是否达到收敛条件。

计算区域P的分割阈值T与对数比差异图DI的预分割阈值T₁的差值，并将该差值与设定的参数10^-4比较进行比较：

若差值小于10^-4则认为阈值收敛，执行步骤(8)；

否则，将区域P的分割阈值T的值赋给对数比差异图预分割阈值T₁，将新的分割图SI的值赋给初步分割图SI₁，返回步骤(5)。

步骤8，对分割图SI进行中值滤波，得到变化检测结果图。

对分割图SI进行3*3中值滤波，采用3*3窗口对二值图像进行滑窗，用窗口内元素中值替代窗口中心元素值，消除点状噪声对变化检测结果的影响，得到最终的SAR图像变化检测图像。

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

1.实验条件与结果

本发明的实验是在计算机配置为core i7 3.40GHZ，内存16G，WINDOWS 7系统和计算机软件配置为Matlab R2016a环境下进行的。

2.实验内容：

实验1，分别应用最大类间方差法Otsu、基于局部邻域模糊C均值聚类法FLICM和本发明方法对两时相Bern地区SAR图像进行变化检测，实验数据及结果如图2所示，其中：

图2(a)、图2(b)是ERS-2获取的1999年4月与1999年5月的两幅SAR图像，影像大小为301×301；

图2(c)是变化检测结果参考图，白色区域代表真实变化，该图中真实的发生变化像元数目为1155，占全部区域的1.27％，满足变化区域占整体区域小部分的情况；

图2(d)是最大类间法变化检测结果；

图2(e)是用现有FLICM方法的变化检测结果；

图2(f)是本发明方法的变化检测结果。

由图2可以看到，最大类间方差法产生了过多虚警，将大量非变化区域划分到变化区域，导致变化检测结果基本失效。FLICM方法检测结果略优于最大类间方差法，但还是有相当一部分非变化区域元素被划分到变化类中。本发明提出的方法通过迭代使用最大类间方差法遍寻全图查找到变化元素与非变化元素相对均匀的区域，然后根据图像自相似性原理将该区域的分割阈值应用到全局，解决了这一问题，变化检测结果图与给定参考图基本一致。

实验2，分别应用最大类间方差法Otsu、基于局部邻域模糊C均值聚类法FLICM和本发明方法对两时相Ottawa地区SAR图像进行变化检测，实验数据及结果如图3所示，其中：

图3(a)、图3(b)是RADARSAT-SAR卫星获取的1997年5月与1997年8月的两幅SAR图像，影像大小为290×350；

图3(c)是变化检测结果参考图，该图中总共有16049个真实变化像素点，占到总体的15.81％，大于10％，不属于变化区域占小比例的情况，将用该数据集验证本发明方法对于不是小面积变化区域情况的适用性；

图3(d)是最大类间方差法变化检测结果；

图3(e)是FLICM方法的变化检测结果；

图3(f)是本发明方法的变化检测结果。

由图3可以看到，对于不是变化区域占小比例的情况，应用最大类间方差法能获得不错的变化检测结果，变化区域基本能被检测到，但是结果图中仍存在大量的杂噪斑点。FLICM方法检测结果优于最大类间方差法，与参考结果基本一致，但依然存在不少虚警点。本发明方法的检测结果与FLICM方法检测结果大致相似，但是值得注意的是在图左下角区域本发明方法检测结果并未出现FLICM方法检测结果类似的条带虚警点。

3.实验结果评价：

为验证变化检测结果质量，本发明选取漏检数、虚警数、总错误检测数、Kappa系数和执行时间作为性能指标参数，评价本发明方法的检测精度和效率。

各评价参数表征如下：

漏检数Fc：变化类像素被判定为非变化类像素个数；

虚警数Fu：非变化类像素被判定为变化类像素个数；

总错误检测数OE：OE＝漏检数Fc+虚警数Fu；

Kappa系数计算公式为：

式中，Pr₀是观察一致率，Pr_c是期望一致率，N是图像总像素数目，其余符号表征意义如表1所述。

表1变化检测结果标识

Kappa系数是对正确检测像素及由漏检和虚警造成的错误检测像素的综合考虑，是关键性的评估指标。Fc、Fu、OE参数越小，Kappa系数越接近1表明检测准确率越高。

在同等实验条件下，执行时间是衡量算法执行效率的重要指标，执行时间越短，算法执行效率越高。

应用上述评价指标对两组实验的结果进行评价，如表2和表3所示，其中：

两种已有方法和本发明方法检测精度的比较结果表2所示：

表2本发明方法与两种已有方法检测精度比较

由表2可以看到，最大类间方差法在实验1变化区域仅占整体区域1.27％情况下，总错误检测数为7189，Kappa系数为0.2219，大量的虚警导致检测结果基本失效。在实验2非小面积变化区域情况下则得到了较好的结果，Kappa系数达到0.8102；FLICM方法也无法对小面积的变化区域进行良好的检测，实验1中Kappa系数仅为0.562，实验2中FLICM方法的检测结果则略优于本发明方法；本发明方法在实验1、2中均获得了良好的检测结果。本发明方法的有效性及普适性得以证明。

两种已有方法和本发明方法执行效率的比较结，如果表3所示：

表3本发明方法与两种已有方法执行效率比较

由表3中可以看出，本发明方法相对于Otsu方法牺牲掉了部分执行效率，但是本发明方法在执行效率上仍优于现有的FLICM方法。

综上所述，本发明提出的基于迭代最大类间方差的SAR图像变化检测方法，可以实现快速有效的变化检测，尤其是对于变化区域仅占整体区域小比例时，本发明方法检测效果优于其他两种已有方法。

Claims (4)

1.一种基于迭代最大类间方差的SAR图像变化检测方法，包括如下：

(1)通过合成孔径雷达获取在不同时相拍摄的同一地区的两幅SAR图像，对这两幅SAR图像依次进行辐射校正、几何校正、图像配准的预处理，获得预处理后的适合变化检测使用的两幅SAR图像I₁和I₂；

(2)计算第一幅SAR图像I₁在坐标(i,j)处的灰度值I₁(i,j)和第二幅SAR图像I₂在坐标(i,j)处的灰度值I₂(i,j)的比值，对比值运算的结果再进行一步对数运算，得到坐标(i,j)处的对数比值DI(i,j)：

DI(i,j)＝log(I₁(i,j)/I₂(i,j))，

遍历两幅SAR图像所有的像素点，得到对数比差异图：DI＝{DI(i,j)}；

(3)使用最大类间方差法对对数比差异图DI进行初分割，得到初步分割阈值T₁；

(4)用阈值T₁对对数比差异图DI进行二值分割得到初步分割图SI₁；

(5)设置大小为w*w的窗口对初步分割图SI₁进行滑窗，统计窗口内零元素与非零元素的个数，直至小窗口内零元素与非零元素个数差小于等于3时，停止滑窗，得到满足条件的窗口中心元素(p,q)，其中w*w窗口大小取值为待处理图像大小的0.5％-1.5％；

(6)取对数比差异图DI中中心像素为(p,q)，窗口大小为w*w的区域P，使用最大类间方差法获取区域P的分割阈值T，用阈值T对对数比差异图DI进行二值分割，得到新的分割图SI；

(7)计算区域P的分割阈值T与对数比差异图DI的预分割阈值T₁的差值，若差值小于10^-4则认为阈值收敛，执行步骤(8)，否则，将区域P的分割阈值T的值赋给对数比差异图预分割阈值T₁，将新的分割图SI的值赋给初步分割图SI₁，返回步骤(5)；

(8)对分割图SI进行3*3中值滤波，消除点状噪声对变化检测结果的影响，生成变化检测结果图。

2.如权利要求1所述方法，其中步骤(3)中使用最大类间方差法对对数比差异图DI进行分割，得到初步分割阈值T₁，按如下步骤进行：

3a)设门限t，并根据该门限将包含S个灰度级的对数比差异图DI划分为C₁和C₂两类，其中C₁表示对数比差异图DI中灰度值小于等于t的像素点集合，C₁＝{1,2,...,t}；C₂表示对数比差异图DI中灰度值大于t的像素点集合，C₂＝{t+1,t+2,...,S-1}，0≤t≤S-1；

3b)按照下式，计算C₁和C₂的类间方差σ²(t)：

σ²(t)＝p₁×(μ₁-μ_T)²+p₂×(μ₂-μ_T)²，

其中，p₁表示C₁出现的概率，n_λ表示灰度值为λ的像素点总数，W表示对数比差异图DI的像素点总数，p₂表示C₂出现的概率，μ₁表示C₁的灰度均值，μ₂表示C₂的灰度均值，μ_T表示对数比差异图DI的灰度均值，

3c)按照下式，计算σ²(t)取得最大值时的阈值T：

T＝Arg max{σ²(t)}，

其中，Arg max表示取最大值对应的变量操作；

3d)将阈值T作为对数比差异图DI的预分割阈值。

3.如权利要求1所述方法，其中步骤(4)中用阈值T₁对对数比差异图DI进行二值分割得到初步分割图SI₁，按如下步骤进行：

4a)判断对数比差异图DI中任意像素点的像素值是否大于阈值T₁，若是，则执行4b)；否则，执行4c)；

4b)将对数比差异图DI中像素点的像素值大于阈值T₁的像素点的像素值置为0；

4c)将对数比差异图DI中像素点的像素值小于等于阈值T₁的像素点的像素值置为1；

4d)将对数比差异图DI中所有像素点的像素值置0置1后的像素点组成初步分割图SI₁。

4.如权利要求1所述方法，其中步骤(8)中所述的对分割图SI进行3*3中值滤波，采用3*3窗口对分割图SI进行滑窗，用窗口内元素中值替代窗口中心元素值，消除点状噪声对变化检测结果的影响，得到最终的SAR图像变化检测图像。