基于STFT和EMD的多普勒信号分离算法

杨彦利，马 德，权建峰

(1.天津工业大学光电检测技术与系统重点实验室，天津 300387；2.西安机电信息技术研究所，陕西 西安 710065)

0 引言

多普勒效应揭示了反射波的频率与物体运动的关系[1]。多普勒信号携带了目标信息，通过检测多普勒信号可以实现对目标的探测。多普勒雷达、多普勒引信等就是利用了多普勒信号来探测和感知目标。多普勒信号的可靠检测与有效识别已成为该体制雷达和引信的核心技术。

当波源离目标很远时，可将目标看成点目标，此时傅里叶变换等方法可实现对多普勒频率的估计；而当波源与目标逐渐靠近时，需用多点目标模型来描述多普勒回波信号[2]，尤其是在近场区，需要考虑目标的体效应[3]。文献[4]通过对实测信号的分析，指出目标效应和噪声干扰对弹目交会信号有较大影响。由于目标和波源的相对运动，多普勒信号的幅度和频率具有时变性。在足够低的信噪比下，只依靠幅度特征已不能有效地对目标信号进行检测和提取[5]。文献[6]研究了经验模态分解(EMD)[7]对多普勒信号的提取技术。文献[8]提出了一种基于快速傅里叶变换(FFT)和自适应滤波并结合Hilbert变换来估计雷达多普勒信号频率的方法。基于时频分布，文献[9]对低信噪比下雷达回波信号瞬时频率的估计进行了有益的尝试。在时频分布上，利用峰值检测与追踪有助于信号瞬时频率的估计[10-11]。利用局部极值点，通过设计合适的滤波器可以将复合信号分解成多个窄带的子信号[12]，文献[13]则基于Wigner-Ville分布提出了一种复合信号的盲成分分离算法。

近场区多普勒回波信号受体目标效应的影响而表现为多个多普勒成分的叠加，准确分离和提取这些成分是进行多普勒频率估计的关键。近场区多普勒信号的识别和瞬时频率的估计目前仍是一个有待解决的问题。本文针对此问题，提出了基于短时傅里叶变换(STFT)和EMD的多普勒信号分离算法。

1 多普勒回波信号模型

在远场区可以用点目标模型来描述多普勒回波信号，可将其表示为[3]

ud=adcos(2πfdt+φd)

(1)

式(1)中，ad为多普勒信号幅值，fd为多普勒信号频率，φd为多普勒信号的初始相位。ad与接收机的探测灵敏度SA、发射波波长λ、天线的方向性系数D、天线的方向图F(φ)、目标的雷达截面积σ和波源与目标的距离R等因素有关，用公式表示为[3]：

(2)

由于式(2)中多个参数随时间变化，导致多普勒信号的幅值具有时变性。对于无线电引信而言，由于发射机和接收机处于同一个载体，故fd可表示为：

fd=±2vr/λ

(3)

式(3)中，vr表示发射机和接收机之间的接近速度。在远场区，多普勒频率基本不变[14]。

当弹目距离逐渐靠近时，需要用多点目标模型来描述多普勒信号，即

(4)

不同散射点反射的多普勒信号在幅值和频率上会有差异，导致整个多普勒回波信号在频谱上具有一定的带宽。

在近场区，体目标效应会影响多普勒回波信号。目标的不同形状，甚至是目标不同部位会具有不同的雷达截面积，再考虑到目标的不同部位导致距离R不同，这些因素会造成反射的多普勒信号幅值出现变化。文献[14]也指出近距离时由于体目标效应显著，多普勒信号幅值会发生畦变。图1所示是某实测信号的STFT时频分布，将STFT变换后的幅值沿时间轴累加求和，得到在不同时刻回波信号的幅值变化曲线，如图2所示。由该图可知，在体目标效应的影响下，多普勒幅值总体上是随着弹目距离接近而增大、反之则减小，但具有很大的波动性。李向军等人[15]建立了近炸引信终点弹道多普勒信号仿真平台，从其仿真的信号以及对比的实测信号也可以发现：从时域上看，在接近目标过程中多普勒信号的幅值总体上是增加的，但有明显的波动。

图1 某实测多普勒信号的时频分布Fig.1 Time-frequency distribution of a measured Doppler signal

图2 STFT变换后幅值沿时间轴的累加和Fig.2 Cumulative sum of amplitudes along the time axis after STFT transformation

对于目标的不同散射点，弹目的接近速度会有所差异，而由式(3)描述的多普勒频率可知，在同一时刻的回波信号中存在多个多普勒成分。近场区，式(4)也可以用来描述多普勒信号，不过需要很多点来描述回波信号。通过图1也可看到近场区多普勒回波信号具有一定的带宽，这为识别多普勒信号和估计多普勒频率带来一些挑战。

2 多普勒信号的分离

根据多普勒信号的特点，设计了一个基于STFT和EMD的多普勒信号分离算法，流程图如图3所示。该算法包括如下几个步骤：首先对采集的多普勒回波信号进行STFT分析；然后将STFT处理后得到的系数沿频率轴进行投影；再根据极值法进行滤波处理得到子信号；将子信号进行EMD处理得到窄带的多普勒信号。

众所周知，给定一个时域信号x(t)，其STFT变换可表示为

(5)

式(5)中，h(τ)是窗函数。计算机处理过程中采用离散STFT变换，得到系数X(tn,fn)。

图3 基于STFT和EMD的多普勒信号分离算法流程图Fig.3 Flow chart of Doppler signal separation algorithm based on STFT and EMD

基于Wigner-Ville分布，可以实现对信号不同成分的分离，但Wigner-Ville分布存在交叉项问题。由于STFT变换没有交叉项，再考虑到在很短一个采样周期内，多普勒信号的频率变化较小，因此，这里采用对STFT的变换结果进行投影，这样还可以减少计算量。对X(tn,fn)沿频率轴方向投影，计算得到

(6)

式(6)中，N表示离散点数。对Y(fn)利用极值法滤波，得到窄带子信号，极值法滤波的流程图如图4所示。为了能够更有效地估计瞬时频率，对滤波得到的子信号进行EMD处理，以使子信号满足本征模态函数(IMF)的要求。

图4 极值法滤波的流程图Fig.4 Flow chart of the extreme filtering method

3 测试分析

对所提出的多普勒信号分离算法用Matlab编程进行测试分析，测试信号为图1所示的实测信号。将STFT系数沿频率轴方向投影，得到如图5所示的一维函数。根据该图的一维函数，利用极值法对原始信号进行滤波，再对滤波结果进行EMD处理，得到一系列子信号，如图6所示，其中图6(b)是图6(a)子信号的频谱图。

为了进一步验证所提算法的有效性，用EMD对原始信号进行分解，分解结果如图7所示。对比图6和图7可知，本文算法得到的子信号带宽更窄，这样更有利于估计信号的瞬时频率。

图5 实测多普勒信号的幅频曲线Fig.5 The amplitude-frequency curve of the measured Doppler signal

图6 本文方法得到的子信号Fig.6 Sub-signals obtained by the proposed method

对于窄带的子信号xi(t)，这里使用下式进行瞬时频率估计[10]：

(7)

瞬时频率的估计结果如图8所示。对比图1的时频分布，可知图8的瞬时频率估计结果反映了回波信号中不同多普勒成分随时间的变化情况。

图7 EMD分解后的子信号Fig.7 Sub-signals decomposed by EMD

利用Hilbert变换来估计瞬时频率也是常用的方法，该方法需要构造解析信号，即[6]

(8)

(9)

对式(9)进行微分，即可得到估计的瞬时频率。经过Hilbert变换后估计的瞬时频率如图9所示， 部分数据如表1所示，由该图可知，Hilbert估计的瞬时频率存在很多尖峰。对比图8和图9可知，虽然图8估计的瞬时频率有些波动，但是相对于图9来说，图8的结果更加合理。

图8 本文估计的瞬时频率Fig.8 The instantaneous frequency estimated by the proposed method

图9 Hilbert法估计的瞬时频率Fig.9 The instantaneous frequency estimated by the Hilbert method

表1 部分Hilbert法估计的瞬时频率值

4 结论

本文提出了基于STFT和EMD的多普勒信号分离算法。该算法根据STFT时频分布，通过极值滤波法并结合EMD将回波信号分解成若干个窄带的子信号，实现了对多普勒信号的分离，进而利用窄带的子信号来估计多普勒信号的瞬时频率。实测验证表明，该算法能够将频率成分相近的多普勒回波信号分解成若干个窄带的子信号，有助于实现对近场区多普勒信号的提取和识别。