CN103353989B - 基于先验和融合灰度与纹理特征的sar图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于先验和融合灰度与纹理特征的SAR图像变化检测方法，主要解决用高斯模型不能完全拟合差异图的分布及只利用SAR图像的像素灰度信息所导致的变化检测准确率低的问题。其实现步骤是：(1)读入两幅已配准、已校正的两时相SAR图像；(2)对两幅图像采用小波融合策略构造差异图；(3)对差异图求类别的先验概率；(4)融合差异图的灰度与纹理信息求得观测量似然概率；(5)由类别先验概率与观测量似然概率计算后验概率；(6)由最大后验概率准则将差异图分为变化类与非变化类；(7)重复步骤(3)到步骤(6)直到满足终止条件，输出最终的变化检测结果。本发明方法具有对SAR图像变化检测精度高的优点，可用于提取和获得SAR图像的变化细节信息。

Description

基于先验和融合灰度与纹理特征的SAR图像变化检测方法

技术领域

本发明属于图像处理领域，涉及一种图像变化检测的方法，尤其涉及一种对不同时刻同一地区已配准的SAR图像变化检测的方法，可用于对多时相SAR图像提取与获得地物的变化信息，提高SAR图像变化检测的精确度，从而对地物信息情况进行更精准的监测与评估。

背景技术

随着合成孔径雷达SAR技术的快速发展，其分辨率不断的提高，获取的合成孔径雷达图像具有不受外界天气条件和太阳光照射强度对地物成像的影响等特点，弥补了光学传感器与红外成像的不足，使得SAR图像的应用日益增加，其中SAR图像的变化检测得到了广泛的关注。SAR图像变化检测是对同一地区不同时间的两幅或多幅SAR图像进行比较，分析出图像之间的差异从而得到所需的地物变化信息。它主要应用于自然灾害、城市扩张情况的分析和军事应用等方面。

近年来，人们根据SAR图像成像的特点及其具有的斑点噪声，提出了许多新颖有效的变化检测方法，用于提高SAR图像变化检测的性能，这些方法大致可以分为图像阈值法和图像分类法两大类。从图像分类的角度，经典的图像分类技术有：MRF模型、贝叶斯技术、模糊集理论等，这些理论在变化检测领域通常被使用。其中基于马尔科夫随机场MRF模型的方法用于SAR图像变化检测的主要过程是对差异图进行初始分割得到一幅二值图像，用MRF对二值图像求取先验概率，再用高斯分布求取似然概率，最后用贝叶斯公式求后验概率得到最后的检测结果。该方法的不足之处在于：(1)一般求得的差异图的分布是一个混合分布，高斯模型并不能完全拟合差异图的分布，并且如何从这种混合分布形式中准确推断出每个高斯分量的均值，方差和形状参数，以及它们的权重也是一个复杂的统计推断问题；(2)在整个求解过程中只是用到了差异图的灰度信息，没有充分利用差异图的其他信息，如：纹理特征和区域特征等，由此SAR图像特有的斑点噪声会对检测结果造成很大的影响，提高了像素的错分率，降低了变化检测的精确度。

发明内容

本发明的目的在于针对上述现有方法的不足，提出一种基于先验和融合灰度与纹理特征的SAR图像变化检测方法，以降低SAR图像斑点噪声的影响，提高图像像素分类的精确性。

实现本发明目的的技术方案是：将变化检测问题看作是两个相互独立分量乘积的最大概率问题。首先，由吉布斯分布求取差异图初始分类的类别先验概率，将其看作是第一个分量。其次，由模糊隶属度求取融合差异图的灰度与纹理特征的模糊隶属概率，将其看作是第二个分量，最后，利用贝叶斯公式求后验概率与独立分布准则，有效地检查出不同时刻同一地点SAR图像变化的区域。其具体步骤如下：

实现本发明目的的技术方案包括如下步骤：

(1)读入已配准和已校正的两幅不同时间同一地点的SAR图像I₁和I₂；

(2)用均值比值法构造这两幅SAR图像I₁和I₂的差异图D₁，用对数比值法构造这两幅SAR图像I₁和I₂的差异图D₂，利用小波变换对所述差异图D₁和所述差异图D₂进行融合，得到融合后差异图D，并初始化迭代次数t=1；

(3)对融合后差异图D分别求类别先验概率P(x)和观测量似然概率P(y)；

(4)根据类别先验概率P(x)和观测量似然概率P(y)，利用贝叶斯公式和独立分布原理求取融合后差异图D的后验概率P(x|y)为：

P(x|y)＝P(y|x)P(x)

＝P(y)P(x)，

其中，P(y|x)是条件概率，x是初始分类结果图，y是融合灰度与纹理信息的观测场，并且初始分类结果图x和融合灰度与纹理信息的观测场y是相互独立的；

(5)根据差异图D的后验概率P(x|y)和MAP准则，得到变化检测结果图为：

$\hat{x}=\arg \max P\left(x\right|y)$

$=\arg \max \left\{P\left(x\right)P\left(y\right)\right\},$

(6)判断循环次数t是否达到了最高迭代次数g_max，如果满足t＞g_max则输出最终的变化检测结果图否则，令t＝t+1，返回步骤(3)，进行下一次迭代。

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

1.本发明方法利用融合差异图灰度与纹理信息的模糊隶属概率作为观测量的似然概率，充分利用了差异图的灰度与纹理信息，避免了用高斯模型拟合差异图及参数的估计，有效地简化了算法实现的过程与复杂度，同时提高了图像像素分类的精确性，为求取最终检测结果提高了检测精度；

2.本发明方法将类别的先验概率和融合灰度与特征的模糊隶属概率看作是两个相互独立的分量，由此将贝叶斯公式求取后验概率看作是对两个独立分量相乘积来求取后验概率，从另一角度分析实现了一种新的SAR图像变化检测方法，提高了变化检测的精确度。

附图说明

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

图2是用本发明和现有MRF方法对1999年4月和1999年5月瑞士Bern地区的图像变化检测实验结果对比图；

图3是用本发明和现有MRF方法对1997年5月和1997年8月加拿大Ottawa地区的图像变化检测实验结果对比图；

图4是用本发明和现有MRF方法对2008年6月和2009年6月黄河入海口部分截图的变化检测实验结果对比图。

具体实施方式

下面结合附图对本发明的具体实施步骤和效果做进一步的详细描述：

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

步骤1，读入两幅已配准、已校正的两幅不同时间同一地点的SAR图像I₁和I₂。

在本发明的实施例中，读入两幅分别在1999年4月和1999年5月通过ERS-2获得的瑞士Bern城市水灾发生前后的SAR图像I₁和I₂，两幅图像的尺寸均为301×301像素，灰度级为256，实际变化像素个数为1155。

步骤2，用均值比值法构造这两幅SAR图像I₁和I₂的差异图D₁，用对数比值法构造这两幅SAR图像I₁和I₂的差异图D₂，利用小波变换对所述差异图D₁和所述差异图D₂进行融合，得到融合后差异图D。

在本发明的实施例中，在利用小波变换对所述差异图D₁和所述差异图D₂进行融合时采用的是三层小波分解。

步骤3，对融合后差异图D求类别先验概率P(x)。

3a)对融合后差异图D用K-means方法进行初始分割，得到初始分类结果图x；

3b)根据吉布斯分布求初始分类结果图x的先验概率：P(x)＝{P(x_i),i＝1,...,L}，L是分类结果图x的总像素点个数，P(x_i)为分类结果图x中第i个像素点x_i的先验概率，表示为：

$P\left(x_i\right)=\frac{\exp (-u\left(x_i\right))}{\underset{x_i\∈x}{\mathrm{\Σ}}\exp (-u\left(x_i\right))},$

其中，u(x_i)是第i个像素点x_i的能量函数，且u(x_i)＝-βΣδ(x_i,x_j)-1， $\mathrm{\δ}(x_i,x_j)=\left\{\begin{array}{cc}1& x_i=x_j\\ 0& x_i\≠x_j\end{array}\right.,$ β是平滑参数，x_j是第i个像素点x_i的邻域集N_i中的第j个像素点，且x_j∈N_i。

在本发明的实施例中，平滑参数β=0.5，求得的先验概率P(x)是一个M×N×k(301×301×2)像素的矩阵，邻域集Ni采用的是邻域窗口大小为(3×3)像素。

步骤4，对融合后差异图D求观测量似然概率P(y)。

4a)根据融合后差异图D进行非下采样小波变换，求取融合后差异图D的纹理特征矩阵G；

在本发明的实施例中，对融合后差异图D提取纹理特征的方法有灰度共生矩阵，小波变换，非下采样小波变换等，本发明方法采用的是非下采样小波变换来提取差异图的纹理特征矩阵G，在非下采样小波变换的方法中采用L=3层小波分解，同时采用区域窗口大小为(15×15)像素，生成的纹理特征矩阵G为(301×301×d)像素，d为特征矩阵G的维数，d=(3×L+1)，L=3。

4b)根据纹理特征矩阵G求取融合后差异图D的特征模糊隶属度矩阵U₁，U₁＝{u_ab,a＝1,...,k,b＝1,...,S}，k为融合后差异图D的分类类别数，S是融合后差异图D的总像素点个数，u_ab是像素点b隶属于第a类的特征模糊隶属度，同时由融合后差异图D中像素点的灰度值求取灰度模糊隶属度矩阵U₂，U₂＝{u_rh,r＝1,...,k,h＝1,...,S}，k为融合后差异图D的分类类别数，u_rh是像素点h隶属于第r类的灰度模糊隶属度；

在本发明实施例中，融合后差异图D的分类类别数k=2，融合后差异图D的总像素点个数S＝301×301，求得的特征模糊隶属度矩阵U₁为(301×301×2)像素，且其每个像素的取值在0～1之间，求得的灰度模糊隶属度矩阵U₂为(301×301×2)像素，且其每个像素的取值在0～1之间；

4c)根据特征模糊隶属度矩阵U₁和灰度模糊隶属度矩阵U₂，得到观测量的似然概率P(y)为：

$P\left(y\right)=\left\{\begin{array}{ccc}\max (U_1,U_2),& \mathrm{if}& k=1\\ \min \left(U_1{,U}_2\right),& \mathrm{if}& k=2\end{array}\right.,$

其中，k是融合后差异图D的分类数，y是融合灰度与纹理信息的观测场。

在本发明的实施例中，由于变化检测过程是将差异图D分为两类，所以k＝1,2，观测量的似然概率P(y)为(301×301×2)像素，且其大小在0～1之间。

步骤5，根据类别先验概率P(x)和观测量似然概率P(y)，利用贝叶斯公式和独立分布原理求取融合后差异图D的后验概率P(x|y)为：

P(x|y)＝P(y|x)P(x)

＝P(y)P(x)，

其中，P(y|x)是条件概率，x是初始分类结果图，y是融合灰度与纹理信息的观测场，并且初始分类结果图x和融合灰度与纹理信息的观测场y是相互独立的；

在本发明的实施例中，利用独立分布原理，将条件概率P(y|x)表示为：P(y|x)＝P(y)，由此推出后验概率P(x|y)＝P(y)P(x)，后验概率P(x|y)为(301×301×2)像素，且其大小在0～1之间。

步骤6，根据差异图D的后验概率P(x|y)和MAP准则，得到变化检测结果图为：

$\hat{x}=\arg \max P\left(x\right|y)$

$=\arg \max \left\{P\left(x\right)P\left(y\right)\right\};$

在本发明的实施例中，得到变化检测结果图的大小为(301×301)像素。

步骤7，判断循环次数t是否达到了最高迭代次数g_max，如果满足t＞g_max则输出最终的变化检测结果图否则，令t＝t+1，返回步骤(3)，进行下一次迭代。

在本发明的实施例中，采用的最大迭代次数g_max＝10。

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

本发明的对比实验为经典的马尔科夫随机场MRF的变化检测方法，图像检测结果以多时相SAR图像进行对比。

1.实验条件：

本发明以三组多时相SAR图像进行实验，一组是分别在1999年4月和1999年5月通过ERS-2获得的瑞士Bern城市的水灾发生前后的SAR图像，两幅图像的尺寸均为301×301像素，灰度级为256，实际变化像素个数为1155。一组是分别在1997年5月和1997年8月通过Radarsat-1SAR获得的加拿大Ottawa地区的图像，两幅图像的尺寸均为350×290像素，灰度级为256，实际变化像素个数为16049。一组是分别在2008年6月和2009年6月由Rasarsat-2SAR获得的反应我国黄河入海口地域的变化情况中截取的典型变化区域，两幅图像的尺寸均为257×289像素，灰度级为256。对以上三组图像分别用马尔科夫随机场MRF方法和本发明方法进行变化检测的实现。

2.实验内容与结果：

实验1，是用本发明方法和MRF方法对第一组瑞士Bern城市水灾发生前后的SAR图像进行变化检测实验，结果如图2，其中，图2(a)为Bern1999年4月的原始图像，图2(b)为Bern1999年5月的原始图像，图2(c)为实际变化检测参考图，图2(d)为采用对比实验MRF方法得到的变化检测结果，图2(e)为采用本发明方法得到的变化检测结果。从图2(d)、图2(e)可以看出：本发明方法与MRF方法相比，不仅减少了伪变化信息，而且检测出较详细的边界信息，降低了错检数，使得到的变化检测结果图更接近参考图。

实验2，是用本发明方法和MRF方法对第二组加拿大Ottawa地区水灾发生前后的SAR图像进行变化检测实验，结果如图3，其中，图3(a)为Ottawa1997年5月的原始图像，图3(b)为Ottawa1997年8月的原始图像，图3(c)为实际变化检测参考图。图3(d)为采用对比实验MRF方法得到的变化检测结果，图3(e)为采用本发明方法得到的变化检测结果。从图3(d)、图3(e)可以看出：本发明方法与MRF方法相比，本发明方法有效地抑制了SAR图像斑点噪声的影响，使变化检测结果含有较少的杂点，对变化区域的细节部分能更好地实现，得到较好的检测结果，提高了变化检测结果的精确度。

实验3，是用本发明方法和MRF方法对第三组反应我国黄河入海口地域的SAR图像进行变化检测实验，结果如图4，其中，图4(a)为黄河入海口某一地区2008年6月的原始图像，图4(b)为黄河入海口某一地区2009年6月的原始图像，图4(c)为实际变化检测参考图。图4(d)为采用对比实验MRF方法得到的变化检测结果，图4(e)为采用本发明方法得到的变化检测结果。从图3(d)、图3(e)可以看出：本发明方法与MRF方法相比，本发明方法提高了变化检测的精确度，能稳定的对SAR图像进行变化检测，证明了本发明方法的有效性和稳定性。

实验4，是用本发明方法和MRF方法对第一组瑞士Bern城市水灾发生前后的SAR图像进行变化检测实验结果的评价，结果如表1。

表1Bern地区实验结果

从表1中数据可以看出：本发明方法与MRF方法相比，错检数FA减少了1911个像素点，漏检数MA增加了148个像素点，但总错误数OE减少了1763个像素点，Kappa系数增加了0.3562，检测正确率PCC增加了1.95%。

实验5，是用本发明方法和MRF方法对第二组加拿大Ottawa地区水灾发生前后的SAR图像进行变化检测实验结果的评价，结果如表2。

表2Ottawa地区实验结果

从表2中数据可以看出：本发明方法与MRF方法相比，错检数FA减少了1911个像素点，漏检数MA增加了624个像素点，但总错误数OE减少了2195个像素点，Kappa系数增加了0.0697，检测正确率PCC增加了2.16%。

实验6，是用本发明方法和MRF方法对第三组反应我国黄河入海口地域的SAR图像进行变化检测实验结果的评价，结果如表3。

表3黄河入海口实验结果

从表3中数据可以看出：本发明方法与MRF方法相比，错检数FA减少了2977个像素点，漏检数MA增加了576个像素点，但总错误数OE减少了1401个像素点，Kappa系数增加了0.0513，检测正确率PCC增加了1.89%。

综上，本发明方法将类别的先验概率和融合灰度与特征的模糊隶属概率看作是两个相互独立的分量，将变化检测问题看作是两个相互独立分量乘积的最大概率问题，有效地简化了算法实现的过程与复杂度，同时该方法利用融合差异图的灰度和纹理特征的模糊隶属概率作为观测量的似然概率，避免了用高斯模型拟合混合分布以及参数的估计，降低了SAR图像斑点噪声的影响，降低了伪变化信息，有效地提高了变化检测结果的精确度。

Claims (2)

1.一种基于先验和融合灰度与纹理特征的SAR图像变化检测方法，包括如下步骤：

(1)读入已配准和已校正的两幅不同时间同一地点的SAR图像I₁和I₂；

(2)用均值比值法构造这两幅SAR图像I₁和I₂的差异图D₁，用对数比值法构造这两幅SAR图像I₁和I₂的差异图D₂，利用小波变换对所述差异图D₁和所述差异图D₂进行融合，得到融合后差异图D，并初始化迭代次数t＝1；

(3)对融合后差异图D分别求类别先验概率P(x)和观测量似然概率P(y)，其中求观测量似然概率P(y)的步骤如下：

3.1)根据融合后差异图D进行非下采样小波变换，求取融合后差异图D的纹理特征矩阵G；

3.2)根据纹理特征矩阵G求取融合后差异图D的特征模糊隶属度矩阵U₁，U₁＝{u_ab,a＝1,...,k,b＝1,...,S}，k为融合后差异图D的分类类别数，S是融合后差异图D的总像素点个数，u_ab是像素点b隶属于第a类的特征模糊隶属度，同时由融合后差异图D中像素点的灰度值求取灰度模糊隶属度矩阵U₂，U₂＝{u_rh,r＝1,...,k,h＝1,...,S}，k为融合后差异图D的分类类别数，u_rh是像素点h隶属于第r类的灰度模糊隶属度；

3.3)根据特征模糊隶属度矩阵U₁和灰度模糊隶属度矩阵U₂，得到观测量的似然概率P(y)为：

$P\left(y\right)=\left\{\begin{array}{ccc}max(U_1,U_2),& if& k=1\\ min(U_1,U_2),& if& k=2\end{array}\right.$

其中，k是融合后差异图D的分类数，y是融合灰度与纹理信息的观测场；

(4)根据类别先验概率P(x)和观测量似然概率P(y)，利用贝叶斯公式和独立分布原理求取融合后差异图D的后验概率P(x|y)为：

P(x|y)＝P(y|x)P(x)

，

＝P(y)P(x)

其中，P(y|x)是条件概率，x是初始分类结果图，y是融合灰度与纹理信息的观测场，并且初始分类结果图x和融合灰度与纹理信息的观测场y是相互独立的；

(5)根据差异图D的后验概率P(x|y)和MAP准则，得到变化检测结果图为：

$\begin{array}{c}\hat{x}=\mathrm{argmax}P\left(x\right|y)\\ =\arg \max \left\{P\left(x\right)P\left(y\right)\right\}\end{array},$

(6)判断循环次数t是否达到了最高迭代次数g_max，如果满足t>g_max则输出最终的变化检测结果图否则，令t＝t+1，返回步骤(3)，进行下一次迭代。

2.根据权利要求1所述的SAR图像变化检测方法，其中步骤(3)所述的对融合后差异图D求类别先验概率P(x)，按照如下步骤实现：

3a)对融合后差异图D用K-means方法进行初始分割，得到初始分类结果图x；

3b)根据吉布斯分布求初始分类结果图x的先验概率P(x)＝{P(x_i),i＝1,...,L}，L是分类结果图x的总像素点个数，P(x_i)为分类结果图x中第i个像素点x_i的先验概率，表示为：

$P\left(x_i\right)=\frac{\exp (-u(x_i\left)\right)}{\underset{x_i\∈X}{\Σ}\exp (-u(x_i\left)\right)},$

其中，u(x_i)是第i个像素点x_i的能量函数，且u(x_i)＝-β∑δ(x_i,x_j)-1， $\mathrm{\δ}(x_i,x_j)=\left\{\begin{array}{cc}1& x_i=x_j\\ 0& x_i\≠x_j\end{array}\right.,$ β是平滑参数，x_j是第i个像素点x_i的邻域集N_i中的第j个像素点，且x_j∈N_i。