目前，国内许多****针对多目标微多普勒特征分离进行了研究，其中，文献[4]指出，多目标的时频图应该是“多目标多散射点”微多普勒的叠加。文献[5]利用Viterbi算法构造时频滤波器，经逆时频变换提取出各信号分量，但该方法只适用于时频交叠程度较弱的多分量微动信号。文献[6]根据对微多普勒点迹不同的相交情况进行关联配对，实现了微多普勒曲线的分离提取，但由于仅利用了斜率和夹角来判断下一时刻点迹，不能很好地解决交叉区域的干扰问题。文献[7]经自适应视野聚类处理得到目标微多普勒的支撑域，通过对食物浓度序列进行匹配选择实现了两旋转目标的分离，但该方法容易误将弱时频分量随背景噪声去除。
本文以典型的无翼弹头类目标为例，针对窄带雷达信号体制下的中段多目标微多普勒特征分离问题展开研究。根据弹道目标多普勒的变化规律，将多条多普勒曲线的提取问题转化为求多条最短路径问题。首先利用骨架处理方法得到时频节点矩阵，再利用改进的拍卖算法迭代提取出多条最短路径，最后利用小波分析的消趋特性对平动分量进行补偿，实现了多分量微多普勒的分离与提取。
1 进动目标多普勒分析 为了保持在大气层外飞行的稳定性和提高命中精度，弹头在中段一般维持进动姿态，本文采用文献[8]描述的滑动型散射中心模型。进动多目标微动模型，如图 1所示。 图中：O-wuv为雷达观测坐标系；Oi-xiyizi为弹头i的滑动坐标系，Oi为弹头i的质心，zi轴为弹头i对应的进动轴的方向，yi轴为zi轴顺时针旋转90°的指向，xi轴符合右手螺旋准则；θi和ωi分别为弹头i的进动角及进动角频率；βi为Si与弹头i自身的旋转轴之间的夹角；ri为弹头i的底面半径；hi为弹头i的质心Oi到底面的距离。
 |
| 图 1 弹道进动多目标微动模型 Fig. 1 Ballistic multi-targets precession model |
| 图选项 |
实际情况中必须考虑散射点的遮蔽问题，雷达的入射电磁波与目标i底面交于一点ui，另一点已被遮蔽，则目标i的顶点Ai、底面交点ui在t时刻处的微距离分别为[8]
 | (1) |
式中：(αi,βi)为观测雷达到弹头i的视线(Si)在Oi-xiyizi中的视角，αi为在Oixiyi平面的投影与xi轴的夹角; aAi、bAi、aui、bui和ci均为调制系数，且只与弹道目标尺寸、进动参数有关，不难推导出
 | (2) |
式中：Hi为弹头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为hi与弹头旋转对称轴的夹角。由于无翼弹头模型为旋转对称目标，可以暂时不考虑弹头的自旋运动对电磁波散射特性的影响[9]。由式(1)和式(2)可知，顶点Ai符合正弦调制规律，而底面交点ui由于受到相位项cos(2ωit)的调制，其对应的微距离近似为多个正弦分量的叠加，且与目标的结构参数以及锥旋频率有关。
由于中段运动只存在地球引力的作用，因此弹道是相对平稳的，可先利用观测时间内某个脉冲测得的速度进行粗补偿[10]，剩余的平动距离分量可近似为3阶多项式：
 | (3) |
式中：r0、v、a1和a2分别为粗补偿后单个雷达波束内包含的目标群在t时刻的径向初始距离、速度、加速度和2阶加速度。
假设雷达发射工作频率为fc的单频信号，经目标散射后得到回波s(t)。通过基带变换后，s(t)变换为sb(t)。对sb(t)的相位项进行求导，得到包含i个目标的多普勒为
 | (4) |
式中：fmi(t)、ftr(t)分别为微多普勒和平动多普勒，i=N+；∪ (·)表示包含关系；c为光速。
在较短脉冲积累时间内可以认为，粗补偿后的弹道目标回波不存在多普勒模糊[10]，且其多普勒带宽小于雷达可探测多普勒带宽，即多普勒曲线在时频面中只出现漂移，不存在折叠现象[11]。
2 基于拍卖算法的多目标分离 2.1 进动周期估计 观察式(4)可知，多目标的瞬时多普勒特征在时频域中表现为一系列总体服从同一分布趋势的交叠曲线，多目标分离的实质就是将各目标包含的各散射点对应的特征曲线单独分离出来。滑动型散射中心产生的微多普勒并不符合正弦模型，不适宜使用参数化提取方法，而传统的一阶条件矩法和峰值位置法均无法处理存在交叉的多分量信号。
待分离曲线可分为两类：一类对应锥顶散射点Ai，另一类对应滑动散射点ui。在相对较短的脉冲积累时间内，同一子目标的ωi近似保持不变，若忽略缓慢变化的平动多普勒，则对于第1类散射点，其单位采样时间间隔内频率的变化可近似表示为
 | (5) |
式中：T为采样时间；Ti为目标i的进动周期；Bs为信号多普勒带宽；Fs为脉冲重复频率；N和M分别为时间采样数和频率采样数；$\left\lceil {} \right\rceil $表示向上取整。
同理，对于第2类散射点，其单位采样时间间隔内频率的变化也近似满足
 | (6) |
由式(5)和式(6)可知，无论哪类散射点，其单位采样时间间隔内频率的变化都在一个与Ti有关的阈值范围之内，这可以作为各曲线上路径选择的依据。本文采用文献[12]中提出的改进自相关法求取多目标回波信号中子目标的个数I不同子目标对应的进动周期Ti，该方法克服了传统自相关法对先验信息要求较高的缺陷，且计算量增大倍数很小，具体步骤见文献[12]。由于模型中每个弹头目标包含2个散射点，时频域中共包含2I条多普勒曲线。
2.2 多目标分离的最短路径描述 雷达录取多目标回波信号后，首先进行预处理。对原始回波进行奇异值分解(SVD)去噪，然后利用短时傅里叶变换(STFT)得到目标回波的时频图，记初始时频矩阵为SM×N。对SM×N进行高斯平滑滤波和二值化后，时频信息只包含0和1的2种元素，利用骨架提取[13]方法进行细化，有效地对多普勒旁瓣进行抑制。最后利用消刺处理对骨架进行修剪，删去伪骨架及孤立点[14]，记预处理后的时频矩阵为SM×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)]，记nman+1l为该路径对应的长度。若确定了各条多普勒曲线的起始节点和最终节点，则多分量瞬时多普勒提取可描述为在带权图D=〈E,G,A〉中求多条最短路径的问题，E、G和A分别为D中的节点集合、路径集合及路径长度集合。
2.3 改进的最短路径拍卖算法 最短路径问题是一个常见的组合优化问题，经典方法包括Dijkstra算法、Kruskal算法和拍卖算法等。其中，拍卖算法适用于大规模稀疏网络最短路径问题的求解，并且便于程序化。该算法最早由Bertsekas和 Castanon[15]提出，通过模拟现实中的拍卖过程来解决指派问题。本文将拍卖算法引入有向图D=〈E,G,A〉中最短路径的求解问题，首先考虑单起点单终点模型。
设P=[(m1,1),(m2,2),…,(mk,k)]为D中的一条路径，若其中包含的节点各不相同，则称P为初等路，节点(mk,k)为该路径的终点，路径长度L可表示为
 | (7) |
该算法在迭代过程中，维护路P并通过扩展和收缩2种操作来修改P。若节点(mk+1,k+1)不在P上，且存在路径[(mk,k),(mk+1,k+1)]，通过(mk+1,k+1)对P的扩展操作将路P变为
 | (8) |
若P包含起点之外的节点，则对P的收缩操作将路P变为
 | (9) |
对D中每个节点(m,n)引入变量pnm，称其为该节点的价格。定义所有节点的价格组成的价格向量为p。算法维护价格向量p和路P，使其满足如下性质：
 | (10) |
 | (11) |
式(10)和式(11)称为互补松弛(CS)条件。可以证明，若二元组(P,p)满足CS条件，则路P上从节点(m1,1)到任意节点(mk,k)的部分是从(m1,1)到(mk,k)的最短路经，而pkmk-p1m1是相应的最短距离。这是因为由式(2)可知，pkmk-p1m1是P上从(m1,1)到(mk,k)的部分的长度，且任意连接这2个节点的长度至少为pkmk-p1m1。
算法迭代进行，将满足CS条件的二元组(P,p)变换为另一组满足CS条件的二元组。每次迭代时，加入新节点扩展P或删除末端节点收缩P。在收缩操作时，末端节点的价格是严格增加的。当路仅包含起点时会出现退化的情况，此时路只能扩展，或者保持不变的同时p1m1严格增加。
传统拍卖算法只能解决单起点单终点的最短路径问题，而多目标多普勒分离属于多起点多终点问题，本文考虑将此类问题转化为若干个单起点多终点问题。解决多条最短路径问题时，可以共用一个价格向量，因为无论选择哪个起点，CS条件总是满足的,具体步骤如下。
步骤 1?为了消除骨架提取处理在两端产生的失真影响，更准确地确定起始节点和最终节点的位置，适当舍弃SM×N首尾处的部分时频信息，记截取的时频矩阵为SM×N。操作时应注意，尽量在不同多普勒曲线的起始节点和最终节点处于分散状态的时刻进行截取。
步骤 2?对所有的起始节点，比较其在SM×N中对应位置处的能量强度，选择强度最大的点(m1,1)为第1条路P1的起点，使所有的最终节点都成为该路的末端节点一次。路径长度nman+1l定义为
 | (12) |
式中:λ为比例因子，阈值Δ的表达式为
 | (13) |
式中:Tmax为求出的各子目标进动周期中的最大值，这是为了兼顾频率变化率较低的节点。
初始解设为P1=[(m1,1)]，pnm=W(m,n)，W为SM×N中各时频点的归一化能量强度矩阵。设(mk,k)为P1的末端节点，若满足
 | (14) |
转到步骤3，否则转到步骤4。
步骤 3 令
 | (15) |
同时若(mk,k)≠(m1,1)，则收缩P1，进行下一轮迭代。
步骤 4 加入节点(m′,k+1)来扩展P1，其中
 | (16) |
若(m′,k+1)是最终节点之一，P1即为第1条最短路径，转步骤5；否则进行下一轮迭代。
步骤 5 将路P1的起始节点和最终节点置零，转步骤2。当SM×N中起始节点和最终节点依次配对完成2I次后，算法终止，得到P1,P2,…,P2I。
由上述步骤可以看出，在提取多条最短路径时，预处理可以有效降低噪声干扰，同时大大减少了节点数量，提高了运算效率，Δ保证了路径的连续性，W考虑了瞬时多普勒时频点的峰值性。通过以上操作，就可以较为准确地提取出多目标回波信号不同散射点对应的多普勒分量。
3 基于小波分析的平动补偿 小波分析法利用小波基函数来逼近信号，其基函数是一组由有限长或快速衰减的小波函数经缩放和平移而生成的函数序列。基于多分辨率分析理论的Mallat算法[16]是一种最常用的小波基函数构建方法，该算法利用尺度函数和小波函数按照塔式结构对信号实现快速分解和重构，其过程如图 2所示。
 |
| 图 2 多尺度小波分析 Fig. 2 Multi-scale wavelet analysis |
| 图选项 |
图 2中，子空间序列Vm、Wm为函数空间L2(R)上的一个多尺度分析，满足
 | (17) |
对于任意函数f(t)∈V0，可以在下一级尺度空间V1和小波空间W1上进行分解：
 | (18) |
式中：p1 f(t)为逼近部分；q1 f(t)为细节部分。将逼近部分进一步分解，如此重复就可得到任意尺度上的逼近部分和细节部分。经m次分解，得到
 | (19) |
式中：pm f(t)为函数f(t)的低频全局信息，pm f(t)对应信号f(t)在分辨率m下的低频趋势项，随着小波分解被单独分离出来；$\mathop \sum \limits_{j = 1}^m $qj f(t)为逐次分解中分离得到的从V0到Vm-1各个尺度上f(t)的相应局部细节信息。
对提取出的各多普勒分量进行一维多尺度小波分解，得到的趋势项即为平动分量，由于目标群中各目标具有相同的平动参数，因此对得到的趋势项进行平均可提高平动分量的估计精度。因此，对多普勒分量进行消除趋势项处理就是平动补偿的过程。
小波基函数的选取对使用小波分析法消除趋势项的准确性有很大影响，其数学特性并不能直接反应处理数据的能力，一般通过比较处理信号的实际效果来判定其适用性。文献[17]构造了消趋误差指数并比较了不同小波基函数的消趋能力，指出sym10的消趋误差最小。因此，本文采用sym10小波基函数实现平动补偿。综上所述，本文方法的流程如图 3所示。
 |
| 图 3 算法流程图 Fig. 3 Algorithm flowchart |
| 图选项 |
4 仿真分析 设雷达发射工作频率fc=10 GHz的单频信号，脉冲重复频率Fs=1 000 Hz，观测时间为3 s。空间中存在2个进动锥体目标Target 1和Target 2。Target 1的参数设置为H1=2.3 m,h1=0.7 m,r1=0.7 m,(α1,β1)=(50.6°,40.9°),θ1=10°,ω1=2π rad/s; Target 2的参数设置为H2=2.1 m,h2=0.5 m,r2=0.5 m,(α2,β2)=(52.4°,38.8°),θ2=12°,ω2=4π rad/s。根据弹道特性，平动参数设置为r0=200 km，v=-2.5 m/s，a1=1.5 m/s2，a2=-0.5 m/s2。
利用MATLAB中的awgn函数向原始信号中加入白噪声，设信噪比SNR=5 dB，图 4为经SVD去噪处理后的初始时频图，背景噪声得到了一定地抑制。图 5为预处理后截取的时频骨架。利用改进自相关法得到Tmax=1.08 s及I=2，代入式(13)可得阈值Δ=3。图 6为采用本文算法得到的4条最短路径，对应4个散射点的多普勒曲线。
 |
| 图 4 初始时频图 Fig. 4 Original time-frequency pattern |
| 图选项 |
 |
| 图 5 截取的时频骨架 Fig. 5 Selective time-frequency skeleton |
| 图选项 |
 |
| 图 6 提取的最短路径 Fig. 6 Extracted shortest paths |
| 图选项 |
图 7为利用sym10小波基函数对4条曲线进行10层分解得到的平均趋势，图 8为消趋并进行平滑处理后的结果。为了验证本文算法的分离效果，以文献[5]中的Viterbi算法作对比分析。图 9为利用Viterbi算法对平动补偿后时频图中微多普勒曲线的提取结果，可以看出该算法易出现频率跳变及错误关联，无法准确地估计各频率分量；而本文方法较好地实现了各微多普勒曲线的分离与提取，克服了交叉区域路径选择易出错的问题。这是由于Viterbi算法实质上主要依据能量的强弱来抽取时频信息，不能兼顾变化程度不同的微多普勒曲线，本文方法不仅考虑能量强度的因素，更进一步利用了多普勒曲线的变化规律。
 |
| 图 7 小波分解得到的平均趋势 Fig. 7 Average trend from wavelet decomposition |
| 图选项 |
 |
| 图 8 消趋后平滑结果 Fig. 8 Result after eliminating trend and smoothing |
| 图选项 |
 |
| 图 9 文献[5] Viterbi算法提取结果 Fig. 9 Extracted result using Viterbi algorithm inRef. [5] |
| 图选项 |
进一步分析本文方法的抗噪性能，在不同信噪比条件下进行100次蒙特卡罗仿真，用相对误差(RR)来衡量提取曲线与理论曲线之间的偏离程度。相对误差的具体表达式为
 | (20) |
式中：L为总仿真次数；C^l(n)为第l次仿真提取的某条微多普勒曲线；C(n)为理论计算得到的该微多普勒曲线。对分离提取出的各条微多普勒曲线的相对误差进行平均，不同方法对应的平均相对误差如图 10所示。可以看出，本方法的平均相对误差明显小于文献[5]方法；随着信噪比的增加，文献[5]方法的平均相对误差下降加快，说明Viterbi算法性能受噪声影响较大；当SNR>20 dB时，本文方法的平均相对误差稳定在0.8左右，这是由于骨架提取固有缺陷的限制。
 |
| 图 10 不同信噪比条件下本文方法与文献[5] Viterbi算法性能分析 Fig. 10 Performance analysis by proposed method andViterbi algorithm in Ref. [5] under different SNRs |
| 图选项 |
5 结 论 本文通过分析进动目标微多普勒的变化规律，提出了一种基于拍卖算法和小波分析的弹道多目标分离方法。
1) 拍卖算法有效地提取出了各多普勒曲线对应的最短路径，该算法对进动周期的估计精度要求不高，因为Tmax仅作为阈值Δ的一个参考参数，路径长度的定义才是影响算法性能的关键；由于构造Δ时认为平动多普勒变化缓慢，忽略了平动的影响，因此本文方法要求粗补偿达到一定精度，使得剩余平动参数的取值在合理的范围内。
2) 利用小波分解的消趋特性，可以在缺乏具体平动模型先验信息的情况下，有效对多普勒曲线的平动分量进行补偿，当小波基函数选定时，分解的层数越高，消趋误差越小，但计算量越大。
3) 若散射点不在整个观测时间段内可见，将会导致提取的骨架出现断裂、不连续，从而得到的最短路径容易出现错误关联，这是由于本文方法依赖于多普勒曲线的连续性变化规律，此时可考虑对数据实现分段提取后再进行融合处理。
4) 本文较好地实现了多目标微多普勒的分离与提取，该方法不仅适用于滑动目标，而且适用于振动、旋转等其他微动形式的弹道目标。如何利用获取的微多普勒信息进行目标分辨与识别，是下一步研究工作的重点。
参考文献
| [1] | DAVID L R. Ballistic missile defense[J].Journal of Electronic Defense, 2006, 29(1): 46–52. |
| [2] | 冯德军, 徐乐涛, 艾小锋. 空间复杂目标群的雷达目标识别技术[J].现代防御技术, 2015, 43(4): 1–6.FENG D J, XU L T, AI X F. Radar recognition technique for complex target groups in space[J].Modern Defence Technology, 2015, 43(4): 1–6.(in Chinese) |
| [3] | CHEN V C.Analysis of radar micro-Doppler signature with time-frequency transform[C]//Proceedings of IEEE Workshop on Statistical Signal and Array Processing.Piscataway,NJ:IEEE Press,2000:463-466. |
| [4] | 胡晓伟, 童宁宁, 胡国平, 等. 基于微多普勒的弹道多目标分离方法[J].系统工程与电子技术, 2015, 37(8): 1734–1740.HU X W, TONG N N, HU G P, et al. Multi-ballistic targets resolution based on micro-Doppler[J].Systems Engineering and Electronics, 2015, 37(8): 1734–1740.(in Chinese) |
| [5] | LI P, WANG D C, WANG L. Separation of micro-Doppler signals based on time frequency filter and Viterbi algorithm[J].Signal, Image and Video Processing, 2013, 7(3): 593–605.DOI:10.1007/s11760-011-0263-3 |
| [6] | 赵盟盟, 张群, 罗迎, 等. 点迹-曲线关联算法的旋转对称群目标分辨[J].空军工程大学学报(自然科学版), 2015, 16(2): 43–48.ZHAO M M, ZHANG Q, LUO Y, et al. Distinguishing of rotationally symmetric group targets based on plot-curve association algorithm[J].Journal of Air Force Engineering University (Natural Science Edition), 2015, 16(2): 43–48.(in Chinese) |
| [7] | 李靖卿, 冯存前, 张栋. 基于自适应视野聚类匹配的多目标分离与提取[J].系统工程与电子技术, 2015, 37(9): 1974–1979.LI J Q, FENG C Q, ZHANG D. Multi-target separation and extraction based on adaptive vision cluster matching[J].Systems Engineering and Electronics, 2015, 37(9): 1974–1979.(in Chinese) |
| [8] | MA L, LIU J, WANG T, et al. Micro-Doppler character of sliding-type scattering center on rotationally symmetric target[J].Science China (Information Sciences), 2011, 54(9): 1957–1967.DOI:10.1007/s11432-011-4254-3 |
| [9] | 雷腾, 刘进忙, 杨少春, 等. 基于三站一维距离像融合的弹道目标特征提取方法研究[J].宇航学报, 2011, 32(2): 228–234.LEI T, LIU J M, YANG S C, et al. Study on feature extraction method of ballistic target based on three-station range profiles[J].Journal of Astronautics, 2011, 32(2): 228–234.(in Chinese) |
| [10] | 贺思三, 赵会宁, 张永顺. 基于延迟共轭相乘的弹道目标平动补偿[J].雷达学报, 2014, 3(5): 505–510.HE S S, ZHAO H N, ZHANG Y S. Translational motion compensation for ballistic targets based on delayed conjugated multiplication[J].Journal of Radars, 2014, 3(5): 505–510.(in Chinese) |
| [11] | 胡晓伟, 童宁宁, 董会旭, 等. 弹道中段群目标平动补偿与分离方法[J].电子与信息学报, 2015, 37(2): 291–296.HU X W, TONG N N, DONG H X, et al. Translation compensation and resolution of multi-ballistic targets in midcourse[J].Journal of Electronics & Information Technology, 2015, 37(2): 291–296.(in Chinese) |
| [12] | 肖立, 周剑雄, 何峻, 等. 弹道中段目标进动周期估计的改进自相关法[J].航空学报, 2010, 31(4): 812–818.XIAO L, ZHOU J X, HE J, et al. Improved autocorrelation method for precession period estimation of ballistic target in midcourse[J].Acta Aeronautica et Astronautica Sinica, 2010, 31(4): 812–818.(in Chinese) |
| [13] | 罗迎, 柏又青, 张群, 等. 弹道目标平动补偿与微多普勒特征提取方法[J].电子与信息学报, 2012, 34(3): 602–608.LUO Y, BAI Y Q, ZHANG Q, et al. Translational motion compensation and micro-Doppler feature extraction of ballistic targets[J].Journal of Electronics & Information Technology, 2012, 34(3): 602–608.(in Chinese) |
| [14] | 赵永嘉, 戴树岭. 基于图像骨架和贪婪算法的无人机航路规划[J].北京航空航天大学学报, 2010, 36(4): 474–477.ZHAO Y J, DAI S L. Unmanned aircraft vehicle path planning based on image skeleton and greedy algorithm[J].Journal of Beijing University of Aeronautics and Astronautics, 2010, 36(4): 474–477.(in Chinese) |
| [15] | BERTSEKAS D P, CASTANON D A. The auction algorithm for the transportation problem[J].Annals of Operations Research, 1989, 20(1): 67–96.DOI:10.1007/BF02216923 |
| [16] | STEPHANE G M. A theory for multi-resolution signal decomposition:The wavelet representation[J].IEEE Transactions on Pattern Analysis and Machine Intelligence, 1989, 11(7): 674–693.DOI:10.1109/34.192463 |
| [17] | 吴志成, 王重阳, 任爱君. 消除信号趋势项时小波基优选方法研究[J].北京理工大学学报, 2013, 33(8): 811–814.WU Z C, WANG C Y, REN A J. Optimal selection of wavelet base functions for eliminating signal trend based on wavelet analysis[J].Transactions of Beijing Institute of Technology, 2013, 33(8): 811–814.(in Chinese) |
