基于匈牙利算法的多分量微动信号提取方法

王 欢 李 峰 宋思盛 姜 洋

(1.西安电子工程研究所 西安 710100;2.空军工程大学 西安 710077)

0 引言

在雷达探测弹道目标场景下，弹道导弹在飞行过程中除了平动飞行以外，同时会绕自身对称轴做自旋运动并有可能绕某一定向轴作锥旋运动。这些精细的运动会对雷达回波产生除平动多普勒以外的附加频率调制。这一现象被称作微多普勒效应[1]。微多普勒为目标识别[2]和雷达成像[3]的研究提供了有效信息。在过去的十年中，窄带雷达弹道目标微动特征提取已得到了深入的研究[4]。

在实际应用中，雷达回波采样的缺损和微动信号在时频域交叠已经成为阻碍微动特征有效提取的两个突出问题。在回波受到强干扰期间，采样的数据不得不被丢弃以免影响后续的数据处理，当雷达出现故障和异常值时也会影响回波信号的有效采集。尽管基于稀疏重构的方案为机动目标和稳定运动目标提供了缺损数据恢复的解决途径，比如贝叶斯压缩感知[5]和自适应稀疏分数模糊函数[6]，然而由于微动信号的瞬时频率会随时间快速变化，上述方法不能有效地应用于微动目标。文献[7]利用一种非参数联合时频字典从损坏的微动信号中重构出聚焦良好的时频分布，但是没有提供针对具有多散射点微动目标信号的提取方法。多散射点微动目标的信号提取实质上可以当作多分量微动信号分解问题处理，一般情况下，基于多分量微动信号分解的参数化方法都会提前设定微动的瞬时多普勒频率多项式曲线[8]或者正弦曲线[9]模型，但是弹道目标微动信号有可能呈现出更为复杂的调制特性(例如滑动散射点的回波信号)，需要设置多个参数来拟合调制特性，从而增加了算法的复杂度。而大多数非参数化的方法，例如变分模态分解[10]、结合短时傅里叶变换的同步挤压变换[11]、稀疏短时傅里叶变换[12]，都无法处理微动信号在时频域出现交叠的情况，然而弹道目标的多分量微动信号在时频域常常出现交叠的情况。文献[13]提出一种岭迹重排的方法提取时频域的瞬时频率曲线，并结合本征线性调频分量分解(ICCD)算法[14]有效分解了在时频域分解重叠的微动信号，但是在时频分布因采样缺损出现散焦和缺损的情况下，岭迹重排算法无法有效运行；短时变分模态分解(STVMD)[15]算法可以在瞬时频率不规则调制和低信噪比情况下正确分解具有交叠状态的微动信号，但是当微动信号在时频域完全重叠时，分解效果不佳。

为了解决上述问题，本文提出一种采样缺失条件下的多分量微动信号分解方法，首先将采样缺失微动信号的时频分布重构问题转换为一个稀疏正则化问题，并采用软阈值迭代算法(ISTA)求解；其次提取重构时频分布的极大值点，将极大值点关联为瞬时频率轨迹的问题构造为指派问题；最后采用卡尔曼滤波和匈牙利算法求解，得到每个散射点的瞬时频率轨迹，实现了雷达采样缺失条件下对多分量微动信号的分解。

1 信号回波模型

在高频电磁条件下，金属目标可以由点散射模型表示[7]，回波信号缺损条件下的微动目标窄带雷达回波信号表示为

(1)

2 微动信号分解

2.1 重构时频表征

s(t)是一个典型的非平稳信号。时频分析是一种在时频域表征非平稳信号能量的有效方法，非平稳信号可以在时频域上表征出频域中无法表示的特征。短时傅里叶变换[16]具备线性和没有交叉项等诸多优势，因此本文采用短时傅里叶变换分析回波信号。短时傅里叶变换可以表示为

(2)

其中TF0(t,f)代表时频表征；f代表瞬时频率；gσ代表长度为σ的窗函数，用来保证加窗信号足够平稳，窗函数有多种形式，本文使用高斯窗，它可以很容易地用其他形式重写，比如矩形窗和汉明窗等。为了方便描述重构算法，结合式(1)，将式(2)转换为矩阵的形式[12]

(3)

其中H=diag{W,…,W}；WN×N是傅里叶变换矩阵；N=Mms,其中ms=1是滑动窗口的步数；Q=[P1…Pm…PM]T是滑动的窗函数，其中Pm=diag{pm(1),…pm(n),…pm(N)},n∈[1,N]是对角矩阵；B=diag{b(1),…b(n),…b(N)},b(n)∈{0,1}是表示信号缺损的对角矩阵；没有缺损的理想回波信号为s0=[s0(1),…s0(n),…s0(N)]T，(·)T代表转置。根据式(3)，采样缺损时的信号为

s=BGf+n

(4)

其中G是HQ的伪逆矩阵；n=[η(1),…η(n),…η(N)]T是加性噪声。为了从幅度起伏和采样缺损的回波s中得到理想的时频表征f，令s=[s(1),…s(n),…s(N)]T，式(4)可以被描述成一个优化问题

(5)

(6)

其中‖·‖1表示1范数运算符。L1正则化模型可以被认为是一个套索回归问题，可以采用ISTA[17]求解。ISTA是一种近端梯度下降法，邻近算子为

(7)

(8)

soft(z,λ1)=sign(z)max{0,|z|-λ1}

(9)

由图1所示，通过软阈值迭代法，得到了缺失采样条件下的稀疏、去噪且锐化的时频谱分布。

图1 缺失采样信号和重构的时频分布

2.2 瞬时频率轨迹的提取

将重构的时频谱分布f取模并重排为由N个时刻频谱组成的时频图TF(t,f)。TF(t,f)中极大值点的位置坐标为qt,g(t)=(t,fg(t))，其中t=t0,…,tN-1，fg(t)代表第t时刻的第g(t)个极值点在频谱上的坐标，t时刻的频谱总共有G(t)个极值点。极大值点的坐标qt,g(t)由式(10)得到[18]

(10)

其中α代表阈值参数。为了避免交叉点附近极值点对后续算法的影响，采用式(11)消除交叉点的极值点

|fg(t)-fg_other(t)|≥df

(11)

其中fg_other(t)代表第t个频谱除了g(t)以外的其他极大值点在频谱上的位置，将不满足式(11)的极大值点剔除，得到qt,g(t)，其在时频平面的位置由图2所示，在本文中df=N/20。

图2 提取的极大值点在时频平面的位置

受到Simple online and realtime tracking (Sort)算法[19]的启发，引入卡尔曼滤波[20]和匈牙利算法[21]相结合的多目标跟踪算法将时频平面的极大值点关联成每个散射点的瞬时频率轨迹。选取G(tm)=P时刻的极大值点qtm,k作为卡尔曼滤波的初始测量值，瞬时频率的位置和速度由线性状态空间描述。为了简化分析，选用匀速(CV)模型，第k个散射点瞬时频率轨迹在tm+1时刻的状态向量为

(12)

(13)

(14)

其中ygk=0代表第g(tm+1)个测量值与第k测量值不匹配；ygk=1代表g(tm+1)个测量值与第k测量值匹配。指派问题的核心问题是如何构造效率矩阵，如式(14)所示的指派问题，其效率矩阵Cm+1可以通过高斯核密度函数构造。

(15)

图3 各散射点的瞬时频率轨迹

3 仿真实验

为了证明本文所提算法的有效性，通过计算机仿真进行验证，并对仿真结果进行分析。设雷达的载频为10GHz，脉冲重复频率为512Hz，观测时间为1s，在加性高斯白噪声环境下，信噪比为10dB，短时傅里叶变换的高斯窗函数长度σ=53。一般情况下，参数λ1的值越大，ISTA算法生成的时频分布就越稀疏；参数α的值越大，从时频图中提取的极大值点越稀疏；参数L的值越大，ICCD算法时频滤波器组的宽度越大。根据经验，当参数值λ1设置在0.01～0.04之间，α设置在0.4～0.6之间，L设置在5～15之间，λ2设置在4 ～ 6之间时算法可以有效运行。在本文中参数λ、α、L、λ2分别设置为0.01、0.5、5、5。

将本文算法应用于基于滑动散射点的微动信号分解中。当弹道目标在空间上进动时，一些散射点的位置随着目标的运动而滑动，称为“滑动散射点”[22]。设弹道目标的锥旋转频率为0.5Hz，锥底半径为0.5m，锥旋角为10°,从锥顶到圆锥质量中心的距离为2m，从圆锥质量中心到圆锥底部中心的距离为0.5m，雷达入射角度为45°。数据随机缺损40%，基于滑动散射模型的缺损采样信号的时频分布如图4(a)所示，可见时频分布是破损且散焦的；结合ICCD算法[22]，本文算法分解缺损采样条件下的滑动散射点信号分量如图4(b)、图4(c)、图4(d)所示。结果表明，该方法能够有效地分解出各滑动散射点的微动信号分量。

图4 基于滑动散射模型的微动信号分解

最后，采用参数δ来评价两种方法的估计精度

(25)

表1列出了两种方法在采样缺损30%和40%条件下分解理想散射目标回波和滑动散射目标回波的精度，可以看出：一方面，采样缺损率越大，两种方法分解得到的微动信号精度越低，这是因为随着采样缺损率的增加，时频图的散焦程度更加严重，导致了两种算法的性能下降；另一方面，本文方法比STVMD算法分解信号的精度平均提高了1dB左右，这是因为当某个散射点信号与其他散射点信号在时频域交叠时，STVMD算法可能会将该散射点的信号能量分解到其他散射点的回波中，导致出现交叠时刻信号丢失的情况。

表1 不同方法分解微动信号的精度

4 结束语

本文提出了一种微动信号采样缺损和回波时频域交叠情况下的多分量微动目标回波分解新方法。首先在构建缺损采样多分量微动信号稀疏正则化模型的基础上，采用ISTA算法对时频分布进行重构；其次结合卡尔曼滤波和匈牙利算法将重构时频分布极大值关联为瞬时频率轨迹；最后仿真结果表明本文相对于STVMD算法方法分解得到的微动信号精度更高。应该指出的是，本文方法作为时频分布的一种后处理方法，其性能受限于时频分布的分辨率，因此需要选择合适的时频分辨率以保证算法的有效运行；其次，与许多已有的多分量信号分解方法相同，本文方法也需要预先设定待分解信号中信号分量的个数。在下一步工作中，我们将重点寻找更有效的字典矩阵以重构时频分布，并改进现有算法使其在不需要预先设定信号分量个数的前提下有效运行。