弹道导弹中段是弹道中时间最长的一段，可以提供较长的识别与拦截窗口，因此被认为是导弹防御的最佳时机^[1]。然而，在中段飞行过程中，所有的弹头、诱饵以及发射碎片只受重力作用，均以相同的速度沿着轨道附近运动^[2]，构成了复杂的突防目标群。中段目标识别的关键就在于从目标群中分辨出真弹头，而基于微多普勒特征的应用为其提供了新的突破口^[3]。

目前，国内许多学者针对多目标微多普勒特征分离进行了研究，其中，文献^[4]指出，多目标的时频图应该是“多目标多散射点”微多普勒的叠加。文献^[5]利用Viterbi算法构造时频滤波器，经逆时频变换提取出各信号分量，但该方法只适用于时频交叠程度较弱的多分量微动信号。文献^[6]根据对微多普勒点迹不同的相交情况进行关联配对，实现了微多普勒曲线的分离提取，但由于仅利用了斜率和夹角来判断下一时刻点迹，不能很好地解决交叉区域的干扰问题。文献^[7]经自适应视野聚类处理得到目标微多普勒的支撑域，通过对食物浓度序列进行匹配选择实现了两旋转目标的分离，但该方法容易误将弱时频分量随背景噪声去除。

本文以典型的无翼弹头类目标为例，针对窄带雷达信号体制下的中段多目标微多普勒特征分离问题展开研究。根据弹道目标多普勒的变化规律，将多条多普勒曲线的提取问题转化为求多条最短路径问题。首先利用骨架处理方法得到时频节点矩阵，再利用改进的拍卖算法迭代提取出多条最短路径，最后利用小波分析的消趋特性对平动分量进行补偿，实现了多分量微多普勒的分离与提取。

1.   进动目标多普勒分析

为了保持在大气层外飞行的稳定性和提高命中精度，弹头在中段一般维持进动姿态，本文采用文献^[8]描述的滑动型散射中心模型。进动多目标微动模型，如图 1所示。 图中：O-wuv为雷达观测坐标系；O_i-x_iy_iz_i为弹头i的滑动坐标系，O_i为弹头i的质心，z_i轴为弹头i对应的进动轴的方向，y_i轴为z_i轴顺时针旋转90°的指向，x_i轴符合右手螺旋准则；θ_i和ω_i分别为弹头i的进动角及进动角频率；β_i为S_i与弹头i自身的旋转轴之间的夹角；r_i为弹头i的底面半径；h_i为弹头i的质心O_i到底面的距离。

图  1  弹道进动多目标微动模型

Figure  1.  Ballistic multi-targets precession model

下载: 全尺寸图片 幻灯片

实际情况中必须考虑散射点的遮蔽问题，雷达的入射电磁波与目标i底面交于一点u_i，另一点已被遮蔽，则目标i的顶点A_i、底面交点u_i在t时刻处的微距离分别为^[8]

式中：(α_i,β_i)为观测雷达到弹头i的视线(S_i)在O_i-x_iy_iz_i中的视角，α_i为在O_ix_iy_i平面的投影与x_i轴的夹角; a_Ai、b_Ai、a_ui、b_ui和c_i均为调制系数，且只与弹道目标尺寸、进动参数有关，不难推导出

式中：H_i为弹头i的高; 底面半径χ_1i、 χ_2i、 χ_3i均为与β_i、θ_i有关的函数，满足χ_1i(β_i,θ_i)=cos β_icos θ_i，χ_2i(β_i,θ_i)=sin β_isin θ_i，χ_3i(β_i,θ_i)=sin β_icos θ_i；δ_i为h_i与弹头旋转对称轴的夹角。由于无翼弹头模型为旋转对称目标，可以暂时不考虑弹头的自旋运动对电磁波散射特性的影响^[9]。由式(1)和式(2)可知，顶点A_i符合正弦调制规律，而底面交点u_i由于受到相位项cos(2ω_it)的调制，其对应的微距离近似为多个正弦分量的叠加，且与目标的结构参数以及锥旋频率有关。

由于中段运动只存在地球引力的作用，因此弹道是相对平稳的，可先利用观测时间内某个脉冲测得的速度进行粗补偿^[10]，剩余的平动距离分量可近似为3阶多项式：

式中：r₀、v、a₁和a₂分别为粗补偿后单个雷达波束内包含的目标群在t时刻的径向初始距离、速度、加速度和2阶加速度。

假设雷达发射工作频率为f_c的单频信号，经目标散射后得到回波s(t)。通过基带变换后，s(t)变换为s_b(t)。对s_b(t)的相位项进行求导，得到包含i个目标的多普勒为

式中：f_mi(t)、f_tr(t)分别为微多普勒和平动多普勒，i=N⁺；∪ (·)表示包含关系；c为光速。

在较短脉冲积累时间内可以认为，粗补偿后的弹道目标回波不存在多普勒模糊^[10]，且其多普勒带宽小于雷达可探测多普勒带宽，即多普勒曲线在时频面中只出现漂移，不存在折叠现象^[11]。

2.   基于拍卖算法的多目标分离

2.1   进动周期估计

观察式(4)可知，多目标的瞬时多普勒特征在时频域中表现为一系列总体服从同一分布趋势的交叠曲线，多目标分离的实质就是将各目标包含的各散射点对应的特征曲线单独分离出来。滑动型散射中心产生的微多普勒并不符合正弦模型，不适宜使用参数化提取方法，而传统的一阶条件矩法和峰值位置法均无法处理存在交叉的多分量信号。

待分离曲线可分为两类：一类对应锥顶散射点A_i，另一类对应滑动散射点u_i。在相对较短的脉冲积累时间内，同一子目标的ω_i近似保持不变，若忽略缓慢变化的平动多普勒，则对于第1类散射点，其单位采样时间间隔内频率的变化可近似表示为

式中：T为采样时间；T_i为目标i的进动周期；B_s为信号多普勒带宽；F_s为脉冲重复频率；N和M分别为时间采样数和频率采样数；$\left\lceil {} \right\rceil$表示向上取整。

同理，对于第2类散射点，其单位采样时间间隔内频率的变化也近似满足

由式(5)和式(6)可知，无论哪类散射点，其单位采样时间间隔内频率的变化都在一个与T_i有关的阈值范围之内，这可以作为各曲线上路径选择的依据。本文采用文献^[12]中提出的改进自相关法求取多目标回波信号中子目标的个数I不同子目标对应的进动周期T_i，该方法克服了传统自相关法对先验信息要求较高的缺陷，且计算量增大倍数很小，具体步骤见文献^[12]。由于模型中每个弹头目标包含2个散射点，时频域中共包含2I条多普勒曲线。

2.2   多目标分离的最短路径描述

雷达录取多目标回波信号后，首先进行预处理。对原始回波进行奇异值分解(SVD)去噪，然后利用短时傅里叶变换(STFT)得到目标回波的时频图，记初始时频矩阵为S_M×N。对S_M×N进行高斯平滑滤波和二值化后，时频信息只包含0和1的2种元素，利用骨架提取^[13]方法进行细化，有效地对多普勒旁瓣进行抑制。最后利用消刺处理对骨架进行修剪，删去伪骨架及孤立点^[14]，记预处理后的时频矩阵为S_M×N。

将值为1的时频点看作节点，记作(m,n)，m和n分别为频率位置和时间位置，且满足{m≤M,n≤N|m,n∈N⁺}，则每个时刻可能有多个节点，其中n=1时刻的节点称为起始节点，n=N时刻的节点称为最终节点。若(l,n+1)为n+1时刻的某个节点，从(m,n)指向(l,n+1)的时频轨迹可以看作一条有向路径，记作[(m,n),(l,n+1)]，记_n^ma_n+1^l为该路径对应的长度。若确定了各条多普勒曲线的起始节点和最终节点，则多分量瞬时多普勒提取可描述为在带权图D=〈E,G,A〉中求多条最短路径的问题，E、G和A分别为D中的节点集合、路径集合及路径长度集合。

2.3   改进的最短路径拍卖算法

最短路径问题是一个常见的组合优化问题，经典方法包括Dijkstra算法、Kruskal算法和拍卖算法等。其中，拍卖算法适用于大规模稀疏网络最短路径问题的求解，并且便于程序化。该算法最早由Bertsekas和 Castanon^[15]提出，通过模拟现实中的拍卖过程来解决指派问题。本文将拍卖算法引入有向图D=〈E,G,A〉中最短路径的求解问题，首先考虑单起点单终点模型。

设P=[(m₁,1),(m₂,2),…,(m_k,k)]为D中的一条路径，若其中包含的节点各不相同，则称P为初等路，节点(m_k,k)为该路径的终点，路径长度L可表示为

该算法在迭代过程中，维护路P并通过扩展和收缩2种操作来修改P。若节点(m_k+1,k+1)不在P上，且存在路径[(m_k,k),(m_k+1,k+1)]，通过(m_k+1,k+1)对P的扩展操作将路P变为

若P包含起点之外的节点，则对P的收缩操作将路P变为

对D中每个节点(m,n)引入变量p_n^m，称其为该节点的价格。定义所有节点的价格组成的价格向量为p。算法维护价格向量p和路P，使其满足如下性质：

式(10)和式(11)称为互补松弛(CS)条件。可以证明，若二元组(P,p)满足CS条件，则路P上从节点(m₁,1)到任意节点(m_k,k)的部分是从(m₁,1)到(m_k,k)的最短路经，而p_k^m_k-p₁^m₁是相应的最短距离。这是因为由式(2)可知，p_k^m_k-p₁^m₁是P上从(m₁,1)到(m_k,k)的部分的长度，且任意连接这2个节点的长度至少为p_k^m_k-p₁^m₁。

算法迭代进行，将满足CS条件的二元组(P,p)变换为另一组满足CS条件的二元组。每次迭代时，加入新节点扩展P或删除末端节点收缩P。在收缩操作时，末端节点的价格是严格增加的。当路仅包含起点时会出现退化的情况，此时路只能扩展，或者保持不变的同时p₁^m₁严格增加。

传统拍卖算法只能解决单起点单终点的最短路径问题，而多目标多普勒分离属于多起点多终点问题，本文考虑将此类问题转化为若干个单起点多终点问题。解决多条最短路径问题时，可以共用一个价格向量，因为无论选择哪个起点，CS条件总是满足的,具体步骤如下。

步骤 1 为了消除骨架提取处理在两端产生的失真影响，更准确地确定起始节点和最终节点的位置，适当舍弃S_M×N首尾处的部分时频信息，记截取的时频矩阵为S_M×N。操作时应注意，尽量在不同多普勒曲线的起始节点和最终节点处于分散状态的时刻进行截取。

步骤 2 对所有的起始节点，比较其在S_M×N中对应位置处的能量强度，选择强度最大的点(m₁,1)为第1条路P₁的起点，使所有的最终节点都成为该路的末端节点一次。路径长度_n^ma_n+1^l定义为

式中:λ为比例因子，阈值Δ的表达式为

式中:T_max为求出的各子目标进动周期中的最大值，这是为了兼顾频率变化率较低的节点。

初始解设为P₁=[(m₁,1)]，p_n^m=W(m,n)，W为S_M×N中各时频点的归一化能量强度矩阵。设(m_k,k)为P₁的末端节点，若满足

转到步骤3，否则转到步骤4。

步骤 3 令

同时若(m_k,k)≠(m₁,1)，则收缩P₁，进行下一轮迭代。

步骤 4 加入节点(m′,k+1)来扩展P₁，其中

若(m′,k+1)是最终节点之一，P₁即为第1条最短路径，转步骤5；否则进行下一轮迭代。

步骤 5 将路P₁的起始节点和最终节点置零，转步骤2。当S_M×N中起始节点和最终节点依次配对完成2I次后，算法终止，得到P₁,P₂,…,P_2I。

由上述步骤可以看出，在提取多条最短路径时，预处理可以有效降低噪声干扰，同时大大减少了节点数量，提高了运算效率，Δ保证了路径的连续性，W考虑了瞬时多普勒时频点的峰值性。通过以上操作，就可以较为准确地提取出多目标回波信号不同散射点对应的多普勒分量。

3.   基于小波分析的平动补偿

小波分析法利用小波基函数来逼近信号，其基函数是一组由有限长或快速衰减的小波函数经缩放和平移而生成的函数序列。基于多分辨率分析理论的Mallat算法^[16]是一种最常用的小波基函数构建方法，该算法利用尺度函数和小波函数按照塔式结构对信号实现快速分解和重构，其过程如图 2所示。

图  2  多尺度小波分析

Figure  2.  Multi-scale wavelet analysis

下载: 全尺寸图片 幻灯片

图 2中，子空间序列V_m、W_m为函数空间L²(R)上的一个多尺度分析，满足

对于任意函数f(t)∈V₀，可以在下一级尺度空间V₁和小波空间W₁上进行分解：

式中：p₁ f(t)为逼近部分；q₁ f(t)为细节部分。将逼近部分进一步分解，如此重复就可得到任意尺度上的逼近部分和细节部分。经m次分解，得到

式中：p_m f(t)为函数f(t)的低频全局信息，p_m f(t)对应信号f(t)在分辨率m下的低频趋势项，随着小波分解被单独分离出来；$\mathop \sum \limits_{j = 1}^m$q_j f(t)为逐次分解中分离得到的从V₀到V_m-1各个尺度上f(t)的相应局部细节信息。

对提取出的各多普勒分量进行一维多尺度小波分解，得到的趋势项即为平动分量，由于目标群中各目标具有相同的平动参数，因此对得到的趋势项进行平均可提高平动分量的估计精度。因此，对多普勒分量进行消除趋势项处理就是平动补偿的过程。

小波基函数的选取对使用小波分析法消除趋势项的准确性有很大影响，其数学特性并不能直接反应处理数据的能力，一般通过比较处理信号的实际效果来判定其适用性。文献^[17]构造了消趋误差指数并比较了不同小波基函数的消趋能力，指出sym10的消趋误差最小。因此，本文采用sym10小波基函数实现平动补偿。综上所述，本文方法的流程如图 3所示。

图  3  算法流程图

Figure  3.  Algorithm flowchart

下载: 全尺寸图片 幻灯片

4.   仿真分析

设雷达发射工作频率f_c=10 GHz的单频信号，脉冲重复频率F_s=1 000 Hz，观测时间为3 s。空间中存在2个进动锥体目标Target 1和Target 2。Target 1的参数设置为H₁=2.3 m,h₁=0.7 m,r₁=0.7 m,(α₁,β₁)=(50.6°,40.9°),θ₁=10°,ω₁=2π rad/s; Target 2的参数设置为H₂=2.1 m,h₂=0.5 m,r₂=0.5 m,(α₂,β₂)=(52.4°,38.8°),θ₂=12°,ω₂=4π rad/s。根据弹道特性，平动参数设置为r₀=200 km，v=-2.5 m/s，a₁=1.5 m/s²，a₂=-0.5 m/s²。

利用MATLAB中的awgn函数向原始信号中加入白噪声，设信噪比SNR=5 dB，图 4为经SVD去噪处理后的初始时频图，背景噪声得到了一定地抑制。图 5为预处理后截取的时频骨架。利用改进自相关法得到T_max=1.08 s及I=2，代入式(13)可得阈值Δ=3。图 6为采用本文算法得到的4条最短路径，对应4个散射点的多普勒曲线。

图  4  初始时频图

Figure  4.  Original time-frequency pattern

下载: 全尺寸图片 幻灯片

图  5  截取的时频骨架

Figure  5.  Selective time-frequency skeleton

下载: 全尺寸图片 幻灯片

图  6  提取的最短路径

Figure  6.  Extracted shortest paths

下载: 全尺寸图片 幻灯片

图 7为利用sym10小波基函数对4条曲线进行10层分解得到的平均趋势，图 8为消趋并进行平滑处理后的结果。为了验证本文算法的分离效果，以文献^[5]中的Viterbi算法作对比分析。图 9为利用Viterbi算法对平动补偿后时频图中微多普勒曲线的提取结果，可以看出该算法易出现频率跳变及错误关联，无法准确地估计各频率分量；而本文方法较好地实现了各微多普勒曲线的分离与提取，克服了交叉区域路径选择易出错的问题。这是由于Viterbi算法实质上主要依据能量的强弱来抽取时频信息，不能兼顾变化程度不同的微多普勒曲线，本文方法不仅考虑能量强度的因素，更进一步利用了多普勒曲线的变化规律。

图  7  小波分解得到的平均趋势

Figure  7.  Average trend from wavelet decomposition

下载: 全尺寸图片 幻灯片

图  8  消趋后平滑结果

Figure  8.  Result after eliminating trend and smoothing

下载: 全尺寸图片 幻灯片

图  9  文献^[5] Viterbi算法提取结果

Figure  9.  Extracted result using Viterbi algorithm inRef. ^[5]

下载: 全尺寸图片 幻灯片

进一步分析本文方法的抗噪性能，在不同信噪比条件下进行100次蒙特卡罗仿真，用相对误差(R_R)来衡量提取曲线与理论曲线之间的偏离程度。相对误差的具体表达式为

式中：L为总仿真次数；C^_l(n)为第l次仿真提取的某条微多普勒曲线；C(n)为理论计算得到的该微多普勒曲线。对分离提取出的各条微多普勒曲线的相对误差进行平均，不同方法对应的平均相对误差如图 10所示。可以看出，本方法的平均相对误差明显小于文献^[5]方法；随着信噪比的增加，文献^[5]方法的平均相对误差下降加快，说明Viterbi算法性能受噪声影响较大；当SNR>20 dB时，本文方法的平均相对误差稳定在0.8左右，这是由于骨架提取固有缺陷的限制。

图  10  不同信噪比条件下本文方法与文献^[5] Viterbi算法性能分析

Figure  10.  Performance analysis by proposed method andViterbi algorithm in Ref. ^[5] under different SNRs

下载: 全尺寸图片 幻灯片

5.   结 论

本文通过分析进动目标微多普勒的变化规律，提出了一种基于拍卖算法和小波分析的弹道多目标分离方法。

1) 拍卖算法有效地提取出了各多普勒曲线对应的最短路径，该算法对进动周期的估计精度要求不高，因为T_max仅作为阈值Δ的一个参考参数，路径长度的定义才是影响算法性能的关键；由于构造Δ时认为平动多普勒变化缓慢，忽略了平动的影响，因此本文方法要求粗补偿达到一定精度，使得剩余平动参数的取值在合理的范围内。

2) 利用小波分解的消趋特性，可以在缺乏具体平动模型先验信息的情况下，有效对多普勒曲线的平动分量进行补偿，当小波基函数选定时，分解的层数越高，消趋误差越小，但计算量越大。

3) 若散射点不在整个观测时间段内可见，将会导致提取的骨架出现断裂、不连续，从而得到的最短路径容易出现错误关联，这是由于本文方法依赖于多普勒曲线的连续性变化规律，此时可考虑对数据实现分段提取后再进行融合处理。

4) 本文较好地实现了多目标微多普勒的分离与提取，该方法不仅适用于滑动目标，而且适用于振动、旋转等其他微动形式的弹道目标。如何利用获取的微多普勒信息进行目标分辨与识别，是下一步研究工作的重点。