基于短时分数阶傅里叶变换的海面微动目标检测和特征提取方法

摘要

本发明涉及一种基于短时分数阶傅里叶变换(STFRFT)的海面微动目标检测和特征提取方法，属于雷达信号处理和检测技术领域。包括以下步骤：1)海尖峰识别，将海杂波分为海尖峰序列和不含海尖峰的海杂波背景序列；2)海杂波数据筛选，选取最小平均功率对应的海杂波背景序列作为待检测数据；3)FRFT域微动目标检测，以FRFT域信号幅值作为检测统计量并与门限比较；4)最佳FRFT域滤波，通过窄带带通滤波器，提取多分量微动信号；5)设置最优时间窗长度，在STFRFT域估计微动特征。该方法能自适应抑制海杂波，改善信杂比，有效分离和提取多分量微动信号，为海面弱目标检测和特征提取提供了新的途径，具有推广应用价值。

说明书

一、技术领域

本发明属于雷达信号处理和检测技术领域，特别涉及对海探测雷达的动目标检测技术。

二、背景技术

海杂波中微弱目标尤其是“低(低掠射角)、慢(静止或慢速运动)、小(目标尺寸小)” 目标的检测技术始终是雷达信号处理领域的难题，不仅具有理论重要性，而且在军用和民用 上均占有非常重要的地位。传统的目标检测方法都是基于统计理论的，将海杂波视为随机过 程，但目标模型和杂波模型均呈多样化发展趋势，尤其是杂波分布模型，在复杂环境中往往 不成立或不完全成立，这就使经典目标检测方法由于模型失配而不能取得预期的检测结果。。 Haykin等从混沌和分形等角度出发，采用关联维数和盒维数等特征量从海杂波中检测微弱目 标，但当信杂比(Signal-to-Clutter Ratio，SCR)持续降低时，检测结果不尽人意，而且很难实 现对运动目标的检测。

海杂波中微弱目标的共同点是，由于雷达分辨低、距离远、背景强等因素，无论是在时 域还是在频域，目标分辨单元中的信杂(噪)比都很低，海杂波会淹没微弱目标信号，能量 积累效果差，这些都严重降低了检测性能。与地物杂波不同，海面的粗糙程度要远远高于地 面，并且海面不断地运动起伏，幅值分布复杂，在例如低掠射角以及高海情的条件下，大量 的杂波尖峰还会造成严重虚警。研究表明，强海杂波背景下，杂波幅度概率密度函数表现出 较长的“拖尾”现象，称为海尖峰效应。此时的多普勒谱包括由“慢变信号”引起的Bragg散射， 以及“快变信号”产生的非Bragg散射(海尖峰，白浪等破碎波)，表现为随机分布在不同距离、 不同角度上零星运动或静止的目标，具有强烈的回波起伏特性，严重影响了雷达对小目标的 探测性能。非Bragg散射导致Bragg谱展宽，体现出明显的时变和非平稳特性，它的统计特 性和多普勒频率随时间变化，因此海杂波的频谱是时间和频率两个变量的函数，应从时间和 频率两方面对海杂波数据进行分析和处理。

当海面环境较为平稳时(通常低于2级海况)，海面目标与雷达的相对平动对回波的多普 勒频率影响最大。然而，在高海况时，海上目标随海面颠簸导致姿态变化，引起雷达回波功 率调制效应，不仅存在平动，舰船还绕参考点作三轴转动(偏航运动、俯仰运动、横滚运动)， 导致散射点的多普勒频率随时间非线性变化。因此，新的检测算法应能较好地处理时变、非 平稳和非线性信号。近年来，微多普勒理论的研究已成为信号处理领域的一个新的技术热点， 将目标或目标部件除质心平动以外的振动、转动和加速运动等微小运动统称为微动，而由目 标微动所引起的多普勒频移称为微多普勒。目标的微动特征反映了目标的精细运动和几何结 构对电磁散射的综合调制特征，反映了多普勒变化特性，为雷达目标检测和特征提取提供了 新的途径。目标距离的周期性将引起目标多普勒频率的周期性，目标姿态变化将对回波产生 调制特性，反映出目标的瞬时速度变化特性。而海面目标，如舰船目标，其回波多普勒谱实 质是平动和微动共同作用的结果。因此，微多普勒非常适于分析海杂波的非Bragg散射特性， 以及高海况条件下的目标检测。

时频分析方法是研究非平稳信号的有力工具，能够给出特定时间和特定频率范围的能量 分布，描述了非平稳信号的频率随时间的变化过程。因此，时频分布被广泛用于分析微多普 勒的特性，并进行微动特征提取与参数估计。第一类时频表示方法是线性时频表示，例如短 时傅里叶变换(Short-Time Fourier Transform，STFT)，无交叉项，但受不确定原理的制约，时 频分辨率较低；另一种时频表示方法是非线性时频表示，例如Winger-Ville分布(Wigner-Ville  Distribution，WVD)，在多目标存在的情况下，运算量大且交叉项严重影响微动信号的检测。 由于微动信号在一定程度上可近似为调幅-线性调频信号(Amplitude Modulation-Linear  Frequency Modulation，AM-LFM)，或者近似分段的AM-LFM信号，定义为

$s\left(t\right)=\underset{i}{\Sigma }{A}_i\left(t\right)\exp (j2\pi f_{i}t+\mathrm{j\pi }{k}_i{t}^2),\left|t\right|\le T$

式中，A_i(t)为第i个微动信号幅度，f和k分别为中心频率和调频率。

分数阶傅里叶变换(Fractional Fourier Transform，FRFT)对LFM有良好的能量聚集性和 检测性能，无交叉项的干扰，但FRFT缺少时域定位的功能。在FRFT中加入滑动的短时窗 函数，得到短时FRFT(Short-Time FRFT，STFRFT)，通过窗函数的滑动完成整个时间上的信 号局部性质分析，可得到在任意时刻的该段信号的频率变化，极大地扩展了FRFT的应用范 围。因此，采用STFRFT处理微动信号，能够提高强海杂波背景下的微动目标检测性能。

三、发明内容

1.要解决的技术问题

本发明的目的在于利用海面目标的微动特征，提供一种基于短时分数阶傅里叶变换的海 面微动目标检测和特征提取方法，其中要解决的技术问题包括：

(1)海杂波的非Bragg散射，如海尖峰，在高海况时尤为明显，使得目标回波被海杂波覆 盖，SCR低，容易造成虚警；

(2)传统的检测方法未利用目标的微动特征，丢失有用信息；

(3)较为合理地设置STFRFT的时间窗长度；

(4)较为合理地设置最佳FRFT域滤波器参数，更好地提取多分量微动信号特征。

2.技术方案

本发明所述的基于短时分数阶傅里叶变换的海面微动目标检测和特征提取方法，其特征 在于包括以下技术措施：

步骤一、海尖峰识别：在接收端，将接收并经过放大和解调处理后得到的同一距离单元 内的雷达回波数据根据尖峰幅度门限、尖峰宽度和尖峰间隔约束条件进行海尖峰判断，将海 杂波数据分为海尖峰序列和不含海尖峰的海杂波背景序列；

步骤二、海杂波数据筛选：分别计算各个海杂波背景序列的平均功率水平，选取最小平 均功率对应的海杂波背景序列作为待检测数据；

步骤三、FRFT域微动目标检测：对待检测数据进行不同变换阶数下的离散FRFT (DFRFT)，取信号的幅值作为检测统计量，与给定虚警概率下的检测门限进行比较，如果检 测统计量高于门限值，判决为存在微动目标信号，否则判决为没有微动目标信号，同时将超 过门限最大峰值所对应的FRFT域作为最佳FRFT域；

步骤四、最佳FRFT域滤波：将最佳FRFT域数据通过构造的窄带带通滤波器，滤除除 第一微动信号分量外的绝大部分信号能量，对滤波后的最佳FRFT域数据进行逆DFRFT运算， 获得第一微动信号数据，将原始雷达回波数据与第一微动信号数据相消，获得去除第一微动 信号数据；

步骤五、STFRFT域微动特征估计：在有限的观测时间范围内设置时间窗长度，对第一 微动信号数据进行最佳变换阶数下的STFRFT运算，得到第一微动信号的瞬时频率估计；

步骤六、对去除第一微动信号的雷达数据继续进行步骤三至步骤五的运算，直到检测不 出微动信号为止。

3.有益效果

对比现有技术，本技术方案所述的基于短时分数阶傅里叶变换的海面微动目标检测和特 征提取方法，有益效果在于：

(1)该方法可以在抑制海尖峰的同时保存更多的信号能量，改善信杂比；

(2)该方法对海杂波与微动目标具有良好的区分能力；

(3)该方法能有效分离和提取多分量微动信号；

(4)该方法具有在强海杂波中检测微动目标和估计微动参数的能力。

四、附图说明

说明书附图1是基于短时分数阶傅里叶变换的海面微动目标检测和特征提取方法的系统 结构图，图2是基于短时分数阶傅里叶变换的海面微动目标检测和特征提取方法的实现流程 图。

五、具体实施方式

以下结合说明书附图1对本发明作进一步详细描述。参照说明书附图1，本发明的处理 流程分以下步骤：

(1)海尖峰识别；

(2)海杂波数据筛选；

(3)FRFT域微动目标检测；

(4)最佳FRFT域滤波；

(5)STFRFT域微动特征估计。

以下结合说明书附图2对本发明作进一步详细描述。参照说明书附图2，本发明的具体 实施方式分以下步骤：

(1)将经过放大和解调处理后得到的同一距离单元内的雷达回波数据送入存储装置1进 行预处理，得到N个采样点的输入信号序列c(i)，i＝1，2，...，N。

(2)将输入信号序列c(i)与给定的海尖峰门限T_s、最小尖峰宽度W_min和最小尖峰间隔I_min一同送入比较器2，进行海尖峰判断，将海杂波数据分为海尖峰序列和不含海尖峰的海杂波 背景序列，判定准则如下

$\left(\begin{array}{c}\left|c\left(i\right)\right|\ge T_s\\ W_s\ge W_{\min }\\ I_s\ge I_{\min }\end{array}\right)$

式中，在实际工程中，L设为5，W_min设为0.1s，I_min取值为0.5s(关于T_s、 W_min和I_min的选取参见文献：M.Greco，P.Stinco，F.Gini，“Identification and analysis of sea radar  clutter spikes，”IET Radar，Sonar and Navigation，vol.4，no.2，pp.239-250，2010.)。

(3)计算装置3接收比较器2的输出结果，并将海尖峰序列置零，仅保留不含海尖峰的 海杂波背景序列。

(4)计算装置4接收计算装置3的输出结果，分别计算各个海杂波背景序列的平均功率 水平，选取最小平均功率对应的海杂波背景序列作为待检测数据，计算方法如下

$i_0=\underset{i}{\arg \min }\left(\frac{1}{M_i}\underset{j=1}{\overset{M_i}{\Sigma }}{\left|x_i\left(j\right)\right|}^2\right),$i＝1，2，...，n，j＝1，2，...，M_i

式中，x_i为第i个海杂波背景序列，M_i为序列长度，i₀为最小平均功率对应的海杂波背景序 列序号。

(5)计算装置5接收计算装置4输出的海杂波背景序列采用H.M.Ozaktas等人提出 的FRFT的分解算法，完成不同变换阶数下的DFRFT运算，基于如下公式

$X_p\left(\frac{m}{2\mathrm{\Delta x}}\right)=A_{\alpha }{e}^{j\frac{1}{2}{\left(\frac{m}{2\mathrm{\Delta x}}\right)}^2(\cot \alpha -\csc \alpha )}\underset{n=-N}{\overset{N}{\Sigma }}\left[x\left(\frac{n}{2\mathrm{\Delta x}}{e}^{j\frac{1}{2}{\left(\frac{m}{2\mathrm{\Delta x}}\right)}^2(\cot \alpha -\csc \alpha )}\right)\right]e^{j\frac{1}{2}{\left(\frac{m-n}{2\mathrm{\Delta x}}\right)}^2\csc \alpha }$

式中，α为旋转角度，α＝pπ/2，p为变换阶数，Δx为信号的带宽，N＝(Δx)²。

(6)将计算装置5的输出幅值作为检测统计量，与给定虚警概率下的检测门限一同送入 比较器6，进行比较，如果检测统计量低于检测门限，判决为没有微动目标信号，结束算法， 若检测统计量高于检测门限，则判决为存在微动目标信号，并将最大峰值(第一微动信号分 量)对应的FRFT域坐标送入存储装置7中，即

$(p_1,u_1)=\underset{p,u}{\arg \min }\left|\mathrm{DFRFT}(p,u)\right|$

将最佳变换域数据送入存储装置8中。

(7)计算装置9根据存储装置7的输出结果，构造最佳FRFT域窄带带通滤波器H(p₁，u)

式中，W1为尖峰宽度，$W_1=\frac{2\pi }{T_{i_0}\csc (p_1\pi /2)},$为观测时长，$T_{i_0}=\frac{M_{i_0}}{f_s},$f_s为采样频率。

(8)将计算装置8和计算装置9的输出结果一同送入计算装置10，进行点乘运算

$X_{p_1}^{\prime }=X_{p_1}·H(p_1,u)$

滤除除第一微动信号分量外的绝大部分信号能量。

(9)计算装置11接收计算装置10的输出结果，对滤波后的最佳FRFT域数据进行-p₁阶DFRFT运算(逆DFRFT)，获得第一微动信号数据。

(10)计算装置12接收存储装置7的最佳变换阶数p₁，估计第一微动信号的调频率，估 计方法为

${\hat{k}}_1=-\cot (p_1\pi /2)/S^2$

式中，$S=\sqrt{T_{i_0}/f_s}$为尺度因子。

(11)计算装置13根据计算装置12输出的调频率估计值，STFRFT运算采用的时间窗 长度为

$N_w=\left(\begin{array}{cc}{f}_s\sqrt{1/{\hat{k}}_1},& \sqrt{1/{\hat{k}}_1}\le T_{i_0}\\ M_{i_0},& \sqrt{1/{\hat{k}}_1}>T_{i_0}\end{array}\right)$

使其既能保证参数估计精度，又能保证较高的时间分辨率。

(12)计算装置14接收计算装置12和计算装置13的输出结果，计算在步骤(11)的时 间窗长度下的第一微动信号最佳变换阶数p′₁

$p_1^{\prime }=\frac{2\mathrm{arccot}\left[{\hat{k}}_1{\left(S^{\prime }\right)}^2\right]}{\pi }+2$

式中，$S^{\prime }=\sqrt{N_w/f_s^2}.$

(13)计算装置15接收计算装置13和计算装置14的输出结果，根据时间窗长度N_w， 对第一微动信号数据进行p′₁阶的STFRFT运算，基于如下公式

${\mathrm{STFRFT}}_{p_1^{\prime }}(t,u)={\int }_{-\infty }^{+\infty }x\left(\tau \right)g(\tau -g)K_{{\alpha }_1^{\prime }}(\tau ,u)\mathrm{d\tau }$

式中，为核函数，$K_{{\alpha }_1^{\prime }}(t,u)=\left(\begin{array}{cc}{A}_{{\alpha }_1^{\prime }}\exp \left\{j\right[\frac{1}{2}{t}^2\cot {\alpha }_1^{\prime }-\mathrm{ut}\csc {\alpha }_1^{\prime }+\frac{1}{2}{u}^2\cot {\alpha }_1^{\prime }\left]\right\}& {\alpha }_1^{\prime }\ne \mathrm{n\pi }\\ \delta [u-{(-1)}^{n}t]& {\alpha }_1^{\prime }=\mathrm{n\pi }\end{array}\right),$g(τ)为窗函数，一般取如下标准高斯窗函数，使STFRFT的聚集性达到最佳

$g\left(\tau \right)=\frac{1}{\sqrt{2\pi }}\exp (-\frac{{\tau }^2}{2}),\left|\tau \right|\le N_w/f_s$

(14)计算装置16接收计算装置15的输出结果，估计第一微动信号的瞬时频率，估计 方法为

${\hat{f}}_1\left(t\right)=u\left(t\right)\csc (p_1^{\prime }\pi /2)/S^{\prime }$

(15)将计算装置3的输出结果和计算装置11的输出结果一同送入计算装置17，作减 法运算，得到去除第一微动信号的雷达数据，并将输出结果送入计算装置5，继续进行步骤 (5)～步骤(14)的运算，直到检测不出微动信号为止。