CN111028256A - 一种基于改进活动轮廓模型的sar图像目标检测方法 - Google Patents

Abstract

本发明公开了一种提出了一种基于改进活动轮廓模型的SAR图像目标检测方法，该方法基于Gamma分布假设，利用Bhattacharyya距离度量推导了基于局部区域统计信息的活动轮廓模型的能量泛函，利用全局最小化框架，实现所提能量泛函的全局能量最小化。本发明方法能够获得精确定位的目标轮廓，从而达到检测目标的目的，并且具有很高的计算效率。

Description

一种基于改进活动轮廓模型的SAR图像目标检测方法

技术领域

本发明属于雷达目标检测技术领域，特别涉及一种在合成孔径雷达SAR图像中快速有效分割目标的方法。

背景技术

图像分割技术在合成孔径雷达SAR图像的自动分析和解译中起着重要作用。作为一种典型的图像分割技术，活动轮廓模型ACM的思想在于：将用于分割目标轮廓的曲线视为可变形的，利用能量泛函最小化引导轮廓线不断向目标边界演化。由于ACM具有对噪声的鲁棒性，并且具有对边界和区域自然表示的优点，得到的分割结果无需进行额外的区域合并或边界连接等后续的处理，目前ACM已广泛应用于光学图像和医学影像中的目标边界定位。然而，由于SAR图像中存在乘性相干斑噪声的影响，导致背景与目标之间的边界变得模糊，并且导致均匀区域的像素强度发生波动。

针对SAR图像分割问题提出的ACM模型，按照可利用的图像特征，可大致分为以下几类：基于像素相似性、基于区域统计信息以及基于能量泛函凸优化的活动轮廓模型等。其中，区域尺度可调RSF活动轮廓模型利用轮廓内部区域和轮廓外部区域的像素强度信息，具有一定的抗噪性，但是由于能量泛函的非凸性，RSF模型的分割结果对于水平集函数的初始化很敏感，并且在轮廓演化的过程中，其能量泛函容易陷入局部极小值。为克服局部极小值的缺点，将RSF的能量泛函引入全局最小化活动轮廓框架中，定义了凸的RSF模型，记为全局最小化区域尺度可调GRSF模型，但是RSF模型和GRSF模型存在相似的缺陷：两者的能量泛函中的区域拟合函数都是根据欧几里得距离计算的，虽然该距离可以用于衡量像素点与轮廓内外区域之间的相似程度，但是由于SAR图像中乘性相干斑噪声的存在，欧氏距离不适用于SAR图像，无法准确反映区域能量的变化趋势。

发明内容

针对上述已有技术的不足，本发明提出了一种基于改进活动轮廓模型的SAR图像目标检测方法，该方法基于Gamma分布假设，推导了基于局部区域统计信息的活动轮廓模型的能量泛函，利用全局最小化框架，实现所提能量泛函的全局能量最小化，能够获得精确定位的目标轮廓，从而达到检测目标的目的。

为达到上述目的，本发明采用如下技术方案予以实现。

一种基于改进活动轮廓模型的SAR图像目标检测方法，包括如下步骤：

步骤1，输入待检测的SAR图像，设其像素强度为I(x)，图像尺寸为H×W像素；

步骤2，根据待检测的SAR图像，计算边缘检测算子g(x)；

步骤3，对各项参数进行初始化处理：设置最大迭代次数为N，初始迭代次数为n＝0，迭代终止条件ε，并根据待检测SAR图像分别设置时间步长τ和参数θ，将待检测的SAR图像I(x)利用最大值进行归一化作为初始的水平集函数，即

设置初始的正则化辅助变量v⁽⁰⁾＝u⁽⁰⁾，初始的对偶变量p⁽⁰⁾＝(p1⁽⁰⁾,p2⁽⁰⁾)＝0，其中p1和p2分别为沿着列方向和沿着行方向的梯度值。

步骤4，根据能量泛函最小化准则进行目标轮廓曲线演化，按照以下子步骤迭代更新得到目标轮廓曲线u：

4a)更新u函数:

4a1)计算div(p⁽ⁿ⁾)，其中div(·)表示散度运算；

4a2)根据计算得到的div(p⁽ⁿ⁾)，计算p⁽ⁿ⁺¹⁾；

4a3)根据计算得到的p⁽ⁿ⁺¹⁾，计算u⁽ⁿ⁺¹⁾＝v⁽ⁿ⁾-θdiv(p⁽ⁿ⁺¹⁾)_；

4b)更新v函数:

4b1)分别计算轮廓内外区域的能量值e₁和e₂；

4b2)根据u⁽ⁿ⁺¹⁾、e₁和e₂，计算v⁽ⁿ⁺¹⁾＝min{max{u⁽ⁿ⁺¹⁾-θ(λ₁e₁-λ₂e₂),0},1}；

4c)更新迭代次数n＝n+1；

4d)判断是否满足终止条件||u⁽ⁿ⁺¹⁾-u⁽ⁿ⁾||₂≤ε(||·||₂表示取模值运算)或者n＞N，若满足，则终止迭代，输出目标轮廓曲线演化结果u；否则，继续执行步骤4a)～步骤4d)之间的操作。

步骤5，根据目标轮廓曲线演化结果u，得到目标检测结果。

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

1.本发明方法采用了GMAC框架，将归一化的图像作为初始的水平集函数，因此最终的目标分割结果的精度不受轮廓初始位置的影响。

2.本发明方法对相干斑噪声具有较强的鲁棒性，能够有效地缓解轮廓演化过程中的局部极小值问题。

附图说明

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

图2(a)是SAR图像处理方法对地面目标的轮廓演化结果；

图2(b)是RSF模型处理方法对地面目标的轮廓演化结果；

图2(c)是GRSF模型处理方法对地面目标的轮廓演化结果；

图2(d)本发明对地面目标的轮廓演化结果；

图3是RSF方法对海面目标的轮廓演化结果；

图4(a)是GRSF方法对海面目标的轮廓演化结果；

图4(b)是GRSF方法对海面目标的目标检测结果；

图5(a)是本发明方法对海面目标的轮廓演化结果

图5(b)是本发明方法对海面目标的目标检测结果。

具体实施方式

下面结合附图对本发明的实施步骤和效果作进一步的说明。

本发明基于改进活动轮廓模型的SAR图像目标检测，提出了一种基于局部统计信息的改进活动轮廓模型，即

式中，λ₁，λ₂，μ＞0和v＞0均为权值系数，e₁和e₂分别为轮廓曲线Ω内外区域的能量函数，第三项

表示轮廓长度项；第四项

为正则项_，用于调整水平集函数以获得稳定的水平集演化结果_。H(φ(x))为Heaviside函数_，可表示为

其中φ(x)为满足以下条件的水平集函数：

根据变分原理

对能量泛函E_modifiedRSF进行全变分，得到如下的水平集函数演化方程：

设n为迭代次数，τ为时间步长，根据梯度下降流算法，则水平集函数的迭代过程可由下式表示：

根据上式中的迭代方程，可对水平集函数迭代求解，完成SAR图像分割。为了消除最终的分割结果对初始轮廓位置的严重依赖性，采用快速全局最小化活动模型框架，限制水平集函数的取值范围为0≤φ(x)≤1，在下文中用新变量u进行代替φ，利用如下的方程式实现能量最小化：

其中，||·||₂表示取模值运算，θ＞0是一个很小的值，v(x)表示正则化辅助变量，g(x)表示边缘检测算子，它是关于图像梯度的单调递减函数。

参照图1，对于待检测的SAR图像，设其像素强度为I(x)，图像尺寸为H×W像素，需要按照如下步骤实现能量泛函E(u,v,e₁,e₂)的总体能量最小化：

步骤1，利用待检测的SAR图像I(x)，按照如下公式计算边缘检测算子g(x)：

其中，r为一个极小值，实验中设定r＝10^-5，

表示二维梯度运算。

步骤2，对各项参数进行初始化处理：利用待检测的SAR图像I(x)设置初始的水平集函数为

设置初始的正则化辅助变量v⁽⁰⁾＝u⁽⁰⁾，初始的对偶变量p⁽⁰⁾＝(p1⁽⁰⁾,p2⁽⁰⁾)＝0，其中p1和p2分别为沿着列方向和沿着行方向的梯度值，设置最大迭代次数为N，初始迭代次数为n＝0，迭代终止条件ε，并根据测试SAR图像分别设置时间步长τ和参数θ。

步骤3，根据能量泛函最小化准则进行目标轮廓曲线演化，按照以下子步骤迭代更新得到最终的水平集函数u，即目标轮廓曲线u：

3a)更新水平集函数u:

3a1)利用如下公式计算div(p⁽ⁿ⁾)：

式中，p＝(p1,p2)是满足下列等式的对偶变量：

3a2)将div(p⁽ⁿ⁾)代入如下公式，计算第n+1次迭代过程中对偶变量p的值p⁽ⁿ⁺¹⁾；

其中，v⁽ⁿ⁾表示第n次迭代过程中的正则化辅助变量；

3a3)根据计算得到的p⁽ⁿ⁺¹⁾，计算第n+1次迭代过程中u函数的值，即

u⁽ⁿ⁺¹⁾＝v⁽ⁿ⁾-θdiv(p⁽ⁿ⁺¹⁾)；

3b)更新正则化辅助变量v:

3b1)基于Gamma分布假设和Bhattacharyya距离度量，分别计算轮廓内外区域的能量值e₁和e₂；

具体地，两个Gamma分布之间的Bhattacharyya距离度量的推导过程如下：

假设轮廓曲线内外以像素点x为中心的两个局部区域内的像素强度分别满足如下的Gamma分布函数：

其中，α_i表示形状参数，β_i表示逆尺度参数。根据Bhattacharyya距离度量，比较轮廓内外的局部区域Gamma分布之间的统计特性差异。对于Gamma分布函数p₁和p₂，两者之间的Bhattacharyya距离具体可用下面的公式表示：

其中，

ln(·)表示对数运算。根据推导得到的Bhattacharyya距离B(p₁,p₂)，利用下式分别计算轮廓内外区域的能量值e₁和e₂：

其中，Γ(·)表示标准的伽马函数，(α_x,β_x)表示以像素点x为中心的局部区域内像素服从的Gamma分布参数对，(α_in,β_in)和(α_out,β_out)分别表示轮廓内部局部区域和轮廓外部局部区域对应的Gamma分布参数对，

具体地，(α_x,β_x)可由下式计算得到：

α_x＝μ²/σ²，β_x＝α_x/μ

其中，μ和σ²分别表示局部区域内的像素强度均值和方差。(α_in,β_in)和(α_out,β_out)分别可由下式计算得到：

其中，

和

分别表示轮廓内部局部区域的像素强度均值和方差，

和

分别表示轮廓外部局部区域的像素强度均值和方差，式中H(φ)和H(-φ)分别表示轮廓内部和轮廓外部的水平集函数取值。

3b2)根据计算得到的u⁽ⁿ⁺¹⁾、e₁和e₂，更新第n+1次迭代过程中v函数的值v⁽ⁿ⁺¹⁾：v⁽ⁿ⁺¹⁾＝min{max{u⁽ⁿ⁺¹⁾-θ(λ₁e₁-λ₂e₂),0},1}

3c)更新迭代次数n＝n+1；

3d)判断是否满足终止条件||u⁽ⁿ⁺¹⁾-u⁽ⁿ⁾||₂≤ε或者n＞N，若满足，则终止迭代，输出目标轮廓曲线演化结果u；否则，继续执行步骤3a)～步骤3d)之间的操作。

步骤4，根据目标轮廓曲线演化得到的水平集函数u，对于任一像素点x，x∈x，若u(x)≥0.5，则将其视为目标像素点；否则为杂波像素点。在对所有像素点进行判别之后，得到目标检测结果。

本发明的效果可以通过以下仿真实验具体说明：

1.仿真实验条件

仿真实验运行平台为MATLAB R2014a，Intel(R)Core(TM)i7-4790 CPU@3.6GHz，内存8GB。

2.仿真实验内容

仿真1：如图2(a)所示，对T-72主战坦克、BTR-70装甲运兵车、BMP-2步兵战车三种不同类型目标的SAR图像进行处理，图像分辨率为0.3m×0.3m，每幅SAR图像的尺寸均为128×128像素，场景中包含了目标及其阴影区域和背景杂波。为验证本发明方法的有效性，设置最大迭代次数N＝100，迭代终止条件ε＝10-⁴，时间步长_τ＝0.05和参数θ＝1，采用RSF模型和GRSF模型作为对比。对比方法的目标分割结果如图2所示。从图2(b)的目标轮廓曲线演化结果可以看出，RSF模型能够准确定位目标边界，但其轮廓曲线在演化的过程中卡在了背景区域里，不能实现进一步的演化，因此无法获得准确的目标分割结果。如图2(c)所示，由于GRSF模型对相干斑噪声的敏感性，该模型不能准确定位目标边界，并且存在局部极小值。本发明方法的轮廓曲线演化结果如图2(d)所示，可以看出，对于三种不同类型的地面目标，本发明方法分割所得的目标边界清晰，并且均能很好地贴合真实目标边界。针对这组SAR图像，本发明方法耗时约为0.2秒，与GRSF模型的耗时相当，而RSF模型平均耗时约为17秒。

仿真2：采用本发明方法对海面舰船目标SAR图像进行处理，距离和方位分辨率均为3m，SAR图像尺寸为1000×800像素，场景中存在四艘舰船，舰船目标的附近存在细长形状的防波堤。图3给出了RSF模型的轮廓曲线演化结果，可以看出尽管该模型能够定位目标边界，但背景中存在无法继续演化的曲线，因此不能得到对应的目标检测结果。GRSF模型获得的目标轮廓曲线如图4(a)所示，可以看出GRSF模型虽然能够检测到目标，但由于相干斑噪声敏感性，目标与杂波之间的边界稍微有些模糊，GRSF模型的目标检测结果如图4(b)所示。本发明方法获得的目标轮廓曲线和目标检测结果分别如图5(a)和5(b)所示，从直观的视觉效果来看，所提方法的目标边界定位准确性优于两种对比模型。

为定量地评估目标检测性能，定义如下的目标检测率和虚警率：

目标检测率DR＝N_dt/N_tt，其中，N_dt为检测到的目标像素点数，N_tt为目标像素总数；虚警率FAR＝N_dc/N_tc，其中，N_dc为杂波像素中被误判为目标像素的虚警数，N_tc为杂波像素总数。

所有对比方法的定量分析结果如表1所示，可以看出，本发明方法具有最佳的检测性能，并且耗时最少。

表1对比方法的目标检测性能分析

对比方法	DR(％)	FAR(％)	耗时(s)
				本发明方法	90.72	0.02	6.8
GRSF模型	67.27	0.05	16.5
				双参数CFAR	62.29	0.21	63.3

Claims (3)

1.一种基于改进活动轮廓模型的SAR图像目标检测方法，其特征在于，包括如下步骤：

步骤1，输入待检测的SAR图像，设其像素强度为I(x)，图像尺寸为H×W像素；

步骤2，利用待检测的SAR图像，按照如下公式计算边缘检测算子g(x)：

其中，r为一个极小值，

表示二维梯度运算；

步骤3，对各项参数进行初始化处理：设置最大迭代次数为N，初始迭代次数为n＝0，迭代终止条件ε，并根据待检测SAR图像分别设置时间步长τ和参数θ，将待检测的SAR图像I(x)利用最大值进行归一化，并作为初始的水平集函数u⁽⁰⁾，即

设置初始的正则化辅助变量v⁽⁰⁾＝u⁽⁰⁾，初始的对偶变量p⁽⁰⁾＝(p1⁽⁰⁾,p2⁽⁰⁾)＝0，其中p1和p2分别为沿着列方向和沿着行方向的梯度值；

步骤4，根据能量泛函最小化准则进行目标轮廓曲线演化，迭代更新水平集函数u和正则化辅助变量v，得到最终的目标轮廓曲线u；

步骤5，根据目标轮廓曲线演化结果u，得到目标检测结果。

2.根据权利要求1所述的方法，其中步骤4根据能量泛函最小化准则进行目标轮廓曲线演化具体按照以下方式进行：

4a)更新水平集函数u：

4a1)利用如下公式计算div(p⁽ⁿ⁾)：

式中，div(·)表示散度运算，其中p⁽ⁿ⁾为迭代n次后的对偶变量，i和j分别表示行列索引，p＝(p1,p2)是满足下列等式的对偶变量：

4a2)利用div(p⁽ⁿ⁾)计算第n+1次迭代过程中对偶变量p的值p⁽ⁿ⁺¹⁾：

其中v⁽ⁿ⁾表示第n次迭代过程中的正则化辅助变量；

4a4)根据计算得到的p⁽ⁿ⁺¹⁾，计算第n+1次迭代过程中u函数的值，即

u⁽ⁿ⁺¹⁾＝v⁽ⁿ⁾-θdiv(p⁽ⁿ⁺¹⁾)；

4b)更新正则化辅助变量v:

4b1)基于Gamma分布假设和Bhattacharyya距离度量，分别利用下式分别计算轮廓内外区域的能量值e₁和e₂：

其中，Γ(·)表示标准的伽马函数，ln(·)表示对数运算，(α_x,β_x)表示以像素点x为中心的局部区域内像素服从的Gamma分布参数对，(α_in,β_in)和(α_out,β_out)分别表示轮廓内部局部区域和轮廓外部局部区域对应的Gamma分布参数对，

具体地，(α_x,β_x)可由下式计算得到：

α_x＝μ²/σ²，β_x＝α_x/μ

其中，μ和σ²分别表示局部区域内的像素强度均值和方差；(α_in,β_in)和(α_out,β_out)分别可由下式计算得到：

其中，

和

分别表示轮廓内部局部区域的像素强度均值和方差，

和

分别表示轮廓外部局部区域的像素强度均值和方差，式中H(φ)和H(-φ)分别表示轮廓内部和轮廓外部的水平集函数取值；

4b2)根据计算得到的u⁽ⁿ⁺¹⁾、e₁和e₂，更新第n+1次迭代过程中v函数的值v⁽ⁿ⁺¹⁾：

v⁽ⁿ⁺¹⁾＝min{max{u⁽ⁿ⁺¹⁾-θ(λ₁e₁-λ₂e₂),0},1}

其中λ₁，λ₂为权值系数；

4c)更新迭代次数n＝n+1；

4d)判断是否满足终止条件||u⁽ⁿ⁺¹⁾-u⁽ⁿ⁾||₂≤ε或者n＞N，若满足，则终止迭代，输出目标轮廓曲线演化结果u；否则，继续执行步骤4a)～步骤4d)之间的操作。

3.根据权利要求1所述的方法，其中步骤5根据目标轮廓曲线演化结果u得到目标检测结果具体按以下方式进行：

对于任一像素点x，x∈x，若u(x)≥0.5，则将其视为目标像素点；否则视为杂波像素点；在对所有像素点进行判别之后，得到目标检测结果。

Images