CN106842181A - 基于经验模态分解的空间锥体目标微动特征提取方法 - Google Patents

Abstract

本发明公开了一种基于经验模态分解的空间锥体目标微动特征提取方法，主要解决现有技术特征提取易失败的问题。其方案是：1.根据窄带线性调频信号模型计算一个脉冲重复周期内发射信号序列；2.根据发射信号和接收到目标回波信号，构造脉压信号矩阵，并根据该矩阵，构造锥体目标的多普勒回波信号；3.根据计算出的锥体目标的多普勒回波信号，利用经验模态分解得到多个本征模态函数；4.根据得到的多个本征模态函数，重构锥体目标散射中心的多普勒信号；5.根据重构出的散射中心信号，构造锥体目标散射中心时频图；6.从时频图中提取出目标的微动特征。本发明在保证准确提取微动特征的同时，提高了特征提取的效率，可用于目标识别。

Description

基于经验模态分解的空间锥体目标微动特征提取方法

技术领域

本发明属于雷达技术领域，特别涉及一种空间锥体目标的微动特征提取方法，可用于目标识别。

背景技术

空间锥体目标或部件在其运行轨道上运动的同时通常伴随着旋转、振动等一些微小的运动，这些微小的运动称做微动。运动目标通常伴有微动，例如人在行走时四肢的摆动、人体呼吸和心跳引起的胸腔振动、车辆在前进时车轮的转动、直升机旋翼的旋转、天线的旋转、发动机工作时的振动等。微动现象首先在相干激光雷达中被观测到，并被应用于测振。2000年，美国海军实验室的V.C.Chen教授将其引入到微波雷达的研究中，并将这种由微动引起的时变多普勒频移定义为雷达中的微多普勒。微多普勒反应了目标微动结构的动力学特征，如果能够加以利用，则可广泛应用到目标识别和分类中。目前微动特征已经成功应用于空中目标、地面车辆、行人、弹道导弹等目标的分类和识别中。

如前所述，若想要获得空间锥体的微动特征，须准确提取其微多普勒信息。但是在实际工作环境中，接收到的窄带锥体目标回波通信号通常是由主体上多个散射中心的回波信号的叠加组成。由于信号在时域和频率相互叠加和干扰，为了提取锥体目标微动特征，有必要分离和提取单分量微多普勒信号。目前，基于窄带微动特征提取方法主要有以下两种：

1.基于时频图的微动特征提取方法。这种方法首先将时域信号转化至时频域，从时频图中分离出目标上各个散射中心的微多普勒，其代表方法是基于多目标追踪的微动特征提取方法，此类方法缺点是提取出的微多普勒具有较低的时频分辨率。

2.参数化模型的微动特征提取方法。该方法的思想是采用特定的参数模型对信号进行描述，这个过程是利用一些目标的先验信息估计出模型参数，其代表方法是基于时变自回归模型的微动特征提取方法。

由于雷达接收到的空间锥体目标的回波信号是锥体上各个散射中心回波信号的叠加，从而回波信号的时频图中各个散射中心的时频曲线相互交叉和重叠，所以难以对各个散射中心的时频曲线在时频图中进行直接提取。基于时频图的微动特征提取方法都需解决频点关联问题，但是频点关联存在两方面缺点，一方面，频点关联易在不同散射中心的时频曲线的交叉点产生错误，导致微动特征提取失败；另一方面，时频关联会导致计算复杂度增加，降低了微动特征提取的效率。

发明内容

本发明的目的在于对上述现有技术的不足，提出一种基于经验模态分解的空间锥体目标微动特征提取方法，以避免由频点关联问题导致的微动特征提取失败问题，提高特征提取效率。

本发明的技术是：利用经验模态分解EMD算法，对空间锥体目标的回波信号在时域中进行经验模态EMD分解，从而得到多个本征模态函数IMF，利用本征模态函数IMF重构各个散射中心的多普勒信号，再将这些信号从时域转化到时频域以提取各个散射中心的微动特征，避免了由频点关联导致的微动特征提取失败的问题，减少计算时间，提高特征提取效率。其实现步骤包括如下：

1、一种基于经验模态分解的空间锥体目标微动特征提取方法，包括：

(1)根据窄带线性调频信号模型，得出一个脉冲重复周期内发射信号：sf＝[sf₁sf₂ ... sf_d ... sf_N]，其中，sf_d代表第d个信号，d∈[1,N]，N表示一个脉冲重复周期内的采样次数；

(2)根据接收到的空间锥体目标回波信号：sr＝{s₁,s₂,...,s_i,...,s_M}，计算s(t)与s_i脉冲压缩后的信号：其中，表示卷积运算，s_i代表第i个脉冲重复周期的回波序列，i∈[1,M]，M表示接收到的脉冲的个数；

(3)将M个压缩后的信号合成一个脉压信号矩阵：sm＝[sm₁ sm₂ … sm_i … sm_M]，提取sm每列中幅度最大的信号x_i，构造锥体目标的多普勒回波信号:x＝{x₁,x₂,...,x_i,...,x_M}；

(4)对锥体目标的多普勒回波信号x进行经验模态分解，得到L个本征模态函数IMF和剩余信号r_L：

其中，IMF_n为第n个本征模态函数，n∈[1,L]；

(5)利用分解得到的L个本征模态函数IMF重构锥体目标的锥顶散射中心多普勒信号s_A和锥底散射中心多普勒信号s_B：

(6)对(5)重构出的散射中心多普勒信号分别进行短时傅里叶变换，得到锥体目标锥顶散射中心时频图P_A和锥底散射中心的时频图P_B：

P_A＝STFT(s_A)，P_B＝STFT(s_B)，

其中，STFT(·)表示短时傅里叶变换；

(7)分别从锥顶散射中心时频图P_A和锥底散射中心的时频图P_B中提取锥顶散射中心的时频曲线f₁和锥底散射中心的时频曲线f₂，完成锥体目标的微动特征提取。

本发明具有如下优点：

1)本发明有效地利用了经验模态分解对接收信号中不同尺度的波动进行分解，该分解过程通过筛选实现，其作用类似于具有不同频率通带的滤波器组，由于这种通带滤波器不需要预先设定好参数，而是基于数据自身特征进行自适应地调整，因此可以自适应将锥体目标多个散射点的叠加信号分离开来；同时由于本发明利用经验模态分解在时域完成对不同散射中心信号的分离，避免了由频点关联问题导致的微动特征提取失败问题，提高了特征提取效率。

2)本发明使用本征模态函数对锥顶和锥底散射中心信号进行重构，利用第一本征模态函数重构锥顶散射中心信号，利用除第一本征模态函数外的其余本征模态函数的和函数重构锥底散射中心信号，由于第一本征模态函数占居接收信号的绝大部分能量，而剩余的本征模态函数占居接收信号的很少一部分能量，因此可以保证重构出的各个散射中心的信号具有足够的能量，以体现出各个散射中心完整清晰的微动特征。

附图说明

图1为本发明的实现流程图；

图2为本发明提取空间锥体进动目标微动的特征的结果图；

图3为本发明提取空间锥体章动目标微动特征的结果图。

具体实施方式

下面结合附图对本发明作进一步说明：

本发明所针对的一收发同置雷达，雷达发射载频为f_c，脉冲宽度为T_p，调频率为γ的窄带线性调频信号，并在锥体目标前进方向下方进行观测。

参照图1，为本发明的实现步骤如下：

步骤1：根据窄带线性调频信号模型计算发射信号sf。

1a)窄带线性调频信号模型如下：

其中，rect(·)表示单位矩形信号，exp(·)表示指数运算，j表示虚数单位，t为采样时刻，γ为调频率，f_c为发射载频，T_p为脉冲宽度；

1b)根据窄带线性调频信号模型,计算一个脉冲重复周期内发射信号：sf＝[sf₁sf ₂... sf_d ... sf_N]，sf_d代表第d个信号，d∈[1,N]，N表示一个脉冲重复周期内的采样次数。

步骤2：获取脉冲压缩信号矩阵。

2a)雷达接收空间锥体目标回波信号：sr＝{s₁,s₂,...,s_i,...,s_M}，其中，s_i代表接收到i个脉冲，i∈[1,M]，M表示接受到的脉冲的个数。

2b)对空间锥体目标回波信号sr和雷达发射信号sf进行脉冲压缩：其中，sm_i是第i个脉冲进行脉冲压缩后的信号。

2a)将M个压缩后的信号合成一个脉冲压缩信号矩阵：sm＝[sm₁ sm₂ … sm_i …sm_M]。

步骤3：利用步骤2中脉压信号矩阵sm，构造锥体目标多普勒回波信号x。

3a)提取脉压信号矩阵sm每列中幅度最大的信号x_i：

3b)构造接收到锥体目标的多普勒回波信号：x＝{x₁,x₂,...,x_i,...,x_M}。

步骤4：对步骤3中锥体目标的多普勒回波信号x进行经验模态分解。

4a)定义临时信号temp＝x；

4b)对临时信号temp进行逐点搜索，记录下临时信号temp的所有极值点；

4c)对于临时信号temp的极小值点和极大值点进行插值，得到临时信号temp的下包络e_min和上包络e_max；

4d)计算临时信号temp的包络均值：

4e)从临时信号temp中减去包络均值，得到新的信号x_new＝x-e_av，并更新临时信号为：temp＝x_new；

重复步骤4b)至4e)，直到新的信号x_new成为本征模函数，得到第一个本征模函数IMF₁＝x_new和剩余信号r₁＝x-IMF₁；

4f)对剩余信号r₁重复4a)到4e)直至剩余信号满足r_L≤SD,其中，SD为事先设定的阈值，最终将多普勒回波信号x分解为剩余信号r_L以及L个本征模函数IMF_n，n∈[1,L]；

步骤5：根据步骤4中的本征模函数，重构出锥体目标锥顶散射中心的信号s_A和锥底散射中心的信号s_B：

步骤6：利用步骤5的结果，得到散射中心的时频图。

6a)对锥体目标的锥顶散射中心信号s_A进行短时傅里叶变换，得到锥顶散射中心的时频图：P_A＝STFT(s_A)，其中，STFT(·)表示短时傅里叶变换；

6b)对锥体目标的锥底散射中心信号s_B进行短时傅里叶变换，得到锥底散射中心的时频图：P_B＝STFT(s_B)。

步骤7：提取时频图中散射中心的时频曲线。

7a)从锥顶散射中心时频图P_A中取每列的最大值f₁(i)，得到由M个最大值组成的序{f₁(1),f₁(2),...,f₁(i),...,f₁(M)}，把该组最大值序列按顺序进行连接，得到锥顶散射中心的时频曲线f₁；

7b)从锥底散射中心时频图P_B中取每列的最大值f₂(i)，得到由M个最大值组成序列：{f₂(1),f₂(2),...,f₂(i),...,f₂(M)}，将这组最大值序列按顺序进行连接，得到锥底散射中心的时频曲线f₂。

至此完成对空间锥体目标的时频曲线的提取。

本发明效果通过以下实验进一步说明。

实验一

1.1)实验参数

使用一收发同置窄带雷达发射线性调频信号。

设雷达参数为：载频10GHZ，带宽1MHZ，脉冲重复频率500HZ，脉冲宽度5μm。

设锥体目标参数为：锥体高1m，底面半径r＝0.25m。

设雷达视线俯仰角60°，方位角90°，雷达驻留时间2s，锥体目标初始姿态角为:[0°90° 0°]，锥体进动参数为：自旋频率0.8HZ，锥旋频率3HZ，锥旋角10°。

1.2)仿真内容

在上述实验条件下用本发明方法提取锥体目标进动状态下锥顶和锥底的时频曲线，结果如图2，其中：

图2(a)是原始接收信号的时频图。

图2(b)是重构的锥底散射中心信号的时频图。

图2(c)是重构的锥顶散射中心信号的时频图。

图2(d)是提取出的锥顶散射中心和锥底散射中心的时频曲线。

由图2可见，在允许一定误差范围内，本发明在目标锥体进动状态下可以较准确提取时频曲线。

实验二

2.1)实验参数：

使用一收发同置窄带雷达发射线性调频信号。

设雷达参数：载频10GHZ，带宽1MHZ，脉冲重复频率500HZ，脉冲宽度5μm。

设锥体目标参数：锥体高1m，底面半径r＝0.25m。

设雷达视线俯仰角60°，方位角90°，雷达驻留时间2s。锥体目标初始姿态角:[0°90° 0°]，锥体进动参数为：自旋频率0.8HZ，锥旋频率3HZ，锥旋角10°，摆动频率0.3HZ，摆动幅度5°。

2.2)仿真内容

在上述2.1)实验条件下，用本发明方法提取锥体目标章动状态下锥顶和锥底的时频曲线，结果如图3，其中：

图3(a)是原始接收信号的时频图，

图3(b)是重构的锥底散射中心信号的时频图，

图3(c)是重构的锥顶散射中心信号的时频图，

图3(d)是提取出的锥顶散射中心和锥底散射中心的时频曲线。

从图3可见，在允许一定误差范围内，本发明在目标锥体章动状态下可以较准确提取时频曲线。

以上描述仅是本发明的一个具体实例，并未构成对本发明的任何限制，显然对于本领域的专业人员来说，在了解了本发明内容和原理后，都可能在不背离本发明原理、结构的情况下，进行形式和细节上的各种修正和改变，但是这些基于本发明思想的修正和改变仍在本发明的权利要求保护范围之内。

Claims (4)

1.一种基于经验模态分解的空间锥体目标微动特征提取方法，包括：

(1)根据窄带线性调频信号模型，得出一个脉冲重复周期内发射信号：sf＝[sf₁sf₂...sf_d...sf_N]，其中，sf_d代表第d个信号，d∈[1,N]，N表示一个脉冲重复周期内的采样次数；

(2)根据接收到的空间锥体目标回波信号：sr＝{s₁,s₂,...,s_i,...,s_M}，计算s(t)与s_i脉冲压缩后的信号：其中，表示卷积运算，s_i代表第i个脉冲重复周期的回波序列，i∈[1,M]，M表示接收到的脉冲的个数；

(3)将M个压缩后的信号合成一个脉压信号矩阵：sm＝[sm₁ sm₂ ... sm_i ... sm_M]，提取sm每列中幅度最大的信号x_i，构造锥体目标的多普勒回波信号：x＝{x₁,x₂,...,x_i,...,x_M}；

(4)对锥体目标的多普勒回波信号x进行经验模态分解，得到L个本征模态函数IMF和剩余信号r_L：

$x=\underset{n=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{IMF}}_n+r_L$

其中，IMF_n为第n个本征模态函数，n∈[1,L]；

(5)利用分解得到的L个本征模态函数IMF重构锥体目标的锥顶散射中心多普勒信号s_A和锥底散射中心多普勒信号s_B：

s_A＝IMF₁，

(6)对(5)重构出的散射中心多普勒信号分别进行短时傅里叶变换，得到锥体目标锥顶散射中心时频图P_A和锥底散射中心的时频图P_B：

P_A＝STFT(s_A)，P_B＝STFT(s_B)，

其中，STFT(·)表示短时傅里叶变换；

(7)分别从锥顶散射中心时频图P_A和锥底散射中心的时频图P_B中提取锥顶散射中心的时频曲线f₁和锥底散射中心的时频曲线f₂，完成锥体目标的微动特征提取。

2.根据权利要求1所述的方法，其中步骤(1)中的窄带线性调频信号模型，表示如下：

$s\left(t\right)=rect\left(\frac{t}{T_p}\right)\exp \left(2\mathrm{\π}j\right(f_{c}t+0.5{\mathrm{\γt}}^2\left)\right)$

其中，rect(·)表示单位矩形信号，T_p为脉冲宽度，exp(·)表示指数运算，f_c是雷达发射机的载频，j表示虚数单位，γ为线性调频信号的调频率，t为采样时刻。

3.根据权利要求1所述的方法，其中步骤(7)中的从锥顶散射中心时频图P_A提取锥顶散射中心的时频曲线f₁，是从锥顶散射中心时频图P_A中取每列的最大值f₁(i)，得到由M个最大值组成序列：{f₁(1),f₁(2),...,f₁(i),...,f₁(M)}，把这些最大值按顺序连接得到锥顶散射中心的时频曲线f₁，其中i∈[1,M]，M为接收到的脉冲的个数。

4.根据权利要求1所述的方法，其中步骤(7)中的从锥顶散射中心时频图P_B提取的锥顶散射中心的时频曲线f₂，是从锥底散射中心时频图P_B中取每列的最大值f₂(i)，得到由M个最大值组成序列：{f₂(1),f₂(2),...,f₂(i),...,f₂(M)}，将些最大值按顺序连接得到锥顶散射中心的时频曲线f₂，i∈[1,M]，M为接收到的脉冲的个数。