CN103456015B - 一种基于最优分数域Gabor谱特征的SAR目标检测方法 - Google Patents

Abstract

本发明公开了一种基于最优分数域Gabor谱特征的SAR目标检测方法，属于图像处理技术领域，其包括如下步骤：读入SAR图像信号；将所述SAR图像信号按行、列方向展开为和的信号；对两个方向的信号设计最优窗函数，并做GT；对两个方向得到的空间‑频率谱做FrFT；对空间‑频率谱进行能量衰减梯度特征提取；将两个方向的特征表示对应空间位置做乘积，从而得到原始SAR图像的特征空间，完成检测。本发明将FrFT与GT结合，因为比常规算法多出了分数域参数，使得算法抗干扰性更强，且能够检测出微弱目标区域，比常规算检测算法精度更高，也适用于各种场景，具有较好的通用性。

Description

一种基于最优分数域Gabor谱特征的SAR目标检测方法

技术领域

本发明涉及图像处理技术领域，具体涉及一种基于最优分数域Gabor谱特征的SAR目标检测方法，可用于SAR目标的检测，其基于分数阶傅里叶变换(Fractional FourierTransform,FrFT)与时频分析方法(Time Frequency Analysis,TFA)，提取图像的分数域Gabor谱特征，从而高精度检测出目标。

背景技术

SAR具有全天时、全天候、穿透能力强、分辨率高等特点，已被广泛运用于军事侦察、地形测绘、地质勘探、海洋应用、农林监测、水资源利用和科学研究等领域。SAR发展迅速，SAR图像系统的信息获取技术明显超前于SAR图像成像后的信息处理技术，因此各个国家都对SAR ATR(Automation Target Recognition,ATR)展开了大量研究工作。SAR目标检测是SAR ATR的关键环节。

FrFT实际上就是Fourier变换的广义形式。FrFT的概念其实最早在1929年，就由N.Wiener提出，但是并没有受到太多的瞩目。1980年Namias从特征值和特征函数的角度，以纯数学的方式重新提出它的概念，并揭示了FRFT的几个特性以后，FrFT变换逐渐受到各国学者的重视。1993年Almeida指出分数阶傅里叶变换可以理解为时频平面的旋转。1996年Ozaktas等提出了一种计算量与FFT相当的离散算法后，开始出现大量相关研究文章。目前FrFT已广泛应用于雷达、声纳、通信、信息安全等众多领域。

时频分析的思想始于上世纪四十年代年D.Gabor提出的Gabor变换。常用的时频分析方法包括短时傅里叶变换(Short Time Fourier Transform，STFT)、Gabor变换(GaborTransform，GT)、小波变换、Winger-Vile分布(Wigner-Vile Distributioan，WVD)、S变换等。Gabor变换是通过时域窗函数将时间信号进行截断，对截断后的局部信号做Fourier变换，并表示为时间和频率的联合分布，来揭示信号的时变谱特征，满足瞬时信号分析的要求，更适用于非平稳信号的研究。对图像信号进行时频分析实际上就是进行图像信号的空间-频率分析。

近年来兴起的分数域时频分析方法正是在FrFT的基础上发展起来的。与传统时频分析方法相比，在某些条件下，分数域时频分析能够取得更好的效果，这是因为分数域时频分析方法比传统时频分析多了一个分数域参量。目前，分数域时频分析方法已经被广泛用在科学研究与工程应用的很多相关领域，如相位恢复、信号重构、信号检测、参数估计、分数域滤波器设计、神经网络、图像数字水印、图像加密、油气勘测等。

由于SAR图像的成像机理，所以其不能很好的描述轮廓和细节；其次，是SAR目标对方位角敏感，因此，不同目标及不同照射方位下所形成的SAR图像，其背景散射强度与目标散射强度也不相同；而且，SAR图像存在特有的斑点噪声，且在成像过程中辐射特征以及几何特征发生畸变，常规的图像检测方法不适用于SAR目标检测当中，因此需要一个高精度，抗噪性强的检测方法。

发明内容

针对上述现有技术，本发明的目的在于如何提供一种基于最优分数域Gabor谱特征的SAR目标检测方法，其旨在抑制SAR图像相干斑噪声，提升传统Gabor变换时频的聚集性，改善时频分布，并能更精确的完成SAR目标检测。

为达到上述目的，本发明采用如下技术方案：

一种基于最优分数域Gabor谱特征的SAR目标检测方法，其特征在于，包括如下步骤：

步骤一:读入m×n SAR图像信号；

步骤二:将所述SAR图像信号按行、列方向展开为1×MN和MN×1的信号；

步骤三：对两个方向的信号设计最优窗函数，并分别对它们做GT；

步骤四：对两个方向得到的空间-频率谱做FrFT，并搜索最优变换阶次，改善空间-频率分布；

步骤五：对经过步骤三和步骤四优化后的空间-频率谱进行能量衰减梯度特征提取，分别得到两个方向的特征表示；

步骤六：将两个方向的特征表示对应空间位置做乘积，从而得到原始SAR图像的特征空间，完成检测。

所述步骤三具体又分解为如下步骤：

步骤3.1：利用FrFT的旋转性，并结合广义时频带宽积准则设计出最优窗，对行方向信号进行GT后，得到MN×MN的空间频率谱；

步骤3.2：对于列方向也进行与步骤3.1相同的处理。

所述步骤四具体包括如下步骤：

步骤4.1：以经过步骤3处理过的行方向的空间-频率谱为基础，固定空间位置，抽取对应位置以频率表示的MN道信号，分别对每一道信号进行不同阶次的FrFT，并通过能量峰值法对每一道信号进行最优阶次的FrFT；

步骤4.2：对于列方向也进行与步骤4.1相同的处理。

所述步骤五具体包括如下步骤：

步骤5.1：以行方向最优空间-频率谱为基础，固定空间位置，抽取对应位置以频率表示的MN道信号，找到每一道信号峰值Max的位置，记为f_max，在f_max之后找到总能量的65％对应的位置，记为f₀；在f_max之后找到总能量的85％对应的位置，记为f₁；计算f₀,f₁,f_max的拟合斜率得到1×MN的行方向空间-频率谱特征表示。

步骤5.2：对于列方向也进行与步骤5.1相同的处理。

步骤5.3：对行方向1×MN按对应位置重置为M×N，得到原始SAR图像的行方向的特征表示；对行方向MN×1按对应位置重置为M×N，得到原始SAR图像的列方向的特征表示；

与现有技术相比，本发明具有以下有益效果：

一、本发明采用分数域时频分析方法，并成功将其应用于SAR图像处理，对非平稳信号和信号瞬态分析提供了一种新的思路；本发明抗噪性强、检测精度高，时频特征剖面加大了目标区域和背景区域的对比度，同时抑制了相干斑噪声，克服了常规检测方法对噪声敏感的缺点。

二、本发明将FrFT引入窗函数的设计，利用分数阶旋转性，并搜索最优旋转因子，找到满足GTBP的最优阶次，完成窗函数设计，显著提升空间-频率聚集性。

三、本发明将FrFT与GT结合，因为比常规算法多出了分数域参数，使得算法抗干扰性更强，且能够检测出微弱目标区域，比常规算检测算法精度更高，也适用于各种场景，具有较好的通用性。

附图说明

图1为本发明方法流程示意图；

图2为FrFT的对时频面的旋转性示意图；

图3为最优旋转因子搜索示意图；

图4为能量衰减梯度示意图。

具体实施方式

下面将结合附图及具体实施方式对本发明作进一步的描述。

一种基于最优分数域Gabor谱特征的SAR目标检测方法，该方法包括以下步骤：

步骤一：读入M×N维SAR图像信号；

步骤二：将读入的图像信号按行、列方向展开为1×MN和MN×1的信号；

步骤三：对两个方向的信号设计最优窗函数，并分别对它们做GT，使得GT的空间-频率聚集性得到显著提升，从而抑制杂波。

步骤四：对两个方向得到的空间-频率谱做FrFT，并搜索最优变换阶次，改善空间-频率分布，从而对微弱目标敏感，提高检测率；

步骤五：对经过步骤三和步骤4优化后的空间-频率谱进行能量衰减梯度特征提取，分别得到两个方向的特征表示；

步骤六：将两个方向的特征表示对应空间位置做乘积，从而得到原始SAR图像的特征空间，完成检测。

所述步骤三具体包括如下步骤：

步骤3.1：对行方向信号进行GT后，得到MN×MN的空间-频率谱。利用FrFT的旋转性，并结合广义时频带宽积准则(Generalized Time-Frequency band product,GTBP)，使得空间宽和频率宽都达到最小，从而设计出最优窗，使得GT拥有很好的空频聚集性。

如图2所示FrFT对时频面的旋转，相当于把信号的时频轴在时频面上逆时针旋转α＝pπ/2，p为变换阶次。

步骤3.1.1：改变旋转阶次使得其在分数域上持续时间最短，从而拥有最大的模值，该阶次既为最优旋转因子。

如图3所示，图3(a)所示为信号的原始时频图，图3(b)所示，当时频图的长轴恰好和分数域频率轴重合，此时时频图的短轴将与时间轴重合时，则信号持续时间最短，时频图的模值最大，据此可以通过遍历不同阶次FrFT，来确定最优旋转因子。

步骤3.1.2：求得最优旋转因子后，计算出最优窗函数；

步骤3.2：对于列方向也进行与步骤3.1相同的处理。

如图4所示，对理论chirp信号分别用最优窗进行GT和普通窗进行GT，时频聚集性上得到了显著的提升，能够很好的抑制因GT本身而产生的杂波，为高精度检测打下基础。

所述步骤四具体包括如下步骤：

步骤4.1：以经过步骤3处理过的行方向的空间-频率谱为基础。固定空间位置，抽取对应位置以频率表示的MN道信号，分别对每一道信号进行不同阶次的FrFT，并通过能量峰值法对对每个空间位置的频率信号都进行最优FrFT，最优阶由下式确定：

P_opt＝arg max(Energy_xp)

其中，Energy_xp表示不同阶次下最优窗Gabor变换得到的空间-频率谱所抽取的每个空间位置的频率信号的能量值，从而改善时频分布，使得本方法对微弱目标更敏感，提高了检测精度。

步骤4.2：对于列方向也进行与步骤4.1相同的处理。

对两个方向得到的空间-频率谱做FrFT，并搜索最优变换阶次，改善空间-频率分布，从而对微弱目标敏感，提高检测率。

所述步骤五具体包括如下步骤：

通过步骤三与步骤四的处理，我们得到了经过最优窗设计和最优FrFT后的最优空间-频率谱。

步骤5.1：以行方向最优空间-频率谱为基础，固定空间位置，抽取对应位置以频率表示的MN道信号，找到每一道信号峰值Max的位置，记为f_max，在f_max之后找到总能量的65％对应的位置，记为f₀；在f_max之后找到总能量的85％对应的位置，记为f₁；计算f₀,f₁,f_max的拟合斜率，实现了空间-频率谱的降维，得到1×MN的行方向空间-频率谱特征表示；

步骤5.2：对于列方向也进行与步骤4.1相同的处理；

步骤5.3：对行方向1×MN按对应位置重置为M×N，得到原始SAR图像的行方向的特征表示；对行方向MN×1按对应位置重置为M×N，得到原始SAR图像的列方向的特征表示。

步骤六：将两个方向的特征表示对应空间位置做乘积，从而得到原始SAR图像的特征空间，完成检测。

传统FT是一种全局性变换，能够得到信号的整体频谱，因此，FT不能用来有效处理非平稳信号。为了克服FT在信号处理中的局限性，更好的分析信号的局部特性，提出了FrFT，时间信号x(t)的p阶FrFT变换定义为：

$X_p\left(u\right)={\∫}_{-\mathrm{\∞}}^{+\mathrm{\∞}}x\left(\mathrm{\τ}\right)K_p\left(u,\mathrm{\τ}\right)dt$

其中，p为FrFT的阶次，其周期为4，α＝pπ/2为分数域坐标轴u相对时间轴t逆时针旋转的角度，当p＝1时实际上就是FT，其中K_p(t,u)表示FrFT的核函数，可以表示为：

$K_p\left(u,\mathrm{\τ}\right)=\left\{\begin{array}{cc}\sqrt{\frac{1-j\cot \mathrm{\α}}{2\mathrm{\π}}}\exp \left(j\frac{{\mathrm{\τ}}^2+u^2}{2}\cot \mathrm{\α}-\frac{j\mathrm{\τ}u}{\sin \mathrm{\α}}\right)& \mathrm{\α}\≠n\mathrm{\π}\\ \mathrm{\δ}\left(\mathrm{\τ}-u\right)& \mathrm{\α}=2n\mathrm{\π}\\ \mathrm{\δ}\left(\mathrm{\τ}+u\right)& \mathrm{\α}=\left(2n\±1\right)\mathrm{\π}\end{array}\right.$

下面给出FrFT的时移性和频移性：

${\mathrm{FrFT}}_p\left(x\right(t-\mathrm{\τ}\left)\right)=e^{{\mathrm{j\π\τ}}^{2}sin\mathrm{\φ}cos\mathrm{\φ}}{x}_p(u-\mathrm{\τ}cos\mathrm{\φ})e^{-j2\mathrm{\π}u\mathrm{\τ}\sin \mathrm{\φ}}$

${\mathrm{FrFT}}_p\left(e^{j2\mathrm{\π}tf}x\right(t\left)\right)=e^{-{\mathrm{j\πf}}^{2}sin\mathrm{\φ}cos\mathrm{\φ}}{x}_p(u-fsin\mathrm{\φ})e^{j2\mathrm{\π}ufcos\mathrm{\φ}}$

FrFT的还有一个重要性质是分数阶核函数的可分离性：

K_p1,p2(τ,η,u,v)＝K_p1,(τ,u)K_p2(η,v)

其中K_p1,p2(τ,η,u,v)为2D-FrFT的核函数，τ,η为空间域变量，u,v表示分数域变量，K_p1,(τ,u),K_p2(η,v)是对不同空间轴进行的不同阶次的FrFT。

FrFT变换的旋转性在本发明中可视为空间-频率坐标轴逆时针旋转角度φ，R表示逆时针旋转，FrFT与Gabor变换的关系如式:

${\mathrm{GT}}_x(t,f)=R_{\mathrm{\φ}}\left\{{\mathrm{GT}}_{x_p}(t,f)\right\}$

Gabor变换是通过时域窗函数将时间信号进行截断，对截断后的局部信号做Fourier变换，并表示为时间和频率的联合分布，来揭示信号的时变谱特征，满足瞬时信号分析的要求，更适用于非平稳信号的研究。对图像信号进行时频分析实际上就是进行图像信号的空间-频率分析。因为相比其他时频分析方法运算速度更快，所以本方法采用Gabor变换来处理原始图像，其定义为：

${\mathrm{GT}}_x(t,\mathrm{\ω})=\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}x\left(\mathrm{\τ}\right)h(\mathrm{\τ}-t)\exp (-j\mathrm{\ω}\mathrm{\τ})d\mathrm{\τ}$

其中，t,ω表示时间、频率分量，h(τ-t)表示窗函数，GT实际上就是加窗傅里叶变换，且窗函数是高斯窗函数。GT后信号被映射到一个2D时频谱，即时间和频率的联合表示，对图像信号而言即空间和频率的联合表示，因为包含了位置信息，因此时频分析方法可以对信号的局部特性进行分析。

本文将从信号能量守恒的角度，提出一种简单、快速的搜索算法，根据FrFT的Parseval定理，FrFT需满足以下关系，

$\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}|x\left(t\right){|}^{2}dt=\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}|x_p\left(u\right){|}^{2}du=C$

其中C表示常数，|x_p(u)|在分数域的积分为一个定值，如果该信号在分数域上持续得越短，则该信号模的峰值越大。

本方法将FrFT与GT相结合，改善时频谱分布。用FrFT的核函数替换傅里叶变换核exp(-jωτ)得到FrGT变换的定义：

${\mathrm{GT}}_p(t,u)=\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}x\left(\mathrm{\τ}\right)h(\mathrm{\τ}-t)K_p(\mathrm{\τ},u)d\mathrm{\τ}$

把2D图像信号的FrGT变换看作一维Gabor变换的扩展，变换后的信号包含4个变量(x,y,u,ν)，时频谱变为4维，得到其表达式为：

${\mathrm{GT}}_{p1,p2}(x,y,u,v)=\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}f(\mathrm{\τ},\mathrm{\η})K_{p1,p2}(\mathrm{\τ},\mathrm{\η},u,v)h(\mathrm{\τ}-x,\mathrm{\η}-y)d\mathrm{\τ}d\mathrm{\η}$

其中f(τ,η)为SAR图像信号，τ,η为空间域变量，u,v表示分数域变量，为2D-FrFT的分数阶核函数，h(τ-x,η-y)为2D高斯窗函数。利用核函数性质，且高斯窗函数也是可分离的，并交换积分次序：

$\begin{array}{cc}\begin{array}{c}{\mathrm{GT}}_{p1,p2}\left(x,y,u,v\right)=\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}f\left(\mathrm{\τ},\mathrm{\η}\right)K_{p1,}\left(\mathrm{\τ},u\right)K_{p2}\left(\mathrm{\η},v\right)\frac{1}{2\mathrm{\π}}\exp \left(-\[\frac{{\left(\mathrm{\τ}-x\right)}^2}{2{\mathrm{\σ}}^2}+\frac{{\left(\mathrm{\η}-y\right)}^2}{2{\mathrm{\σ}}^2}\]\right)d\mathrm{\τ}d\mathrm{\η}\\ =\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}\[\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}f\left(\mathrm{\τ},\mathrm{\η}\right)K_{p2}\left(\mathrm{\η},v\right)\frac{1}{\sqrt{2\mathrm{\π}}}\exp \left(-\frac{{\left(\mathrm{\η}-y\right)}^2}{2}\right)d\mathrm{\η}\]K_{p1}\left(\mathrm{\τ},u\right)\frac{1}{\sqrt{2\mathrm{\π}}}\exp \left(-\frac{{\left(\mathrm{\τ}-x\right)}^2}{2}\right)d\mathrm{\τ}\\ =\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}\[\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}f\left(\mathrm{\τ},\mathrm{\η}\right)K_{p2}\left(\mathrm{\η},v\right)h\left(\mathrm{\η}-y\right)d\mathrm{\η}\]K_{p1}\left(\mathrm{\τ},u\right)h\left(\mathrm{\τ}-x\right)d\mathrm{\τ}\end{array}& (\end{array}$

由上式可知，进行2维信号的FrGT时，可以先沿着一个方向进行GT，然后再沿着另一个方向进行，再对两个方向的计算结果求与，最终完成2D-FrGT，对M×N维SAR图像信号，其离散形式可表示为：

$\begin{array}{c}{\mathrm{GT}}_{p1,p2}\left(m,n,\frac{k}{{\mathrm{KT}}_x},\frac{q}{{\mathrm{QT}}_y}\right)\\ =\underset{x=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{k=0}{\overset{K-1}{\mathrm{\Σ}}}\[\underset{y=0}{\overset{N-1}{\Σ}}\underset{q=0}{\overset{Q-1}{\Σ}}f(x,y)K_{p2}(\mathrm{\η},\frac{q}{{\mathrm{QT}}_y})h(y-n)\]K_{p1}(x,\frac{k}{{\mathrm{KT}}_x})h(x-m)\end{array}$

式中K,Q为分数域采样点数，Tx,Ty为采样间隔，M,N为图像尺寸。

以上所述仅为本发明的实施例，并非因此限制本发明的专利范围，凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其他相关的技术领域，均同理包括在本发明的专利保护范围内。

Claims (1)

1.一种基于最优分数域Gabor谱特征的SAR目标检测方法，包括如下步骤：

步骤一:读入m×n SAR图像信号；

步骤二:将所述SAR图像信号按行、列方向展开为1×MN和MN×1的信号；

步骤三：对两个方向的信号设计最优窗函数，并分别对它们做GT；

步骤四：对两个方向得到的空间-频率谱做FrFT，并搜索最优变换阶次，改善空间-频率分布；

步骤五：对经过步骤三和步骤四优化后的空间-频率谱进行能量衰减梯度特征提取，分别得到两个方向的特征表示；

步骤六：将两个方向的特征表示对应空间位置做乘积，从而得到原始SAR图像的特征空间，完成检测；

所述步骤三具体又分解为如下步骤：

步骤3.1：利用FrFT的旋转性，并结合广义时频带宽积准则设计出最优窗,对行方向信号进行GT后，得到MN×MN的空间-频率谱；

步骤3.2：对于列方向也进行与步骤3.1相同的处理；

所述步骤四具体包括如下步骤：

步骤4.1：以经过步骤3处理过的行方向的空间-频率谱为基础，固定空间位置，抽取对应位置以频率表示的MN道信号，分别对每一道信号进行不同阶次的FrFT，并通过能量峰值法对每一道信号进行最优阶次的FrFT；

步骤4.2：对于列方向也进行与步骤4.1相同的处理；

所述步骤五具体包括如下步骤：

步骤5.1：以行方向最优空间-频率谱为基础，固定空间位置，抽取对应位置以频率表示的MN道信号，找到每一道信号峰值Max的位置，记为f_max，在f_max之后找到总能量的65％对应的位置，记为f₀；在f_max之后找到总能量的85％对应的位置，记为f₁；计算f₀,f₁,f_max的拟合斜率得到1×MN的行方向空间-频率谱特征表示；

步骤5.2：对于列方向也进行与步骤5.1相同的处理；

步骤5.3：对行方向1×MN按对应位置重置为M×N，得到原始SAR图像的行方向的特征表示；对列方向MN×1按对应位置重置为M×N，得到原始SAR图像的列方向的特征表示。