CN101493520B - 一种基于二维Gamma分布的SAR图像变化检测方法 - Google Patents

Abstract

一种基于二维Gamma分布的SAR图像变化检测方法，包括以下步骤：根据输入的SAR待测图像和参考图像数据，利用矩估计法估计二维Gamma分布的参数；根据Neyman-Pearson准则构造似然比统计量，在二维Gamma分布的基础上利用图像数据相关性进行杂波抑制，得到杂波抑制后的图像；对杂波抑制后的图像进行CFAR归一化，并设定全局阈值将图像二值化，得到初步的检测结果；对检测后得到的二值图像进行形态学处理、计数滤波及目标聚类，进一步消除孤立虚警点，得到最终的检测结果；本发明在低虚警率的基础上达到了较高的检测率，适合在各种杂波环境，特别是强杂波环境下检测人造目标。

Description

一种基于二维Gamma分布的SAR图像变化检测方法

技术领域

本发明属于SAR图像处理领域，涉及一种基于二维Gamma分布的SAR图像变化检测方法。

背景技术

合成孔径雷达(Synthetic Aperture Radar，SAR)本身是一种利用微波感知的主动式传感器，它不受天气、光照等条件限制，可对感兴趣目标进行全天候、全天时的侦察；此外，由于微波具有一定的穿透探测能力，可以探测到一定深度地表下的或是其他经过伪装或掩盖的目标，这使得SAR在目标探测方面具有巨大的应用潜力。

在杂波较强的环境中，单一的CFAR目标检测技术受到了较大的制约。随着高分辨、短周期重复观测等技术走向成熟，变化检测技术得到了快速的发展，应用于目标检测的变化检测技术也成为了重要的研究方向。为了在大场景中检测出人造目标，像元级的直接比较检测法是比较可行且可靠的一类方法，即直接比较同一位置上不同时相的像元特征值来检测变化，通常先通过数学变换的方式产生差异图像，再对差异图像进行阈值化处理，从中提取变化区域。

在这类变化检测方法中，较常用的基于一维分布模型的图像差值法或比值法虽然简单易行，但是无法充分利用不同时相图像数据的相关性；而在利用二维分布模型的检测方法中，基于二维Gaussian分布的杂波抑制变化检测法虽然利用图像相关性进行杂波抑制，但是SAR图像一般并不是简单地服从Gaussian分布，因此这种方法存在着分布模型的拟合精度问题；而基于二维Gamma分布的相关系数变化检测法，虽然利用了拟合精度较高的分布模型，但是仅利用相关系数作为判决依据，难以设定合理的检测阈值。

发明内容

本发明的要解决的技术问题是：克服现有技术的不足，提供一种基于二维Gamma分布的SAR图像变化检测方法，该方法利用拟合精度较高的二维Gamma分布模型改善算法的检测性能；根据Newman-Pearson准则构造似然比统计量抑制杂波，减小背景图像强散射区域中虚警目标的产生；同时，利用CFAR对抑制杂波后的图像进行归一化，阈值设置简单快捷。

本发明的解决其技术问题所采用的技术方案为：一种基于二维Gamma分布的SAR图像变化检测方法，其特征在于以二维Gamma分布为基础，构造似然比统计量抑制杂波，并通过CFAR归一化设置检测阈值并进行形态学处理，从而在大场景中检测出人造目标，其具体步骤如下：

(1)根据输入的SAR待测图像和参考图像数据，利用矩估计法估计二维Gamma分布的参数；

(2)根据Neyman-Pearson准则构造似然比统计量，在二维Gamma分布的基础上利用图像数据相关性进行杂波抑制，得到杂波抑制后的图像；

(3)对杂波抑制后的图像进行CFAR归一化，并设定全局阈值将图像二值化，得到初步的检测结果；

(4)对检测后得到的二值图像进行形态学处理、计数滤波及目标聚类，进一步消除孤立虚警点，得到最终的检测结果；

所述的步骤(4)中，为进一步消孤立虚警点，首先对二值化的检测结果进行腐蚀膨胀操作；再根据目标尺寸大小与SAR图像分辨率进行计数滤波；最后再进行目标聚类，得到最终的检测结果；

所述的步骤(2)中，采用如下方式进行杂波抑制：

${(1-\frac{\mathrm{\&zeta;}}{z_1})}^{q-1}\&CenterDot;\frac{f_q({\mathrm{cz}}_1{z}_2-\mathrm{c\&zeta;}{z}_2)}{f_q\left({\mathrm{cz}}_1{z}_2\right)}=\left\{\begin{array}{cc}>\mathrm{\&lambda;}& \&RightArrow;H_1\\ <\mathrm{\&lambda;}& \&RightArrow;H_0\end{array}\right.$

其中z₁和z₂分别为待测图像和参考图像数据，参数c与二维.Gamma分布参数有关，而参数ζ为根据图像强度值数量级调节的比例系数，典型取值为：

$\mathrm{\&zeta;}=\frac{(E\left(z_1\right)+E\left(z_2\right))}{2}.$

本发明的原理是：根据乘性噪声模型，SAR多视强度图像服从一维Gamma分布。由于同一场景不同时相的SAR图像间具有一定的相关性，因此选择二维Gamma分布作为待测图像和参考图像数据的联合分布模型。在此基础上，根据Neyman-Pearson准则构造似然比统计量，即等效地利用了待测图像和参考图像的相关性进行了杂波抑制。对杂波抑制后的图像进行CFAR归一化，并设定全局阈值，则可实现变化检测的目的。

本发明与现有技术相比的优点在于：较之于基于一维分布的SAR图像差值法或比值法，本发明更加充分地利用了不同时相图像数据之间的相关性；较之于基于二维Gaussian分布的杂波抑制变化检测法，本发明采用了更加适合SAR图像的二维Gamma分布，通过拟合精度较高的分布模型来改善检测性能；最后，较之于基于二维Gamma分布的相关系数变化检测法，本发明并不只是简单地利用相关系数作为判决依据，而是利用相关性抑制杂波，同时CFAR归一化也使本发明的阈值设定更加合理、方便、快捷。

附图说明

图1为本发明的SAR图像变化检测方法流程图；

图2为CFAR归一化的模板示意图；

图3为形态学处理、计数滤波和目标聚类的流程图。

具体实施方式

下面结合附图及具体实施方式详细介绍本发明。

如图1所示，本发明的基于二维Gamma分布的SAR图像变化检测方法的具体实施步骤如下：

(1)、根据输入的SAR待测图像和参考图像数据，利用矩估计法估计二维Gamma分布的参数。对于多维向量x＝(x₁，…，x_d)^T，若任意x_i的边缘分布都服从一维Gamma分布，则认为向量x服从多维Gamma分布。但是服从这一条件的分布族是非常庞大的，因此本发明采用S.barlev和P.bernardoff的定义方式(矩母函数和拉普拉斯变换)进行进一步限制，定义如下：

假设P(z)为关于多维向量Z＝(z₁，…，z_d)的多项式，若对于任意的

都可表示为Az_j+B，其中A和B为关于其它z_i(i≠j)的多项式，则多项式P(z)被称为是仿射的。若随机向量x＝(x₁，…，x_d)^T的矩母函数(也称为拉普拉斯变换)定义如式(1)，则称该随机向量服从形状参数为q、尺度参数为P的多维Gamma分布，记为x～G(q，P)。

${\mathrm{\&psi;}}_{G(q,P)}=E\left(e^{-{\mathrm{\&Sigma;}}_{i=1}^d{x}_i{z}_i}\right)={\left[P\left(z\right)\right]}^{-q}---\left(1\right)$

其中q≥0，P为一仿射多项式。当d＝2时，对应的是二维的情况，此时相应的矩母函数为：

Ψ(z₁，z₂)＝(1+p₁z₁+p₂z₂+p₁₂z₁z₂)^-q    (2)

其中p₁＞0，p₂＞0，p₁₂＞0，p₁p₂-p₁₂≥0，由此可得二维Gamma分布的概率密度函数表达式为：

$f_{2D}\left(x\right)=\exp (-\frac{p_2{x}_1+p_1{x}_2}{p_{12}})\frac{x_1^{q-1}{x}_2^{q-1}}{p_{12}^q\mathrm{\&Gamma;}\left(q\right)}{f}_q(\frac{p_1{p}_2-p_{12}}{p_{12}^2}\&CenterDot;x_1{x}_2)---\left(3\right)$

其中函数f_q(z)与合流超几何函数相关，定义如式(4)：

$f_q\left(z\right)=\underset{k=0}{\overset{\&infin;}{\mathrm{\&Sigma;}}}\frac{z^k}{k!\mathrm{\&Gamma;}(q+k)}---\left(4\right)$

对矩母函数进行泰勒级数展开，有：

$\mathrm{\&psi;}(z_1,z_2)=\underset{k,l\&GreaterEqual;0}{\mathrm{\&Sigma;}}\frac{{(-1)}^{k+l}}{k!l!}E\left[x_1^k{x}_2^l\right]z_1^l{z}_2^l---\left(5\right)$

对式(5)作微分，则可得到二维Gamma分布的各阶矩。因此，可利用矩估计法对二维Gamma分布进行参数估计。构造

可得统计量x的矩为：

$E\left(x\right)=[m_1,m_2,m_1^2(1+q^{-1}),m_2^2(1+q^{-1}),m_1{m}_2(1+{\mathrm{rq}}^{-1})]---\left(6\right)$

用(f₁，f₂，f₃，f₄，f₅)表示E(x)，则可推出如下的关系式：

$m_1=f_1,m_2=f_2,r=\frac{f_5-f_1{f}_2}{\sqrt{\left(f_3{-f}_1^2\right)(f_4-f_2^2)}}---\left(7\right)$

由此，可进一步推出参数p₁、p₂和p₁₂的估计式如下：

$\left\{\begin{array}{c}\\ {\hat{p}}_1=E\left(x_1\right)\\ {\hat{p}}_2=E\left(x_2\right)\\ {\hat{p}}_{12}=E\left(x_1\right)E\left(x_2\right)\&CenterDot;(1-\frac{\mathrm{cov}(x_1,x_2)}{\sqrt{\mathrm{var}\left(x_1\right)}\sqrt{\mathrm{var}\left(x_2\right)}})\end{array}\right.---\left(8\right)$

式中E(·)表示数学期望，var(·)表示方差，cov(·，·)表示协方差，因此参数p₁、p₂反映了单个变量幅度的大小，而p₁₂则反映了两变量间相关性程度。

(2)、根据Neyman-Pearson准则构造似然比统计量，在二维Gamma分布的基础上利用图像数据相关性进行杂波抑制，得到杂波抑制后的图像。对于不同时刻获得的SAR图像数据，假定目标、杂波和噪声模型表示如下：

$s=\left[\begin{array}{c}{s}_1\\ s_2\end{array}\right]$ $c=\left[\begin{array}{c}{c}_1\\ c_2\end{array}\right]$ $n=\left[\begin{array}{c}{n}_1\\ n_2\end{array}\right]---\left(9\right)$

其中下标分别对应着两个不同时刻获得的回波信号。同时，对于有无变化发生作如下假设：

H₀：z＝q  (no change)

                            (10)

H₁：z＝s+q  (change)

式中z表示SAR图像数据，q＝c+n表示背景杂波和噪声。即只有背景杂波和噪声，而没有目标时，表示没有变化发生；当有目标存在时，表示有变化发生。为了在固定的虚警率下达到最大的检测概率，根据N-P准则，采用似然比检验的方法，构造统计量如下：

$\mathrm{\&Lambda;}\left(z\right)=\frac{P\left(z\right|H_1)}{P\left(z\right|H_0)}---\left(11\right)$

其中P表示图像数据z分别在H₀和H₁下的条件概率。此时，根据二维Gamma分布，可得不同假设条件下相应的概率密度函数表达式为：

$P\left(z\right|H_0)=\exp (-\frac{p_2{z}_1+p_1{z}_2}{p_{12}})\frac{z_1^{q-1}{z}_2^{q-1}}{p_{12}^q\mathrm{\&Gamma;}\left(q\right)}{f}_q\left({\mathrm{cz}}_1{z}_2\right)---\left(12\right)$

$P\left(z\right|H_1)=\exp [-\frac{p_2(z_1-s_1)+p_1(z_2-s_2)}{p_{12}}]\&CenterDot;\frac{{(z_1-s_1)}^{q-1}{(z_2-s_2)}^{q-1}}{p_{12}^q\mathrm{\&Gamma;}\left(q\right)}---\left(13\right)$

$\&times;f_q\left[c(z_1-s_1)(z_2-s_2)\right]$

其中

$f_q\left(z\right)=\underset{k=0}{\overset{\&infin;}{\mathrm{\&Sigma;}}}\frac{z^k}{k!\mathrm{\&Gamma;}(q+k)},$ $c=\frac{p_1{p}_2-p_{12}}{p_{12}^2}---\left(14\right)$

而p₁、p₂和p₁₂是待估计的二维Gamma分布的参数，q则为已知的SAR图像视数。去除一些与变量z无关的项，式(11)可进一步等效为：

${\left[\frac{(z_1-s_1)(z_2-s_2)}{z_1{z}_2}\right]}^{q-1}\&CenterDot;\frac{f_q\left(c(z_1-s_1)(z_2-s_2)\right)}{f_q\left({\mathrm{cz}}_1{z}_2\right)}=\left\{\begin{array}{cc}>\mathrm{\&lambda;}& \&RightArrow;H_1\\ <\mathrm{\&lambda;}& \&RightArrow;H_0\end{array}\right.---\left(15\right)$

此时，虽然目标信号s是未知的，但是对于目标与场景的关系来说，不外乎两种情况：目标存在与目标不存在。因此可以简单地假设目标信号模型为：

$s=\left[\begin{array}{c}1\\ 0\end{array}\right]$ 或 $s=\left[\begin{array}{c}0\\ 1\end{array}\right]---\left(16\right)$

为了便于说明，不妨假设目标只出现在待测图像中而非参考图像中，取式(16)中的前者，即s＝[1 0]^T。但由于图像数据z与目标信号s直接进行加减处理，考虑到数量级可能的差异，可将目标信号模型乘以一个相应的系数，如下：

$s=\left[\begin{array}{c}\mathrm{\&zeta;}\\ 0\end{array}\right]\mathrm{\&eta;}\left(E\right\{z_1\}+E\{z_2\left\}\right)\&CenterDot;\left[\begin{array}{c}1\\ 0\end{array}\right]---\left(17\right)$

其中η为比例系数，可适当调节，典型值可取0.5。将目标信号代入式(15)中，可简化为：

${(1-\frac{\mathrm{\&zeta;}}{z_1})}^{q-1}\&CenterDot;\frac{f_q({\mathrm{cz}}_1{z}_2-\mathrm{c\&zeta;}{z}_2)}{f_q\left({\mathrm{cz}}_1{z}_2\right)}=\left\{\begin{array}{cc}>\mathrm{\&lambda;}& \&RightArrow;H_1\\ <\mathrm{\&lambda;}& \&RightArrow;H_0\end{array}\right.---\left(18\right)$

其中z₁和z₂分别为待测图像和参考图像数据，参数c与二维Gamma分布参数有关，而参数ζ为根据图像强度值数量级调节的比例系数，典型取值为：

$\mathrm{\&zeta;}=\frac{(E\left(z_1\right)+E\left(z_2\right))}{2},$ $c=\frac{p_1{p}_2-p_{12}}{p_{12}^2}---\left(19\right)$

式中E(·)表示数学期望。

根据式(18)等号左端的处理方式，利用参考图像数据对待测图像进行遍历，则可得到杂波抑制后的结果。

(3)、对杂波抑制后的图像进行CFAR归一化，并设定全局阈值将图像二值化，得到初步的检测结果。为了进一步确定全局阈值λ，本发明利用图2的CFAR检测模板对杂波抑制后的图像进行归一化，从而达到恒虚警的目的。以每一个待检测的像素为中心，按照式(20)归一化并进行判决：

$\frac{x-E}{\mathrm{\&sigma;}}=\left\{\begin{array}{cc}>\mathrm{\&lambda;}& \&RightArrow;t\arg \mathrm{et}\\ <\mathrm{\&lambda;}& \&RightArrow;\mathrm{clutter}\end{array}\right.---\left(20\right)$

其中x为杂波抑制后的图像数据，E和σ分别为空心滑窗内样本数据的均值和标准差。通过CFAR归一化处理可以方便快捷的设置全局阈值，得到初步的检测结果。

(4)、对检测后得到的二值图像进行形态学处理、计数滤波及目标聚类，进一步消除孤立虚警点，得到最终的检测结果。在SAR图像中，由于乘性斑点噪声的影响，在初步检测得到的二值图像中，往往会存在一些孤立的虚警点。同时，在高分辨的SAR图像中，目标呈分布式，一个目标往往占有多个分布单元；而目标本身散射特性的起伏使得对应于同一目标的像素点往往不能形成连通区域，造成目标的误判。因此本发明按照附图3所示的流程图通过一些形态学的方法对检测结果进行后处理。首先进行形态学处理，此处采用腐蚀膨胀操作，填补连通区域、消除一些孤立的点目标；接着根据目标尺寸大小与SAR图像的分辨率，进行计数滤波，进一步消除虚警点；最后再进行目标聚类，若任意两不连通区域距离小于门限，则认为其属于同一个目标，进行目标合并，从而得到最终的检测结果。

本发明说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。

尽管为说明目的公开了本发明的最佳实施例和附图，但是本领域的技术人员可以理解：在不脱离本发明及所附的权利要求的精神和范围内，各种替换、变化和修改都是可能的。因此，本发明所保护的技术方案不应局限于最佳实施例和附图所公开的内容。

Claims (3)

1.一种基于二维Gamma分布的SAR图像变化检测方法，其特征在于：包括以下几个步骤：

(1)根据输入的SAR待测图像和参考图像数据，利用矩估计法估计二维Gamma分布的参数；

(2)根据Neyman-Pearson准则构造似然比统计量，在二维Gamma分布的基础上利用图像数据相关性进行杂波抑制，得到杂波抑制后的图像；

(3)对杂波抑制后的图像进行CFAR归一化，并设定全局阈值将图像二值化，得到初步的检测结果；

(4)对检测后得到的二值图像进行形态学处理、计数滤波及目标聚类，进一步消除孤立虚警点，得到最终的检测结果；

所述的步骤(4)中，为进一步消除孤立虚警点，首先对二值化的检测结果进行腐蚀膨胀操作；再根据目标尺寸大小与SAR图像分辨率进行计数滤波；最后再进行目标聚类，得到最终的检测结果；

所述的步骤(2)中，采用如下方式进行杂波抑制：

H₀：z＝q

H₁：z＝s+q

式中z表示SAR图像数据，q＝c+n表示背景杂波和噪声；

其中z₁和z₂分别为待测图像和参考图像数据，参数c与二维Gamma分布参数有关，而参数ζ为根据图像强度值数量级调节的比例系数，典型取值为：

式中E(·)表示数学期望；p₁、p₂和p₁₂是待估计的二维Gamma分布的参数，q则为已知的SAR图像视数，λ为全局阈值。

2.根据权利要求1所述的基于二维Gamma分布的SAR图像变化检测方法，其特征在于：所述的步骤(1)的二维Gamma分布采用的是S.barlev和P.bernardoff通过矩母函数和拉普拉斯变换的方式定义的Gamma分布函数族，二维Gamma分布的概率密度函数表达式为：

其中q为已知的SAR图像视数，函数f_q(z)与合流超几何函数相关，表达式为：

参数估计采用矩估计法，估计式为：

式中E(·)表示数学期望，var(·)表示方差，cov(·，·)表示协方差，因此参数p₁、p₂反映了单个变量幅度的大小，而p₁₂则反映了两变量间相关性程度。

3.根据权利要求1所述的基于二维Gamma分布的SAR图像变化检测方法，其特征在于：所述的步骤(3)中，为进一步明确阈值λ，采用空心滑窗对抑制杂波后的数据进行CFAR归一化：

其中x为杂波抑制后的图像数据，E和σ分别为空心滑窗内样本数据的均值和标准差；设定全局阈值λ对CFAR归一化后的数据作判决，得到二值化的检测结果。 

Images