CN108446637B - 基于立体图模型的sar图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于立体图模型的SAR图像变化检测方法，其实现过程是：首先对单时刻SAR图像进行高斯滤波处理，而后提取该滤波图像中的局部最大值点作为变化检测研究的关键点；根据关键点的位置属性，在两时SAR图像上构建立体图模型，进而获取两时关键点间的多跨度邻域及空间信息；依据此类信息，结合对数比和均值比算子对两时图像上关键点的变化情况进行度量；最后利用K‑means聚类算法对变化度量结果进行聚类分析，以分离变化区域和非变化区域。由于多跨度邻域及空间信息的引入，使得本发明方法在有效抑制相干斑噪声的同时能进一步提高变化检测方法的准确度。

Description

基于立体图模型的SAR图像变化检测方法

技术领域

本发明隶属于图像处理领域，其中涉及合成孔径雷达(Synthetic ApertureRadar,SAR)图像的变化检测，可以广泛应用到自然灾害监测、土地覆盖与利用、资源和环境监测、城市发展与规划以及军事侦察和打击效果评估等诸多领域。

背景技术

SAR图像变化检测是一种利用不同时刻、相同地理区域的两时SAR图像来获取地表变化信息的重要技术。而由于SAR微波成像的特殊机制，其可在任何大气和阳光条件下进行主动成像，具有全天候、全天时工作的特点；这使得SAR图像在自然灾害、国民经济和军事等应用领域的需求越来越多，因此开发具有高精度、高效率的SAR图像变化检测方法成为目前国内外学者研究的热点。

目前使用最为广泛SAR变化检测方法是学者Rignot提出的基于均值比(meanratio detector,MRD)的变化检测方法，该方法通过考虑两时图上对应像素点的局部均值来反映该点的变化情况，从而降低相干斑噪声对变化检测结果的影响。然而，如果该区域的局部纹理发生变化的同时均值保持不变，则MRD可能无法检测到该区域的变化信息。学者Bazi等人提出利用对数比(log ratio detector,LRD)方法在变化检测过程中减少乘性相干斑噪声的影响，以提高变化检测方法的抗噪性。然而，对数运算的特征在于增强低强度区域的像素值的同时削弱高强度区域的像素值，这样则压缩了差异图像的对比度，进而降低了图像的可分性。最近学者Pham提出使用关键点来进行变化检测研究，这些关键点被作为图的顶点，并通过图模型中的边进行连接，而后利用依据邻接顶点的变化信息来求解中心顶点的变化程度。与传统的基于像素的方法相比，该Graph-based方法有效降低了SAR图像变化检测的冗余度的同时提高了检测方法抗噪性。然而，该方法在变化区域边缘处的表征无法令人满意，往往会出现检测过界的现象。

发明内容

本发明的目的是针对现有变化检测方法所存在的问题，提出了一种基于立体图模型的SAR图像变化检测方法，以提高对图像空间信息的利用，从而提升变化检测方法的抗噪能力及检测精度。

为了达到上述目的，本发明所采用的实施步骤如下：

(1)输入两时SAR图像：

输入经过配准及几何校正的同一地区不同时刻的两幅SAR图像X₁和X₂，其获取时间分别为t₁和t₂。

(2)高斯滤波：

对SAR图像X₁进行高斯滤波处理，得到滤波后的SAR图像X'₁。

(3)关键点提取：

利用局部窗口在滤波后的SAR图像X'₁中漫游，提取窗口内强度最大的像素点作为研究的关键点，该提取过程可由如下公式表示：

式中p表示坐标为(i,j)的像素点，X'₁(p)为SAR图像X'₁中坐标(i,j)处的强度值；

是大小为ω×ω的局部邻域，q为该局部邻域

中的像素点，X'₁(q)为邻域像素点q的强度值。S_ω(X'₁)为要提取的局部最大值关键点集。此外，值得注意的是滤波图像X'₁仅用于特征点的提取，并不参与后面的立体图模型构建和变化度量研究。

(4)立体图模型构建：

4a)将所提取的关键点集S_ω(X'₁)作为图的顶点V，记关键点集内总点数量为N。计算图中各顶点间的欧式距离，而后依据K近邻规则，相互连接K个距离最近的邻接顶点(p_n,p_k)，从而在SAR图像X₁上建立平面图G₁。其中参数n＝(1,2,...,N)，k＝(1,2,...,K)。

4b)依据平面图G₁顶点分布及拓扑结构，在SAR图像X₂上建立具有与平面图G₁相同结构的平面图G₂；这两个平面图除顶点强度大小不同外，其顶点分布与拓扑结构完全相同。

4c)记

为SAR图像X₁上图G₁的顶点，

为SAR图像X₂上图G₂的顶点；在连接两平面图中坐标相同的对应顶点

的基础上，依据K近邻规则，连接不同平面图内的K个距离最近的邻接顶点

和

以建立完整的立体图模型。

定义立体图

其表达公式如下：

式中

为顶点

在图G₁上的K近邻域，

为该邻域内的顶点；

为顶点

在图G₂上的K近邻域，

为该邻域内的顶点；E_t为立体图时间跨度上的边，E_s为立体图空间跨度上的边，而E_st为立体图时空跨度上的边，上述E_t、E_s、E_st这三种边能够较好的反映两时SAR图像内各顶点在时间跨度、空间跨度及时空跨度上的邻域及空间信息；w(p_n,p_k)为空间距离权重，其主要反映的是顶点间的位置关系，D(p_n,p_k)为顶点p_n和p_k间的欧式距离，D_max为K个邻域点p_k中距离顶点p_n最远的距离，其主要用于对距离D(p_n,p_k)进行标准化，以确保空间距离权重w(p_n,p_k)的值在合理的区间内。

(5)变化度量：

结合立体图模型中各边关系，可利用如下公式对两时SAR图像上的关键点变化情况进行度量：

CM(p_n)＝CM_ts-CM_s+CM_t

式中CM_t为顶点间的时间相关性，CM_s为顶点间的空间相关性，CM_ts为顶点间的时空相关性。

其中：

式中

为图像X₁中顶点

附近ω_p×ω_p大小的块均值，

为顶点

附近的块均值；

为

处顶点强度的对数值，

与之类似。

(6)K-means聚类：

利用K-means聚类方法，对(5)中变化度量结果进行聚类，从而将各关键点分离为非变化类和变化类，获得变化检测结果图。

(7)输出变化检测结果。

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

1)本发明由于利用局部最大值点作为关键点进行变化检测研究，从而极大的降低了变化检测过程的冗余度。

2)本发明由于利用立体图模型获取的多跨度邻域信息，而后结合均值比和对数比算子来引导变化度量的进行，使得变化度量结果更加平滑和准确。

实验结果表明，本发明与现有的变化检测方法(MRD、LRD和Graph-based)相比，具有更高的检测精确度和更强的抗相干斑噪声能力。

附图说明

图1是本发明方法的实现流程图。

图2是本发明方法中立体图模型的结构图。

图3是本发明应用于Ottawa地区两时SAR图像的变化检测结果图；

图4是对发明应用于Bern地区两时SAR图像的变化检测结果图；

具体实施方式

参照图1所示，基于立体图模型的SAR图像变化检测方法的具体实施包括以下步骤：

步骤1.输入两时SAR图像：输入经过配准及几何校正的同一地区不同时刻的两幅SAR图像X₁和X₂。

步骤2.高斯滤波：设置高斯滤波尺度参数σ＝0.5，而后对SAR图像X₁进行高斯滤波处理，得到滤波后的SAR图像X'₁。

步骤3.关键点提取：利用3×3大小的局部窗口在滤波后的SAR图像X'₁中漫游，提取窗口内强度最大的像素点作为研究的关键点。

步骤4.立体图模型构建：将所提取的关键点集S_ω(X'₁)作为图的顶点V，并求取各顶点间的欧式距离；而后依据K近邻规则，相互连接K个距离最近的邻接顶点(p_n,p_k)，从而在SAR图像X₁上建立平面图G₁。依据平面图G₁顶点分布及拓扑结构，在SAR图像X₂上建立具有与平面图G₁相同结构的平面图G₂。在连接两平面图中坐标相同的对应顶点

的基础上，依据K近邻规则，连接不同平面图内的K个距离最近的邻接顶点

和

以建立完整的立体图模型，如参照图2所示。

步骤5.变化度量：依据立体图模型中三种不同类型的边，计算顶点间的时间相关性、空间相关性和时空相关性，结合三种相关性对各顶点的变化情况进行度量。

步骤6.K-means聚类：利用K-means聚类方法对步骤5中变化度量结果进行聚类，从而将各关键点分离为非变化类和变化类，获得变化检测结果图。

步骤7.输出变化检测结果。

通过以下两组真实SAR图像实验进一步说明本发明的有效性。

1.实验设置：

实验环境：MATLAB R2014，Intel(R)Core(TM)i7-4790CPU3.60GHz，Window7旗舰版。

实验数据：图3为Radarsat-1SAR系统获取的Ottawa地区的SAR图像，分辨率为10m，图像大小为290×350像素；图4为ERS-2SAR系统获取的Bern地区的SAR图像，分辨率为30m，图像大小为301×301像素。

实验参数：高斯滤波尺度参数σ＝0.5，关键点的搜索窗口大小ω＝3，均值块的大小ω_p＝3，邻域点个数K＝10。

实验对比方法：分别采用了MRD方法、LRD方法和Graph-based方法。

实验结果评价指标：分别选择漏检数、误检数、正确率及Kappa系数对实验结果进行客观评价。

2.实验结果与分析：

实验1分别应用MRD方法、LRD方法、Graph-based方法和本发明方法对Ottawa地区的两时SAR图像进行变化检测。实验结果如图3所示，其中图3(a)为Ottawa地区前一时刻洪水灾害发生时的SAR图像，图3(b)Ottawa地区后一时刻洪水退去时的SAR图像，图3(c)是人工判读的变化检测参考图，图3(d)是提取的关键点分布图，图3(e)～(h)分别对应着MRD方法、LRD方法、Graph-based方法和本发明方法的变化检测结果。

实验2分别四种变化检测方法对Ottawa地区的两时SAR图像进行研究。实验结果如图4所示，其中图4(a)是为Bern地区未发生洪水灾害时的SAR图像，图4(b)是为Bern地区洪水灾害发生时的SAR图像，图4(c)是人工判读的变化检测参考图，图4(d)是提取的关键点分布图，图4(e)～(h)分别对应着MRD方法、LRD方法、Graph-based方法和本发明方法的变化检测结果。

从图3、图4和表1可以看到，本发明在检测精确度及抗相干斑噪声的能力上均优于其他变化检测方法。这是由于本发明利用的立体图模型能更全面的获取多跨度领域和空间信息，且将对数比算子和均值比算子进行了有机结合，进而更好的反映非变化区域像素和变化区域像素间的区别，使得变化检测的准确度明显高于其他方法。

表1 Ottawa及Bern地区SAR图像变化检测结果

Claims (1)

1.一种基于立体图模型的SAR图像变化检测方法，包括以下步骤：

(1)输入两时SAR图像：

输入经过配准及几何校正的同一地区不同时刻的两幅SAR图像X₁和X₂，其获取时间分别为t₁和t₂；

(2)高斯滤波：

对SAR图像X₁进行高斯滤波处理，得到滤波后的SAR图像X′₁；

(3)关键点提取：

利用局部窗口在滤波后的SAR图像X′₁中漫游，提取窗口内强度最大的像素点作为研究的关键点，该提取过程可由如下公式表示：

式中p表示坐标为(i,j)的像素点，X′₁(p)为SAR图像X′₁中坐标(i,j)处的强度值；

是大小为ω×ω的局部邻域，q为该局部邻域

中的像素点，X′₁(q)为邻域像素点q的强度值；S_ω(X′₁)为要提取的局部最大值关键点集；此外，值得注意的是滤波图像X′₁仅用于特征点的提取，并不参与后面的立体图模型构建和变化度量研究；

(4)立体图模型构建：

4a)将所提取的关键点集S_ω(X′₁)作为图的顶点V，记关键点集内总点数量为N；计算图中各顶点间的欧式距离，而后依据K近邻规则，相互连接K个距离最近的邻接顶点(p_n,p_k)，从而在SAR图像X₁上建立平面图G₁；其中参数n＝(1,2,...,N)，k＝(1,2,...,K)；

4b)依据平面图G₁顶点分布及拓扑结构，在SAR图像X₂上建立具有与平面图G₁相同结构的平面图G₂；这两个平面图除顶点强度大小不同外，其顶点分布与拓扑结构完全相同；

4c)记

为SAR图像X₁上图G₁的顶点，

为SAR图像X₂上图G₂的顶点；在连接两平面图中坐标相同的对应顶点

的基础上，依据K近邻规则，连接不同平面图内的K个距离最近的邻接顶点

和

以建立完整的立体图模型；

定义立体图

其表达公式如下：

式中

为顶点

在图G₁上的K近邻域，

为该邻域内的顶点；

为顶点

在图G₂上的K近邻域，

为该邻域内的顶点；E_t为立体图时间跨度上的边，E_s为立体图空间跨度上的边，而E_st为立体图时空跨度上的边，上述E_t、E_s、E_st这三种边能够较好的反映两时SAR图像内各顶点在时间跨度、空间跨度及时空跨度上的邻域及空间信息；w(p_n,p_k)为空间距离权重，其主要反映的是顶点间的位置关系，D(p_n,p_k)为顶点p_n和p_k间的欧式距离，D_max为K个邻域点p_k中距离顶点p_n最远的距离，其主要用于对距离D(p_n,p_k)进行标准化，以确保空间距离权重w(p_n,p_k)的值在合理的区间内；

(5)变化度量：

结合立体图模型中各边关系，可利用如下公式对两时SAR图像上的关键点变化情况进行度量：

CM(p_n)＝CM_ts-CM_s+CM_t

式中CM_t为顶点间的时间相关性，CM_s为顶点间的空间相关性，CM_ts为顶点间的时空相关性；

其中：

式中

为图像X₁中顶点

附近ω_p×ω_p大小的块均值，

为顶点

附近的块均值；

为

处顶点强度的对数值，

与之类似；

(6)K-means聚类：

利用K-means聚类方法，对(5)中变化度量结果进行聚类，从而将各关键点分离为非变化类和变化类，获得变化检测结果图；

(7)输出变化检测结果。

Images