随着航空航天技术的不断发展，具有高速、高机动以及超远程特性的目标不断涌现，例如飞行速度可达20马赫(约合6.8 km/s)的弹道导弹以及速度可达25马赫(约合8.5 km/s)的空天飞机等，对此类目标进行探测和跟踪是当今雷达技术领域的热点和难点之一。

对于来自高速运动目标的雷达回波信号，在对其进行积累的过程中会遇到由于目标高速运动带来的跨距离单元(cross range unit, ARU)和由于目标在雷达径向方向上的加速度导致的多普勒频率扩散(Doppler frequency migration, DFM)问题^[1]。ARU与DFM的出现削弱了雷达对目标进行能量积累的增益。长时间的相干积累是提高雷达回波能量积累增益的有效手段，对于装备了相控阵天线和使用了数字波束形成技术的雷达，长时间的相干积累成为可能^[2]。对于远距离的高速运动目标，其在雷达视线方向的运动难以等效为匀速运动模型，这使得研究者不得不考虑在长的相干积累间隔(coherent processing interval, CPI)内非线性运动分量带来的影响^[3]。为了实现对高速运动目标雷达回波信号的有效积累，常规方法利用线性调频信号(linear frequency modulated signal, LFM)的性质实现脉冲内目标能量的积累，在脉冲间则是通过相干或是非相干的方式实现能量积累。

对于非相干积累方法，不需要雷达系统严格的相干性，容易实施，但是由于其不能充分利用回波的相位信息，使得非相干积累方法的积累增益较低。Hough变换是一种常见的非相干积累方法，被频繁应用于图像处理中，能够很好地定位噪声平面的直线^[4-6]。Radon变换也是一种常见的非相干积累方法，它与Hough变换的区别在于前者是直线参数变换的离散形式，后者则是直线参数变换的连续形式，Radon变换常被应用在雷达信号处理和其他领域^[7]。

相干积累方法同时利用雷达回波信号中的幅度和相位信息，因此，相比于非相干积累方法，相干积累方法可以获得更高的积累增益。拉东傅里叶变换(Radon Fourier transform, RFT)方法是一种相干积累方法，该方法基于运动目标的径向速度与ARU、多普勒频率的关系进行实现，通过联合搜索雷达与目标的径向距离和速度实现长时间的相干积累^[8-11]。RFT引入了一个多普勒滤波器，该滤波器可以抵消速度带来的多普勒相位项的拉东变换，从而实现一个CPI内的相干积累。然而，该方法存在算法复杂度过大、盲速旁瓣和参数设计复杂等缺点。

轴旋转目标检测方法通过旋转回波数据的位置坐标来校正目标回波的线性ARU，并且可以通过搜索旋转角度估计目标的径向速度，然后通过慢时间方向的傅里叶变换实现目标回波的能量积累^[12]。

以上所提相干积累方法仅仅考虑目标速度对ARU的线性影响关系，针对该关系进行校正或者速度的搜索补偿实现目标信号的能量积累。上述方法对于目标径向加速度带来影响是通过限制脉冲的积累数量加以限制的，对于存在径向加速度的目标有一定的性能损失，并且无法充分利用所有回波脉冲。

针对二阶ARU和DFM的校正需要考虑到径向加速度带来的影响，对于径向加速度带来的非线性ARU，可以从距离频率域的补偿或者时域的积累方式考虑。

拉东分数阶傅里叶变换(Radon fractional Fourier transform, RFFT)方法不仅可以通过选择一个合适的角度，消除DFM的影响，而且还可以不受ARU的影响, 实现长时间的相干积累^[13]。RFFT的内核函数包含1个2次指数形式，这使得LFM信号在RFFT域中表现为一个脉冲。RFFT是一个线性变换方法，不会存在交叉项的影响，可以适用于多目标的积累。与RFT相比，RFFT的积累时间由于不受距离弯曲的限制，可以得到很大程度上的扩展。RFFT由于是三维的参数搜索过程，因此计算复杂度很高。

广义拉东傅里叶变换(generalized Radon Fourier transform, GRFT)定义了一个多参数的多普勒滤波器，通过联合搜索参数的方式实现在参数域的相干积累^[14]，RFT是GRFT的一个特殊例子。GRFT是一个线性变换，因此对多目标进行GRFT不会产生交叉项。但是，由于GRFT的多维参数搜索需要对每一个参数的范围区间进行遍历，这就给GRFT造成了很大的计算负担。

综上所述，现有的高速目标检测算法或是未能考虑目标非线性运动带来的影响，或是时间复杂度高，存在对目标进行检测的干扰项。

高速运动目标的检测及其运动参数的估计需要同时解决ARU、DFM和速度模糊问题。为了实现对高速运动目标信号的有效积累，提高实时性要求，获取高精度的目标运动参数，本文提出一种基于运动参数估计的高速目标检测方法(motion parameter estimation moving target detection, MPE-MTD)。首先通过分段相干积累的方式抑制分段下的ARU现象，获取每个分段的相干积累结果；其次，通过恒虚警(constant false alarm rate, CFAR)检测的方法获取段内目标峰值点的位置信息，联合各分段的峰值点位置信息估计目标运动方程参数；再次，采用一维的可变步长的参数搜索方式获取精确的目标径向加速度，构造多普勒滤波器补偿原回波数据中的非线性分量；再次，通过Keystone变换校正线性ARU，利用估计获得的速度确定模糊数搜索范围，构造滤波器对模糊项进行补偿；最后，通过方位向的快速傅里叶变换(fast Fourier transform, FFT)实现对目标信号有效的相干积累。

1 回波信号模型 1.1 高速运动目标雷达回波信号模型

对于来自目标的回波，其回波信号可以建模为

$ \begin{align*} S\left(t_{\mathrm{r}}, t_{\mathrm{a}}\right)= & A_{0} \sum\limits_{i=1}^{N} \operatorname{rect}\left(t_{\mathrm{r}}-\frac{2 R_{i}\left(t_{\mathrm{a}}\right)}{c}-\frac{t_{\mathrm{p}}}{2} ; t_{\mathrm{p}}\right) . \\ & \exp \left(\mathrm{j} {\rm{\mathsf{π}}} K_{\mathrm{r}}\left(t_{\mathrm{r}}-\frac{2 R_{i}\left(t_{\mathrm{a}}\right)}{c}-\frac{t_{\mathrm{p}}}{2}\right)^{2}\right) . \\ & \exp \left(-\mathrm{j} 4 {\rm{\mathsf{π}}} f_{\mathrm{c}} \frac{R_{i}\left(t_{\mathrm{a}}\right)}{c}\right), \end{align*} $

(1)

其中：$N$表示目标个数，$t_{\mathrm{r}}$代表距离向快时间，$t_{\mathrm{a}}$代表方位向慢时间，$A_{0}$表示回波幅度，$K_{\mathrm{r}}$表示调频率，$c$代表光速，$f_{\mathrm{c}}$代表载频，$T_{\mathrm{p}}$代表脉冲宽度，$R_{i}\left(t_{\mathrm{a}}\right)$表示雷达与目标的瞬时斜距，$\operatorname{rect}\left(t ; t_{\mathrm{p}}\right)$表示矩形函数，其可以表达为

$ \operatorname{rect}\left(t ; t_{\mathrm{p}}\right)= \begin{cases}1, & -\frac{t_{\mathrm{p}}}{2} \leqslant t \leqslant \frac{t_{\mathrm{p}}}{2}, \\ 0, & t<-\frac{t_{\mathrm{p}}}{2} \text { 或 } t>\frac{t_{\mathrm{p}}}{2} .\end{cases} $

(2)

目标实际的运动分量可以被划分到径向和切向上，但在实际中只有径向方向的运动分量会对回波造成影响。由于目标不一定沿着雷达径向方向运动，并且速度也可能随时发生改变，因此其在径向方向的运动应等效为非匀速运动模型。目标沿着雷达径向方向在一个CPI内的运动可以根据其初始斜距、相对径向速度和相对的径向加速度建模为一个多元函数^[15]，表达为

$ R_{i}\left(t_{\mathrm{a}}\right)=R_{0 i}+v_{i} \cdot t_{\mathrm{a}}+\alpha_{i} \cdot t_{\mathrm{a}}^{2}, $

(3)

其中：$v_{i}$表示第$i$个目标的径向速度，$\alpha_{i}$表示第$i$个目标的径向加速度，$R_{0 i}$表示第$i$个目标与雷达的初始斜距。对回波数据进行脉冲压缩之后，信号在距离-频率-方位-时间域的表达式近似为

$ \begin{align*} & S_{\mathrm{hr}}\left(f_{\mathrm{r}}, t_{\mathrm{a}}\right)=A_{1} \operatorname{rect}\left(f_{\mathrm{r}}, t_{\mathrm{a}}\right) \sum\limits_{i=1}^{N} \cdot \\ & \quad \exp \left(-\mathrm{j} 4 {\rm{\mathsf{π}}} \frac{f_{\mathrm{r}}}{c}\left(R_{0 i}+v_{i} \cdot t_{\mathrm{a}}+\alpha_{i} \cdot t_{\mathrm{a}}^{2}\right)\right). \\ & \quad \exp \left(-\mathrm{j} 4 {\rm{\mathsf{π}}} \frac{f_{\mathrm{c}}}{c}\left(R_{0 i}+v_{i} \cdot t_{\mathrm{a}}+\alpha_{i} \cdot t_{\mathrm{a}}^{2}\right)\right) , \end{align*} $

(4)

其中：$f_{\mathrm{r}}$表示回波的距离向频率。将其转换回距离-时间-方位时间域下，其表达式为

$ \begin{array}{r} S_{\mathrm{hr}}\left(t_{\mathrm{r}}, t_{\mathrm{a}}\right)=A_{2} \sum\limits_{i=1}^{N} \operatorname{sinc}\left(B\left(t_{\mathrm{r}}-\frac{2 R_{i}\left(t_{\mathrm{a}}\right)}{c}\right)\right) \cdot \\ \exp \left(-\mathrm{j} 4 {\rm{\mathsf{π}}} \frac{f_{\mathrm{c}}}{c}\left(R_{0 i}+v_{i} \cdot t_{\mathrm{a}}+\alpha_{i} \cdot t_{\mathrm{a}}^{2}\right)\right), \end{array} $

(5)

其中：B表示LFM信号的带宽。

1.2 ARU与DFM现象

在式(5)中可以看到sinc包络的回波延时是关于慢时间$t_{\mathrm{a}}$的函数，即对于不同的方位向脉冲，sinc函数的峰值是变化的，其在距离时间轴的位置为$2 R_{i}\left(t_{\mathrm{a}}\right) / c$。

根据傅里叶变换的性质，观察式(4)可知第1个指数项$\exp \left(-\mathrm{j} 4 {\rm{\mathsf{π}}} \frac{f_{\mathrm{r}}}{c}\left(R_{0 i}+v_{i} \cdot t_{\mathrm{a}}+\alpha_{i} \cdot t_{\mathrm{a}}^{2}\right)\right)$中的$f_{\mathrm{r}}$与$t_{\mathrm{a}}$的耦合关系正是造成ARU的原因。因此，为了校正ARU，可以考虑对$f_{\mathrm{r}}$与$t_{\mathrm{a}}$进行解耦合操作。

式(4)中的第2个指数项$\exp \left(-\mathrm{j} 4 {\rm{\mathsf{π}}} \frac{f_{\mathrm{c}}}{c}\left(R_{0 i}+\right.\right. \left.v_{i} \cdot t_{\mathrm{a}}+\alpha_{i} \cdot t_{\mathrm{a}}^{2}\right)$)与常规匀速运动目标模型的不同之处在于运动模型中的径向加速度引入了多余的2次相位$\exp \left(-\mathrm{j} 4 {\rm{\mathsf{π}}} \frac{f_{\mathrm{c}}}{c} \alpha_{i} \cdot t_{\mathrm{a}}^{2}\right)$。此时对方位向进行FFT会导致多普勒频率主瓣的展宽，进一步恶化能量积累的效果^[16-17]。

1.3 速度模糊现象

通过分析多普勒频率可以获取目标的径向速度，对于高速目标来说，雷达系统的PRF选择将会影响目标径向速度的获取结果。PRF选择过小会导致速度模糊现象的产生，相当于方位向信号的采样率小于目标速度产生的多普勒频率，导致目标的多普勒频率发生折叠。根据式(5), 单个目标的多普勒频率项为

$ f_{\mathrm{d}}=\exp \left(-\mathrm{j} 4 {\rm{\mathsf{π}}} \frac{f_{\mathrm{c}}}{c} v_{\text {real }} t_{\mathrm{a}}\right), $

(6)

其中：$v_{\text {real }}$为通过多普勒频率分析实际获取的目标径向速度。它与目标的真实径向速度$v_{i}$的关系可以表示为

$ v_{\text {real }}=v_{i}+N_{\mathrm{am}} v_{\mathrm{am}}, $

(7)

其中：N_am为速度模糊数；v_am为当前PRF下不模糊的最大运动径向速度，其与PRF的关系可以表达为

$ v_{\mathrm{am}}=\frac{\mathrm{PRF} \cdot \lambda}{2} . $

(8)

因此，式(6)可以改写为

$ f_{\mathrm{d}}=\exp \left(-\mathrm{j} 4 {\rm{\mathsf{π}}} \frac{f_{\mathrm{c}}}{c} v_{i} t_{\mathrm{a}}\right) \exp \left(-\mathrm{j} 4 {\rm{\mathsf{π}}} \frac{f_{\mathrm{c}}}{c} N_{\mathrm{am}} v_{\mathrm{am}} t_{\mathrm{a}}\right) . $

(9)

从式(9)可以看到，针对DFM现象进行补偿时必须考虑模糊项的影响。

2 运动参数估计 2.1 分段相干积累

由于目标径向上速度和加速度的存在，无法直接对脉冲压缩后的信号进行相干积累。为得到目标运动方程的估计参数，首先采用分段相干积累的方式对回波数据矩阵进行处理，在每个分段下通过限制分段脉冲数来限制回波脉冲的跨距离单元数，从而在一定程度上抑制ARU。

根据式(5)，可知对于单个目标，不同脉冲下的跨距离单元数与方位向时间t_a的关系为

$ n_{\mathrm{ARU}}\left(t_{\mathrm{a}}\right)=\frac{2\left(v_{i} \cdot t_{\mathrm{a}}+\alpha_{\max } \cdot t_{\mathrm{a}}^{2}\right)}{\Delta t_{\mathrm{r}} \cdot c}, $

(10)

其中：Δt_r表示快时间采样间隔。假设某一分段的起始时间为t_s0，分段结束时间为t_s1。于是该分段下的距离单元跨越数量可以表示为

$ \begin{align*} \Delta n_{\mathrm{ARU}} & =\frac{2\left(v_{i} \cdot\left(t_{\mathrm{s} 1}-t_{\mathrm{s} 0}\right)+\alpha \cdot\left(t_{\mathrm{s} 1}^{2}-t_{\mathrm{s} 0}^{2}\right)\right)}{\Delta t_{\mathrm{r}} \cdot c} \\ & \approx \frac{2\left(v_{i} \cdot\left(t_{\mathrm{s} 1}-t_{\mathrm{s} 0}\right)\right)}{\Delta t_{\mathrm{r}} \cdot c} \\ & =\frac{2\left(v_{i} \cdot \mathrm{PRI} \cdot n_{\mathrm{s}}\right)}{\Delta t_{\mathrm{r}} \cdot c}, \end{align*} $

(11)

其中：PRI表示脉冲重复间隔，$n_{\mathrm{s}}$表示分段脉冲数。为提高相干积累的效率，适当地选择分段脉冲数，控制分段内的距离单元跨越数在一个距离单元以内，可以提高脉冲的相干积累效果。于是，每个分段下的脉冲数选择可以根据下式进行

$ \begin{align*} & \Delta n_{\mathrm{ARU}} \leqslant 1 \\ = & >\frac{2\left(v_{i} \cdot \mathrm{PRI} \cdot n_{\mathrm{s}}\right)}{\Delta t_{\mathrm{r}} \cdot c} \leqslant 1 \\ = & n_{\mathrm{s}} \leqslant \frac{\Delta t_{\mathrm{r}} \cdot c}{2 v_{\max } \cdot \mathrm{PRI}}, \end{align*} $

(12)

其中：$v_{\text {max }}$表示实际应用中所感兴趣的径向速度的最大值。在进行目标检测时，目标径向上的速度与加速度不一定接近最大值且方向不一定相同。因此，分段内可以容忍一定的ARU并且适当扩大分段的脉冲数选择。

按照抑制ARU的条件选取的分段下脉冲数越多，分段相干积累的增益就越大，越有利于后续对积累结果进行峰值提取。然而，分段下的脉冲数选择过多会导致分段数过少，此时会影响后续的参数估计精度。这种情况下可以通过延长整个相干积累时间，或者使分段下脉冲数和分段数的选取达到平衡进行解决，适当地提升雷达系统的PRF也能在满足抑制ARU的条件下提升分段的脉冲数。

2.2 二维CFAR

在分段相干积累后，需要通过CFAR对积累的峰值位置进行获取，以便后续结合分段相干积累结果估计目标的运动参数。

CFAR的原理是在保持一定虚警率的条件下获取判决门限，对接收机输出进行判决，确定是否有目标存在。二维CFAR处理窗口可以分为待检单元、保护单元和参考单元。其中待检测单元是判断目标是否存在的单元，保护单元是考虑到目标的能量扩散效应而设置的，参考单元则是用来计算噪声电平、获取判决门限的。判决门限的设置如下

$ \beta=N_{\mathrm{u}} \cdot\left(p_{\mathrm{fa}}^{-\frac{1}{N_{\mathrm{u}}}}-1\right), $

(13)

其中：β即为最终的判决门限值，N_u为参考单元数，p_fa为虚警率。保护单元的大小需要根据具体情况进行设定。设置保护单元窗口的宽度，需要考虑分段内的ARU和脉冲积累后的能量宽度，设置保护单元窗口的高度则需要考虑DFM带来的多普勒频率展宽。保护单元窗口宽度w_pt和高度h_pt的设置如下

$ \begin{gather*} w_{\mathrm{pt}}=\Delta n_{\mathrm{ARU}}+\frac{2}{B \cdot \Delta t_{\mathrm{r}}}, \\ h_{\mathrm{pt}}=n_{\mathrm{s}}\left(n_{\mathrm{s}}-1\right) \mathrm{PRI}^{2} \cdot \frac{4 f_{\mathrm{c}} \alpha_{\max }}{c} . \end{gather*} $

(14)

由于CFAR只能以一定的概率从分段下检测到积累峰值，因此，当信噪比(signal-to-noise ratio, SNR)很低时，CFAR能够提取到的峰值结果就会变少，此时增加分段数获取尽可能多的样本点将有利于后续的运动参数估计精度。

2.3 K-means聚类

CFAR仅仅是保证在恒定虚警率下进行目标检测，并不能保证不产生虚警点。对于存在虚警点和能量扩散的检测点时，后续目标运动方程参数估计将会引入误差，影响高阶多普勒项的补偿结果，也会对Keystone变换的效果造成影响。因此，Keystone变换只适合对存在线性运动分量的目标回波进行补偿。利用K-means聚类方法可以获得聚类后的能量中心，利用获得的能量中心进行运动参数的估计，提高运动参数的估计精度。

常规的K-means聚类算法通常采用随机化初始k值的策略，后续通过迭代的方式使得损失函数最小以达到一个最优的分类数。由于初始k值的选取决定了迭代数和最终的分类效果，因此结合二维CFAR的检测点特征，本文首先对检测图像做一个最大池化操作，相当于将一定范围内目标点以最大值进行取代。最大池化的窗口选择应当与二维CFAR保护窗口的设置一致，这样可以保证池化窗口包含整个能量扩散范围。

在获得初始的k值后，计算样本点与每一个聚类中心的距离，将距离最近的点归为一类，而当样本点与质心的最小距离超过某一门限值r_t/2时将样本点自身作为一个新的聚类中心，而当质心间的距离小于r_t/2时将质心合并，r_t的定义如下

$ r_{\mathrm{t}}=\sqrt{h_{\mathrm{pt}}^{2}+w_{\mathrm{pt}}^{2}}. $

(15)

对于每一类别，利用下式求解其质心坐标($\operatorname{cen} X_{i}, \operatorname{cen} Y_{j}$)

$ \operatorname{cen} X_{i}=\frac{\sum\limits_{n=1}^{N} x_{i n}}{N}, \quad \operatorname{cen} Y_{i}=\frac{\sum\limits_{n=1}^{N} y_{i n}}{N}, $

(16)

每次迭代直至质心坐标不再发生变化。

K-means算法的执行流程如下：

1) 池化操作减小检测散射点数，并且将池化窗口下幅度最大的检测点作为初始聚类中心，记为$\boldsymbol{\mu}_{1}^{(0)}, \boldsymbol{\mu}_{2}^{(0)}, \cdots, \boldsymbol{\mu}_{k}^{(0)}$，数量作为初始$k$值；

2) 计算聚类中心到彼此的距离$r_{i j}(i, j$表示第$i$个聚类中心到第$j$个聚类中心的距离)，合并距离$r_{i j}<r_{\mathrm{t}} / 2$的聚类中心

$ r_{i j}=\left\|\boldsymbol{\mu}_{i}-\boldsymbol{\mu}_{j}\right\| ; $

3) 当样本点到自己聚类中心的距离$d_{i}> r_{\rm t} / 2$时，自己作为新的聚类中心存在；

4) 对于每个类$j$，重新计算该类中心

$ \operatorname{cen} X_{i}=\frac{\sum\limits_{i=1}^{N_{j}} x_{i j}}{N}, \quad \operatorname{cen} Y_{j}=\frac{\sum\limits_{i=1}^{N_{j}} y_{i j}}{N} ; $

5) 定义距离损失函数：$R\left(\boldsymbol{\mu}_{i}, \boldsymbol{p}_{j}\right)= \min \sum_{j=1}^{M}\left\|\boldsymbol{p}_{j}-\boldsymbol{\mu}_{i}\right\|^{2}$；

6) 令$n_{\mathrm{it}}=0, 1, 2, \cdots$为迭代次数，重复如下过程直至$R\left(\boldsymbol{\mu}_{i}, \boldsymbol{p}_{j}\right)$收敛；

7) 对于检测点$\boldsymbol{p}_{j}$，将其分配至最近的能量中心$\boldsymbol{\mu}_{i}$，其中$i$满足

$ i=\operatorname{argmin}\left\|\boldsymbol{p}_{j}-\boldsymbol{\mu}_{i}\right\|^{2} ; $

8) 当样本点到所在聚类中心的距离$d_{i}> r_{t} / 2$时，自己作为新的聚类中心存在；

9) 对于每个类j，重新计算该类中心。

2.4 最小二乘法估计运动方程参数

对于下式所示的目标函数$h_{\theta}\left(x_{1}, x_{2}, \cdots, x_{n}\right)$

$ h_{\theta}\left(x_{1}, x_{2}, \cdots, x_{n}\right)=\theta_{0}+\theta_{1} x_{1}+\cdots+\theta_{n} x_{n}, $

(17)

其中：$\left(x_{1}, x_{2}, \cdots, x_{n}\right)$表示样本值，$\left(\theta_{1}, \theta_{2}, \cdots\right.$，$\theta_{n}$)表示模型参数。式(17)可以表达为矩阵形式

$ H_{\theta}(X)=X \boldsymbol{\theta}, $

(18)

为使得样本更加拟合目标函数，定义一个损失函数$\boldsymbol{J}(\boldsymbol{\theta})$

$ \boldsymbol{J}(\boldsymbol{\theta})=\frac{1}{2}(\boldsymbol{X} \boldsymbol{\theta}-\boldsymbol{Y})^{\mathrm{T}}(\boldsymbol{X} \boldsymbol{\theta}-\boldsymbol{Y}), $

(19)

其中：Y是样本的输出向量。当模型参数很好地拟合目标函数时，式(19)应该取得最小值，对其求偏导数，并令导数为0，可以得到

$ \boldsymbol{\theta}=\left(\boldsymbol{X}^{\mathrm{T}} \boldsymbol{X}\right)^{-1} \boldsymbol{X}^{\mathrm{T}} \boldsymbol{Y}, $

(20)

根据式(20)即可利用样本值求解目标函数模型参数的估计值。在雷达回波数据矩阵中，距离跨越单元数是关于方位向时间t_a的2次函数，因此可以将目标函数建模为

$ y(x)=a_{0}+a_{1} x+a_{2} x^{2}, $

(21)

其中：$x$即为检测样本点的横坐标。于是可以将$\boldsymbol{X}$矩阵建模为

$ \boldsymbol{X}=\left(\begin{array}{ccc} 1 & x_{1} & x_{1}^{2} \\ \vdots & & \vdots \\ 1 & x_{n} & x_{n}^{2} \end{array}\right), $

(22)

Y矩阵可以表示为

$ \boldsymbol{Y}=\left(\begin{array}{c} y_{1} \\ \vdots \\ y_{n} \end{array}\right). $

(23)

将式(22)和式(23)代入式(20)中即可求解运动方程参数。在利用样本点进行运动方程参数估计的过程中，目标的峰值轨迹上样本点越多，通过最小二乘法拟合的目标函数结果就越准确。因此，整个分段数量的选择不应过小。

由于最小二乘法无法分辨出样本点属于哪一个目标，因此，本文联合分段结果和距离向坐标排序对整个参数估计算法进行设计，通过误差判决的方式使得通过最小二乘法能够获得多个目标的样本点分类，并且获得多个目标的运动参数。

多目标参数估计算法的执行步骤如下：

1) 首先将聚类获得的聚类中心作为样本点的输入，按照距离向坐标排序；

2) 根据实验参数设置分段之间的最大偏移量d_max

$ d_{\max }=\Delta n_{\mathrm{ARU}} ; $

3) 从第1个分段开始加入样本点作为目标i峰值轨迹上的能量中心；

4) 计算当前分段样本点和上一分段样本点的距离向偏移量$d$，如果$d>d_{\text {max }}$则跳过该样本点，否则执行下一步；

5) 计算当前分段中第k个样本点加入后目标i轨迹的参数θ

$ \boldsymbol{\theta}=\left(\boldsymbol{X}^{\mathrm{T}} \boldsymbol{X}\right)^{-1} \boldsymbol{X}^{\mathrm{T}} \boldsymbol{Y} ; $

6) 根据运动方程参数计算轨迹曲线上对应距离坐标x_ij(表示第i个目标轨迹上的第j个点)的理论方位向坐标y_pre

$ y_{\mathrm{pre}}=\theta_{0}+\theta_{1} x_{i j}+\theta_{2} x_{i j}^{2} ; $

7) 计算当前分段中第$k$个点的平均误差$\operatorname{err}_{k}$

$ \operatorname{err}_{k}=\frac{1}{N_{i}} \sum\limits_{j=1}^{N_{i}}\left|y_{\mathrm{pre}}-y_{i j}\right|; $

8) 选取该分段下$\mathrm{err}_{k}$最小的样本点加入到当前目标轨迹下，若当前分段下所选样本点的误差过大，则跳过该分段进入下一个分段。回到步骤4)，直到所有分段迭代结束；

9) 当目标$i$轨迹上的样本点数量低于分段数的$1 / 4$时，将目标$i$轨迹上的样本点作为不属于任意目标的样本点剔除，并排除该目标的存在。否则样本点属于存在的目标$i$；

10) 回到步骤3)开始迭代未完成分类的样本点，直到所有样本点都完成归类。

2.5 多目标信号分离

获得目标的运动参数后需要进行速度补偿以及高阶多普勒项的补偿才能实现有效的能量积累。然而，对于多目标的场景，由于目标之间的运动参数并不相同，难以统一地完成多个目标的速度补偿和高阶多普勒项的补偿。在获得每个目标的运动参数后，可以通过设置矩形窗的位置和宽度获得分离的目标信号，矩形窗的宽度根据目标信号的能量峰值区间进行设定。

对每一脉冲$x=1, 2, \cdots, N_{\mathrm{a}}, y$的取值由式(21)决定，得到的$\left(x_{i n}, y_{i n}\right)$即为第$n$个目标在第$i$个脉冲下的能量中心坐标，于是目标的能量峰值区间可以表示为

$ \left(y_{i n}-\frac{f_{\mathrm{s}}}{B}, y_{i n}+\frac{f_{\mathrm{s}}}{B}\right) . $

(24)

3 ARU与DFM的补偿与校正 3.1 高阶多普勒项补偿

在对不同目标的信号进行分离之后，首先应当依次补偿各个目标的高阶多普勒分量。补偿高阶多普勒分量的关键在于获得精确的径向加速度。

在获取目标运动参数的估计值后，由于径向加速度的估计值与理论值之间存在较大的误差，因此不能直接使用估计的径向加速度项进行补偿。

注意到式(4)可以写为

$ \begin{align*} S_{\mathrm{hr}}\left(f_{\mathrm{r}}, t_{\mathrm{a}}\right)= & A_{1} \operatorname{rect}\left(f_{\mathrm{r}}, t_{\mathrm{a}}\right) \sum\limits_{i=1}^{N} \exp \left(-\mathrm{j} 4 {\rm{\mathsf{π}}} \frac{f_{\mathrm{r}}+f_{\mathrm{c}}}{c} R_{0 i}\right) \cdot \\ & \exp \left(-\mathrm{j} 4 {\rm{\mathsf{π}}} \frac{f_{\mathrm{r}}+f_{\mathrm{c}}}{c} v_{i} \cdot t_{\mathrm{a}}\right) \cdot \\ & \exp \left(-\mathrm{j} 4 {\rm{\mathsf{π}}} \frac{f_{\mathrm{r}}+f_{\mathrm{c}}}{c} \alpha_{i} \cdot t_{\mathrm{a}}^{2}\right) . \end{align*} $

(25)

利用估计得到的径向速度，构建下式所示的滤波器补偿式(25)中的速度项

$ h p_{v}\left(f_{\mathrm{r}}, t_{\mathrm{a}}\right)=\exp \left(\mathrm{j} 4 {\rm{\mathsf{π}}} \frac{f_{\mathrm{r}}+f_{\mathrm{c}}}{c} v_{\text {pre }} \cdot t_{\mathrm{a}}\right), $

(26)

其中：v_pre表示通过参数估计获得的目标径向速度。补偿后可以得到

$ \begin{align*} S_{0 v}\left(f_{\mathrm{r}}, t_{\mathrm{a}}\right)= & A_{1} \operatorname{rect}\left(f_{\mathrm{r}}, t_{\mathrm{a}}\right) \sum\limits_{i=1}^{N} \exp \left(-\mathrm{j} 4 {\rm{\mathsf{π}}} \frac{f_{\mathrm{r}}+f_{\mathrm{c}}}{c} R_{0 i}\right) . \\ & \exp \left(-\mathrm{j} 4 {\rm{\mathsf{π}}} \frac{f_{\mathrm{r}}+f_{\mathrm{c}}}{c} \alpha_{i} \cdot t_{\mathrm{a}}^{2}\right) . \end{align*} $

(27)

假设径向加速度的理论值和估计值的误差是δ_α，径向加速度的搜索范围设置如下

$ \left(\alpha_{\mathrm{pre}}-\delta_{\alpha}, \alpha_{\mathrm{pre}}+\delta_{\alpha}\right) , $

(28)

其中：α_pre表示通过参数估计获得的目标径向加速度。构造如下所示的滤波器

$ h p_{\alpha}\left(f_{\mathrm{r}}, t_{\mathrm{a}} ; \alpha_{\text {search }}\right)=\exp \left(\mathrm{j} 4 {\rm{\mathsf{π}}} \frac{f_{\mathrm{r}}+f_{\mathrm{c}}}{c} \alpha_{\text {search }} \cdot t_{\mathrm{a}}^{2}\right). $

(29)

径向加速度的搜索值应该满足

$ \alpha_{\text {search }}=\operatorname{argmax}\left\{\operatorname{abs}\left[\operatorname{FFT}\left(\operatorname{IFFT}\left(S_{0 v} \cdot h p_{\alpha}, f_{\mathrm{r}}\right), t_{\mathrm{a}}\right)\right]\right\}, $

(30)

其中：$\operatorname{IFFT}\left({ }^{*}, f_{\mathrm{r}}\right)$表示沿着距离向频率轴做$\operatorname{IFFT}, \operatorname{FFT}\left({ }^{*}, t_{\mathrm{a}}\right)$表示沿着方位向时间轴做FFT。

3.2 Keystone变换

式(4)中$t_{\mathrm{a}}$与$f_{\mathrm{r}}$的耦合项$\exp \left(-\mathrm{j} 4 {\rm{\mathsf{π}}} \frac{f_{\mathrm{r}}}{c} v_{i} \cdot t_{\mathrm{a}}\right)$是导致回波信号在距离时间和方位时间域下产生ARU的原因，因此只要消除该耦合项即可校正ARU。

Keystone变换即令$t_{\mathrm{a}}=\frac{f_{\mathrm{c}}}{f_{\mathrm{c}}+f_{\mathrm{r}}} t_{\mathrm{a}}^{\prime}$，则耦合指数项变为$\mathrm{e}^{-\mathrm{j} 4 {\rm{\mathsf{π}}} \frac{f_{\mathrm{c}}}{c} v_{i} \cdot t_{\mathrm{a}}^{\prime}}$，实现对耦合项的消除，达到校正ARU的目的^[18-22]。

利用chirp-z变换可以快速实现Keystone变换。Keystone的变量代换相当于对信号$\exp \left(-\mathrm{j} 4 {\rm{\mathsf{π}}} \frac{f_{\mathrm{r}}+f_{\mathrm{c}}}{c} v_{i} \cdot t_{\mathrm{a}}\right)$进行尺度缩放的操作。利用DFT的性质，信号$g(u t)$的傅里叶变换为

$ \frac{1}{|u|} G\left(\frac{w}{u}\right) . $

(31)

只要能够获得一个信号的变尺度傅里叶变换，就可以对其进行逆傅里叶变换操作，获得尺度变换后的信号。通过如下所示的变DFT即可获得信号g(ut)的变尺度傅里叶变换

$ G\left(\frac{k}{u}\right)=\sum\limits_{m=0}^{M-1} g(m) \mathrm{e}^{-\mathrm{j} \frac{2 {\rm{\mathsf{π}}} m k}{u M}} . $

(32)

之后, 利用IFFT进行逆傅里叶变换的操作，可以减小算法的复杂度。然而变DFT依旧会带来较大的计算量，采用chirp-z的方法可以利用2次FFT实现Keystone变换。

chirp-z变换将做如下的形式变换

$ \begin{gather*} X(k)=\sum\limits_{n=0}^{N-1} x(n) A^{-n} W^{n k} \Rightarrow \\ X(k)=W^{\frac{k^{2}}{2}} \sum\limits_{n=0}^{N-1}\left[x(n) A^{-n} W^{\frac{n^{2}}{2}}\right] W^{-\frac{(k-n)^{2}}{2}}. \end{gather*} $

(33)

变换之后的累加符号及其后的表达式，可以看作是2个函数线性卷积的结果。线性卷积可以利用圆周卷积进行运算，再利用FFT实现圆周卷积即可。将式(33)变换之前的形式对比式(32)可知变量之间的关系为

$ \begin{align*} & A=1, \\ & W=\mathrm{e}^{-\mathrm{j}_{u M}^{2 {\rm{\mathsf{π}}}}} . \end{align*} $

(34)

3.3 解速度模糊

由式(9)可知，需要构造下式进行模糊数的搜索。

$ h\left(N_{\mathrm{am}}\right)=\exp \left(\mathrm{j} 4 {\rm{\mathsf{π}}} N_{\mathrm{am}} v_{\mathrm{am}} \frac{f_{\mathrm{c}}}{c} t_{\mathrm{a}}^{\prime} \cdot \frac{f_{\mathrm{c}}}{f_{\mathrm{c}}+f_{\mathrm{r}}}\right) . $

(35)

假设径向速度的理论值和估计值的误差为$\delta_{v}$，考虑到$\delta_{v}$可能使得径向速度的估计值向上或者向下跨越，因此令模糊数的搜索范围为

$ \left(\left\lfloor\frac{v_{\mathrm{pre}}}{v_{\mathrm{am}}}\right\rfloor-\left\lceil\frac{\delta_{v}}{v_{\mathrm{am}}}\right\rceil, \left\lceil\frac{v_{\mathrm{pre}}}{v_{\mathrm{am}}}\right\rceil+\left\lceil\frac{\delta_{v}}{v_{\mathrm{am}}}\right\rceil\right), $

(36)

其中: *表示向下取整，*表示向上取整。

4 数值实验

为了更清楚地说明本文所提算法的处理步骤，将算法执行流程展示如图 1所示。

雷达的仿真参数分别设置为$f_{\mathrm{c}}=3 \mathrm{GHz}$，LFM信号带宽$B=10 \mathrm{MHz}$，LFM信号脉冲宽度$T_{\mathrm{p}}= 10 \mu \mathrm{~s}$，采样频率$f_{\mathrm{s}}=12 \mathrm{MHz}$，脉冲重复频率PRF $=2000 \mathrm{~Hz}$，积累脉冲数$N_{\mathrm{a}}=1024$，光速$c=3 \times$$10^{8} \mathrm{~m} / \mathrm{s}$，目标A的斜距$R_{\mathrm{A}}=37494 \mathrm{~m}$，径向速度$v_{\mathrm{A}}=-2830 \mathrm{~m} / \mathrm{s}$，径向加速度$\alpha_{\mathrm{A}}=300 \mathrm{~m} / \mathrm{s}^{2}$，目标B的斜距$R_{\mathrm{B}}=56241 \mathrm{~m}$，径向速度$v_{\mathrm{B}}=1310 \mathrm{~m} / \mathrm{s}$，径向加速度$\alpha_{\mathrm{B}}=100 \mathrm{~m} / \mathrm{s}^{2}$。

在-30 dB的SNR条件下对2个目标进行处理。首先对回波矩阵进行分段相干积累，选取的分段数为20。分段相干积累效果如图 2所示。

从图 2(a)和2(b)的对比可以看出通过分段相干积累可以获得一定的SNR增益，这将有利于后续的二维CFAR获取积累峰值位置。当所设置的分段数无法很好兼顾到不同运动参数的目标时，CFAR会检测到能量的扩散点。由于需要根据能量中心的位置进行运动参数的估计，根据扩散点的特性，使用K-means聚类方法获取扩散点的质心。

在进入K-means算法的迭代之前，首先会通过池化进行初始k值的选取，k值选取合理可以减小算法的迭代次数，使得结果收敛得更快，目标A的聚类结果如图 3所示。

获取的目标能量中心的位置将有利于后续的运动参数估计，利用最小二乘法估计目标的运动参数需要将能量中心坐标作为样本点代入式(20)，图 4(b)展示了通过运动参数估计获取的拟合峰值轨迹曲线。

从图 4中可以看出运动参数估计得到的峰值轨迹能够很好地拟合当前目标的峰值轨迹，但是运动参数的估计是存在误差的，表 1展示了2个目标的估计值。

将表 1结果与实验所设参数进行对比，可以看出在运动参数的估计结果中径向速度和斜距的估计误差较小，而径向加速度的估计误差较大。通过蒙特卡洛实验给出不同SNR下的斜距、径向速度和径向加速度估计结果的均方根误差(root mean square error, RMSE)，结果如图 5所示。

由图 5可知斜距和径向速度的RMSE相对较小，而径向加速度的RMSE相对较大。因此，可以利用估计得到的速度来补偿掉大部分的由于目标高速产生的ARU影响，为后续径向加速度精确值的搜索提供保障。

由于目标的运动参数不同，因此在进一步的处理之前需要再次利用估计得到的运动方程对得到的每个目标加窗进行信号的分离，通过式(24)设置矩形窗的宽度获得的分离信号如图 6所示。

由图 6可以看出，目标A和B的信号能量主体已经被成功分离，接下来使用估计得到的径向速度构造滤波器对相应的目标信号进行补偿，对补偿之后的数据矩阵进行径向加速度的搜索。当搜索参数接近目标值时，补偿后的信号进行方位向FFT将出现峰值。当然，由于速度并非精确的补偿，径向加速度搜索结果的精度也会受到影响。考虑通过变步长的方式进行径向加速度的搜索，先设置一个较大的搜索步长，目标A与目标B的径向加速度搜索结果如图 7所示。

由图 7可知，径向加速度搜索值越接近理论值时峰值越高。在当前的径向加速度下，通过小范围的径向加速度搜索与速度估计值的联合补偿，可以获得更精确的径向加速度。径向加速度的精确搜索结果如图 8所示。

在不同的SNR下通过蒙特卡洛实验获得精确搜索径向加速度的RMSE，结果如图 9所示。

由图 9可以看出，此时的径向加速度误差较小。在获得精确的径向加速度后，构造式(29)所示的滤波器对分离的目标信号进行补偿，通过使用Keystone变换和估计的径向速度搜索速度模糊数并进行补偿, 可以很好地校正ARU。最终的校正和积累结果如图 10所示。

由图 10(a)和10(b)可以看出，目标A和目标B的ARU完全校正了，图 10(c)~10(e)展示了良好的相干积累效果，这说明本文所提算法能够很好地解决高速目标的ARU和DFM问题，提高目标检测的效果。目标A和目标B的SNR提升见表 2。

由表 2可以看出, 积累前后目标处的SNR有较大的提升，说明本文算法能够有效实现目标信号的能量积累。最终的目标运动参数获取结果如表 3所示。

通过蒙特卡洛实验, 将所提方法与GRFT、RFT和传统的运动目标检测(moving target detection, MTD)方法进行比较。为方便观察，设置GRFT的斜距搜索次数固定为10，并能保证取到目标处的斜距值，相关算法的检测概率和运行时间结果如图 11所示。

可以看到，所提算法MPE-MTD相比于GRFT的SNR要求是更高的，但是算法的运行时间大幅度降低了。相比于RFT和传统的MTD方法，所提算法对SNR的要求更低，这是因为RFT未能消除DFM现象，MTD也无法对ARU和DFM进行校正补偿。因此，RFT和MTD也获取不了高精度的目标运动参数。

为方便比较算法的计算量，将实数加法和实数减法当作运算量相同的一次运算，省略运算量中的常数、常数系数与低阶运算量表示。对于RFT与GRFT方法，分别以N_i表示对斜距的搜索次数，以N_j表示对径向速度的搜索次数，以N_k表示对径向加速度的搜索次数。由图 1可知，本文所提算法的处理步骤可以拆解为不同模块。分段相干积累可以等效为针对整个回波数据做一个方位向FFT，其运算量为N_aN_rlog₂(N_a)。二维CFAR的运行时间主要花费在门限值计算上，令$N_{\mathrm{ref}}$表示所设参考单元的数目，则其运算量表示为$N_{\mathrm{ref}} N_{\mathrm{a}} N_{\mathrm{r}}$。K-means聚类算法的运算量主要和迭代数、检测获得的样本点数和聚类中心的选取数量有关，令$N_{\mathrm{it}}$表示迭代的次数，$N_{\mathrm{kd}}$表示聚类中心的数量，$N_{\mathrm{d}}$表示检测样本点的数量，则运算量表示为$N_{\mathrm{it}} N_{\mathrm{kd}} N_{\mathrm{d}}$。多目标参数估计的运算量受到目标数和分段数的影响，在求解运动方程的过程中主要涉及矩阵乘法和矩阵求逆，令$N_{\mathrm{t}}$表示目标的数量，$N_{\mathrm{s}}$表示分段数，则运算量表示为$N_{\mathrm{t}} N_{\mathrm{s}}^{3}$。径向加速度搜索由于是一维的参数搜索过程，因此可以使用变步长的方法降低运算量，此处主要涉及多次的FFT和IFFT，令加速度$\alpha$的搜索次数为$N_{\alpha}$，则运算量为$N_{\mathrm{t}} N_{\alpha} N_{\mathrm{a}} N_{\mathrm{r}} \log _{2}\left(N_{\mathrm{a}}\right)+ N_{\mathrm{t}} N_{\alpha} N_{\mathrm{a}} N_{\mathrm{r}} \log _{2}\left(N_{\mathrm{r}}\right)$。通过chirp-z变换实现的Keystone变换也是主要涉及FFT操作，此部分的运算量为$N_{\mathrm{a}} N_{\mathrm{r}} \log _{2}\left(N_{\mathrm{a}}\right)$。模糊补偿部分主要受搜索次数影响，主要操作为补偿和FFT，令$N_{\mathrm{Nam}}$表示模糊数的搜索次数，则运算量主要为$N_{\mathrm{t}} N_{\mathrm{Nam}} N_{\mathrm{a}} N_{\mathrm{r}} \log _{2}\left(N_{\mathrm{a}}\right)$。最后的相干积累通过方位向的FFT来实现，因此运算量为$N_{\mathrm{a}} N_{\mathrm{r}} \log _{2}\left(N_{\mathrm{a}}\right)$。算法总的运算量与RFT、GRFT相比可见表 4。

对于高速目标而言，$N_{i}$的数量级在3以上，$N_{j}$的数量级是$3, N_{k}$的数量级是2，这对于GRFT和RFT这类算法带来了很大的计算负担。本文所提MPE-MTD的运算量实际上可以看作主要受到FFT次数的影响。FFT操作带来的时间消耗主要花费在径向加速度搜索和模糊数搜索上，即主要计算量为$N_{\mathrm{t}}\left[N_{\alpha} \log _{2}\left(N_{\mathrm{a}}\right)+\right. \left.N_{\alpha} \log _{2}\left(N_{\mathrm{r}}\right)+N_{\text {Nam }} \log _{2}\left(N_{\mathrm{a}}\right)\right] N_{\mathrm{a}} N_{\mathrm{r}}$。由于已经获取到目标估计的运动参数，因此可以在确定范围内对其进行一维搜索。$\log _{2}\left(N_{\mathrm{a}}\right) 、\log _{2}\left(N_{\mathrm{r}}\right) 、N_{\mathrm{ref}} 、N_{\mathrm{t}} 、N_{\mathrm{it}} 、N_{\mathrm{s}} 、N_{\alpha}$和$N_{\mathrm{Nam}}$的数量级都是1，而$N_{\mathrm{kd}}$与$N_{\mathrm{d}}$的数量级大致在2左右，相比于GRFT与RFT这2种参数搜索方法已经有了较大提升。

5 结论

本文所提方法MPE-MTD在进行分段相干积累的基础上通过二维CFAR处理获得各分段下的峰值位置信息，联合利用各分段下的峰值位置实现目标运动方程参数的估计，搜索获得目标精确的径向加速度，分离各运动目标信号并补偿各个目标径向加速度引起的DFM现象，再通过Keystone变换和速度模糊项的补偿完成目标ARU的校正，实现可靠的相干积累，仿真实验验证了所提方法的有效性。