CN108564603B - 基于改进的高斯混合模型的sar图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于改进的高斯混合模型的SAR图像变化检测方法，具体包括如下步骤：步骤1，对SAR图像进行滤波预处理；步骤2，利用log‑ratio运算符获取差异图像X；步骤3，划分差异图像X的特定区域；步骤4，在高斯模型中加入平移参数m获得改进的高斯模型，利用改进的高斯模型建立统计模型对差异图像X的特定区域的直方图进行拟合，并基于GKI算法最小化代价函数获得最终的初始阈值T*；步骤5，建立能量方程；步骤6，寻找步骤5建立的能量方程的一个局部极小值，获得对差异图像X的分割结果；步骤7，对步骤6得到的分割结果的变化类区域标记为“1”，非变化类标记为“0”，输出最终的变化检测结果。

Description

基于改进的高斯混合模型的SAR图像变化检测方法

技术领域

本发明属于遥感图像处理技术领域，涉及基于改进的高斯混合模型的SAR图像变化检测方法。

背景技术

SAR图像变化检测技术的目的是通过分析在两个不同时期获取的待检测区域的图像来识别土地覆盖的变化。由于这一技术是从大量应用背景派生而来的，因此广泛引起了图像处理领域学者们的研究兴趣。对于合成孔径雷达图像而言，由于受到斑点噪声的影响，变化检测任务表现出更多的困难。但是，SAR图像具有全天候、全天时的特性，这使得SAR图像的变化检测是监测环境变化，如自然资源、土地、城市环境等，尤其是在自然灾害评估分析方面已得到广泛地应用。

SAR图像变化检测旨在确定不同时间获取的同一场景土地覆盖的变化情况。为了区分变化类和非变化类区域，在对两类像素统计建模的基础上，阈值算法是一类行之有效的非监督变化检测算法。然而这类算法存在的最大障碍就是选择的统计模型直接决定了分类精度，这成为限制此类算法的最大因素。因此，如何提高对差异图像的统计建模质量是SAR图像变化检测技术中的一项关键内容。

差异图像的概率统计分布属性为变化区域的分析和检测提供了关键依据。然而SAR的受到斑点噪声、几何校正和辐射校正等多种因素的影响，导致在差异图像的直方图上，变化类和非变化类像素相互重叠，使得基于统计模型的变化检测算法困难重重。在大量的文献中，研究者们尝试了各种统计模型对差异图像进行建模，期待获得良好的阈值分割结果。张小华等人对高斯模型和广义高斯模型进行了对比，结果表明后者对于变化类的分布情况不能够准确描述。此外，Ban等人对比了对数正态分布、广义高斯模型、Nakagami分布和Weibull分布对变化检测结果的影响，认为对数正态分布和Nakagami分布模型对城市地区变化类建模能够获得更好的结果。可以看出，统计模型的方法在拟合差异图像变化类和非变化类分布特性，处理SAR图像变化检测问题时，虽然有效但不具有广泛适用性。由于基于概率统计模型的阈值算法是一种非监督的变化检测技术。Bruzzone等人提出使用EM算法自适应地选择阈值，使得后验概率误差最小化。然而，这种方法获得的变化检测结果的精度依然取决于给定的统计模型与差异图像直方图的拟合度，也就是说，统计模型对变化类或非变化类区域描述越准确获得的变化检测结果正确率越高。除此此外，通过最小化代价函数，GKI阈值选择算法能够确定最佳阈值，然而，这种方法也存在着同样的问题。

发明内容

本发明的目的是提供基于改进的高斯混合模型的SAR图像变化检测方法，能够通过加入平移参数m使得差异图像的直方图与给定的统计模型尽可能逼近，从而通过代价函数最小化快速确定最优初始阈值，为后续的图像分割算法提供有利的初始条件。

本发明所采用的技术方案是，基于改进的高斯混合模型的SAR图像变化检测方法，具体包括如下步骤：

步骤1，选取窗口大小为3×3的Lee滤波器分别对两幅原始的获取于不同时间的SAR图像X₁，X₂进行滤波预处理；

步骤2，利用log-ratio运算符获取差异图像X；

步骤3，利用高斯模型对差异图像X的变化类和非变化类像素建立统计模型，然后基于GKI算法获取该统计模型的初始阈值T₀，并根据初始阈值T₀划分差异图像X的特定区域；

步骤4，在步骤3划分的特定区域内，在高斯模型中加入平移参数m获得改进的高斯模型，利用改进的高斯模型建立统计模型对差异图像X的特定区域的直方图进行拟合，并基于GKI算法最小化代价函数获得最终的初始阈值T*；

步骤5，基于步骤4获得的最终的初始阈值T*建立能量方程；

步骤6，采用α-expansion算法通过迭代运算来寻找步骤5建立的能量方程的一个局部极小值，获得对差异图像X的分割结果；

步骤7，对步骤6得到的分割结果的变化类区域标记为“1”，非变化类标记为“0”，输出最终的变化检测结果。

本发明的特点还在于，

步骤2中，采用log-ratio运算符获取的差异图像X为：

其中，X₁，X₂分别为获取于不同时间的SAR图像，log为log运算符。

步骤3的具体过程如下：

步骤3.1，利用高斯对差异图像X的变化类和非变化类像素建立统计模型，建立的统计模型如下公式(2)、(3)所示：

其中，Ω＝{变化类，非变化类}，P_ω为属于ω类的像素的先验概率；μ_ω为属于ω类的像素的均值；

为属于ω类的像素的方差，T表示统计模型的阈值；p(x_p|μ_ω,σ_ω,T)为基于P_ω、μ_ω、σ_ω、T获得的像素灰度值为x_p的概率；

步骤3.2，基于GKI算法获取步骤3.1中统计模型的初始阈值T⁰，如下公式(4)、(5)所示：

其中，h(x_p)代表真实的差异图像X的特定区域的直方图；p(x_p|μ_ω,σ_ω)为基于μ_ω、σ_ω获得的像素灰度值为x_p的概率；

步骤3.3，基于步骤3.2获得的初始阈值T⁰确定差异图像X的特定区域的下限R_low和上限R_up；

首先，通过比较直方图和差异图像X中非变化类像素的最大值来定义特定区域的下限R_low，如下公式(6)所示：

通过如下公式(7)确定特定区域的上限R_up：

其中，ST为平滑参数；

步骤3.4，根据步骤3.3获得的特定区域的下限R_low和上限R_up划分差异图像X上的特定区域。

步骤4的具体过程如下：

步骤4.1，在高斯模型中加入平移参数m，获得改进的高斯模型，如下公式(8)、(9)、(10)所示：

步骤4.2，根据如下公式(11)获得平移量m，获得平移参数m的过程，即为步骤4.1获得的改进的高斯模型分布曲线逐渐逼近特定分布区域的直方图的过程，该过程即完成了对差异图像X的特定区域的直方图进行拟合：

其中，

步骤4.3，基于GKI算法最小化代价函数获得最终的初始阈值T*，如下公式(12)、(13)所示：

步骤5的具体过程如下：

利用步骤4获得的最终的初始阈值T*对差异图像X进行阈值分割获得初始标签L⁰，基于初始标签L⁰建立的能量方程如下：

其中，N表示像素的集合。

本发明的有益效果是，本发明提供的基于改进的高斯混合模型的SAR图像变化检测方法，能够提高能量方程的局部极小值搜索效率；此外本发明对于复杂分布的SAR图像，仍然能够有效的检测变化区域，可用于对环境变化的自动检测。实验证明，本发明对SAR图像可以准确的检测变化区域。

附图说明

图1是本发明基于改进的高斯混合模型的SAR图像变化检测方法中改进后的高斯模型获取的结果图；

图2(a-1)、图2(a-2)、图2(a-3)表示黄河河口-农田的原始SAR的图像；

图2(b-1)、图2(b-2)、图2(b-3)表示黄河河口-河流的原始SAR的图像；

图2(c-1)、图2(c-2)、图2(c-3)表示黄河河口-海岸线的原始SAR的图像；

图2(d-1)、图2(d-2)表示较大区域的原始SAR的图像；

图3(a)表示是采用本发明基于改进的高斯混合模型的SAR图像变化检测方法对伯尔尼区域进行实验得到的变化区域的检测结果图；

图3(b)表示是采用本发明基于改进的高斯混合模型的SAR图像变化检测方法对渥太华区域进行实验得到的变化区域的检测结果图；

图3(c)表示是采用本发明基于改进的高斯混合模型的SAR图像变化检测方法对石门水库区域进行实验得到的变化区域的检测结果图；

图3(d)表示是采用本发明基于改进的高斯混合模型的SAR图像变化检测方法对黄河河口-农田区域进行实验得到的变化区域的检测结果图。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细说明。

本发明基于改进的高斯混合模型的SAR图像变化检测方法，具体包括如下步骤：

步骤1，选取窗口大小为3×3的Lee滤波器分别对两幅原始的获取于不同时间的SAR图像X₁，X₂进行滤波预处理；

步骤2，利用log-ratio运算符获取差异图像X；

采用log-ratio运算符获取的差异图像X为：

其中，X₁，X₂分别为获取于不同时间的SAR图像，log为log运算符。

步骤3，利用高斯模型对差异图像X的变化类和非变化类像素建立统计模型，然后基于GKI算法获取该统计模型的初始阈值T₀，并根据初始阈值T₀划分差异图像X的特定区域；步骤3的具体过程为：

步骤3.1，利用高斯对差异图像X的变化类和非变化类像素建立统计模型，建立的统计模型如下公式(2)、(3)所示：

其中，Ω＝{c，u}，c为变化类，u代表非变化类。，P_ω为属于ω类的像素的先验概率；μ_ω为属于ω类的像素的均值；

为属于ω类的像素的方差，T表示统计模型的阈值；p(x_p|μ_ω,σ_ω,T)为基于P_ω、μ_ω、σ_ω、T获得的像素灰度值为x_p的概率；

步骤3.2，基于GKI算法获取步骤3.1中统计模型的初始阈值T⁰，如下公式(4)、(5)所示：

其中，GL为最大灰度值，等于255。h(x_p)代表真实的差异图像X的特定区域的直方图；p(x_p|μ_ω,σ_ω)为基于μ_ω、σ_ω获得的像素灰度值为x_p的概率；

步骤3.3，基于步骤3.2获得的初始阈值T⁰确定差异图像X的特定区域的下限R_low和上限R_up；

首先，通过比较直方图和差异图像X中非变化类像素的最大值来定义特定区域的下限R_low，如下公式(6)所示：

这样的好处在于可以确保特定区域位于直方图的右侧。

通过如下公式(7)确定特定区域的上限R_up：

其中，ST为平滑参数；t为灰度值偏移变量。

步骤3.4，根据步骤3.3获得的特定区域的下限R_low和上限R_up划分差异图像X上的特定区域。

步骤4，在步骤3划分的特定区域内，在高斯模型中加入平移参数m获得改进的高斯模型，利用改进的高斯模型建立统计模型对差异图像X的特定区域的直方图进行拟合，并基于GKI算法最小化代价函数获得最终的初始阈值T*；

步骤4的具体过程如下：

步骤4.1，在高斯模型中加入平移参数m，获得改进的高斯模型，如下公式(8)、(9)、(10)所示：

通过参数m的变化使目标区域面积减少，达到最佳的拟合效果；

步骤4.2，根据如下公式(11)获得平移量m，获得平移参数m的过程，即为步骤4.1获得的改进的高斯模型分布曲线逐渐逼近特定分布区域的直方图的过程，该过程即完成了对差异图像X的特定区域的直方图进行拟合：

其中，

h(·)代表基于改进后的高斯模型与特定区域拟合后的真实的直方图；h(·)的括号中的参量可根据不同的直方图填入不同的参数；p_k＝p_m(x_p)；p(x_p|μ_ω,σ_ω,m,T⁰)为基于μ_ω、σ_ω、m、T⁰获得的像素灰度值为x_p的概率；k为灰度值偏移变量。

步骤4.3，基于GKI算法最小化代价函数获得最终的初始阈值T*，如下公式(12)、(13)所示：

其中，Ω＝{c，u}，c为变化类，u代表非变化类。在公式(13)中，P_ω为属于ω类的像素的先验概率；μ_ω为属于ω类的像素的均值；

为属于ω类的像素的方差，T表示统计模型的阈值；p(x_p|μ_ω,σ_ω,T)为基于P_ω、μ_ω、σ_ω、T获得的像素灰度值为x_p的概率；GL为最大灰度值，等于255。R_c＝{x_p|L_p＝c}R_u＝{x_p|L_p＝u}分别表示属于变化类和非变化类的像素集合。

步骤5，基于步骤4获得的最终的初始阈值T*建立能量方程；

利用步骤4获得的最终的初始阈值T*对差异图像X进行阈值分割获得初始标签L⁰，基于初始标签L⁰建立的能量方程，本发明中差异图像X能量函数的数据项可以通过下式描述：

D_p(L_p)＝(θ_ω-x_p)²x_p∈R_ω,ω∈Ω (14)；

数据项使用分段常数模型，其中θ_ω代表分段常数模型的参数。一个不连续保持函数作为平滑项V_p,q，定义为：

V_p,q(L_p,L_q)＝|L_p-L_q| (15)；

这里，如果L_p＝L_q，则V_p,q(L_p,L_q)值为零，否则，该值是一个常数C。

由于SAR图像非线性和复杂分布的特性，分段常数模型不能直接适用于表达差异图像中变化类和非变化类像素的统计分布。因此，新算法使用核函数K(i,j)＝Ψ(i)^TΨ(j)来提高数据在特征空间中的可分性。根据Mercer定理，特征空间中点积形式可以被用来构造图像能量方程的数据项。本发明使用了在图像处理领域的应用最为广泛的径向基函数(RBF)—高斯核函数K(i,j)＝exp(-(i-j)²/(2σ²))。因此，建立的能量方程描述如下：

这里

(Ψ(θ_ω)-Ψ(x_p))²

＝(Ψ(θ_ω)-Ψ(x_p))^T(Ψ(θ_ω)-Ψ(x_p))

＝K(θ_ω,θ_ω)+K(x_p,x_p)-2K(θ_ω,x_p)

Ν代表对像素的集合。

步骤6，采用α-expansion算法通过迭代运算来寻找步骤5建立的能量方程的一个局部极小值，获得对差异图像X的分割结果；给定初始标签L⁰，通过图切算法找到标签L使能量方程中的E最小化获得变化检测的最终结果，其迭代过程包含了更新标签和更新参数两个步骤。

步骤7，对步骤6得到的分割结果的变化类区域标记为“1”，非变化类标记为“0”，输出最终的变化检测结果。

实施例

实验所用微机CPU为Intel Pentium4 3.0GHz内存1GB，编程平台是Matlab 7.0.1。实验中采用的图像数据为自黄河河口实验数据集的不同的地理位置选择的三个具有代表性的区域，分别是黄河河口-农田，黄河河口-河流和黄河河口-海岸线。该数据集来自Radarsat-2传感器在C-波段HH极化“带状地图”模式下获得的一组山东东营黄河河口数据。空间分辨率为8m×8m。黄河河口-农田数据集显示一块农田内陆，属于变化区域的部分相对较大且有规律，大小为306×291像素。黄河河口-河流数据集给出了一块内陆水域的部分区域，大小为291×444像素。黄河河口-海岸线数据集所示的变化区域是海岸线附近的海洋表面，大小为450×280像素。较大区域数据集选择了一块面积大小为1272×1149像素的区域。

2.实验内容

首先，对原始的SAR图像进行预处理，产生差异图像。其次，使用改进的高斯混合模型对差异图像中的变化类和非变化类像素建立统计模型，利用GKI算法获得初始标签。再次，建立能量函数，采用分段常数分割模型作为能量函数的数据项，采用一个不连续保持函数作为平滑项。最后，利用α-expansion算法通过迭代运算来寻找到一个局部极小值，获得最终的变化检测结果。

通过对变化检测结果参考图和实验得到的变化检测结果图对比，来评价本发明的效果。

3.实验结果

图1是利用函数J(T|m)值选取变化示例，图2是SAR图像数据集，图3是本发明实验得到的变化检测结果图。

图2(a-1)、图2(a-2)、图2(a-3)表示黄河河口-农田的原始SAR的图像；图2(b-1)、图2(b-2)、图2(b-3)表示黄河河口-河流的原始SAR的图像；图2(c-1)、图2(c-2)、图2(c-3)表示黄河河口-海岸线的原始SAR的图像；图2(d-1)、图2(d-2)表示较大区域的原始SAR的图像。图3(a)表示伯尔尼区域的检测结果图；图3(b)表示渥太华区域的检测结果图；图3(c)表示石门水库区域的检测结果图；图3(d)表示黄河河口-农田区域的检测结果图。

从图2可以看出，SAR图像背景中斑点噪声的影响，且变化区域形状上也存在很大差异。通过比较图2和图3可以看出，图2中的变化区域都可以被准确的检测出来。

Claims (5)

1.基于改进的高斯混合模型的SAR图像变化检测方法，其特征在于：具体包括如下步骤：

步骤1，选取窗口大小为3×3的Lee滤波器分别对两幅原始的获取于不同时间的SAR图像X₁，X₂进行滤波预处理；

步骤2，利用log-ratio运算符获取差异图像X；

步骤3，利用高斯模型对差异图像X的变化类和非变化类像素建立统计模型，然后基于GKI算法获取该统计模型的初始阈值T₀，并根据初始阈值T₀划分差异图像X的特定区域；

步骤4，在步骤3划分的特定区域内，在高斯模型中加入平移参数m获得改进的高斯模型，利用改进的高斯模型建立统计模型对差异图像X的特定区域的直方图进行拟合，并基于GKI算法最小化代价函数获得最终的初始阈值T*；

步骤5，基于步骤4获得的最终的初始阈值T*建立能量方程；

步骤6，采用α-expansion算法通过迭代运算来寻找步骤5建立的能量方程的一个局部极小值，获得对差异图像X的分割结果；

步骤7，对步骤6得到的分割结果的变化类区域标记为“1”，非变化类标记为“0”，输出最终的变化检测结果。

2.根据权利要求1所述的基于改进的高斯混合模型的SAR图像变化检测方法，其特征在于：

步骤2中，采用log-ratio运算符获取的差异图像X为：

其中，X₁，X₂分别为获取于不同时间的SAR图像，log为log运算符。

3.根据权利要求1所述的基于改进的高斯混合模型的SAR图像变化检测方法，其特征在于：所述步骤3的具体过程如下：

步骤3.1，利用高斯对差异图像X的变化类和非变化类像素建立统计模型，建立的统计模型如下公式(2)、(3)所示：

其中，Ω＝{变化类，非变化类}，P_ω为属于ω类的像素的先验概率；μ_ω为属于ω类的像素的均值；

为属于ω类的像素的方差，T表示统计模型的阈值；p(x_p|μ_ω,σ_ω,T)为基于μ_ω、σ_ω、T获得的像素灰度值为x_p的概率；

步骤3.2，基于GKI算法获取步骤3.1中统计模型的初始阈值T⁰，如下公式(4)、(5)所示：

其中，h(x_p)代表真实的差异图像X的特定区域的直方图；p(x_p|μ_ω,σ_ω)为基于μ_ω、σ_ω获得的像素灰度值为x_p的概率；GL为最大灰度值；

步骤3.3，基于步骤3.2获得的初始阈值T⁰确定差异图像X的特定区域的下限R_low和上限R_up；

首先，通过比较直方图和差异图像X中非变化类像素的最大值来定义特定区域的下限R_low，如下公式(6)所示：

通过如下公式(7)确定特定区域的上限R_up：

其中，ST为平滑参数；

步骤3.4，根据步骤3.3获得的特定区域的下限R_low和上限R_up划分差异图像X上的特定区域。

4.根据权利要求3所述的基于改进的高斯混合模型的SAR图像变化检测方法，其特征在于：步骤4的具体过程如下：

步骤4.1，在高斯模型中加入平移参数m，获得改进的高斯模型，如下公式(8)、(9)、(10)所示：

其中，p(x_p|μ_ω,σ_ω,m,T⁰)为基于μ_ω、σ_ω、m、T⁰获得的像素灰度值为x_p的概率，p_k＝p_m(x_p)，k为灰度值偏移变量；

步骤4.2，根据如下公式(11)获得平移量m，获得平移参数m的过程，即为步骤4.1获得的改进的高斯模型分布曲线逐渐逼近特定分布区域的直方图的过程，该过程即完成了对差异图像X的特定区域的直方图进行拟合：

其中，

步骤4.3，基于GKI算法最小化代价函数获得最终的初始阈值T*，如下公式(12)、(13)所示：

5.根据权利要求1所述的基于改进的高斯混合模型的SAR图像变化检测方法，其特征在于：步骤5的具体过程如下：

利用步骤4获得的最终的初始阈值T*对差异图像X进行阈值分割获得初始标签L⁰，基于初始标签L⁰建立的能量方程如下：

其中，N表示像素的集合，

如果L_p＝L_q，则V_p,q(L_p,L_q)值为零，否则，该值是一个常数C；

Ω＝{c，u}，c为变化类，u代表非变化类。

Images