CN103500453B - 基于伽玛分布和邻域信息的sar图像显著性区域检测方法 - Google Patents

Abstract

本发明公开了一种基于伽玛分布和邻域信息的SAR图像显著性区域检测方法，主要解决斑点噪声下现有算法不能稳定有效检测SAR图像显著性区域的问题。其实现步骤是：(1)利用8个灰度变化值检测窗口得到像素点的变化值和变化方向；(2)在不同尺度下构建像素点的变化值和变化方向的二维直方图对，计算邻域变化不一致性度量，并确定像素点邻域的显著性和显著性尺度；(3)通过斑点噪声的伽玛分布模型估计像素点在显著性尺度邻域的灰度直方图，计算显著性尺度的局部显著性度量以及像素点的显著性度量；(4)通过迭代方法得到稳定的显著性区域坐标及其区域半径。本发明减少了斑点噪声的影响，检测的稳定性以及有效性均得到提高，可用于SAR图像配准、变化检测和目标识别。

Description

基于伽玛分布和邻域信息的SAR图像显著性区域检测方法

技术领域

本发明属于图像处理技术领域，涉及SAR图像显著性区域检测，可用作SAR图像配准、变化检测和目标识别。

背景技术

合成孔径雷达(syntheticapertureradar,SAR)系统具有全天时，全天候，穿透性等特点，是一种重要的遥感数据来源。随着SAR图像在军事和民用领域应用需求的增长，对SAR图像进行自动分析和解译的需求也越来越多。利用SAR图像的局部特征对图像内容进行描述，不但可以减少图像分析的复杂度，还可以为SAR图像的匹配、变化检测、目标识别与分类打下基础。因此，有效的、准确的局部图像特征提取技术可以提高SAR图像的自动解译性能。

通过对人类视觉系统的研究，人们发现早期视觉首先感知突出于相邻区域的局部特征。这种局部特征具有某种形式的空间不连续性，被称为显著性区域。人类视觉系统将该特征作为基本元素对场景进行进一步的分析和理解。因此，在图像的分析和理解中，显著性区域也被当作一种有效的局部特征用于描述图像内容。针对显著性区域，Kadir和Brady提出了一种可靠的检测算法(Kadir.T,Brady.M.Saliency,ScaleandImageDescription.InternationalJournalofComputerVision,45(2),83-105,2001)。该算法步骤为：首先利用灰度直方图计算像素点不同尺度邻域的信息熵，并将其定义为不同尺度的局部显著性度量；对比这些熵值，判断该像素点的显著性，将局部最大熵的尺度定义为显著性尺度；利用像素点的显著性尺度邻域与其相邻尺度邻域的灰度直方图计算它们之间的不一致性，并将其定义为显著性尺度下的不一致性度量；最后，联合显著性尺度的两种度量给出像素点在该尺度的显著性度量。相比其它众多局部区域检测方法，该显著性区域检测算法只依赖于图像内容本身，从而避免了对先验信息的依赖。并且，该显著性区域检测算法利用特定大小的图像区域信息作为区域提取依据，它的结果更为鲁棒。但研究发现，将该方法应用于SAR图像时，由于SAR系统成像时产生的大量斑点噪声对灰度直方图的影响，得到的像素点的显著性判断及显著性度量并不准确，从而影响了显著性区域的位置和尺度的判定。

发明内容

本发明的目的在于克服上述已有问题的缺点，提出一种基于伽玛分布和邻域信息的SAR图像显著性区域检测方法，以提高像素点的显著性判断及显著性度量的准确性，避免出现对区域显著性及其坐标和尺度的误判，为后续图像分析和理解奠定良好的基础。

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

(1)输入一幅I×J大小的SAR图像，利用8个不同的3×3灰度变化值检测窗口，计算该图像的每个像素点在8个方向上的灰度变化值；

(2)将每个像素点在8个方向上灰度变化值的最大值，记录为该像素点的变化值，并将灰度变化值的最大值产生方向，记录为该像素点的变化方向，其中变化方向标记为1到8的整数；

(3)将坐标为第i行第j列的像素点(i,j)的变化值和变化方向分别表示为xc_i,j和xd_i,j，1≤i≤I，1≤j≤J；

(4)根据每个像素点(i,j)的变化值和变化方向，1≤i≤I，1≤j≤J，计算初始显著性矩阵Y_T：

4a)计算像素点(i,j)的每个尺度r的二维直方图对{ch_i,j(xc,xd,r),h_i,j(xc,xd,r)}，2≤r≤r_max：

${\mathrm{ch}}_{i,j}(\mathrm{xc},\mathrm{xd},r)=\underset{(m,n)\∈R_{i,j}^{r-1}}{\mathrm{\Σ}}{\mathrm{\δ}}_{{\mathrm{xc},\mathrm{xc}}_{m,n}}\×{\mathrm{\δ}}_{{\mathrm{xd},\mathrm{xd}}_{m,n}}\×\frac{\mathrm{xc}}{{\mathrm{xc}}_{\max }},$

$\begin{array}{c}{h}_{i,j}(\mathrm{xc},\mathrm{xd},r)={\mathrm{ch}}_{i,j}(\mathrm{xc},\mathrm{xd},r)+\underset{(m,n)\∈(R_{i,j}^{r-1},R_{i,j}^r]}{\mathrm{\Σ}}{\mathrm{\δ}}_{\mathrm{xc},{\mathrm{xc}}_{m,n}}\×{\mathrm{\δ}}_{{\mathrm{xd},\mathrm{xd}}_{m,n}}\\ \×\frac{(g_{\min }+g\×r-g)(2\×{g\×g}_{\min }+g^2)}{g_{\min }^2\×g\×(2\×g_{\min }+2\×g\×r-g)}\×\frac{\mathrm{xc}}{{\mathrm{xc}}_{\max }}\end{array},$

其中，ch_i,j(xc,xd,r)为像素点(i,j)在尺度r下对应的尺度r-1邻域内像素点的变化值和变化方向二维直方图，h_i,j(xc,xd,r)为像素点(i,j)的尺度r邻域内像素点的变化值和变化方向二维直方图，xc包含图像上像素点的变化值的所有取值，xd包含像素点的8个变化方向标记，xc_max为整个图像上最大的灰度变化值，xc_m,n为像素点(m,n)的变化值，xd_m,n为像素点(m,n)的变化方向，表示像素点(m,n)包含于距离像素点(i,j)不超过g_min+g×(r-1)-g个像素点数的圆形邻域内，表示像素点(m,n)包含于距离像素点(i,j)大于g_min+g×(r-1)-g个像素点数又不超过g_min+g×r-g个像素点数的环形邻域内，g_min是尺度1邻域半径的像素点数，g是相邻尺度邻域半径增量的像素点数，为关于xc和xc_m,n的克罗内克函数， ${\mathrm{\δ}}_{\mathrm{xc},{\mathrm{xc}}_{m,n}}=\left\{\begin{array}{cc}1,& {\mathrm{xc}}_{m,n}=\mathrm{xc}\\ 0,& {\mathrm{xc}}_{m,n}\≠\mathrm{xc}\end{array}\right.,{\mathrm{\δ}}_{\mathrm{xd},{\mathrm{xd}}_{m,n}}$ 为关于xd和xd_m,n的克罗内克函数， ${\mathrm{\δ}}_{\mathrm{xd},{\mathrm{xd}}_{m,n}}=\left\{\begin{array}{cc}1,& {\mathrm{xd}}_{m,n}=\mathrm{xd}\\ 0,& {\mathrm{xd}}_{m,n}\≠\mathrm{xd}\end{array}\right.;$

4b)利用4a)的(r_max-1)个二维直方图对，计算像素点(i,j)的每个尺度r的邻域变化不一致性度量2≤r≤r_max：

$\begin{array}{c}{W}_{i,j}^r=\underset{\mathrm{xc}=0}{\overset{{\mathrm{xc}}_{\max }}{\mathrm{\Σ}}}\underset{\mathrm{xd}=1}{\overset{8}{\mathrm{\Σ}}}|{\mathrm{ch}}_{i,j}(\mathrm{xc},\mathrm{xd},r)/g_{\min }^2-h_{i,j}(\mathrm{xc},\mathrm{xd},r)/(g_{\min }^2+2g\×g_{\min }+g^2)|\\ \×g_{\min }^2/{(g_{\min }+g\×r-g)}^2\end{array};$

4c)从所有求得的(r_max-1)个邻域变化不一致性度量中找出最大的邻域变化不一致性度量,如果最大的邻域变化不一致性度量对应的尺度为2或r_max，则像素点(i,j)不具有显著性区域，不再定义该像素点的显著性度量；否则，定义最大的邻域变化不一致性度量对应的尺度为像素点(i,j)的显著性尺度，并用r′_i,j表示；

4d)利用斑点噪声的伽玛分布,估计像素点(i,j)在显著性尺度r′_i,j邻域的灰度直方图

${\mathrm{hc}}_{i,j}^{r_{i,j}^{\text{\′}}}\left(y\right)=\underset{x\∈D}{\mathrm{\Σ}}\frac{x}{y^2}\frac{{\mathrm{ENL}}_{i,j}^{r_{i,j}^{\′}{\mathrm{ENL}}_{i,j}^{r_{i,j}^{\′}}}{\left(\frac{x}{y}\right)}^{{\mathrm{ENL}}_{i,j}^{r_{i,j}^{\′}}-1}}{\mathrm{\Γ}\left({\mathrm{ENL}}_{i,j}^{r_{i,j}^{\′}}\right)}\exp \left({-\mathrm{ENL}}_{i,j}^{r_{i,j}^{\′}}\left(\frac{x}{y}\right)\right)h_{i,j}^{r_{i,j}^{\′}}\left(x\right),$

其中，y包含从0到255的所有整数灰度值，x是灰度值，x∈D表示x在灰度值范围D内取值，D包含从0到255的所有整数灰度值，是像素点(i,j)在显著性尺度r′_i,j邻域的灰度直方图，Γ(·)是伽玛函数，是像素点(i,j)的显著性尺度r′_i,j邻域的等效视数；

4e)利用4d)中估计得到的显著性尺度r′_i,j邻域的灰度直方图计算像素点(i,j)的显著性尺度r′_i,j的局部显著性度量

4f)利用显著性尺度r′_i,j的局部显著性度量和邻域变化不一致性度量计算像素点(i,j)的显著性度量S_i,j，并将行向量(i,j,g_min+g×(r′_i,j-1),S_i,j)加入初始显著性矩阵Y_T；

(5)选取初始显著性矩阵Y_T中的前ε%个最大的显著性度量对应的行构建显著性矩阵Y_T′，0＜ε≤100，再通过迭代方法得到稳定显著性矩阵Y_S，提取稳定显著性矩阵Y_S每一行的前两个元素作为圆心的行列坐标以及第3个元素作为圆半径的像素数，在SAR图像中画出相应的圆形显著性区域。

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

(1)本发明由于采用受斑点噪声影响小的像素点的变化值及变化方向构建具有像素点变化幅度信息的二维直方图对，并在计算邻域变化不一致性度量时保持不同相邻尺度区域的面积比例不变，使得邻域变化不一致性度量更加有效的描述了邻域变化的不一致性，从而提高了区域显著性判断以及显著性尺度判断的准确性；

(2)本发明基于斑点噪声的伽玛分布估计显著性尺度邻域的灰度直方图，减少了斑点噪声对直方图统计的影响，提高了局部显著性度量的准确性。

仿真结果表明，本发明与现有的Kadir和Brady显著性区域检测方法相比，有效的消除了斑点噪声的影响，增加了显著性区域的检测准确性，并且提高有效显著性区域在检出的显著性区域中的比例。

附图说明

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

图2是本发明中8个3×3灰度变化值检测窗口与实测SAR图像进行相关操作后得到的像素点的变化值图和变化方向图；

图3是本发明应用于含有显著性目标的实测SAR图像的显著性区域检测结果图；

图4是本发明应用于两幅不同比例的城市实测SAR图像的显著性区域检测结果图以及显著性区域匹配图；

图5是本发明应用于两幅不同时刻的机场实测SAR图像的显著性区域检测结果图以及显著性区域匹配图。

具体实施方式

参照图1，本发明的具体实施步骤如下：

步骤1.输入SAR图像，计算图像的每个像素点在8个方向上灰度变化值。

输入一幅I×J大小的SAR图像，利用8个不同的3×3灰度变化值检测窗口，计算该图像的每个像素点(i,j)在8个方向上的灰度变化值，1≤i≤I，1≤j≤J：

v_i,j＝x_i+1,j-1+2×x_i+1,j+x_i+1,j+1-x_i-1,j-1-2×x_i-1,j-x_i-1,j+1

v′_i,j＝x_i-1,j-1+2×x_i-1,j+x_i-1,j+1-x_i+1,j-1-2×x_i+1,j-x_i+1,j+1

h_i,j＝x_i-1,j-1+2×x_i,j-1+x_i+1,j-1-x_i-1,j+1-2×x_i,j+1-x_i+1,j+1

h′_i,j＝x_i-1,j+1+2×x_i,j+1+x_i+1,j+1-x_i-1,j-1-2×x_i,j-1-x_i+1,j-1，

r_i,j＝x_i,j-1+2×x_i+1,j-1+x_i+1,j-x_i-1,j-2×x_i-1,j+1-x_i,j+1

r′_i，j＝x_i-1,j+2×x_i-1,j+1+x_i,j+1-x_i,j-1-2×x_i+1,j-1-x_i+1,j

l_i,j＝x_i,j+1+2×x_i+1,j+1+x_i+1,j-x_i-1,j-2×x_i-1,j-1-x_i,j-1

l′_i,j＝x_i-1,j+2×x_i-1,j-1+x_i,j-1-x_i,j+1-2×x_i+1,j+1-x_i+1,j

其中，v_i,j为像素点(i,j)的垂直向上灰度变化值，v′_i,j为像素点(i,j)的垂直向下灰度变化值，h_i,j为像素点(i,j)的水平向右灰度变化值，h′_i,j为像素点(i,j)的水平向左灰度变化值，r_i,j为像素点(i,j)的左下至右上灰度变化值，r′_i，j为像素点(i,j)的右上至左下灰度变化值，l_i,j为像素点(i,j)的右下至左上灰度变化值，l′_i,j为像素点(i,j)的左上至右下灰度变化值，x_i-1,j-1为像素点(i-1,j-1)的灰度值，x_i-1,j为像素点(i-1,j)的灰度值，x_i-1,j+1为像素点(i-1,j+1)的灰度值，x_i,j-1为像素点(i,j-1)的灰度值，x_i,j+1为像素点(i,j+1)的灰度值，x_i+1,j-1为像素点(i+1,j-1)的灰度值，x_i+1,j为像素点(i+1,j)的灰度值，x_i+1,j+1为像素点(i+1,j+1)的灰度值。

步骤2.将每个像素点在8个方向上灰度变化值的最大值，记录为该像素点的变化值，并将灰度变化值的最大值产生方向，记录为该像素点的变化方向，其中变化方向标记为1到8的整数；

步骤3.将坐标为第i行第j列的像素点(i,j)的变化值和变化方向分别表示为xc_i,j和xd_i,j，1≤i≤I，1≤j≤J。

步骤4.根据每个像素点(i,j)的变化值和变化方向，1≤i≤I，1≤j≤J，计算初始显著性矩阵Y_T。

4a)设定参数：最大尺度为r_max、尺度1邻域半径的像素点数为g_min、相邻尺度邻域半径增量的像素点数为g、像素点(i,j)的初始坐标为i＝1，j＝1；

4b)计算像素点(i,j)的每个尺度r的二维直方图对{ch_i,j(xc,xd,r),h_i,j(xc,xd,r)}，2≤r≤r_max：

${\mathrm{ch}}_{i,j}(\mathrm{xc},\mathrm{xd},r)=\underset{(m,n)\∈R_{i,j}^{r-1}}{\mathrm{\Σ}}{\mathrm{\δ}}_{{\mathrm{xc},\mathrm{xc}}_{m,n}}\×{\mathrm{\δ}}_{{\mathrm{xd},\mathrm{xd}}_{m,n}}\×\frac{\mathrm{xc}}{{\mathrm{xc}}_{\max }},$

$\begin{array}{c}{h}_{i,j}(\mathrm{xc},\mathrm{xd},r)={\mathrm{ch}}_{i,j}(\mathrm{xc},\mathrm{xd},r)+\underset{(m,n)\∈(R_{i,j}^{r-1},R_{i,j}^r]}{\mathrm{\Σ}}{\mathrm{\δ}}_{\mathrm{xc},{\mathrm{xc}}_{m,n}}\×{\mathrm{\δ}}_{{\mathrm{xd},\mathrm{xd}}_{m,n}}\\ \×\frac{(g_{\min }+g\×r-g)(2\×{g\×g}_{\min }+g^2)}{g_{\min }^2\×g\×(2\×g_{\min }+2\×g\×r-g)}\×\frac{\mathrm{xc}}{{\mathrm{xc}}_{\max }}\end{array},$

其中，ch_i,j(xc,xd,r)为像素点(i,j)在尺度r下对应的尺度r-1邻域内像素点的变化值和变化方向二维直方图，

h_i,j(xc,xd,r)为像素点(i,j)的尺度r邻域内像素点的变化值和变化方向二维直方图，

xc包含图像上像素点的变化值的所有取值，

xd包含像素点的8个变化方向标记，

xc_max为整个图像上最大的灰度变化值，

xc_m,n为像素点(m,n)的变化值，

xd_m,n为像素点(m,n)的变化方向，

表示像素点(m,n)包含于距离像素点(i,j)不超过g_min+g×(r-1)-g个像素点数的圆形邻域内，

表示像素点(m,n)包含于距离像素点(i,j)大于g_min+g×(r-1)-g个像素点数又不超过g_min+g×r-g个像素点数的环形邻域内，

为关于xc和xc_m,n的克罗内克函数， ${\mathrm{\δ}}_{\mathrm{xc},{\mathrm{xc}}_{m,n}}=\left\{\begin{array}{cc}1,& {\mathrm{xc}}_{m,n}=\mathrm{xc}\\ 0,& {\mathrm{xc}}_{m,n}\≠\mathrm{xc}\end{array}\right.,$

为关于xd和xd_m,n的克罗内克函数， ${\mathrm{\δ}}_{\mathrm{xd},{\mathrm{xd}}_{m,n}}=\left\{\begin{array}{cc}1,& {\mathrm{xd}}_{m,n}=\mathrm{xd}\\ 0,& {\mathrm{xd}}_{m,n}\≠\mathrm{xd}\end{array}\right.;$

4c)利用4b)的(r_max-1)个二维直方图对，计算像素点(i,j)的每个尺度r的邻域变化不一致性度量2≤r≤r_max：

$\begin{array}{c}{W}_{i,j}^r=\underset{\mathrm{xc}=0}{\overset{{\mathrm{xc}}_{\max }}{\mathrm{\Σ}}}\underset{\mathrm{xd}=1}{\overset{8}{\mathrm{\Σ}}}|{\mathrm{ch}}_{i,j}(\mathrm{xc},\mathrm{xd},r)/g_{\min }^2-h_{i,j}(\mathrm{xc},\mathrm{xd},r)/(g_{\min }^2+2g\×g_{\min }+g^2)|\\ \×g_{\min }^2/{(g_{\min }+g\×r-g)}^2\end{array};$

4d)从所有求得的(r_max-1)个邻域变化不一致性度量中找出最大的邻域变化不一致性度量,如果最大的邻域变化不一致性度量对应的尺度为2或r_max，则像素点(i,j)不具有显著性区域，不再定义该像素点的显著性度量，并根据像素点坐标跳转至相应步骤：

当i＝I,j＝J时，执行步骤4；

当j≠J时，令j＝j+1，并返回步骤4b)；

当i≠I,j＝J时，令i＝i+1,j＝1，并返回步骤4b)；

否则，定义最大的邻域变化不一致性度量对应的尺度为像素点(i,j)的显著性尺度，并用r′_i,j表示；

4e)利用斑点噪声的伽玛分布，估计像素点(i,j)在显著性尺度r′_i,j邻域的灰度直方图

${\mathrm{hc}}_{i,j}^{r_{i,j}^{\text{\′}}}\left(y\right)=\underset{x\∈D}{\mathrm{\Σ}}\frac{x}{y^2}\frac{{\mathrm{ENL}}_{i,j}^{r_{i,j}^{\′}{\mathrm{ENL}}_{i,j}^{r_{i,j}^{\′}}}{\left(\frac{x}{y}\right)}^{{\mathrm{ENL}}_{i,j}^{r_{i,j}^{\′}}-1}}{\mathrm{\Γ}\left({\mathrm{ENL}}_{i,j}^{r_{i,j}^{\′}}\right)}\exp \left({-\mathrm{ENL}}_{i,j}^{r_{i,j}^{\′}}\left(\frac{x}{y}\right)\right)h_{i,j}^{r_{i,j}^{\′}}\left(x\right),$

其中，y包含从0到255的所有整数灰度值，x是灰度值，x∈D表示x在灰度值范围D内取值，D包含从0到255的所有整数灰度值，是像素点(i,j)在显著性尺度r′_i,j邻域的灰度直方图，Γ(·)是伽玛函数，是像素点(i,j)的显著性尺度r′_i,j邻域的等效视数，按下式计算：

${\mathrm{ENL}}_{i,j}^{r_{i,j}^{\′}}=\frac{\underset{(m,n)\∈R_{i,j}^{r_{i,j}^{\′}}}{\mathrm{mean}}{\left(x_{m,n}\right)}^2}{\underset{(m,n)\∈R_{i,j}^{r_{i,j}^{\′}}}{\mathrm{var}}\left(x_{m,n}\right)},$

其中，mean(·)是求均值，var(·)是求方差，表示像素点(m,n)包含于距离像素点(i,j)不超过g_min+g×(r′_i,j-1)-g个像素点数的圆形邻域内。

4f)利用4e)中估计得到的显著性尺度r′_i,j邻域的灰度直方图计算像素点(i,j)的显著性尺度r′_i,j的局部显著性度量

$H_{i,j}^{r_{i,j}^{\′}}=-\underset{y\∈D}{\mathrm{\Σ}}{\mathrm{hc}}_{i,j}^{r_{i,j}^{\′}}\left(y\right)\log \left({\mathrm{hc}}_{i,j}^{r_{i,j}^{\′}}\left(y\right)\right),$

其中，y是灰度值，y∈D表示y在灰度值范围D内取值，D包含从0到255的所有整数灰度值；

4g)利用显著性尺度r′_i,j的局部显著性度量和邻域变化不一致性度量计算像素点(i,j)的显著性度量S_i,j：

$S_{i,j}=H_{i,j}^{r_{i,j}^{\′}}\×W_{i,j}^{r_{i,j}^{\′}},$

并将行向量(i,j,g_min+g×(r′_i,j-1),S_i,j)加入初始显著性矩阵Y_T，根据像素点坐标跳转至相应步骤：

当i＝I,j＝J时，执行步骤4；

当j≠J时，令j＝j+1，并返回步骤4b)；

当i≠I,j＝J时，令i＝i+1,j＝1，并返回步骤4b)；

步骤5.选取初始显著性矩阵Y_T中的前ε%个最大的显著性度量对应的行构建显著性矩阵Y_T′，0＜ε≤100，根据显著性矩阵Y_T′按如下步骤得到稳定显著性矩阵Y_S：

(5a)设置显著性尺度差阈值sca_max，显著性度量差阈值sal_max，相邻像素点数量阈值n_b，并设稳定显著性矩阵Y_S为空矩阵；

(5b)选择显著性矩阵Y_T′中显著性度量最大的像素点作为候选点，并计算显著性矩阵Y_T′中所有其它像素点坐标到候选点坐标的欧式距离；

(5c)选择欧式距离小于候选点的显著性尺度邻域半径并与候选点的显著性尺度差绝对值小于sca_max和显著性度量差绝对值小于sal_max的像素点，作为候选点的相邻像素点；

(5d)如果相邻像素点数小于n_b，则将候选点所在行从显著性矩阵Y_T′中去掉；否则将候选点所在行加入稳定显著性矩阵Y_S，再将候选点及其所有相邻像素点所在行从显著性矩阵Y_T′中去掉；

(5e)判断显著性矩阵Y_T′是否为空，若为空则停止，并输出稳定显著性矩阵Y_S；否则，回到步骤(5b)。

步骤6.提取稳定显著性矩阵Y_S中每一行的前两个元素作为圆心的行列坐标以及第3个元素作为圆半径的像素数，在SAR图像中画出相应的圆形显著性区域。

本发明效果可以通过以下仿真进一步说明：

1.实验条件和内容

实验仿真环境为：MATLABR2011b,Intel(R)Corei5-3470CPU3.2GHz，Window7专业版。

2.实验内容及结果：

实验1，将8个3×3灰度变化值检测窗口作用于实测SAR图像，根据灰度变化值的最大值得到像素点的变化值图和变化方向图，如图2所示，其中图2(a)是实测SAR图像，图2(b)是变化值图，图2(c)是变化方向图。

由图2(b)可以看出，3×3灰度变化值检测窗口具有一定的降噪作用，明亮的大灰度变化区域和较暗的匀质灰度变化区域受到斑点噪声的影响均不明显，此外图2(b)还提供了像素点的灰度变化幅度情况。图2(c)的大灰度变化区域的变化方向基本没有受到噪声影响，但匀质灰度变化区域的变化方向由于斑点噪声的影响具有随机性，这并不影响邻域变化不一致性度量。

实验2，分别用本发明方法和现有的KadirBrady显著性区域检测算法对含有显著性目标的实测SAR图像进行显著性区域检测，结果如图3，其中图3(a)是含有显著性目标的实测SAR图像，图3(b)是用KadirBrady显著性区域检测算法对图3(a)的显著性区域检测结果，图3(c)是本发明方法对图3(a)的显著性区域检测结果。

由图3(b)可以看出，Kadir和Brady显著性区域检测算法漏检漏左边第二个显著性目标，而在右下方的匀质变化区域误检为显著性区域。此外，该方法检测的显著性尺度大于显著性目标的实际尺度。由图3(c)可以看出，本方法正确的检测出来所有显著性目标，并且显著性尺度与显著性目标的实际大小相符。

实验3，分别用本发明方法和现有的KadirBrady显著性区域检测算法对两幅不同比例的城市实测SAR图像进行显著性区域检测，结果如图4，其中图4(a)是大比例的城市实测SAR图像，图4(b)是小比例的城市实测SAR图像，图4(c)是KadirBrady显著性区域检测算法对图4(a)的显著性区域检测结果，图4(d)是KadirBrady显著性区域检测算法对图4(b)的显著性区域检测结果，图4(e)是本发明方法对图4(a)的显著性区域检测结果，图4(f)是本发明方法对图4(b)的显著性区域检测结果，图4(g)是图4(c)中与图4(d)匹配的显著性区域，图4(h)是图4(d)中与图4(c)匹配的显著性区域，图4(i)是图4(e)中与图4(f)匹配的显著性区域，图4(j)是图4(f)中与图4(e)匹配的显著性区域。

由图4(c)可以看出，部分匀质灰度变化区域被KadirBrady显著性区域检测算法误检成显著性区域，并且其中的显著性区域的显著性尺度偏大；由图4(d)可以看出，虽然图像大小缩小，但KadirBrady显著性区域检测算法检出的显著性区域的显著性尺度与图像4(c)中的显著性区域的显著性尺度相比变化不大；由图4(g)和图4(h)可以看出，在图像大小变化后，KadirBrady显著性区域检测算法检出的能够匹配的显著性区域少；由图4(e)可以看出，本发明方法能够检出大量的显著性区域；由图4(f)可以看出，本发明方法检出的显著性区域的显著性尺度随图像大小的缩小而缩小；由图4(i)和图4(j)可以看出，本发明方法检出显著性区域的匹配对数多。

实验4，分别用本发明方法和现有的KadirBrady显著性区域检测算法对两幅不同时刻的机场实测SAR图像进行显著性区域检测，结果如图5，其中图5(a)是时刻一的机场实测SAR图像，图5(b)是时刻二的机场实测SAR图像，图5(c)是KadirBrady显著性区域检测算法对图5(a)的显著性区域检测结果，图5(d)是KadirBrady显著性区域检测算法对图5(b)的显著性区域检测结果，图5(e)是本发明方法对图5(a)的显著性区域检测结果，图5(f)是本发明方法对图5(b)的显著性区域检测结果，图5(g)是图5(c)中与图5(d)匹配的显著性区域，图5(h)是图5(d)中与图5(c)匹配的显著性区域，图5(i)是图5(e)中与图5(f)匹配的显著性区域，图5(j)是图5(f)中与图5(e)匹配的显著性区域。

由图5(c)可以看出，匀质灰度变化区域被KadirBrady显著性区域检测算法误检成显著性尺度较大的显著性区域；从图5(d)可以看出，在图像视角变化后，KadirBrady显著性区域检测算法在匀质灰度变化区域依然出现大量误检；由图5(g)和图5(h)可以看出，在图像视角变化后，KadirBrady显著性区域检测算法检出的能够匹配的显著性区域少；由图5(e)可以看出，本发明方法能够在灰度变化大的区域检出大量的显著性区域；由图5(f)可以看出，本发明方法检出的显著性区域位置和数量不会随图像视角的变化而变化；由图5(i)和图5(j)可以看出，本发明方法检出显著性区域的匹配对数多。

表1本发明方法与KadirBrady显著性区域检测算法的比较结果

为验证本发明方法和KadirBrady显著性区域检测算法的稳定性和有效性，本发明选取显著性区域检测数量、显著性区域匹配数量以及匹配率，对上述实验的结果进行统计，结果如表1所示。

从表1可见，KadirBrady显著性区域检测算法检出的显著性区域的数量少，并由于较低的匹配率使得显著性区域匹配数量少，所以在图像发生变化后KadirBrady显著性区域检测算法的稳定性和有效性低；相反，本发明方法检出的显著性区域匹配率高，由于检出的显著性区域数量大使得匹配的显著性区域多，可见本发明方法具有较高的稳定性以及较好的有效性。

Claims (3)

1.一种基于伽玛分布和邻域信息的SAR图像显著性区域检测方法，包括以下步骤：

(1)输入一幅I×J大小的SAR图像，利用8个不同的3×3灰度变化值检测窗口，计算该图像的每个像素点在8个方向上的灰度变化值；

(2)将每个像素点在8个方向上灰度变化值的最大值，记录为该像素点的变化值，并将灰度变化值的最大值产生方向，记录为该像素点的变化方向，其中变化方向标记为1到8的整数；

(3)将坐标为第i行第j列的像素点(i,j)的变化值和变化方向分别表示为xc_i,j和xd_i,j，1≤i≤I，1≤j≤J；

(4)根据每个像素点(i,j)的变化值和变化方向，1≤i≤I，1≤j≤J，计算初始显著性矩阵Y_T：

4a)计算像素点(i,j)的每个尺度r的二维直方图对{ch_i,j(xc,xd,r),h_i,j(xc,xd,r)}，2≤r≤r_max：

${\mathrm{ch}}_{i,j}(xc,xd,r)=\underset{(m,n)\∈R_{i,j}^{r-1}}{\Σ}{\mathrm{\δ}}_{xc,{\mathrm{xc}}_{m,n}}\×{\mathrm{\δ}}_{xd,{\mathrm{xd}}_{m,n}}\×\frac{xc}{{\mathrm{xc}}_{max}},$

$\begin{array}{c}{h}_{i,j}\left(xc,xd,r\right)={\mathrm{ch}}_{i,j}\left(xc,xd,r\right)+\underset{\left(m,n\right)\∈\left(R_{i,j}^{r-1},R_{i,j}^r\right)}{\Σ}{\mathrm{\δ}}_{xc,{\mathrm{xc}}_{m,n}}\×{\mathrm{\δ}}_{xd,{\mathrm{xd}}_{m,n}}\\ \×\frac{\left(g_{\min }+g\×r-g\right)\left(2\×g\×g_{\min }+g^2\right)}{g_{\min }^2\×g\×\left(2\×g_{\min }+2\×g\×r-g\right)}\×\frac{xc}{{\mathrm{xc}}_{\max }}\end{array},$

其中，ch_i,j(xc,xd,r)为像素点(i,j)在尺度r下对应的尺度r-1邻域内像素点的变化值和变化方向二维直方图，h_i,j(xc,xd,r)为像素点(i,j)的尺度r邻域内像素点的变化值和变化方向二维直方图，xc包含图像上像素点的变化值的所有取值，xd包含像素点的8个变化方向标记，xc_max为整个图像上最大的灰度变化值，xc_m,n为像素点(m,n)的变化值，xd_m,n为像素点(m,n)的变化方向，表示像素点(m,n)包含于距离像素点(i,j)不超过g_min+g×(r-1)-g个像素点数的圆形邻域内，表示像素点(m,n)包含于距离像素点(i,j)大于g_min+g×(r-1)-g个像素点数又不超过g_min+g×r-g个像素点数的环形邻域内，g_min是尺度1邻域半径的像素点数，g是相邻尺度邻域半径增量的像素点数，为关于xc和xc_m,n的克罗内克函数， ${\mathrm{\δ}}_{xc,{\mathrm{xc}}_{m,n}}=\left\{\begin{array}{cc}1,& {\mathrm{xc}}_{m,n}=xc\\ 0,& {\mathrm{xc}}_{m,n}\≠xc\end{array}\right.,$ 为关于xd和xd_m,n的克罗内克函数， ${\mathrm{\δ}}_{xd,{\mathrm{xd}}_{m,n}}=\left\{\begin{array}{cc}1,& {\mathrm{xd}}_{m,n}=xd\\ 0,& {\mathrm{xd}}_{m,n}\≠xd\end{array}\right.;$

4b)利用4a)的(r_max-1)个二维直方图对，计算像素点(i,j)的每个尺度r的邻域变化不一致性度量2≤r≤r_max：

$\begin{array}{c}{W}_{i,j}^r=\underset{xc=0}{\overset{{\mathrm{xc}}_{\max }}{\Σ}}\underset{xd=1}{\overset{8}{\Σ}}|{\mathrm{ch}}_{i,j}(xc,xd,r)/g_{min}^2-h_{i,j}(xc,xd,r)/(g_{min}^2+2g\×g_{\min }+g^2)|\\ \×g_{\min }^2/{\left(g_{\min }+g\×r-g\right)}^2\end{array};$

4c)从所有求得的(r_max-1)个邻域变化不一致性度量中找出最大的邻域变化不一致性度量,如果最大的邻域变化不一致性度量对应的尺度为2或r_max，则像素点(i,j)不具有显著性区域，不再定义该像素点的显著性度量；否则，定义最大的邻域变化不一致性度量对应的尺度为像素点(i,j)的显著性尺度，并用r′_i,j表示；

4d)利用斑点噪声的伽玛分布,估计像素点(i,j)在显著性尺度r′_i,j邻域的灰度直方图

${\mathrm{hc}}_{i,j}^{r_{i,j}^{\′}}\left(y\right)=\underset{x\∈D}{\Σ}\frac{x}{y^2}\frac{{\mathrm{ENL}}_{i,j}^{r_{i,j}^{\′}{\mathrm{ENL}}_{i,j}^{r_{i,j}^{\′}}}{\left(\frac{x}{y}\right)}^{{\mathrm{ENL}}_{i,j}^{r_{i,j}^{\′}}-1}}{\mathrm{\Γ}\left({\mathrm{ENL}}_{i,j}^{r_{i,j}^{\′}}\right)}\exp \left(-{\mathrm{ENL}}_{i,j}^{r_{i,j}^{\′}}\left(\frac{x}{y}\right)\right)h_{i,j}^{r_{i,j}^{\′}}\left(x\right),$

其中，y包含从0到255的所有整数灰度值，x是灰度值，x∈D表示x在灰度值范围D内取值，D包含从0到255的所有整数灰度值，是像素点(i,j)在显著性尺度r′_i,j邻域的灰度直方图，Γ(·)是伽玛函数，是像素点(i,j)的显著性尺度r′_i,j邻域的等效视数，

${\mathrm{ENL}}_{i,j}^{r_{i,j}^{\′}}=\frac{\underset{\left(m,n\right)\∈R_{i,j}^{r_{i,j}^{\′}}}{mean}{\left(x_{m,n}\right)}^2}{\underset{\left(m,n\right)\∈R_{i,j}^{r_{i,j}^{\′}}}{\mathrm{var}}\left(x_{m,n}\right)},$

其中，mean(·)是求均值，var(·)是求方差，表示像素点(m,n)包含于距离像素点(i,j)不超过g_min+g×(r′_i,j-1)-g个像素点数的圆形邻域内；

4e)利用4d)中估计得到的显著性尺度r′_i,j邻域的灰度直方图计算像素点(i,j)的显著性尺度r′_i,j的局部显著性度量

$H_{i,j}^{r_{i,j}^{\′}}=-\underset{y\∈D}{\Σ}{\mathrm{hc}}_{i,j}^{r_{i,j}^{\′}}\left(y\right)log\left({\mathrm{hc}}_{i,j}^{r_{i,j}^{\′}}\right(y\left)\right),$

其中，y是灰度值，y∈D表示y在灰度值范围D内取值，D包含从0到255的所有整数灰度值；

4f)利用显著性尺度r′_i,j的局部显著性度量和邻域变化不一致性度量计算像素点(i,j)的显著性度量：并将行向量(i,j,g_min+g×(r′_i,j-1),S_i,j)加入初始显著性矩阵Y_T；

(5)选取初始显著性矩阵Y_T中的前ε％个最大的显著性度量对应的行构建显著性矩阵Y_T′，0<ε≤100，再通过迭代方法得到稳定显著性矩阵Y_S，提取稳定显著性矩阵Y_S每一行的前两个元素作为圆心的行列坐标以及第3个元素作为圆半径的像素数，在SAR图像中画出相应的圆形显著性区域。

2.根据权利要求1所述方法，其中所述步骤1)每个像素点在8个方向上的灰度变化值，按下式计算：

v_i,j＝x_i+1,j-1+2×x_i+1,j+x_i+1,j+1-x_i-1,j-1-2×x_i-1,j-x_i-1,j+1

v′_i,j＝x_i-1,j-1+2×x_i-1,j+x_i-1,j+1-x_i+1,j-1-2×x_i+1,j-x_i+1,j+1

h_i,j＝x_i-1,j-1+2×x_i,j-1+x_i+1,j-1-x_i-1,j+1-2×x_i,j+1-x_i+1,j+1

h′_i,j＝x_i-1,j+1+2×x_i,j+1+x_i+1,j+1-x_i-1,j-1-2×x_i,j-1-x_i+1,j-1

，

r_i,j＝x_i,j-1+2×x_i+1,j-1+x_i+1,j-x_i-1,j-2×x_i-1,j+1-x_i,j+1

r′_i,j＝x_i-1,j+2×x_i-1,j+1+x_i,j+1-x_i,j-1-2×x_i+1,j-1-x_i+1,j

l_i,j＝x_i,j+1+2×x_i+1,j+1+x_i+1,j-x_i-1,j-2×x_i-1,j-1-x_i,j-1

l′_i,j＝x_i-1,j+2×x_i-1,j-1+x_i,j-1-x_i,j+1-2×x_i+1,j+1-x_i+1,j

其中，i为像素点所在行，1≤i≤I，j为像素点所在列，1≤j≤J，v_i,j为像素点(i,j)的垂直向上灰度变化值，v′_i,j为像素点(i,j)的垂直向下灰度变化值，h_i,j为像素点(i,j)的水平向右灰度变化值，h′_i,j为像素点(i,j)的水平向左灰度变化值，r_i,j为像素点(i,j)的左下至右上灰度变化值，r′_i,j为像素点(i,j)的右上至左下灰度变化值，l_i,j为像素点(i,j)的右下至左上灰度变化值，l′_i,j为像素点(i,j)的左上至右下灰度变化值，x_i-1,j-1为像素点(i-1,j-1)的灰度值，x_i-1,j为像素点(i-1,j)的灰度值，x_i-1,j+1为像素点(i-1,j+1)的灰度值，x_i,j-1为像素点(i,j-1)的灰度值，x_i,j+1为像素点(i,j+1)的灰度值，x_i+1,j-1为像素点(i+1,j-1)的灰度值，x_i+1,j为像素点(i+1,j)的灰度值，x_i+1,j+1为像素点(i+1,j+1)的灰度值。

3.根据权利要求1所述方法，其中所述步骤5)通过迭代方法得到稳定显著性矩阵Y_S，按如下步骤得到：

(5a)设置显著性尺度差阈值sca_max，显著性度量差阈值sal_max，相邻像素点数量阈值n_b，并设稳定显著性矩阵Y_S为空矩阵；

(5b)选择显著性矩阵Y_T′中显著性度量最大的像素点作为候选点，并计算显著性矩阵Y_T′中所有其它像素点坐标到候选点坐标的欧式距离；

(5c)选择欧式距离小于候选点的显著性尺度邻域半径并与候选点的显著性尺度差绝对值小于sca_max和显著性度量差绝对值小于sal_max的像素点，作为候选点的相邻像素点；

(5d)如果相邻像素点数小于n_b，则将候选点所在行从显著性矩阵Y_T′中去掉；否则将候选点所在行加入稳定显著性矩阵Y_S，再将候选点及其所有相邻像素点所在行从显著性矩阵Y_T′中去掉；

(5e)判断显著性矩阵Y_T′是否为空，若为空则停止，并输出稳定显著性矩阵Y_S；否则，回到步骤(5b)。