CN109697474B - 基于迭代贝叶斯的合成孔径雷达影像变化检测方法 - Google Patents

Abstract

一种基于迭代贝叶斯的合成孔径雷达影像变化检测方法，由构造差异图、局部特征提取、计算似然函数、进行迭代贝叶斯计算四个步骤对合成孔径雷达影像进行变化检测。其中，用对数比操作构造差异图；局部特征通过主成分分析方法提取，并且使用最小描述长度规则来确定特征的维数。本发明针对同一地区、不同时刻的两幅合成孔径雷达影像，引入贝叶斯推理的思想，首先对数比值差异图作为先验信息，然后通过和提取到的主要特征进行整合，得出变化信息的后验概率，使用最大类间方差法对后验概率进行处理来初始化下一次迭代的先验概率。该发明取得了较好的效果，具有较高的Kappa系数、正确率和G‑mean。

Description

基于迭代贝叶斯的合成孔径雷达影像变化检测方法

技术领域

本发明属于遥感影像应用技术领域，具体地涉及到合成孔径雷达影像变化检测方法。

背景技术

合成孔径雷达(SAR)具有全天候、全天时的特点，可方便地获得同一地区不同时刻的图像，于是SAR图像变化检测是它的一个重要的应用领域，它通过对不同时期图像的比较分析，根据图像之间的差异来得到我们所需要的变化信息。SAR图像变化检测技术可以应用于很多方面，例如对地震区域的定位和灾害评估、对农作物生长状况的监测、对土地使用的监测等等。

近年来，国内外学者在应用SAR图像变化检测问题上做了不少工作，但该领域仍然是目前许多专家关注的焦点。变化检测算法的归宿是检测图像发生显著变化的区域，舍弃那些不重要的、低置信度的伪变化。然而，由于合成孔径雷达系统采用合成相干方式成像，所获得的合成孔径雷达图像中包含大量的相干斑信息，这使得在做合成孔径雷达影像处理时，必须要充分考虑相干斑的影响，因此，急需要一种合成孔径雷达影像变化检测方法既要消除相干斑的影响又要保留真实变化信息不受到损失。

发明内容

本发明所要解决的技术问题在于提供一种鲁棒性好、正确率高、分类精度高的基于迭代贝叶斯的合成孔径雷达影像变化检测方法。

解决上述技术问题所采用的技术方案由下述步骤组成：

(1)用对数比操作对同一地区不同时刻的合成孔径雷达影像图I₁和I₂构造对数比值差异图D，作为先验概率p(s)；

(2)通过主成分分析和最小描述长度规则提取差异图D的K维特征，K为正整数；

(3)通过最小误差阈值法对先验概率p(s)进行二值化处理，得到估计的目标点集S_T和背景点集B_T，修改差异图D每维特征的重要性权重系数ω_i，根据已获得的估计的目标点集S_T和背景点集B_T以及每维特征的重要性权重系数ω_i，得到目标似然函数p(x|s)和背景似然函数p(y|b)：

式中，p(x_i|S_T)是第i维特征下目标点集S_T的分布函数，x_i是目标点集中含有第i维特征的像素，p(y_i|B_T)是第i维特征下背景点集B_T的分布函数，y_i是背景点集中含有第i维特征的像素；

(4)根据先验概率p(s)和目标似然函数p(x|s)及背景似然函数p(y|b)，通过贝叶斯公式整合全局的先验信息和多维的局部特征，得到后验概率p(s|x)为：

(5)根据后验概率得到灰度图，通过最大类间方差法确定分割图像的阈值，将灰度图变为包含变化和非变化两部分的二值图；

(6)对后验概率p(s|x)进行滤波得到下一次迭代的先验概率，迭代步骤(3)～步骤(5)，直到达到最大迭代次数T，停止迭代，输出最终变化检测结果图。

作为一种优选的技术方案，所述步骤(2)中主成分分析和最小描述长度规则提取差异图D的K维特征的方法如下：

(a)将差异图D按照3×3的像素进行非重叠分块；

(b)利用所有图像块进行奇异值分解，建立特征空间；

(c)使用最小描述长度规则确定特征空间的维数K：

式中p＝9，k＝1,2,...,p，λ_j是第j个特征值，n为将差异图D进行非重叠分块得到块的个数；

(d)对差异图D的每个像素点，利用周围3×3大小的数据块，向特征空间的K个主要方向进行投影，得到特征；

作为一种优选的技术方案，所述步骤(3)中修改差异图D每维特征的重要性权重系数ω_i的公式为：

式中，

是包含第i维特征的目标点集的平均值，/>

是包含第μ维特征的目标点集的平均值，/>

是包含第i维特征的背景点集的平均值，/>

是包含第μ维特征的背景点集的平均值。

本发明的有益效果如下：

本发明采用基于迭代贝叶斯的方法对同一地区、不同时刻的两幅合成孔径雷达影像进行变化检测，贝叶斯模型利用概率统计知识来整合先验信息和多维的局部特征，再通过迭代对预估计的变化信息进行不断的更新。似然函数中的权重系数反映了每一维特征对于整个联合概率分布函数的贡献，并且随着迭代过程被不断更新。主成分分析方法提取特征一方面利用了邻域像素，增加了特征信息；另一方面，使用最小描述长度规则降低了特征的维数，由于特征矢量互不相关，特征更有利于将变化和未变化的像素区分开。

因此与现有技术相比，本发明具有方法简单、参数选择合理、鲁棒性好、分类精度高等优点，有较好地变化检测效果。

附图说明

图1是本发明的流程图。

图2是实施例1变化前的影像。

图3是实施例1变化后的影像。

图4是实施例1的变化参考图。

图5是PCA+K-means方法对图2、图3进行变化检测的结果图。

图6是本发明对图2、图3进行变化检测的结果图。

图7是实施例2变化前的卡通图。

图8是实施例2变化后的卡通图。

图9是实施例2的变化参考图。

图10是对图7加入2视的散斑噪声的变化前的含噪卡通图。

图11是对图8加入2视的散斑噪声的变化后的含噪卡通图。

图12是PCA+K-means方法对图10、图11进行变化检测的结果图。

图13是本发明方法对图10、图11进行变化检测的结果图。

具体实施方式

下面结合附图和实施例对本发明进进一步详细说明，但本发明不限于下述实施方式。

实施例1

本实施例的图像是1997年5月和1997年8月获得的关于加拿大Ottawa地区的合成孔径雷达影像，如图2、图3，图像大小为352×313，变化参考图如图4，在图4中，白色部分为发生变化的参考区域，黑色部分为未发生变化的参考区域。如图1，本实施例的基于迭代贝叶斯的合成孔径雷达影像变化检测方法，具体操作步骤如下：

(1)用对数比操作对同一地区不同时刻的合成孔径雷达影像图I₁和I₂即图2和图3，构造对数比值差异图D，作为先验概率p(s)；

(2)通过主成分分析和最小描述长度规则提取差异图D的K维特征，步骤如下：

(a)将差异图D按照3×3的像素进行非重叠分块；

(b)利用所有图像块进行奇异值分解，建立特征空间；

(c)使用最小描述长度规则确定特征空间的维数K：

式中p＝9，k＝1,2,...,p，λ_j是第j个特征值，n为将差异图D进行非重叠分块得到块的个数，n＝12168；

(d)对差异图D的每个像素点，利用周围3×3大小的数据块，向特征空间的K＝5个主要方向进行投影，得到特征；

(3)通过最小误差阈值法对先验概率p(s)进行二值化处理，得到估计的目标点集S_T和背景点集B_T，修改差异图D每维特征的重要性权重系数ω_i，根据已获得的估计的目标点集S_T和背景点集B_T以及每维特征的重要性权重系数ω_i，得到目标似然函数p(x|s)和背景似然函数p(y|b)：

/>

p(x_i|S_T)是第i维特征下目标点集S_T的分布函数，x_i是目标点集中含有第i维特征的像素，p(y_i|B_T)是K＝5维特征下背景点集B_T的分布函数，y_i是背景点集中含有第i维特征的像素；

本实施例修改差异图D每个特征的重要性权重系数ω_i的公式为：

式中，

是包含第i维特征的目标点集的平均值，/>

是包含第μ维特征的目标点集的平均值，/>

是包含第i维特征的背景点集的平均值，/>

是包含第μ维特征的背景点集的平均值；

(4)根据先验概率p(s)和目标似然函数p(x|s)及背景似然函数p(y|b)，通过贝叶斯公式整合全局的先验信息和多维的局部特征，得到后验概率p(s|x)为：

(5)根据后验概率得到灰度图，通过最大类间方差法确定分割图像的阈值，将灰度图变为包含变化和非变化两部分的二值图；

(6)对后验概率进行滤波得到下一次迭代的先验概率，迭代步骤(3)～步骤(5)，直到达到最大迭代次数T＝5，停止迭代，输出最终变化检测结果图，如图6。

实施例2

本实施例的图像是给原卡通图即图7和图8添加2视的散斑噪声得到的合成散斑图像即图10和图11，图像大小为512×512。

本实施例基于迭代贝叶斯的合成孔径雷达影像变化检测方法，具体操作步骤如下：

(1)用对数比操作对同一地区不同时刻的合成孔径雷达影像图I₁和I₂即图10和图11，构造对数比值差异图D，作为先验概率p(s)；

(2)通过主成分分析和最小描述长度规则提取差异图D的K维特征，步骤如下：

(a)将差异图D按照3×3的像素进行非重叠分块；

(b)利用所有图像块进行奇异值分解，建立特征空间；

(c)使用最小描述长度规则确定特征空间的维数K：

式中p＝9，k＝1,2,...,p，λ_j是第j个特征值，n为将差异图D进行非重叠分块得到块的个数，n＝28900；

(d)对差异图D的每个像素点，利用周围3×3大小的数据块，向特征空间的K＝6个主要方向进行投影，得到特征；

步骤(3)～(5)与实施例1相同；

(6)对后验概率进行滤波得到下一次迭代的先验概率，迭代步骤(3)～步骤(5)，直到达到最大迭代次数T＝5，停止迭代，输出最终变化检测结果图，如图13。

为了验证本发明的有益效果，发明人进行了2组实验如下：

1、应用PCA+K-means方法对图2和图3进行变化检测实验，该实验结果如图5，该实验与本发明实施例1的变化检测评价指标对比，如表1。

表1

	FA	MA	Kappa	OA	F1	TPR	TNR	G-mean
									PCA+K-means	957	1948	0.8905	0.9736	0.9058	0.8776	0.9899	0.9321
实施例1	1076	1537	0.9029	0.9763	0.9167	0.9034	0.9886	0.9451

2、应用PCA+K-means方法对图10和图11进行变化检测实验，该实验结果如图12，该实验与本发明实施例2的变化检测评价指标对比，如表2。

表2

	FA	MA	Kappa	OA	F1	TPR	TNR	G-mean
									PCA+K-means	842	3664	0.8639	0.9828	0.8731	0.8088	0.9965	0.8978
实施例2	956	3319	0.8724	0.9837	0.8811	0.8268	0.9961	0.9075

实验结论

由表1、表2可见，本发明获得了MA、Kappa、OA、F1、TPR和G-mean的最佳值，F1同时兼顾了分类模型的查准率和查全率，G-mean同时兼顾了分类模型的TPR和TNR，能衡量不平衡数据的模型表现，本发明的F1和G-mean都较高，说明在不平衡的数据模型上，本发明的性能也很好，并且本发明具有较高的Kappa系数值，这意味着本发明在实际的多时相SAR图像中的性能也是合理的。本发明与现有技术相比，具有方法简单、参数选择合理、鲁棒性好、分类精度高等优点，可见，用主成分分析方法提取特征，在迭代过程中不断更新每维特征的重要性权重，使用贝叶斯模型来整合先验信息和多维的局部特征，并且通过迭代进行更新，具有较好地变化检测效果。

Claims (2)

1.一种基于迭代贝叶斯的合成孔径雷达影像变化检测方法，其特征在于由下述步骤组成：

(1)用对数比操作对同一地区不同时刻的合成孔径雷达影像图I₁和I₂构造对数比值差异图D，作为先验概率p(s)；

(2)通过主成分分析和最小描述长度规则提取差异图D的K维特征，K为正整数；

所述主成分分析和最小描述长度规则提取差异图D的K维特征的方法如下：

(a)将差异图D按照3×3的像素进行非重叠分块；

(b)利用所有图像块进行奇异值分解，建立特征空间；

(c)使用最小描述长度规则确定特征空间的维数K：

式中p＝9，k＝1,2,...,p，λ_j是第j个特征值，n为将差异图D进行非重叠分块得到块的个数；

(d)对差异图D的每个像素点，利用周围3×3大小的数据块，向特征空间的K个主要方向进行投影，得到特征；

(3)通过最小误差阈值法对先验概率p(s)进行二值化处理，得到估计的目标点集S_T和背景点集B_T，修改差异图D每维特征的重要性权重系数ω_i，根据已获得的估计的目标点集S_T和背景点集B_T以及每维特征的重要性权重系数ω_i，得到目标似然函数p(x|s)和背景似然函数p(y|b)：

式中，p(x_i|S_T)是第i维特征下目标点集S_T的分布函数，x_i是目标点集中含有第i维特征的像素，p(y_i|B_T)是第i维特征下背景点集B_T的分布函数，y_i是背景点集中含有第i维特征的像素；

(4)根据先验概率p(s)和目标似然函数p(x|s)及背景似然函数p(y|b)，通过贝叶斯公式整合全局的先验信息和多维的局部特征，得到后验概率p(s|x)为：

(5)根据后验概率得到灰度图，通过最大类间方差法确定分割图像的阈值，将灰度图变为包含变化和非变化两部分的二值图；

(6)对后验概率p(s|x)进行滤波得到下一次迭代的先验概率，迭代步骤(3)～步骤(5)，直到达到最大迭代次数T，停止迭代，输出最终变化检测结果图。

2.根据权利要求1所述的一种基于迭代贝叶斯的合成孔径雷达影像变化检测方法，其特征在于所述步骤(3)中修改差异图D每维特征的重要性权重系数ω_i的公式为：

式中，

是包含第i维特征的目标点集的平均值，/>

是包含第μ维特征的目标点集的平均值，/>

是包含第i维特征的背景点集的平均值，/>

是包含第μ维特征的背景点集的平均值。/>

Images