CN102866391A - 基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法 - Google Patents

Abstract

本发明公开的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法属雷达目标探测技术领域，本方法首先利用短时傅里叶变换进行信号初步检测，然后采用二值化方法对初检结果进行处理，在处理中保留信号的相位信息，接着对短时傅里叶逆变换还原后的信号，采用分数阶傅里叶变换进行检测，采用多种方法进行联合处理的优点是可以克服强信号旁瓣对弱信号主瓣的压制现象，提高待检测信号的信噪比，解决其它已有方法在低信噪比下检测信号时出现的虚警概率较大的问题；同时图像对比度方法和逐步消去法的采用，利用了信号的空间和功率强度信息，能够检测调频率不同或相同的多个强弱信号，进一步提高检测概率和计算效率，易于工程实现，值得采用和推广。

Description

基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法

一.技术领域

本发明公开的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法属雷达目标探测技术领域，具体涉及的是一种利用短时傅里叶变换和分数阶傅里叶变换完成对匀加速/匀减速多目标的检测，解决弱目标信号主瓣被强目标旁瓣干扰的问题，提高复杂环境中多目标的检测能力。

二.背景技术

线性调频信号(LFM)的检测问题在非平稳信号处理领域中占有重要地位，它有以下两点原因：(1)LFM信号具有时宽带宽积大，抗干扰性能强，频移不敏感等性质，可作为雷达、声纳、通信等系统的发射调制信号；(2)在雷达、声纳和通信等领域，变速运动目标的反射信号可近似为LFM信号或近似为分段LFM信号。因此，上述系统的性能与LFM信号检测方法的优劣有直接关系。针对单目标环境中的LFM信号检测问题，国内外已有很多研究成果可应用于实际系统；但在雷达、通信及电子对抗系统中，经常需要处理未知参数的多分量线性调频信号。如弹道导弹逼近告警雷达系统时，导弹会释放多个诱饵或“子弹头”，这些多目标相对于雷达具有不同的径向加速度，其回波表现为多个LFM信号的叠加。而在电子对抗系统中，为了对敌方雷达系统实施干扰，经常需要在复杂的电磁环境中对未知参数的多个LFM雷达信号进行检测。多弹头、多诱饵、编队飞行等目标具有在空间维和距离维上分布较近的特点，现有窄带雷达无法分辨，可通过加长积累时间从频率维上进行分辨，这一问题能等效为调频率相同，初始频率不同的多个LFM信号的检测问题。目前国内外学者对多分量LFM信号的分析主要集中于线性变换的短时傅里叶变换(STFT)和基于时频分布的魏格纳(WVD)方法。STFT的基本思想是首先对观测信号进行加窗移位，然后求取加窗信号的傅里叶变换，因此，STFT可用一段时间信号来表示观测信号在这个时段的频谱特性，但受到不确定性定理的约束，时间分辨率和频率分辨率不可能同时提高。窗越宽，时间分辨率越低；窗宽度越窄，时间分辨率会提高，但频率分辨率又会降低。为了提高时频分辨率，WVD分布被广泛研究，它作为一种双线性时频分布，对单分量LFM信号具有很好的能量聚集性，但是在分析多分量LFM信号时会产生严重的交叉项，使得时频分布的能量聚集性也随之下降，尽管目前已提出许多抑制交叉项的有效方法，但是它们是以降低信号的时频聚集性为代价。近年来分数阶傅里叶变换(FRFT)作为一种广义的傅里叶变换具有较强的LFM信号分析能力和良好的噪声抑制能力受到关注。分数傅里叶变换将信号转换到一维参数空间获取LFM信号的调频率估计，然后采用解调频的方法获得其它参数估计值。但对多分量的LFM信号检测会存在强目标信号旁瓣对弱目标信号主瓣的干扰，或者主瓣信号被旁瓣信号干扰，为此有研究者提出利用逐次消去法(CLEAN)进行迭代求解，以获得多个信号的参数信息，但该方法仍存在低信噪比下检测性能低，且运算量较大的问题。

本发明在分析多种方法基础上，提出采用多种方法进行联合处理，可进一步提高检测概率和计算效率。本发明针对线性调频窄带雷达在角度维、距离维无法分辨多个目标的问题，提出通过多普勒频率维对多目标进行检测。在FRFT基础上，提出了一种首先利用STFT进行信号初步检测，提取可能的目标信号，然后采用二值化方法对其进行处理，这里与图像处理中二值化方法不同的是要保留信号的相位信息，然后利用STFT的可逆性进行信号还原，接下来对还原后的信号利用分数阶核函数角度旋转的特点，在最优旋转角时得到信号的最佳能量集聚，同时配合图像对比度方法和CLEAN方法进行逐次滤波检测，多方法的联合使用过程如图一所示。本方法不仅能够检测调频率不同的多个强弱信号，而且也能检测调频率相同的多个强弱信号，计算量较小，易于工程实现。

三.发明内容

本发明的目的是：向社会提供这种基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法。本发明在分析多种方法基础上，提出采用多种方法进行联合处理，可进一步提高检测概率和计算效率。本方法不仅能够检测调频率不同的多个强弱信号，而且也能检测调频率相同的多个强弱信号，具有计算量较小、易于工程实现等优点。

本发明的技术方案是这样的：这种基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法，技术特点在于：所述的该多目标检测方法包括如下步骤：

步骤一、利用线性调频信号雷达，对P个目标的回波信号进行数据采集，得到的第n次多目标回波信号可表示为：

$s(m,n)=\underset{p=1}{\overset{P}{\mathrm{\Σ}}}\left[\exp \left(j2\mathrm{\π}{f}_c(m-{\mathrm{\τ}}_{\mathrm{np}})\right)\exp \left(\mathrm{j\π\μ}{(m-{\mathrm{\τ}}_{\mathrm{np}})}^2\right)\right]---\left(1\right)$

(1)式中m为快时间域采样单元，即目标回波时延单元，n为慢时间域积累时间单元，即积累脉冲数，p为目标个数，f_c为载波频率，单位为赫兹，μ为发射信号调频率，单位为赫兹/秒，τ_np为第n次回波中第p个目标相对于雷达的延迟时间，单位为秒，表示为：

${\mathrm{\τ}}_{\mathrm{np}}=\frac{2[R_{0\mathrm{np}}-v_{0\mathrm{np}}\mathrm{nT}-0.5a_{0\mathrm{np}}{\left(\mathrm{nT}\right)}^2]}{c}---\left(2\right)$

(2)式中R_0np为第n次回波中第p个目标相对于雷达的初始距离，单位为米，v_0np为第n次回波中第p个目标相对于雷达的初始径向速度，单位为米/秒，a_0np为第n次回波中第p个目标相对于雷达的初始径向加速度，单位为米/秒²，这个加速度是由目标机动引起，或是由目标相对雷达视线夹角的变化引起，T为线性调频信号脉冲重复时间，单位为秒，c为电磁波传播速度，等于3×10⁸米/秒。该步骤完成了对雷达回波信号的离散采集，方便后续的数字化处理。

步骤二、(1)式经过下变频和脉压处理后组成的时延-积累时间单元矩阵S_nm可表示为：

$S_{\mathrm{nm}}=\left[\begin{array}{cccc}{Y}_{11}& Y_{12}& ...& Y_{1m}\\ Y_{21}& Y_{22}& ...& Y_{2m}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ Y_{n1}& Y_{n2}& ...& Y_{\mathrm{nm}}\end{array}\right],$ m＝1，2，…M，n＝1，2，…N    (3)

(3)式中T_Ynm代表第n次目标回波经脉压处理后在第m个时延单元的值，单位为秒，即：

$Y_{\mathrm{nm}}=\sqrt{D}\underset{p=1}{\overset{P}{\mathrm{\Σ}}}\exp \left(\mathrm{j\π}{f}_{\mathrm{dnp}}\frac{2R_{0\mathrm{np}}}{c}\right)\sin c\left[\mathrm{\πB}(m-\frac{2R_{0\mathrm{np}}}{c}+\frac{2v_{0\mathrm{np}}}{c}\mathrm{nT}+\frac{f_{\mathrm{dnp}}}{\mathrm{\μ}})\right]\exp \left(\mathrm{j\π}{f}_{\mathrm{dnp}}m\right)---\left(4\right)$

$\·\exp \left(j2\mathrm{\π}{f}_{\mathrm{dnp}}\mathrm{nT}\right)\exp \left[\mathrm{j\π}{k}_{\mathrm{np}}{\left(\mathrm{nT}\right)}^2\right]$

(4)式中D＝BT₀为时宽带宽积，单位为赫兹·秒，B为线性调频信号带宽，近似于接收机带宽，单位为赫兹，T₀为线性调频信号脉冲宽度，单位为秒，为第n次回波中第p个目标速度引起的多普勒频率，单位为赫兹，为第n次回波中第p个目标加速度引起的调频率，单位为赫兹/秒。该步骤完成了雷达回波信号的脉压处理，提高了信号的信噪比，有利于目标在低信噪比环境中的检测。

步骤三、对矩阵S_nm内元素Y_nm按列进行短时傅里叶变换，然后对其绝对值进行恒虚警处理，虚警门限设置为Th₁，单位为相对值，无量纲，这里Th₁取低值，以保证微弱目标信号的信息保留，门限Th₁的取值满足虚警概率为10^-4～10^-3；该步骤不但可以检测匀变速的目标，也可以检测非匀变速的目标，该步骤具体过程分为以下2步：

(a)矩阵S_nm中m列元素分别作短时傅里叶变换，变换结果f_nm可表示为：

$f_{\mathrm{nm}}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{Y}_{\mathrm{nm}}g(n-l)\exp (-j2\mathrm{\πfnT})---\left(5\right)$

该步骤可以提高信号在所分析信号时段内的信噪比，利于弱目标信号的检测。(5)式中g(n)为高斯窗函数，l是窗函数滑动步进，这里取值为窗宽度的四分之一或二分之一。

(b)对(5)式进行恒虚警处理的结果f′_nm可表示为：

该步骤经过归一化的处理后，可克服强目标信号旁瓣对弱目标信号主瓣的压制现象。

步骤四、利用短时傅里叶变换的逆变换性质还原(6)式的相位信息，(6)式还原后的时域形式S′_nm可表示为：

S′_nm＝ISTFT{f′_nm·angle[S_nm]}    (7)

(7)式中ISTFT表示为短时傅里叶变换的逆变换，即可恢复原始信号的相位信息，有利于步骤五的相参处理。

步骤五、对步骤四结果进行归一化的分数阶傅里叶变换，其表达式为：

F′_α(u)＝F_α(S′_nm/max|S′_nm|)    (8)

(8)式中α为分数阶域变换阶次。进行分数阶傅里叶变换有助于对加速或减速目标信号进行相参处理，可以提高被检测目标信号的信噪比，利于目标检测和参数估计。

步骤六、对步骤五结果进行二维搜索，把超过门限Th₂的点作为目标信号，记输出结果为：

(9)式中i为满足(9)式的目标数目，门限Th₂的单位为相对值，无量纲，其取值满足虚警概率为10^-6～10^-5。进行二维搜索后，可以把满足(9)式的所有目标点都找到，可提高目标信号的检测概率。

步骤七、为了减少脉冲噪声或强信号对弱信号造成的虚警影响，对(9)式获得的i个点目标信号采用图像对比度方法来进行评价，图像对比度定义为：

$D_i=\frac{1}{Q}\mathrm{\Σ}{\left|\right[F^{\′\′}({\stackrel{\‾}{u}}_i,{\stackrel{\‾}{\mathrm{\α}}}_i),(a,b)\left]\right|}^2-{\left|\frac{1}{Q}\mathrm{\Σ}\right[F^{\′\′}({\stackrel{\‾}{u}}_i,{\stackrel{\‾}{\mathrm{\α}}}_i),(a,b)\left]\right|}^2---\left(10\right)$

(10)式中D_i表示图像对比度，值越大，代表变换阶次估计值越接近真实值；(a，b)代表空间搜索区域的范围，Q为(a，b)区域内总的点数。采用所述的图像对比度处理利用了目标点周围的空间信息，有助于提高检测概率，降低虚警概率。

步骤八、对步骤七中获得的图像对比度进行降序排列，对超过门限Th₃的目标点进行标记，并把对应的目标点

记录下来，门限Th₃的单位为相对值，无量纲，其取值满足虚警概率为10^-6～10^-5。该步骤可降低目标的虚警概率，提高检测概率。

步骤九、对步骤八中记录下来的点信号进行参数估计，其表达式如下：

${\hat{f}}_i={\stackrel{\‾}{u}}_i\csc {\stackrel{\‾}{\mathrm{\α}}}_i,$ ${\hat{k}}_i=-\cot {\stackrel{\‾}{\mathrm{\α}}}_i---\left(11\right)$

(11)式中

对应目标速度引起的多普勒频率，单位为赫兹，

对应目标加速度引起的调频率，单位为赫兹/秒。该步骤处理后可直接获得目标的运动参数信息，可方便指挥控制员对目标的识别和后续处理。

步骤十、采用逐次消去法去除步骤八中记录的所有点，然后重复步骤三到步骤九，直到所有目标点在分数阶傅里叶变换域的峰值低于预设门限为止，这里选择本次目标峰值的50％为下一次的预设门限值。该步骤处理后有利于对不同加速或减速的目标信号进行检测。

步骤十一、为了降低虚警概率，对获得的所有目标参数进行凝聚处理，凝聚处理后的目标个数和参数即该多目标检测方法检测的多目标个数和参数。所述的该步骤可以进一步提高目标的信杂比，降低虚警概率。

根据以上所述的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法，技术特点还有：所述的该多目标检测方法中：a.所述的步骤三(a)中高斯窗函数的选择：高斯函数中窗越宽，时间分辨率越低；窗越窄，时间分辨率会提高，但频率分辨率又会降低；本发明窗函数选择高斯窗函数表示为：

g(n)＝exp[-π(nt_s)²B_x/T_x]    (12)

(12)式中t_s为信号采样间隔，单位为秒，T_x为被测信号x(nt_s)对应的时域支撑宽度，单位为秒，B_x为被测信号x(nt_s)对应的频域支撑宽度，单位为赫兹，且

T_x，B_x分别为定义为：

$T_x=\frac{\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{\left[{({\mathrm{nt}}_s-{\mathrm{\η}}_n)}^2{\left|x\left(n{t}_s\right)\right|}^2\right]}^{1/2}}{\left|\right|x\left|\right|},$ $B_x=\frac{\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{\left[{(\mathrm{n\Δf}-{\mathrm{\η}}_f)}^2{\left|X\left(\mathrm{n\Δf}\right)\right|}^2\right]}^{1/2}}{\left|\right|X\left|\right|}---\left({13}^{\′}\right)$

这里

x(nt_s)和X(nΔf)是一对傅里叶变换对，Δf为傅里叶变换域的频率分辨率，单位为赫兹，||x||为x(nt_s)的范数，||X||为X(nΔf)的范数，但须注意本发明中用到的信号x(nt_s)是指其公式(3)中的离散信号S_nm。采用高斯函数可以获得目标信号的最高时频分辨率，利于多目标信号的分离。b.所述的步骤三(a)中窗宽度BW的选择：由(12)式知高斯函数最大值下降0.707倍对应的窗宽度为：

$\mathrm{BW}\≈0.5\sqrt{\frac{T_x}{\mathrm{\π}{B}_x}}---\left(14\right)$

从(14)式知窗的宽度与T_x/B_x成正比，T_x/B_x值越大，BW越大，反之亦然，T_x和B_x的取值依据(13)式中被测信号的能量谱密度，功率谱密度及信号在时频域的幅度决定；(13)式中η_n为被测信号所占时间在整个测试时间内的平均分布值，单位为秒，η_f为被测信号所占频率在整个频带上的平均分布值，单位为赫兹，|x(nt_s)|表示信号在时域的幅度，单位为伏特，|X(nΔf)|表示信号在傅里叶变换域的幅度，单位为伏特。该步骤可以根据输入信号的特征，自适应确定高斯函数的窗长，提高了算法的灵活性。c.所述的步骤三、步骤六和步骤八中门限的选择：多次对门限阈值进行选择，以步骤三门限选择为例，这里假设噪声服从高斯分布，其模值服从瑞利(Rayleigh)分布，则有：

$f\left(r\right)=\left\{\begin{array}{cc}\frac{r}{{\mathrm{\&delta;}}^2}\exp (-\frac{r^2}{2{\mathrm{\&delta;}}^2})& r\&GreaterEqual;0\\ 0& r<0\end{array}\right.---\left(15\right)$

(15)式中δ与噪声样本均值mean的关系为：

$\mathrm{\&delta;}=\sqrt{\frac{2}{\mathrm{\&pi;}}}\&CenterDot;\mathrm{mean}---\left(16\right)$

(16)式中δ表示噪声的方差，mean表示噪声的均值。另外，(15)式对应的恒虚警概率可表示为：

$P_{\mathrm{fa}}=\underset{\mathrm{Th}}{\overset{\&infin;}{\&Integral;}}f\left(r\right)\mathrm{dr}=\exp (-\frac{{\mathrm{Th}}^2}{2{\mathrm{\&delta;}}^2})---\left(17\right)$

(17)式中检测门限Th与虚警概率P_fa的关系为：

$\mathrm{Th}=\mathrm{\&delta;}\sqrt{2\ln (1/P_{\mathrm{fa}})}---\left(18\right)$

根据(18)式，步骤三中检测门限Th₁选用的P_fa取值范围10^-4～10^-3，步骤六中检测门限Th₂选用的P_fa取值范围10^-6～10^-5，步骤八中检测门限Th₃选用的P_fa取值范围10^-6～10^-5。

根据以上所述的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法，技术特点还有：所述的该多目标检测方法中关于分数阶傅里叶变换的变换阶次选择方法是：利用初步搜索和精细搜索相互配合完成：首先利用目标部分先验信息对回波信号中线性调频信号的调频率进行分析和估计。其次，在所估计的调频率范围内进行分数阶傅里叶变换运算，分数阶傅里叶变换阶次的搜索范围为[0，2]，步进为Δ，选择为0.1或0.01；接着在分数阶傅里叶变换域的二维平面内搜索最大峰值，利用峰值坐标得到线性调频信号调频率参数对应的分数阶阶次初估值

接着以该参数为中心，搜索范围为

步进为0.1Δ进行下一次估计，依次循环缩小搜索范围使调频率估计精度逐渐增加，同时设置分数阶阶次顺序估计偏差门限，偏差门限选择为0.1Δ，当两次顺序估计偏差低于该门限就停止循环，得到线性调频信号调频率的最优估计解。所述的该步骤可以提高分数阶傅里叶变换算法的运算速度，增强了工程实用性。

根据以上所述的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法，技术特点还有：所述的该多目标检测方法中：a.所述的短时傅里叶变换的短时窗宽度选择原则：首先将待检测的信号代入(13)式，接着(13)式根据被测信号的功率谱密度、幅度计算时域支撑宽度T_x、频域支撑宽度B_x。其次，将计算获得的T_x，B_x代入(14)式计算短时傅里叶变换所需的窗宽度。b.所述的步骤三中窗函数滑动步进l选择取值范围为窗宽度的四分之一到二分之一。c.所述的步骤五中分数阶傅里叶变换的变换阶次α选择取值范围：首先对被测信号在[0，2]的阶次范围内以Δ为0.1或0.01的步进进行初步检测，接着以初步检测获得的结果为中心，依次循环缩小搜索范围，同时设置分数阶阶次顺序估计偏差门限，偏差门限选择为0.1Δ，当两次顺序估计偏差低于该门限就停止循环。d.所述的步骤六中对步骤五结果以分数阶傅里叶域u和变换阶次域α为二维域进行幅值搜索，把幅值超过预设门限的点作为目标信号记录下来，其中门限的设置要符合公式(18)，这里设虚警概率P_fa的范围为10^-6～10^-5。e.所述的步骤七中(10)式的图像对比度D_i选择取值范围为0.6～1。f.所述的步骤九中对步骤八中记录下来的目标点

带入(11)式进行多普勒频率

的估计以及调频率

的估计，然后将代入

中求得目标速度v_i，

代入

中求得目标加速度a_i。

本发明的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法优点有：1.本发明针对线性调频窄带雷达在角度维，距离维无法分辨多个目标的问题，提出通过多普勒频率维对多目标进行检测，通过采用分数阶傅里叶变换的方法提高了多普勒频率的分辨率，从频率维解决了多目标的分离问题；2.本发明还提出首先利用STFT进行信号初步检测，提取可能的目标信号，然后采用二值化方法对其进行处理，这里与图像处理中二值化方法不同的是要保留信号的相位信息，然后利用STFT的可逆性进行信号还原，优点是可以克服强信号旁瓣对弱信号主瓣的压制现象，解决了强弱信号分离较困难的问题；3.本发明在FRFT基础上，对还原后的信号利用分数阶核函数角度旋转的特点，在最优旋转角时得到信号的最佳能量集聚，优点是提高了待检测信号的信噪比，解决了低信噪比下检测信号虚警概率较大的问题；4.本发明同时配合图像对比度方法和逐次消去方法进行逐次滤波检测，优点是利用信号的空间信息，可提高信号的检测概率，降低了噪声对目标信号的影响。5.本发明采用多种方法进行联合处理，不仅能够检测调频率不同的多个强弱信号，而且也能检测调频率相同的多个强弱信号，可进一步提高检测概率和计算效率，使计算量减小，易于工程实现。这种基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法值得采用和推广。

四.附图说明

本发明的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法说明书附图共有十幅：

图1为基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法总体框图；

图2为经短时傅里叶变换结果分析示意图；

图3为经分数阶傅里叶变换结果分析示意图；

图4为本多目标检测方法中经步骤三结果分析示意图；

图5为本多目标检测方法中经步骤六结果分析示意图；

图6为本多目标检测方法中经步骤七结果分析示意图；

图7为本多目标检测方法的理论分析结果示意图；

图8为多目标检测方法中经步骤十结果分析示意图；

图9为本多目标检测方法中目标检测概率与输入幅值比的变化关系示意图；

图10为本多目标检测方法中目标检测概率与加速度差值的变化关系示意图。

在各图中采用了统一标号，即同一物件在各图中用同一标号。在各图中：1.多目标信号采集；2.短时傅里叶变换(STFT)；3.短时傅里叶域二维矩阵；4.短时傅里叶域虚警门限Th₁；5.归一化后保留相位信息的二维矩阵；6.短时傅里叶变换的逆变换(ISTFT)；7.原始信号恢复；8.分数阶傅里叶变换(FRFT)；9.分数阶傅里叶域二维矩阵；10.分数阶傅里叶域虚警门限Th₂；11.图像对比度统计；12.图像对比度虚警门限Th₃；13.目标检测及参数估计；14.逐步消去法；15.目标峰值判决；16.算法结束；17.Z轴：归一化幅度，无量纲；18.傅里叶域多目标检测结果；19.傅里叶域扩散信号和噪声；20.X轴：傅里叶域，量纲：赫兹；21.分数阶傅里叶域多目标检测结果，峰值处为目标2检测结果；22.分数阶傅里叶域扩散信号和噪声；23.X轴：分数阶傅里叶域，量纲：赫兹；24.Y轴：变换阶次域，无量纲；25.短时傅里叶域多目标处理结果，峰值处为能量归一化后的信号；26.短时傅里叶域扩散信号和噪声；27.短时傅里叶域时间坐标轴，量纲：秒；28.分数阶傅里叶域多目标检测结果，峰值处为目标1检测结果；29.图像对比度方法处理后的多目标检测结果，峰值处为目标2检测结果；30.图像对比度方法处理后的多目标检测结果，峰值处为目标1检测结果；31.理论分析后的多目标检测结果，峰值处为目标2检测结果；32.理论分析后的多目标检测结果，峰值处为目标1检测结果；33.逐次消去法处理后的多目标检测结果，峰值处为目标3检测结果；34.Y轴：目标检测概率，无量纲；35.X轴：两目标幅值比，量纲：分贝；36.本发明方法处理后的检测概率曲线；37.直接分数阶傅里叶变换处理后的检测概率曲线；38.X轴：两目标加速度差，量纲：米/秒²。

五.具体实施方式

本发明的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法非限定实施例如下：

实施例一.基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法

该例的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法具体情况由图1～图8联合示出。该例的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法具体实施步骤如下：步骤一，是多目标信号采集，针对LFM脉冲体制雷达，采用线性调频信号雷达，对P个目标的回波信号进行数据采集，得到的第n次多目标回波信号可表示为：

$s(m,n)=\underset{p=1}{\overset{P}{\mathrm{\&Sigma;}}}\left[\exp \left(j2\mathrm{\&pi;}{f}_c(m-{\mathrm{\&tau;}}_{\mathrm{np}})\right)\exp \left(\mathrm{j\&pi;\&mu;}{(m-{\mathrm{\&tau;}}_{\mathrm{np}})}^2\right)\right]---\left(1\right)$

(1)式中m为快时间域采样单元，即目标回波时延单元，n为慢时间域积累时间单元，即积累脉冲数，p为目标个数，f_c为载波频率，单位为赫兹，μ为发射信号调频率，单位为赫兹/秒，τ_np为第n次回波中第p个目标相对于雷达的延迟时间，单位为秒，表示为：

${\mathrm{\&tau;}}_{\mathrm{np}}=\frac{2[R_{0\mathrm{np}}-v_{0\mathrm{np}}\mathrm{nT}-0.5a_{0\mathrm{np}}{\left(\mathrm{nT}\right)}^2]}{c}---\left(2\right)$

(2)式中R_0np为第n次回波中第p个目标相对于雷达的初始距离，单位为米，v_0np为第n次回波中第p个目标相对于雷达的初始径向速度，单位为米/秒，a_0np为第n次回波中第p个目标相对于雷达的初始径向加速度，单位为米/秒²，这个加速度是由目标机动引起，或是由目标相对雷达视线夹角的变化引起，T为线性调频信号脉冲重复时间，单位为秒，c为电磁波传播速度，等于3×10⁸米/秒。该例中采用的数据分别为：雷达工作在S波段，最大探测距离为15km，信号带宽为2MHz，脉冲重复周期为100μs，脉冲宽度为10μs，积累时间约为800ms，采样频率为2MHz。该例假设有三个加速运动的编队目标被探测，目标1初始距离为15km，速度为100m/s，加速度为20m/s²；目标2初始距离为15.1km，速度为150m/s，加速度为20m/s²；目标3初始距离为14.9km，速度为150m/s，加速度为15m/s²。其中目标1和目标2的加速度相同，速度不同；目标2和目标3的速度相同，加速度不同。所设目标1信噪比为-17dB，目标3信噪比为-5dB，目标2信噪比为-2dB。将以上数据按公式(1)、(2)计算该例的三个目标的第n次回波信号s(m，n)和三个目标第n次回波相对于雷达的延迟时间τ_np。这里需注意本发明中雷达的发射信号为窄带LFM信号，其距离分辨率为150米，大于三目标相互间的距离间隔(三目标间的最大距离间隔为100米)，故无法从距离维上分辨。下面应用多普勒频率维进行分辨，通过步骤二：(1)式经过下变频和脉压处理后组成的时延-积累时间单元矩阵S_nm可表示为：

$S_{\mathrm{nm}}=\left[\begin{array}{cccc}{Y}_{11}& Y_{12}& ...& Y_{1m}\\ Y_{21}& Y_{22}& ...& Y_{2m}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ Y_{n1}& Y_{n2}& ...& Y_{\mathrm{nm}}\end{array}\right],$ m＝1，2，…M，n＝1，2，…N    (3)

(3)式中Y_nm代表第n次目标回波经脉压处理后在第m个时延单元的值，单位为秒，即：

$Y_{\mathrm{nm}}=\sqrt{D}\underset{p=1}{\overset{P}{\mathrm{\&Sigma;}}}\exp \left(\mathrm{j\&pi;}{f}_{\mathrm{dnp}}\frac{2R_{0\mathrm{np}}}{c}\right)\sin c\left[\mathrm{\&pi;B}(m-\frac{2R_{0\mathrm{np}}}{c}+\frac{2v_{0\mathrm{np}}}{c}\mathrm{nT}+\frac{f_{\mathrm{dnp}}}{\mathrm{\&mu;}})\right]\exp \left(\mathrm{j\&pi;}{f}_{\mathrm{dnp}}m\right)---\left(4\right)$

$\&CenterDot;\exp \left(j2\mathrm{\&pi;}{f}_{\mathrm{dnp}}\mathrm{nT}\right)\exp \left[\mathrm{j\&pi;}{k}_{\mathrm{np}}{\left(\mathrm{nT}\right)}^2\right]$

(4)式中D＝BT₀为时宽带宽积，单位为赫兹·秒，B为线性调频信号带宽，近似于接收机带宽，单位为赫兹，T₀为线性调频信号脉冲宽度，单位为秒，

为第n次回波中第p个目标速度引起的多普勒频率，单位为赫兹，

为第n次回波中第p个目标加速度引起的调频率，单位为赫兹/秒。按该例的三个目标的各数据代入公式(3)、(4)分别进行计算，得到S_nm时延-积累时间单元矩阵的大小为100×8000，其中为了利用傅里叶的快速计算方法，这里积累时间单元取为8192，即最后的矩阵大小为100×8192，另外，Y_nm表示了8192次回波信号的脉压结果，每次脉压针对100个时延单元进行处理，其获得时宽带宽积为20。步骤三：对矩阵S_nm内元素Y_nm按列进行短时傅里叶变换(STFT)，然后对其绝对值进行恒虚警处理，虚警门限设置为Th₁，单位为相对值，无量纲，这里Th₁取低值，以保证微弱目标信号的信息保留，门限Th₁的取值满足虚警概率为10^-6或10^-5。该步骤具体过程分为以下2步：(a)矩阵S_nm中m列元素的短时傅里叶变换变换结果f_nm可表示为：

$f_{\mathrm{nm}}=\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{Y}_{\mathrm{nm}}g(n-l)\exp (-j2\mathrm{\&pi;fnT})---\left(5\right)$

(5)式中g(n)为高斯窗函数，l是窗函数滑动步进，该例的步骤三中窗函数滑动步进l选择取值范围为窗宽度的四分之一到二分之一，该例里取值为窗宽度的二分之一。该例的步骤三(a)中高斯窗函数的选择：高斯函数中窗越宽，时间分辨率越低；窗越窄，时间分辨率会提高，但频率分辨率又会降低；本发明选择的高斯窗函数表示为：

g(n)＝exp[-π(nt_s)²B_x/T_x]    (12)

(12)式中t_s为信号采样间隔，单位为秒，T_x为被测信号x(nt_s)对应的时域支撑宽度，单位为秒，B_x为被测信号x(nt_s)对应的频域支撑宽度，单位为赫兹，且

T_x，B_x分别为定义为：

$T_x=\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{\left[{({\mathrm{nt}}_s-{\mathrm{\&eta;}}_n)}^2{\left|x\left(n{t}_s\right)\right|}^2\right]}^{1/2}}{\left|\right|x\left|\right|},$ $B_x=\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{\left[{(\mathrm{n\&Delta;f}-{\mathrm{\&eta;}}_f)}^2{\left|X\left(\mathrm{n\&Delta;f}\right)\right|}^2\right]}^{1/2}}{\left|\right|X\left|\right|}---\left(13\right)$

这里

x(nt_s)和X(nΔf)是一对傅里叶变换对，Δf为傅里叶变换域的频率分辨率，单位为赫兹，||x||为x(nt_s)的范数，||X||为X(nΔf)的范数，但须注意本发明中用到的信号x(nt_s)是指其公式(3)中第m个时延单元中n次目标回波经脉压处理后的离散信号S_nm。该例的步骤三(a)中窗宽度BW的选择：由(12)式知高斯函数最大值下降0.707倍对应的窗宽度为：

$\mathrm{BW}\&ap;0.5\sqrt{\frac{T_x}{\mathrm{\&pi;}{B}_x}}---\left(14\right)$

从(14)式知窗的宽度与T_x/B_x成正比，T_x/B_x值越大，BW越大，反之亦然，T_x和B_x的取值依据(13)式中被测信号的能量谱密度，功率谱密度及信号在时频域的幅度决定；(13)式中η_n为被测信号所占时间在整个测试时间内的平均分布值，单位为秒，η_f为被测信号所占频率在整个频带上的平均分布值，单位为赫兹，|x(nt_s)|表示信号在时域的幅度，单位为伏特，|X(nΔf)|表示信号在傅里叶变换域的幅度，单位为伏特。该例的短时傅里叶变换的短时窗宽度选择原则：首先将待检测的信号代入(13)式，接着(13)式根据被测信号的功率谱密度、幅度计算时域支撑宽度T_x、频域支撑宽度B_x。其次，将计算获得的T_x，B_x代入

(14)式计算短时傅里叶变换所需的窗宽度。该例的步骤三中门限的选择为：这里假设噪声服从高斯分布，其模值服从瑞利(Rayleigh)分布，则有：

$f\left(r\right)=\left\{\begin{array}{cc}\frac{r}{{\mathrm{\&delta;}}^2}\exp (-\frac{r^2}{2{\mathrm{\&delta;}}^2})& r\&GreaterEqual;0\\ 0& r<0\end{array}\right.---\left(15\right)$

(15)式中δ与噪声样本均值mean的关系为：

$\mathrm{\&delta;}=\sqrt{\frac{2}{\mathrm{\&pi;}}}\&CenterDot;\mathrm{mean}---\left(16\right)$

(16)式中δ表示噪声的方差，mean表示噪声的均值。另外，(15)式对应的恒虚警概率可表示为：

$P_{\mathrm{fa}}=\underset{\mathrm{Th}}{\overset{\&infin;}{\&Integral;}}f\left(r\right)\mathrm{dr}=\exp (-\frac{{\mathrm{Th}}^2}{2{\mathrm{\&delta;}}^2})---\left(17\right)$

(17)式中检测门限Th与虚警概率P_fa的关系为：

$\mathrm{Th}=\mathrm{\&delta;}\sqrt{2\ln (1/P_{\mathrm{fa}})}---\left(18\right)$

根据(18)式，步骤三中检测门限Th₁选用的P_fa取值范围10^-4～10^-3。

(b)对(5)式进行恒虚警处理的结果f′_nm可表示为：

该例的步骤四.利用短时傅里叶变换的逆变换性质还原(6)式的相位信息，(6)式还原后的时域形式S′_nm可表示为：

S′_nm＝ISTFT{f′_nm·angle[S_nm]}    (7)

(7)式中ISTFT表示为短时傅里叶变换的逆变换。按该例的三目标的各数据代入公式(5)计算矩阵S_nm中m列元素的短时傅里叶变换变换结果f_nm，根据(12)式选择高斯窗函数、根据(14)式选择窗宽度、窗函数滑动步进值，其中STFT窗宽度经式(14)计算为430个采样单元，为了快速傅里叶变换计算方便，这里选取点数为512点，并且窗移动步长选为窗宽的二分之一，即256点。根据(15)-(18)式选择步骤三的检测门限Th₁值为307和虚警概率P_fa值为10^-3。按(6)式计算(7)式短时傅里叶变换的逆变换S′_nm值。该例的步骤五.对步骤四结果进行归一化的分数阶傅里叶变换(FRFT)，其表达式为：

F′_α(u)＝F_α(S′_nm/max|S′_nm|)    (8)

(8)式中(8)式中α为分数阶域变换阶次。该例关于分数阶傅里叶变换的变换阶次选择方法是：利用初步搜索和精细搜索相互配合完成：首先利用目标部分先验信息对回波信号中线性调频信号的调频率进行分析和估计；其次，在所估计的调频率范围内进行分数阶傅里叶变换运算，分数阶傅里叶变换阶次的搜索范围为[0，2]，步进为0.1；接着在分数阶傅里叶变换域的二维平面内搜索最大峰值，利用峰值坐标得到线性调频信号，调频率参数对应的分数阶阶次初估值接着以该参数为中心，搜索范围为

步进为0.01进行下一次估计，依次循环缩小搜索范围使调频率估计精度逐渐增加，同时设置分数阶阶次顺序估计偏差门限，偏差门限选择为0.01，当两次顺序估计偏差低于该门限就停止循环，得到线性调频信号调频率的最优估计解。该例的步骤五中分数阶傅里叶变换的变换阶次α选择取值范围：首先对被测信号在[0，2]的阶次范围内以0.01的步进进行初步检测，接着以初步检测获得的结果为中心，依次循环缩小搜索范围，同时设置分数阶阶次顺序估计偏差门限，偏差门限选择为0.01，当两次顺序估计偏差低于该门限就停止循环。按该例的三目标各参数根据(8)式在FRFT处理步骤中变换阶次以[0 2]粗搜索范围进行步进间隔为0.1的初次搜索，获得的分数阶傅里叶变换阶次为1.1，并以0.1的步进间隔进行二维平面初次搜索，然后在[1.0 1.2]范围内再采用0.01的步进再进行分数阶傅里叶变换对可疑目标信号作精细搜索。该例共进行20次分数阶傅里叶变换后，获得的最佳分数阶阶次为1.02。实施步骤六：对步骤五在分数阶傅里叶变换域的二维平面上，接着进行二维搜索，在二维平面内搜索最大峰值，利用峰值坐标得到线性调频信号，把超过门限Th₂的点作为目标信号，记输出结果为：

(9)式中i为满足(9)式的目标数目，门限Th₂的单位为相对值，无量纲，该例的步骤六中对步骤五结果以分数阶傅里叶域u和变换阶次域α为二维域进行幅值搜索，把幅值超过预设门限的点作为目标信号记录下来，其中门限的设置要符合(18)式，这里设虚警概率P_fa的范围为10^-6值。根据(15)-(18)式选择步骤六的检测门限Th₂值0.89。实施该例的步骤七：为了减少脉冲噪声或强信号对弱信号造成的虚警影响，对(9)式获得的i个点目标信号采用图像对比度方法来进行评价，图像对比度定义为：

$D_i=\frac{1}{Q}\mathrm{\&Sigma;}{\left|\right[F^{\&prime;\&prime;}({\stackrel{\&OverBar;}{u}}_i,{\stackrel{\&OverBar;}{\mathrm{\&alpha;}}}_i),(a,b)\left]\right|}^2-{\left|\frac{1}{Q}\mathrm{\&Sigma;}\right[F^{\&prime;\&prime;}({\stackrel{\&OverBar;}{u}}_i,{\stackrel{\&OverBar;}{\mathrm{\&alpha;}}}_i),(a,b)\left]\right|}^2---\left(10\right)$

(10)式中D_i表示图像对比度，值越大，代表变换阶次估计值越接近真实值；(a，b)代表空间搜索区域的范围，Q为(a，b)区域内总的点数，这里取a为10，b为5。接着实施步骤八：对步骤七中获得的图像对比度进行降序排列，对超过门限Th₃的目标点进行标记，并把对应的目标点

记录下来，门限Th₃的单位为相对值，无量纲，其取值满足虚警概率为10^-6。接着实施该例的步骤九：对步骤八中记录下来的点信号进行参数估计，其表达式如下：

${\hat{f}}_i={\stackrel{\&OverBar;}{u}}_i\csc {\stackrel{\&OverBar;}{\mathrm{\&alpha;}}}_i,$ ${\hat{k}}_i=-\cot {\stackrel{\&OverBar;}{\mathrm{\&alpha;}}}_i---\left(11\right)$

(11)式中

对应目标速度引起的多普勒频率，单位为赫兹，

对应目标加速度引起的调频率，单位为赫兹/秒。按该例的步骤七中(10)式的图像对比度D_i选择取值范围为0.9。对通过(9)式获得该例的三个点目标信号采用图像对比度方法来进行评价，根据(15)-(18)式选择步骤八的检测门限Th₃值0.87。按步骤七中获得的图像对比度进行降序排列，对超过门限Th₃的目标点进行标记，并把对应的目标点

记录下来，接着对步骤八中记录下来的点信号进行参数估计：如目标距离、目标速度、目标加速度等参数进行估计。该例的步骤九中对经过步骤八处理后记录下来的目标点

带入(11)式进行多普勒频率

的估计以及调频率

的估计，然后将

代入

中求得目标速度v_i，

代入

中求得目标加速度a_i。利用式(11)可求的目标1的初始频率为2010Hz，调频率为420Hz/s，目标2的初始频率为3030Hz，调频率为420Hz/s；继续实施该例的步骤十：采用逐次消去法去除步骤八中记录的所有点，具体做法是重复步骤三到步骤九，直到所有目标点在分数阶傅里叶变换域的峰值低于预设门限为止，这里选择本次目标峰值的50％为下一次的预设门限值。当步骤十目标峰值判决，判断目标峰值不大于预设值时结束。实施步骤十一：为了降低虚警概率，对获得的所有目标参数进行凝聚处理，凝聚处理后的目标个数和参数即该多目标检测方法检测的多目标个数和参数。该例的方法检测的多目标个数为3，目标1的速度为101米/秒，加速度为21米/秒²，目标2的速度为152米/秒，加速度为21米/秒²，目标3的速度为152米/秒，加速度为16米/秒²。该例的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法具体按照图1的总体实现框图给出重要步骤的仿真示意结果。图1的总体实现框图为：1为多目标信号采集，2为短时傅里叶变换(STFT)，3为短时傅里叶域二维矩阵，4为短时傅里叶域虚警门限Th₁，5为归一化后保留相位信息的二维矩阵，6为短时傅里叶变换的逆变换(ISTFT)，7为原始信号恢复，8为分数阶傅里叶变换(FRFT)，9为分数阶傅里叶域二维矩阵，10为分数阶傅里叶域虚警门限Th₁，11为图像对比度统计，12为图像对比度虚警门限Th₃，13为目标检测及参数估计，14为逐步消去法，15为目标峰值判决，16为算法结束。图2，图3分别为直接进行短时傅里叶变换与分数阶傅里叶变换的分析结果，可以明显看出图中只显示有一个目标，目标峰值信号对应的分数阶阶次为1.02，其它目标看不到。图4为采用短时傅里叶变换的分析结果，此时门限为307，虚警概率为10-3，虚警较多，不易检测到真正目标的个数。图5为采用步骤六方法进行的分数阶傅里叶变换，这时可以清楚的看到有两个目标，并且由于采用了归一化方法，使得两个目标的幅值基本一致且其峰值对应的分数阶阶次均为1.02，故本方法能够克服调频率相同时强目标信号对弱目标信号干扰的影响。图6为采用图像对比度方法(对比度门限为0.9，虚警概率为10^-6)获得的分析结果，可以看出其二维平面上的能量分布比图5均匀且平滑度要好于图5，易于实现目标的检测。图7为了检测图六中结果是否正确，给出了目标1和目标2信噪比相同时的直接分数阶傅里叶变换结果，对比图6可以看出，两者结果一致，证明了本发明所提出的方法正确性。图8为采用逐次消去法方法后对目标3的检测结果，其峰值对应的分数阶阶次为1.01，与所设值相同。从以上仿真实例图示可以看出，本发明的方法可以实现多个强弱目标信号在调频率相同和调频率不同时的检测。通过该例验证本发明方法检测多个加速目标的有效性。

实施例二.基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法

该例的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法具体情况可用类似图1～图8等联合示出，该例的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法的检测步骤如实施例一中所述的共十一步骤，其检测过程也如图1中所述的共16过程，不再重述。雷达系统为：雷达工作在S波段，最大探测距离为15km，信号带宽为1MHz，脉冲重复周期为100μs，脉冲宽度为10μs，积累时间约为800ms，采样频率为2MHz。该例假设有三个加速运动的编队目标被探测，目标1初始距离为15km，速度为100m/s，加速度为20m/s²；目标2初始距离为15.1km，速度为150m/s，加速度为20m/s²；目标3初始距离为14.9km，速度为150m/s，加速度为15m/s²。其中目标1和目标2的加速度相同，速度不同；目标2和目标3的速度相同，加速度不同。所设信噪比目标1最小，目标3次之，目标2最大。该例的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法与实施例一的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法不同点有：1.计算短时傅里叶变换所需的窗宽度是256。2.所述的步骤三中窗函数滑动步进l选择取值范围为128。3.所述的步骤三中门限的选择：门限的设置要符合公式(18)，门限Th₁值是154，这里设虚警概率P_fa的范围为10^-3。4.所述的步骤五中分数阶傅里叶变换的变换阶次α选择的粗搜索范围为[0 2]，步进为0.1，选择的精搜索范围为[1.0 1.2]，步进为0.01，最后获得的目标1，目标2的分数阶阶次为1.05，目标3的分数阶阶次为1.03。5.所述的步骤六中门限的选择：门限的设置要符合公式(18)，门限Th₂值0.86，这里设虚警概率P_fa的范围为10^-6。6.所述的步骤七中(10)式的图像对比度D_i选择取值范围为0.91。7.所述的步骤八中门限的选择门限的设置要符合公式(18)，门限Th₃值0.89，这里设虚警概率P_fa的范围为10^-6。8.然后将代入

中求得目标速度v_i，

代入

中求得目标加速度a_i即：目标1的速度为101.8米/秒，加速度为21.7米/秒²，目标2的速度为152.4米/秒，加速度为21.7米/秒²，目标3的速度为152.4米/秒，加速度为16.7米/秒²。对图2～图8的显示结果，表明了本发明所述方法的正确性和有效性。该例的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法其余未述的，全同于实施例一中所述的，不再重述。

实施例三.基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法

该例的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法具体情况可用类似图1～图8等联合示出，该例的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法的检测步骤如实施例一中所述的共十一步骤，其检测过程也如图1中所述的共16过程，不再重述。雷达系统为：雷达工作在S波段，最大探测距离为15km，信号带宽为1MHz，脉冲重复周期为100μs，脉冲宽度为10μs，积累时间约为800ms，采样频率为2MHz。该例假设有三个加速运动的编队目标被探测，目标1初始距离为15km，速度为100m/s，加速度为20m/s²；目标2初始距离为15.1km，速度为150m/s，加速度为20m/s²；目标3初始距离为14.9km，速度为150m/s，加速度为15m/s²。其中目标1和目标2的加速度相同，速度不同；目标2和目标3的速度相同，加速度不同。所设信噪比目标1最小，目标3次之，目标2最大。该例的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法与实施例一、实施例二的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法不同点有：1.计算短时傅里叶变换所需的窗宽度1024。2.所述的步骤三中窗函数滑动步进l选择取值范围为窗宽度的四分之一。3.所述的步骤三中门限的选择：门限的设置要符合公式(18)，门限Th₁值650，这里设虚警概率P_fa的范围为10^-4。4.所述的步骤五中分数阶傅里叶变换的变换阶次α选择的粗搜索范围为[0 2]，步进为0.01，选择的精搜索范围为[1.01 1.21]，步进为0.001，最后获得的目标1，目标2的分数阶阶次为1.018，目标3的分数阶阶次为1.008。5.所述的步骤六中门限的选择：门限的设置要符合公式(18)，门限Th₂值0.82，这里设虚警概率P_fa的范围为10^-5。6.所述的步骤七中(10)式的图像对比度D_i选择取值范围为0.8。7.所述的步骤八中门限的选择门限的设置要符合公式(18)，门限Th₃值0.85，这里设虚警概率P_fa的范围为10^-5。8.然后将

代人中求得目标速度v_i，

代入

中求得目标加速度a_i，即：目标1的速度为100.5米/秒，加速度为20.8米/秒²，目标2的速度为151.2米/秒，加速度为20.4米/秒²，目标3的速度为151.7米/秒，加速度为15.6米/秒²。该例的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法其余未述的，全同于实施例一、实施例二中所述的，不再重述。

实施例四.基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法

该例的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法具体情况可用类似图1～图8等联合示出，该例的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法的检测步骤如实施例一中所述的，其检测过程也如图1中所述的，不再重述。雷达系统为：雷达工作在S波段，最大探测距离为15km，信号带宽为1MHz，脉冲重复周期为100μs，脉冲宽度为10μs，积累时间约为800ms，采样频率为2MHz。该例假设有三个加速运动的编队目标被探测，目标1初始距离为15km，速度为100m/s，加速度为20m/s²；目标2初始距离为15.1km，速度为150m/s，加速度为20m/s²；目标3初始距离为14.9km，速度为150m/s，加速度为15m/s²。其中目标1和目标2的加速度相同，速度不同；目标2和目标3的速度相同，加速度不同。所设信噪比目标1最小，目标3次之，目标2最大。该例的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法与实施例一、实施例二、实施例三的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法不同点有：1.计算短时傅里叶变换所需的窗宽度是2048。2.所述的步骤三中窗函数滑动步进l选择取值范围为窗宽度的三分之一。3.所述的步骤三中门限的选择：门限的设置要符合公式(18)，门限Th₁值1208，这里设虚警概率P_fa的范围为0.5×10^-3。4.所述的步骤五中分数阶傅里叶变换的变换阶次α选择的粗搜索范围为[0 2]，步进为0.01，选择的精搜索范围为[1.01 1.21]，步进为0.001，最后获得的目标1，目标2的分数阶阶次为1.016，目标3的分数阶阶次为1.003。5.所述的步骤六中门限的选择：门限的设置要符合公式(18)门限Th₂值0.86，这里设虚警概率P_fa的范围为0.5×10^-6。6.所述的步骤七中(10)式的图像对比度D_i选择取值范围为0.6。7.所述的步骤八中门限的选择门限的设置要符合公式(18)，门限Th₃值0.87，这里设虚警概率P_fa的范围为0.5x10^-6。8.然后将

代入中求得目标速度v_i，

代入

中求得目标加速度a_i，即：目标1的速度为101.2米/秒，加速度为20.4米/秒²，目标2的速度为150.7米/秒，加速度为20.4米/秒²，目标3的速度为150.9米/秒，加速度为15.9米/秒²。该例的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法其余未述的，全同于实施例一、实施例二、实施例三中所述的，不再重述。

实施例五.基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法

该例的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法具体情况可用图1～图9等联合示出，该例的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法的检测步骤如实施例一中所述的，其检测过程也如图1中所述的，不再重述。本例在实施例一～实施例四的基础上，进一步验证本发明方法对多个目标幅值比或信噪比的分辨能力。该例的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法与实施例一～实施例四的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法不同点有：1.雷达工作在L波段，信号带宽为1MHz，脉冲重复周期为600μs，脉冲宽度为60μs，积累时间约为600ms，采样频率为2MHz。2.假设在编队飞行情况下有两个加速运动目标被探测，两目标相距60m，速度均为200m/s，目标1加速度为10m/s²，目标2加速度为20m/s²。3.计算短时傅里叶变换所需的窗宽度64。4.所述的步骤三中窗函数滑动步进l选择取值范围为窗宽度的二分之一。5.所述的步骤三中门限的选择：门限的设置要符合公式(18)门限Th₁值38，这里设虚警概率P_fa的范围为10^-3。6.所述的步骤五中分数阶傅里叶变换的变换阶次α选择的粗搜索范围为[02]，步进为0.01，选择的精搜索范围为[1.000 1.200]，步进为0.001，最后获得的目标1，目标2的分数阶阶次为1.010，目标2的分数阶阶次为1.021。7.所述的步骤六中门限的选择：门限的设置要符合公式(18)，门限Th₂值0.89，这里设虚警概率P_fa的范围为10^-6。8.所述的步骤七中(10)式的图像对比度D_i选择取值范围为0.9。9.所述的步骤八中门限的选择门限的设置要符合公式(18)，门限Th₃值0.84，这里设虚警概率P_fa的范围为10^-6。10.然后将

代入中求得目标速度v_i，

代入

中求得目标加速度a_i。这里以两目标为例进行分析。根据信号带宽可知目标最小距离分辨率为150m，故相距为60m的两目标无法从距离上区分，下面应用本发明方法从FRFT域分析两目标在不同幅值比情况下的目标检测性能。取两目标输入幅值比变化范围为-40dB-40dB，经过200次蒙特卡罗仿真实验，其目标检测概率(这里设虚警率为10^-6)与输入幅值比的变化关系如图九所示。设雷达正常工作时所需的检测概率必须大于80％以上，从图9可以看出，采用FRFT方法对幅值比范围约为[-20dB，20dB]的两目标可以分辨，而采用本发明方法(STFT中高斯窗的宽度为64，门限为38)的幅值比范围可扩展为[-30dB，30dB]，显然本发明方法比传统方法更适应对动态范围大的强弱目标信号进行检测。验证本发明方法对多目标在不同幅值比下的检测性能。该例的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法其余未述的，全同于实施例一～实施例四中所述的，不再重述。

实施例六.基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法

该例的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法具体情况可用图1～图10等联合示出，该例的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法的检测步骤如实施例一中所述的，其检测过程也如图1中所述的，不再重述。本例在实施例一～例五的基础上，进一步验证本发明方法对多个目标加速度的分辨能力。该例的处理结果可用图10示出，该例设值的目标参数，雷达参数与实施例五均相同，该例的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法与实施例一～实施例五的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法不同点有：1.计算短时傅里叶变换所需的窗宽度是32。2.所述的步骤三中窗函数滑动步进l选择取值范围为窗宽度的四分之一。3.所述的步骤三中门限的选择：门限的设置要符合公式(18)，门限Th₁值20，这里设虚警概率P_fa的范围为10^-3。4.所述的步骤五中两目标加速度差值变化范围为-8～8m/s²时，分数阶傅里叶变换的变换阶次α选择的粗搜索范围为[0 2]，步进为0.01，选择的精搜索范围为[0.91 1.10]，步进为0.001，最后获得的目标分数阶阶次的范围为[0.992，1.008]。5.所述的步骤六中门限的选择：门限的设置要符合公式(18)，门限Th₂值0.83，这里设虚警概率P_fa的范围为10^-6。6.所述的步骤七中(10)式的图像对比度D_i选择取值范围为1.0。7.所述的步骤八中门限的选择门限的设置要符合公式(18)，门限Th₃值0.88，这里设虚警概率P_fa的范围为10^-6。8.然后将

代入中求得目标速度v_i，

代入

中求得目标加速度a_i。其中目标信号所加的高斯噪声为-5dB，下面应用本发明方法分析两目标在不同加速度差值下的目标检测性能，取两目标加速度差值变化范围为-8～8m/s²，经过200次蒙特卡罗仿真实验，其目标检测概率与加速度差值的变化关系如图十所示。从图10可以看出，本发明方法和FRFT方法的加速度最小分辨能力基本相同，都为2.5m/s²，故该例中本方法只能对加速度差值大于最小分辨值2.5m/s²的多个目标进行分辨和检测。该例的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法其余未述的，全同于实施例一～实施例五中所述的，不再重述。

Claims (4)

1.一种基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法，特征在于：所述的该多目标检测方法包括如下步骤：

步骤一、利用线性调频信号雷达，对P个目标的回波信号进行数据采集，得到的第n次多目标回波信号可表示为：

$s(m,n)=\underset{p=1}{\overset{P}{\mathrm{\&Sigma;}}}\left[\exp \left(j2\mathrm{\&pi;}{f}_c(m-{\mathrm{\&tau;}}_{\mathrm{np}})\right)\exp \left(\mathrm{j\&pi;\&mu;}{(m-{\mathrm{\&tau;}}_{\mathrm{np}})}^2\right)\right]---\left(1\right)$

(1)式中m为快时间域采样单元，即目标回波时延单元，n为慢时间域积累时间单元，即积累脉冲数，p为目标个数，f_c为载波频率，μ为发射信号调频率，τ_np为第n次回波中第p个目标相对于雷达的延迟时间，表示为：

${\mathrm{\&tau;}}_{\mathrm{np}}=\frac{2[R_{0\mathrm{np}}-v_{0\mathrm{np}}\mathrm{nT}-0.5a_{0\mathrm{np}}{\left(\mathrm{nT}\right)}^2]}{c}---\left(2\right)$

(2)式中R_0np为第n次回波中第p个目标相对于雷达的初始距离，v_0np为第n次回波中第p个目标相对于雷达的初始径向速度，a_0np为第n次回波中第p个目标相对于雷达的初始径向加速度，T为线性调频信号脉冲重复时间，c为电磁波传播速度；

步骤二、(1)式经过下变频和脉压处理后组成的时延-积累时间单元矩阵S_nm可表示为：

$S_{\mathrm{nm}}=\left[\begin{array}{cccc}{Y}_{11}& Y_{12}& ...& Y_{1m}\\ Y_{21}& Y_{22}& ...& Y_{2m}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ Y_{n1}& Y_{n2}& ...& Y_{\mathrm{nm}}\end{array}\right],$ m＝1，2，…M，n＝1，2，…N    (3)

(3)式中Y_nm代表第n次目标回波经脉压处理后在第m个时延单元的值：

$Y_{\mathrm{nm}}=\sqrt{D}\underset{p=1}{\overset{P}{\mathrm{\&Sigma;}}}\exp \left(\mathrm{j\&pi;}{f}_{\mathrm{dnp}}\frac{2R_{0\mathrm{np}}}{c}\right)\sin c\left[\mathrm{\&pi;B}(m-\frac{2R_{0\mathrm{np}}}{c}+\frac{2v_{0\mathrm{np}}}{c}\mathrm{nT}+\frac{f_{\mathrm{dnp}}}{\mathrm{\&mu;}})\right]\exp \left(\mathrm{j\&pi;}{f}_{\mathrm{dnp}}m\right)---\left(4\right)$

$\&CenterDot;\exp \left(j2\mathrm{\&pi;}{f}_{\mathrm{dnp}}\mathrm{nT}\right)\exp \left[\mathrm{j\&pi;}{k}_{\mathrm{np}}{\left(\mathrm{nT}\right)}^2\right]$

(4)式中D＝BT₀为时宽带宽积，B为线性调频信号带宽，T₀为线性调频信号脉冲宽度，

为第n次回波中第p个目标速度引起的多普勒频率，

为第n次回波中第p个目标加速度引起的调频率；

步骤三、对矩阵S_nm内元素Y_nm按列进行短时傅里叶变换，然后对其绝对值进行恒虚警处理，虚警门限设置为Th₁，这里Th₁取低值，以保证微弱目标信号的信息保留，门限Th₁的取值满足虚警概率为10^-4～10^-3；该步骤具体过程分为以下2步：

(a)矩阵S_nm中m列元素分别作短时傅里叶变换，变换结果f_nm可表示为：

$f_{\mathrm{nm}}=\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{Y}_{\mathrm{nm}}g(n-l)\exp (-j2\mathrm{\&pi;fnT})---\left(5\right)$

(5)式中g(n)为高斯窗函数，l是窗函数滑动步进；

(b)对(5)式进行恒虚警处理的结果f′_nm可表示为：

步骤四、利用短时傅里叶变换的逆变换性质还原(6)式的相位信息，(6)式还原后的时域形式S′_nm可表示为：

S′_nm＝ISTFT{f′_nm·angle[S_nm]}    (7)

(7)式中ISTFT表示为短时傅里叶变换的逆变换；

步骤五、对步骤四结果进行归一化的分数阶傅里叶变换，其表达式为：

F′_α(u)＝F_α(S′_nm/max|S′_nm|)    (8)

(8)式中α为分数阶域变换阶次；

步骤六、对步骤五的结果进行二维搜索，把超过门限Th₂的点作为目标信号，其输出结果记为：

(9)式中i为满足(9)式的目标数目，门限Th₂取值满足虚警概率为10^-6～10^-5；

步骤七、为了减少脉冲噪声或强信号对弱信号造成的虚警影响，对(9)式获得的i个点目标信号采用图像对比度方法来进行评价，图像对比度定义为：

$D_i=\frac{1}{Q}\mathrm{\&Sigma;}{\left|\right[F^{\&prime;\&prime;}({\stackrel{\&OverBar;}{u}}_i,{\stackrel{\&OverBar;}{\mathrm{\&alpha;}}}_i),(a,b)\left]\right|}^2-{\left|\frac{1}{Q}\mathrm{\&Sigma;}\right[F^{\&prime;\&prime;}({\stackrel{\&OverBar;}{u}}_i,{\stackrel{\&OverBar;}{\mathrm{\&alpha;}}}_i),(a,b)\left]\right|}^2---\left(10\right)$

(10)式中D_i表示图像对比度，值越大，代表变换阶次估计值越接近真实值；(a，b)代表空间搜索区域的范围，Q为(a，b)区域内总的点数；

步骤八、对步骤七中获得的图像对比度进行降序排列，对超过门限Th₃的目标点进行标记，并把对应的目标点

记录下来，门限Th₃取值满足虚警概率为10^-6～10^-5；

步骤九、对步骤八中记录下来的点信号进行参数估计，其表达式如下：

${\hat{f}}_i={\stackrel{\&OverBar;}{u}}_i\csc {\stackrel{\&OverBar;}{\mathrm{\&alpha;}}}_i,$ ${\hat{k}}_i=-\cot {\stackrel{\&OverBar;}{\mathrm{\&alpha;}}}_i---\left(11\right)$

(11)式中

对应目标速度引起的多普勒频率，

对应目标加速度引起的调频率；

步骤十、采用逐次消去法去除步骤八中记录的所有点，然后重复步骤三到步骤九，直到所有目标点在分数阶傅里叶变换域的峰值低于预设门限为止，这里选择本次目标峰值的50％为下一次的预设门限值；

步骤十一、为了降低虚警概率，对获得的所有目标参数进行凝聚处理，凝聚处理后的目标个数和参数即该多目标检测方法检测到的多目标个数和参数。

2.根据权利要求1所述的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法，特征在于：所述的该多目标检测方法中：

a.所述的步骤三(a)中高斯窗函数的选择：

高斯函数中窗越宽，时间分辨率越低；窗越窄，时间分辨率会提高，但频率分辨率又会降低；本发明选择的高斯窗函数表示为：

g(n)＝exp[-π(nt_s)²B_x/T_x]    (12)

(12)式中t_s为信号采样间隔，T_x为被测信号x(nt_s)对应的时域支撑宽度，B_x为被测信号x(nt_s)对应的频域支撑宽度，且

T_x，B_x分别定义为：

$T_x=\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{\left[{({\mathrm{nt}}_s-{\mathrm{\&eta;}}_n)}^2{\left|x\left(n{t}_s\right)\right|}^2\right]}^{1/2}}{\left|\right|x\left|\right|},$ $B_x=\frac{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{\left[{(\mathrm{n\&Delta;f}-{\mathrm{\&eta;}}_f)}^2{\left|X\left(\mathrm{n\&Delta;f}\right)\right|}^2\right]}^{1/2}}{\left|\right|X\left|\right|}---\left(13\right)$

这里

其中x(nt_s)和X(nΔf)是一对傅里叶变换对，Δf为傅里叶变换域的频率分辨率，||x||为x(nt_s)的范数，||x||为X(nΔf)的范数，但须注意本发明中用到的信号x(nt_s)是指其公式(3)中第m个时延单元中n次目标回波经脉压处理后的离散信号S_nm；

b.所述的步骤三(a)中窗宽度BW的选择：

由(12)式知高斯函数最大值下降0.707倍对应的窗宽度为：

$\mathrm{BW}\&ap;0.5\sqrt{\frac{T_x}{\mathrm{\&pi;}{B}_x}}---\left(14\right)$

从(14)式知窗的宽度与T_x/B_x成正比，T_x/B_x值越大，BW越大，反之亦然，T_x和B_x的取值依据(13)式中被测信号的能量谱密度，功率谱密度及信号在时频域的幅度决定；(13)式中η_n为被测信号所占时间在整个测试时间内的平均分布值，η_f为被测信号所占频率在整个频带上的平均分布值，|x(nt_s)|表示信号在时域的幅度，|X(nΔf)|表示信号在傅里叶变换域的幅度；

c.所述的步骤三、步骤六和步骤八中门限的选择：

多次对门限阈值进行选择，以步骤三门限选择为例，这里假设噪声服从高斯分布，其模值服从瑞利(Rayleigh)分布，则有：

$f\left(r\right)=\left\{\begin{array}{cc}\frac{r}{{\mathrm{\&delta;}}^2}\exp (-\frac{r^2}{2{\mathrm{\&delta;}}^2})& r\&GreaterEqual;0\\ 0& r<0\end{array}\right.---\left(15\right)$

(15)式中δ与噪声样本均值mean的关系为：

$\mathrm{\&delta;}=\sqrt{\frac{2}{\mathrm{\&pi;}}}\&CenterDot;\mathrm{mean}---\left(16\right)$

(16)式中δ表示噪声的方差，mean表示噪声的均值；

另外，(15)式对应的恒虚警概率可表示为：

$P_{\mathrm{fa}}=\underset{\mathrm{Th}}{\overset{\&infin;}{\&Integral;}}f\left(r\right)\mathrm{dr}=\exp (-\frac{{\mathrm{Th}}^2}{2{\mathrm{\&delta;}}^2})---\left(17\right)$

(17)式中检测门限T_h与虚警概率P_fa的关系为：

$\mathrm{Th}=\mathrm{\&delta;}\sqrt{2\ln (1/P_{\mathrm{fa}})}---\left(18\right)$

根据(18)式，步骤三中检测门限Th₁选用的P_fa取值范围为10^-4～10^-3，步骤六中检测门限Th₂选用的P_fa取值范围为10^-6～10^-5，步骤八中检测门限Th₃选用的P_fa取值范围为10^-6～10^-5。

3.根据权利要求2所述的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法，特征在于：所述的该多目标检测方法中关于分数阶傅里叶变换的变换阶次选择方法是：利用初步搜索和精细搜索相互配合完成：首先利用目标部分先验信息对回波信号中线性调频信号的调频率进行分析和估计；其次，在所估计的调频率范围内进行分数阶傅里叶变换运算，分数阶傅里叶变换阶次的搜索范围为[0，2]，步进Δ为0.1或0.01；接着在分数阶傅里叶变换域的二维平面内搜索最大峰值，利用峰值坐标得到线性调频信号调频率参数对应的分数阶阶次初估值接着以该参数为中心，搜索范围为

步进为0.1Δ进行下一次估计，依次循环缩小搜索范围使调频率估计精度逐渐增加，同时设置分数阶阶次顺序，估计偏差门限，偏差门限选择为0.1Δ，当两次顺序估计偏差低于该门限就停止循环，得到线性调频信号调频率的最优估计解。

4.根据权利要求3所述的基于短时傅里叶变换和分数阶傅里叶变换的多目标检测方法，特征在于：所述的该多目标检测方法中：

a.所述的短时傅里叶变换的短时窗宽度选择原则：首先将待检测的信号代入(13)式，接着(13)式根据被测信号的功率谱密度、幅度计算时域支撑宽度T_x，频域支撑宽度B_x；其次，将计算获得的T_x，B_x代入(14)式计算短时傅里叶变换所需的窗宽度；

b.所述的步骤三中窗函数滑动步进l选择取值范围为窗宽度的四分之一到二分之一的范围；

c.所述的步骤五中分数阶傅里叶变换的变换阶次α选择取值范围：首先对被测信号在[0，2]的阶次范围内以Δ为0.1或0.01的步进进行初步检测，接着以初步检测获得的结果为中心，依次循环缩小搜索范围，同时设置分数阶阶次顺序估计偏差门限，偏差门限选择为0.1Δ，当两次顺序估计偏差低于该门限就停止循环；

d.所述的步骤六中对步骤五结果以分数阶傅里叶域u和变换阶次域α为二维域进行幅值搜索，把幅值超过预设门限的点作为目标信号记录下来，其中门限的设置要符合(18)式，这里设虚警概率P_fa的范围为10^-6～10^-5；

e.所述的步骤七中(10)式的图像对比度D_i选择取值范围为0.6～1；

f.所述的步骤九中对经过步骤八处理后记录下来的目标点

带入(11)式进行多普勒频率的估计以及调频率

的估计，然后将

代入

中求得目标速度v_i，

代入

中求得目标加速度a_i。

Images