在现代战争中，随着导弹、歼击机等空中目标的机动性越来越强，雷达检测与追踪它们的难度也日益增加，所以提高雷达的性能就成为了夺取空战胜利中不可或缺的一环。而机载预警雷达将平台架设到飞机上，具有机动性强、覆盖范围广等优点^[1]，既提高了雷达在战场环境的生存能力，也弥补了地面雷达在探测距离上和检测低空飞行目标上的不足。但机载预警雷达的下视工作状态决定了其较地基雷达面临更加难以处理的地杂波问题，这使得它在目标检测和参数估计问题的处理上变得尤为困难^[2]。空时自适应处理(Space-Time Adaptive Processing,STAP)是一种较为可行的机载雷达地杂波抑制方法，但STAP方法在处理雷达目标回波时，假设回波数据在相干处理时间(Coherent Processing Interval,CPI)内的多普勒频率不变。而机载预警雷达在检测机动目标时，其回波会产生多普勒走动^[2](其多普勒频率在CPI内发生变化)，这样也使得原有的STAP方法对目标能量的相参积累性能大大下降，从而导致雷达对机动目标的检测性能下降^{[3, 4]}。因此研究机动目标检测和参数估计具有重要意义。

当机动目标在一个CPI内做匀加速运动时，目标回波可以看作线性调频(Linear Frequency Modulation,LFM)信号^[5]。目前LFM信号检测和参数估计方法主要有最大似然估计^[6](Maximum Likelihood,ML)和时频分析^[7]等方法。虽然最大似然估计精度很高，对参数的估计方差接近于Cramer-Rao下界，但其过大的运算量不利于该方法在实际情况中进行实时处理和工程实现。而时频分析方法要求有足够的脉冲采样点数，否则目标参数的估计精度难以满足要求^{[8, 9, 10, 11, 12]}。当脉冲重复频率恒定时，脉冲数的增加会使得CPI变大，此时目标和杂波的回波信号可能会产生距离走动，使得后续的处理更加困难。文献^[9]利用Radon-傅里叶变换(Radon-Fourier Transform,RFT)校正目标回波的距离走动，但仅适于匀速运动目标。文献^[10]利用Radon-分数阶傅里叶变换(Radon-Fractional Fourier Transform,RFRFT)校正目标回波的距离走动和多普勒走动，但在多天线情况下无法直接应用。文献^[11]在单天线体制下利用Radon-分数阶模糊函数(Radon-Fractional Ambiguity Function,RFRAF)技术能够有效地校正高机动目标(目标加速度变化，即存在急动度)的距离走动和多普勒走动，同样不能直接应用于多天线体制。另外文献^[12]提出的基于重构时间采样的机动目标检测和参数估计算法能够在脉冲数较少的情况下获得精度较高的参数估计，但该方法需要目标方位先验信息已知，当方位信息未知时，该方法无法直接应用。

本文提出了一种基于WVD和重构时间采样的空中机动目标检测和参数估计方法，该方法能够在脉冲数目较少，方位信息未知的情况下，利用重构时间采样和WVD方法对目标回波信号进行有效积累，然后采用变步长3维搜索实现对方位未知机动目标的参数估计。仿真结果验证了所提方法的有效性。

设机载预警雷达为N元均匀线阵，雷达发射脉冲波长为λ，阵元之间的间隔为d=0.5λ。雷达在一个相干处理间隔内发射K个脉冲，x_nk 为第n个阵元在第k个脉冲上的采样值，则在一个相干处理间隔内的待检测单元回波数据是一个N×K的矩阵^[12]：

把X中的元素拉成NK×1的列向量，记作x=vec(X)，其中${\mathop{\rm vec}\nolimits} ( * )$表示矩阵向量化操作算子，这样就得到了一个空时快拍数据，如式(2)：

当待检测单元内只有1个目标时，该单元的一次快拍数据为：

其中，x_s,x_c和x_n分别表示目标、杂波和噪声的采样数据，其中x_s可以表示为：

其中$\otimes$表示Kronecker积，a_s表示目标回波复幅度，$s({\varphi _{\rm{s}}},{\varphi _{\rm{t}}})$为目标空时导向矢量，$s({\varphi _{\rm{s}}})$为空间导向矢量，$s({\varphi _{\rm{t}}})$为时间导向矢量^[10]，可分别表示为：

其中${\varphi _{\rm{s}}} = 2\pi d\cos {\theta _{\rm{s}}}\cos {\gamma _{\rm{s}}}$表示归一化的空间角频率，${\theta _{\rm{s}}}$表示目标的方位角，${\gamma _{\rm{s}}}$表示目标的俯仰角，由式(5)可以看出，匀速目标(即a=0时)和机动目标的归一化空间角频率和空间导向矢量完全相同，而归一化时间角频率和时间导向矢量并不相同。机动目标的归一化时间角频率${\varphi _{\rm{t}}}$由两部分组成：第1部分为目标初速度对应的归一化时间角频率；第2部分为目标加速度对应的归一化时间角频率，会随时间发生改变，为多普勒走动项^[12]。

Wigner-Ville分布是一种最基本、也是应用最广的时频分布之一^[13]。信号x(t)的Wigner-Ville分布为：

式(6)可以看作是信号x(t)的自相关函数${R_{\rm{x}}}(t,\tau ) = x \left(t + \frac{\tau }{2} \right){x^*} \left( t - \frac{\tau}{2} \right)$对t的傅里叶变换。LFM信号时频分布的能量集中在其瞬时频率周围，故WVD对LFM信号有着最理想的时频聚集性^[14]。

重构时间采样是利用空间采样重构时间采样，对各个阵元杂波抑制后的回波数据进行相位补偿后首尾拼接，这相当于增加了单个阵元的脉冲点数进行提高参数估计精度^{[12, 15]}。图 1为4个阵元补偿相位前后时频图，图 1(a)为拼接前的时频效果图，图 1(b)为拼接后的时频效果图。由图中看出，只要对每个阵元补偿其相应的相位，就可以达到增加单个阵元脉冲点数的效果。

图 1(Fig. 1)

图 1 4个阵元补偿相位前后时频图 Fig. 1 The time-frequency distribution of 4 elements’ data

在此为简化推导过程，下面以2个阵元为例讨论将多阵元首尾拼接时每个阵元所需补偿的相位，这两个阵元杂波抑制后的回波信号分别为：

其中，K为脉冲点数，f_p为脉冲重复频率，${f_{\rm{d}}} = \frac{{2v}}{{\lambda \cdot {f_{\rm{p}}}}}$为初始频率，${a_{\rm{d}}} = \frac{{2a}}{{\lambda \cdot \mathop f\nolimits_{\rm{p}}^2 }}$为调频率，v为目标的初始速度，a为目标的加速度，λ为雷达工作波长，$\theta$是目标方位。

当第1个阵元的脉冲数由K增加到2K时，第1个阵元接收到的回波数据可表示为：

其中$\odot$表示Hadamard积，由式(8)可见，只要补偿第2个阵元的相位后再进行拼接就相当于将两个阵元的脉冲点拼接到第1个阵元上，由式(8)得到阵元2相应的相位补偿量应为：

以此类推，第n个阵元对应的补偿相位为：

其中，n=2,3,$\cdots$,N。根据式(10)对N个阵元进行相位补偿后，可以将N个阵元的脉冲数据拼接到1个阵元上，起到了利用空域信息来重构时域信息的效果，进行WVD变换后可使目标能量得到有效积累，用于进行目标参数估计。

由式(10)可知，在重构时间采样时，所需的补偿相位中包含未知的方位角$\theta$和初始频率项f_d以及调频率项a_d。${f_{\rm{d}}} = \frac{{2v}}{{\lambda \cdot {f_{\rm{p}}}}}$中包含未知参数v,${a_{\rm{d}}} = \frac{{2a}}{{\lambda \cdot \mathop f\nolimits_{\rm{p}}^2 }}$中包含未知参数项a，所以无法直接将多个阵元的回波数据(杂波抑制后)进行拼接。在构造代价函数时需进行速度、加速度、方位角3维搜索。代价函数可表示为：

其中x_s,n为第n个阵元杂波抑制后的回波数据，${\varphi _n}$为第n个阵元所需的补偿相位。W_s($\cdot$)为WVD变换算子符号，${\small \sum} {\left\{ \cdot \right\}}$是对矩阵内所有元素相加。当拼接效果最好时，交叉项成分最少，由于WVD的交叉项成分具有震荡特性，且交叉项的平均值为零^[12]，故此时WVD变换后能量密度的绝对值之和最小。将WVD变换后能量密度的绝对值之和的倒数的最大值作为代价函数，参数估计的结果可以通过搜索代价函数的最大值得到。

本文所提方法利用了重构时间采样(其作用等效于增加单个阵元脉冲采样数)和WVD对LFM信号的频率聚集性，二者结合来提高参数估计精度，通过构造代价函数进行3维搜索实现目标参数估计。图 2给出了该方法的实现流程，随后给出了具体操作步骤。

图 2(Fig. 2)

图 2 本文方法实现流程图 Fig. 2 Flow chart of the proposed method

步骤1  杂波抑制。本文采用文献^{[16, 17]}中的方法对地杂波进行抑制，杂波抑制后的数据可表示为：

其中$\widehat {{\rm{ }}R}$为估计杂波协方差矩阵，利用待检测单元邻近的距离单元训练样本进行估计得到的，即：

其中，L为训练样本的个数，x(i)为第i(i=1,2,$\cdots$,L)个训练样本的空时快拍数据。

步骤2  确定初值进行搜索。为降低3维搜索的运算量，先确定参数对中$a,v,\,\theta$初值，首先利用延迟自相关解线性调频算法估计加速度初值，设数据x_s,n延迟$\tau$后的数据矢量为x_ds,n，则瞬时自相关函数为：

对瞬时自相关函数r_n做FFT估计其频率估值f_n，对应的目标加速度的估值为a_n=f_n/2$\tau$。对每一个阵元的接收数据估计加速度并求平均，可得粗略估值作为加速度初值表示为：

另外，本文分别使用了单脉冲测角法和FFT测频法来估计方位角$\theta$和速度v的初值，估计初速度时，先用每个阵元数据估计目标速度初值，最终初值为各阵元数据估计结果的平均值。以这些初值构造搜索范围进行3维搜索。

步骤3  3维搜索中，先重构时间采样后WVD变换构造代价函数。本文方法在3维搜索中，先对个阵元的数据进行相位拼接完成重构时间采样，WVD变换后得到重构数据的WVD能量谱。根据式(11)构造代价函数，在最优匹配时WVD变换后能量密度值的绝对值之和最小，即在变步长搜索中代价函数取得最大值，此时所对应的参数v₀,a₀为该目标参数的最优解。

在本文方法中，为降低算法运算量，3维搜索采用了变步长搜索，先使用较大搜索步长进行粗搜索，来确定后续精搜索的初值，取该初值的左右各两个粗搜步长构造精搜索范围，精搜索使用较小的搜索步长来确保搜索的精确度。

本节将通过仿真实验来验证该方法的有效性。仿真参数设置：天线为阵元间距是d=0.5λ，阵元数为N=16的正侧视理想均匀线阵。载机速度为120 m/s，载机平台高度为10 km，雷达的距离分辨率为20 m，雷达工作波长为λ=0.32 m，相干处理脉冲数为K=64，脉冲重复频率为f_p=1500 Hz，输入信噪比SNR=0 dB，杂噪比CNR=50 dB。设机动目标处于待测单元内，处于方位角90°处，初速度为24.01 m/s，加速度为99.9 m/s²。实验估计参数均方根误差进行了200次蒙特卡罗实验。

图 3(a)为总回波的功率谱，因为信杂比较低，回波信号中杂波信号完全覆盖了目标信号。图 3(b)是用杂波协方差矩阵求逆法来抑制杂波后的功率谱，由于未补偿目标加速度，目标在杂波背景下仍然难以检测。图 3(c)为本文方法得出的功率谱，加速度补偿后，目标能量在多普勒域中得以聚集。图 4中蓝线和红线分别为图 3(b)和图 3(c)中$\cos \psi = 0$(其中$\cos \psi = \cos \theta \cos \gamma$)时的功率谱。由图 4可见，未补偿加速度时，由于目标加速度导致的多普勒走动，使得其回波信号功率谱在多普勒域上有着一定的展宽，不利于对机动目标进行检测。而本文方法利用步骤2得出的加速度初值对目标加速度进行补偿后目标能量得到了更为有效的累积，有利于对机动目标的检测。

图 3(Fig. 3)

图 3 回波信号傅氏谱 Fig. 3 Spectrum of signals

图 4(Fig. 4)

图 4 对加速度项补偿前后功率谱图 Fig. 4 Spectrum before and after compensating the acceleration

图 5表示所估参数的均方根误差随信噪比的变化。本文所使用的3维搜索为变步长搜索，具体搜索参数设置为：加速度的搜索范围为初值的±100 m/s²，粗搜索步长为4.0 m/s²，细搜索步长为0.5 m/s²，速度的搜索范围为初值的±60 m/s，粗搜索步长为2.0 m/s，细搜索步长为0.2 m/s，方位角的搜索范围为初值的±10°，粗搜索步长为1°，细搜索步长为0.2°。图中单阵元法为仅利用1个阵元数据做分数阶傅里叶变换后的参数估计结果，非相干积累法是指对所有阵元做分数阶傅里叶变换后的参数估计结果进行统计平均后得到的均值。图 5(a)为初速度均方根误差随信噪比的变化，图 5(b)为加速度均方根误差随信噪比的变化。结果表明，在补偿加速度并进行重构时间采样后，本文方法比单阵元法和非相参积累法参数估计精度有了显著的提高，且当信噪比较低时本文方法估计效果更好。图 6给出了不同搜索间隔条件下的参数估计性能曲线，由此可见，使用较小的步长可以提高参数估计精度。

图 5(Fig. 5)

图 5 本文所提方法的参数估计性能曲线 Fig. 5 Estimation performance of the parameters with the proposed method

图 6(Fig. 6)

图 6 不同搜索间隔条件下的参数估计性能曲线 Fig. 6 Estimation performance of the parameter with different interval

表 1给出了本文方法中各个步骤所需运算量，其中L为估计协方差矩阵所需样本数，V_cz1,a_cz1,$\theta$_cz1和$\Delta {v_1}$,$\Delta {a_1}$,$\Delta {\theta _1}$分别表示速度、加速度和方位的粗搜范围和步长，V_cz2,a_cz2,$\theta$_cz2和$\Delta {v_2}$,$\Delta {a_2}$,$\Delta {\theta _2}$分别表示速度、加速度和方位的细搜范围和步长。

表 1(Tab. 1)

表 1 各个步骤所需运算量分析Tab. 1 Computation analysis for each steps 本文方法步骤 复乘次数 估计协方差矩阵并求逆抑制杂波 $L\cdot{N^2}\cdot{K^2} + O[{(N·K)^3}] + NK·NK$ 估计参数搜索初值 $N\cdot\frac{K}{2}\cdot\log _2^K + N\cdot(\frac{K}{2}\cdot\log _2^K + K) + N·$ 重构时间采样和WVD变换 $(N - 1)·K\cdot{(NK)^2}$ 3维搜索估计参数 $(\frac{{{V_{cz1}}}}{{\Delta {v_1}}} + 1)(\frac{{{V_{cz2}}}}{{\Delta {v_2}}} + 1)(\frac{{{a_{cz1}}}}{{\Delta {a_1}}} + 1)(\frac{{{a_{cz2}}}}{{\Delta {a_2}}} + 1)(\frac{{{a_{cz1}}}}{{\Delta {a_1}}} + 1)(\frac{{{a_{cz2}}}}{{\Delta {a_2}}} + 1)$

表 1 各个步骤所需运算量分析 Tab. 1 Computation analysis for each steps

表 2给出了本文中3种方法进行1次蒙特卡洛实验的运算时间对比，仿真实验采用Matlab R2014a软件进行仿真，计算机参数为Inter(R) Core(TM) i3-3220处理器；主频：3.30 GHz；内存：8.00 GB。本文方法虽然较之单阵元法和非相干积累法运算量有着增大，但其估计精度有着显著的提高。

表 2(Tab. 2)

表 2 运算时间比较Tab. 2 Comparison of computation time 参数估计方法 运算时间(s) 本文方法 418.17 单阵元 6.06 非相干积累 17.53

表 2 运算时间比较 Tab. 2 Comparison of computation time

5    结论

本文所提出的基于WVD的机动目标检测与参数估计方法能够在目标方位未知，脉冲数目较少的情况下，得到较高精度的参数估计。通过仿真实验表明，该方法利用WVD和重构时间采样对回波信号进行相干积累，进行变步长3维搜索后得到的参数估计精度较传统的单阵元法和非相参积累法更高，仿真实验证明了该方法的有效性。另外，在本文中虽然通过变步长搜索来实现降低算法运算量，但其运算量仍然较大，所以在保证参数估计精度的同时降低运算量是下一步研究和改进的重点。