基于CEMD的旋翼微动目标杂波抑制方法

2018-09-29夏赛强向虎陈文峰杨军陈一畅

航空学报 2018年9期

夏赛强，向虎，陈文峰，杨军，陈一畅

空军预警学院，武汉 430019

当目标中存在旋转、振动、进动和章动等微动形式时，雷达回波将产生除目标刚体部分引起的多普勒频率外的附加频率调制，称为微多普勒现象(Micro-Doppler Effect，MDE)。自2000年Chen[1]提出微多普勒概念以来，对微动目标微多普勒现象的研究迅速发展起来，其中，旋翼旋转作为典型的微动形式，受到国内外学者的广泛研究[2-4]。实际中，旋翼微动目标的背景通常含有杂波，因此需要对杂波进行抑制。目前针对旋翼微动目标的杂波抑制方法大多在时域完成，如文献[5]直接在多普勒域中利用频域滤波的方法对杂波进行了抑制，但是由于该方法采用的是理想滤波器，信号由频域逆变换到时域时会出现拖尾效应。文献[6-7]利用延迟线技术建立动目标显示(Moving Target Indicator,MTI)滤波器，在时域对信号进行对消，达到杂波抑制的效果。文献[8-9]利用杂波相邻距离、方位单元内的相关性来进行杂波对消，即求取时域回波的协方差矩阵，然后利用特征分解和主元分析法获得地杂波的估计值，最后用原始数据减去估计值，达到杂波抑制的效果。

然而，上述时域杂波抑制方法存在回波信号能量损失大的缺点，针对这一问题，本文将回波信号进行时频分析，研究利用复数据经验模态分解(CEMD)进行杂波抑制，从能量的角度选择微动信号的本征模态函数(IMF)分量对微动信号进行重构来抑制旋翼微动目标回波中的地杂波，并与MTI杂波抑制方法进行了对比。理论分析和仿真结果均表明该方法既能够有效保留微动信号，又可以有效抑制杂波。

1 杂波背景下旋翼微动目标信号建模分析

基于窄带雷达，在观测时间较短的情况下，旋翼目标不会出现距离单元走动，若观测时间较长，旋翼目标出现距离单元走动现象，此时对旋翼目标进行运动补偿后，旋翼目标的运动状态可等效为悬停状态[10]。为简化分析，假设目标与雷达处于同一个平面(若不处同一平面，只需要将目标投影到雷达所处平面上)，旋翼目标模型如图1所示，雷达到旋翼目标旋转中心C的距离为RC，假设旋翼微动目标为散射点模型，旋翼某一个叶片上的散射点P到旋转中心的距离为r(0≤r≤l,l为旋翼叶片长度)，点P到雷达的距离为RP，同时点P以角速度ω绕旋转中心C旋转，初始时刻点P的初始旋转角为θ。

设雷达发射信号为窄带线性调频信号[2]

(1)

杂波背景下旋翼叶片上散射点的回波信号可表示为

(2)

式中：K为叶片数目；i为旋翼上的第i个叶片；N为目标单个叶片上的散射点数目；j为叶片上的第j个散射点；σij为散射点的散射系数；c(t)为杂波；n(t)为噪声；c为光速；Rij(tm)为旋翼上散射点到雷达的距离，即[2]

Rij(tm)=RC+rijcos(ωtm+θij)

(3)

杂波c(t)是具有高斯型杂波谱的地杂波，高斯型杂波归一化功率谱密度函数[9]可表示为

(4)

讨论具有高斯型功率谱的杂波谱宽，定义高斯型功率谱的3 dB宽度为杂波的谱宽Sw，即

(5)

从表1中杂波速度起伏展宽值σv的典型值和式(5)来看，地杂波的谱宽较窄，若雷达发射波长λ=0.3 m，那么Sw只有几赫兹到几十赫兹，并且在时频域上，杂波能量和频谱主要集中在零频附近。

式(2)脉冲压缩后的结果可表示为

(6)

表1 σv的典型值Table 1 Typical values of σv

且本文定义的信噪比(SNR)和信杂比(SCR)均为脉压后的信噪比和信杂比。

由多普勒项可以看出，旋翼目标微多普勒相位可表示为

(7)

由式(7)可知，散射点的微多普勒呈现余弦函数形式，而余弦函数的周期与旋翼旋转周期是一致的。且不同叶片上散射点引起的微多普勒曲线有着不同的初相，而同一叶片散射点初相是一致的，差异在于多普勒频率峰值大小的不同。对应式(7)，可知旋翼目标微多普勒频率为

(8)

从式(8)中可以看出，当sin(ωtm+θij)=1且rij=l时，旋翼目标在该散射点的瞬时多普勒频率最大，此时

(9)

假设叶片转动速度ω=6π rad/s，叶片长度l=6 m，雷达发射信号波长λ=0.3 m，此时旋翼目标瞬时微多普勒最大频率fm-Dmax=754 Hz，并且其变化规律受到了正弦的调制。而根据表1中杂波速度起伏值知道地杂波的典型值均小于1 m/s，根据式(5)计算出杂波在时频域的展宽最大值为11.1 Hz，并且分布在零频附近，因此微动目标信号与杂波在时频域上存在展宽以及分布位置上的区别，因此可在时频域上将杂波和微动信号区分开。

2 基于CEMD的杂波抑制方法原理

基于CEMD的旋翼微动目标杂波抑制方法，首先,需要对脉压后的回波进行CEMD处理，得到不同的IMF分量，通过定义的选择方式，选择出高频IMF分量在时域对微动信号进行重构;然后,对重构的信号进行时频分析和过门限处理[12];最终,得到杂波抑制后的微动曲线。处理基本过程如图2所示。

由于杂波和微动目标信号在时频域上存在明显的展宽以及分布位置的区别，对于这种情况，可以利用CEMD方法将杂波和微动目标信号分离开来。

传统的经验模态分解(EMD)可以将复杂的多分量实信号自适应分解为若干个IMF分量的和，对每个IMF分量进行Hilbert变换，可求出瞬时频率和幅值，从而得到信号完整的时频分布[13]。但是对于雷达回波而言，回波都是复信号形式，EMD方法无法对其进行处理，文献[14]针对EMD方法无法处理复信号的缺点，将EMD方法由实信号拓展到复信号，提出了CEMD方法。该方法采用提取零均值的旋转分量来计算信号局部细节，并将信号分解为一系列从快变到慢变的旋转矢量。实质上，CEMD的分解过程就是一系列的滤波过程，将信号中存在的不同频率或趋势逐级分解开来，产生一系列具有不同频率尺度的IMF，并且含高频分量的IMF优先被分离出来，也就是说，CEMD分解的过程就是一种从高频到低频的自适应滤波过程。CEMD分解的示意图如图3所示，能有效地将包含了多个频率分量的信号分解成频率不同的单频信号和一个余量[15]。

由分析可知，同时包含了微动目标回波和杂波的雷达回波信号，如果直接用加窗的方法来分开低频和高频，将会对微动信号产生影响。但是从时频域分析来看，其微动目标回波受到旋转运动的调制，引起的微多普勒频率是从负频率到正频率变化的宽谱，而地杂波通常只在零频附近，存在明显的区别，因此通过CEMD方法可以有效地将杂波和微动目标回波区分开来，然后在时域和时频域进行微动目标信号的重构，从而达到杂波抑制的效果。CEMD方法具体流程[16]如图4所示。

图4中：x(t)为输入信号,记为Si,k(t)；imfi,k-1为模态函数；Pφi,k(t)为在一系列方向φi,k(φi,k=kπ/K, 0≤k≤K)上对imfi,k-1进行投影，并取实部；Sφi,k(t)为将实信号Pφi,k(t)的IMF在φi,k上投影得到的复IMF。

通常，在CEMD分解得到了多个IMF分量之后，通过计算各IMF分量与回波信号的相关系数来判断哪个IMF分量是主杂波，并认为主杂波之后的IMF分量均为杂波分量[15]。基于该方法运算量大的缺点，本文将CEMD分离后的各个模态函数IMF(i)取出，分别求取各个模态函数内信号的能量，并求出各个模态函数占全部模态函数总能量的比重T

(10)

由第1节分析可知，杂波只存在于低频部分，而目标信号不仅仅局限于低频部分，还存在大量的高频分量。进行CEMD处理后，根据IMF(i)从高频到低频的排列方式可对信号进行频率区分，在求出各分量能量比重后，其比重T突变较快的位置都是可能的杂波和目标在频率上的界限，当比重T在高频区域急剧下降或者在低频区域急剧上升的位置即可判定为杂波和目标的界限，在本文方法中，这个界限的经验值通常取IMF(5)。将这个分界IMF作为选择标准来区分微动信号IMF分量和杂波IMF分量，然后将占据能量比重较高的低频杂波IMF分量去除，仅保留具有一定能量的高频微动信号IMF分量对微动目标信号进行重构，即对符合标准的IMF分量求和。重构得到的微动时域信号为

(11)

式中：IMF(j)为选择后的高频模态分量。对重构的微动时域信号进行Gabor变换可以重构出微动目标信号的时频图，Gabor变换属于加窗短时傅立叶变换，是短时Fourier变换中当窗函数取为高斯函数时的一种特殊情况，由于高斯函数满足测不准原理，因此，在所有可能的窗口函数中，高斯窗口函数能得到最好的时频效果。其表达式[17]为

(12)

从上述分析可知，通过分析目标微动信号与杂波的特性和区别，利用CEMD方法分离出代表杂波和微动目标的各个模态函数，再通过微动目标的模态函数重构微动信号，达到抑制杂波、保留微动信号的目的。

3 仿真分析

为了更加清晰得到本文方法对杂波的抑制效果，这里选择3个旋翼叶片顶点和1个旋转中心点作为旋翼散射点，模型如图5所示，其叶片长度为l=6 m，旋翼转速ω=6π rad/s。且本文定义的SCR和SNR均为脉压后的信杂比和信噪比。具体仿真参数如表2所示。

参数PRF/HzTp/μsB/MHzfs/MHzTa/sσv/(m·s-1)数值4 000100120.5120.32

3.1 回波建模仿真

回波建模仿真如图6所示，其中图6(a)为脉压后SNR=5 dB时的回波仿真结果；图6(b)为在图6(a)中添加杂波，使脉压后SCR=0 dB时的回波；图6(c)为图6(b)的距离多普勒表示。

从图6可以看出，在窄带条件下，旋翼目标回波信号脉压后均在相同距离单元内，不会产生距离单元走动，在所有距离单元内都存在杂波和噪声，仿真的地杂波在方位向是相关且是慢起伏的。

3.2 杂波抑制性能仿真与分析

杂波抑制性能仿真如图7所示，其中图7(a)为SCR=0 dB时，旋翼目标回波信号添加杂波后的脉压结果；图7(b)为对图7(a)中目标所在的距离单元进行时频分析的结果；图7(c)为进行CEMD分解后各个IMF所占总能量的比重T之间的关系；图7(d)为杂波抑制后对所有距离单元在时域重构的结果；图7(e)为对图7(d)中微动目标所在的距离单元进行时频分析以及过门限处理的结果；图7(f)为对图7(e)进行过门限处理后的结果。

从图7可以看出，在SCR=0 dB时，添加杂波后对目标所在的距离单元进行时频分析可以清晰完整地得到时频图。其中微动信号的时频曲线呈现余弦变化规律，并且从负频率到正频率占据了很大的谱宽，其值与叶片转速和长度以及雷达工作波长有关，而杂波的时频曲线只是一条在零频附近占据了很小一段谱宽的直线，并且在整个时频域都充满了噪声。得到各个IMF所占总能量的比重后，以IMF(5)为界，可以判断IMF(1)～IMF(4)是旋翼目标信号和噪声的高频IMF，IMF(6)～IMF(11)是包含杂波、噪声和旋翼目标信号的低频部分的低频IMF。在进行CEMD分离后，利用本文方法选择出的高频IMF可以很好地在时域重构出旋翼目标信号，有效地进行了杂波抑制。并且重构出的旋翼目标信号与无杂波情况的旋翼目标信号幅度类似，并且对信噪比有一定的提升，这是因为在进行CEMD后，噪声的低频部分连同杂波可一起被抑制。对重构出的时域旋翼目标信号进行时频分析，可以有效地得到杂波抑制后的旋翼目标信号时频图，但是由于仅取出了高频IMF，导致了旋翼目标信号时频图在低频部分有缺失。使用过门限处理方法和骨架提取方法对旋翼目标信号重构的时频图进行处理可以得到很清晰的骨架图像，能够较为准确地读取微动信息，此时旋翼目标瞬时微多普勒最大频率为759 Hz，与理论值754 Hz比较接近。

3.3 杂波抑制性能对比仿真与分析

为了直观地对本文方法与文献[6]中所提方法(简称方法1)进行杂波抑制性能对比，仿真分为两组进行，第1组是在无噪声情况下进行的，此时SCR=-54 dB；第2组是在噪声背景下进行的，此时SCR=-54 dB，SNR=5 dB。图8和图9分别为第1组仿真和第2组仿真的结果，其中图8(a)和图9(a)为方法1杂波抑制后的时域对消结果，图8(b)和图9(b)为对目标所在的距离单元进行时频分析的结果，图8(c)和图9(c)为本文方法杂波抑制后对全部距离单元在时域进行重构的结果，图8(d)和图9(d)为本文方法杂波抑制后对微动信号在时频域进行重构的结果。

从仿真结果可以看到，经本文方法CEMD杂波抑制后，只存在较少的杂波剩余，并且对信噪比有一定提升，这是由于本文方法在抑制杂波的同时对低频部分的噪声同样具有抑制效果，但是经方法1对杂波进行对消后，信号幅度明显变小且存在明显杂波剩余，并且利用本文方法在时频域重构出的微动信号时频曲线基本保持完整，仅在低频部分出现了少许缺失。而利用方法1进行杂波抑制后，目标所在距离单元的时频分析结果不仅丢失了微动信号时频曲线低频部分的信息，还对不同频率位置的微动曲线有削弱和放大的效果，并且对高频部分的噪声放大效果明显，信噪比恶化严重，这是由于构建的方法1中的滤波器幅频特性导致的。

对比两种杂波抑制方法，在无噪声和有噪声情况下，本文方法均能够很好地抑制杂波，杂波剩余较少，并且能够重构出较为完整的旋翼微动信号，而利用方法1在进行杂波抑制时存在明显杂波剩余情况，并且丢失了大量低频附近的信号，只能重构出部分旋翼微动信号。因此，本文方法在杂波抑制时性能明显优于方法1。