一种雷达目标运动特征提取方法及系统

1.本发明涉及信号处理技术领域，特别是涉及一种雷达目标运动特征提取方法及系统。


背景技术：

2.雷达是信息获取与感知领域中的重要系统，具有全天时、全天候、作用距离远等独特优势，利用雷达系统获取空天、地面、海洋等环境中目标的回波信号，进而通过噪声抑制、成像、特征分析等信号处理工具，获取目标的几何形状、运动状态、目标类型、工作状态等关键信息，雷达目标特征提取能够为目标的辨识与解译提供条件。
3.雷达目标的运动特征是辨识目标属性的关键要素，包括平动、自旋、进动等具体形式，反映了目标的动态特性。研究表明，雷达目标运动特征提取精度不仅依赖于雷达系统本身，还与雷达目标特性以及精细信号处理等方面紧密相关。现有目标运动特征方法主要是针对目标单一运动、低速运动、不存在耦合运动等低复杂度情况下研发设计的，雷达目标运动多样化、高速化、耦合化等新的特点对雷达目标运动特征提取带来诸多问题与挑战。
4.首先，雷达目标运动更加多样，导致现有运动特征提取方法难以兼顾多种复合运动。例如，失去了姿态控制能力而失效的卫星，由于空间碎片撞击或电力供应不足等原因，卫星表现为旋转、翻滚的运动形式，且转动速度会慢慢变小，使得运动具有非均匀性。对于部分空间目标，由于目标在撞击或分离过程产生的横向作用力，还可能产生进动。总的来说，空间目标除轨道运动外还会产生自旋、翻滚、进动、非均匀转动等多种运动形式。其次，雷达目标运动高速化，导致雷达回波信号质量下降。例如，导弹等超音速运动雷达目标的高速运动，导致雷达回波误差加大、同时信号信噪比变低，直接影响后续运动特征提取。最后，雷达目标多散射点导致的运动耦合，导致运动特征提取容易出现误判。例如目标多个散射点做同样的自旋运动，容易导致运动特征提取的旋转周期加倍，进而使得运动特征提取失败。
5.面对上述瓶颈问题，迫切需要发展新的雷达目标运动特征提取方法。


技术实现要素：

6.为了克服现有技术的不足，本发明的目的是提供一种雷达目标运动特征提取方法及系统。
7.为实现上述目的，本发明提供了如下方案：
8.一种雷达目标运动特征提取方法，包括：
9.对雷达一维距离像数据进行预处理，得到目标的相位信号；
10.利用经验模态分解对所述相位信号进行信号分离，得到多个本征模态函数；
11.基于傅里叶分析的频谱分析方法，对所述本征模态函数分别进行运动参数估计，得到雷达目标运动特征。
12.优选地，所述对雷达一维距离像数据进行预处理，得到目标的相位信号，包括：
13.对所述雷达一维距离像x∈cm×n计算绝对值，再对所述绝对值进行逐列求和，得到求和参数其中，x为所述雷达一维距离像，y(n)为所述求和参数，c为复数域，m为一维距离像时间序列长度，n为一维距离像距离单元个数，m为一维距离像时间序列索引，n为一维距离像距离单元索引，n＝1,2,
…
,n，y∈rn；
14.设置阈值设置阈值t1∈r，针对y∈rn，统计所述求和参数的指标集并根据所述指标集提取强散射单元的时序信号；其中γ为所述指标集，γ中的元素为雷达一维距离像x∈cm×n中强散射单元所在的列，k为一维距离像距离单元索引；
15.令y＝arg(x)表示所述时序信号的相位数据；其中y＝arg(x)为提取复数相位的算子，且y∈rm；
16.根据所述时序信号的相位数据计算相位梯度d(m)＝y(m)-y(m+1)；其中m＝1,2,
…
,m-1；
17.根据所述时序信号的相位数据计算缠绕相位差：e(m)＝arctan(sin(y(m))/cos(y(m)))，其中m＝1,2,
…
,m-1；
18.令z(1)＝y(1)；z(1)为初始化数据；
19.根据公式z(m)＝y(m)+e(m-1)对缠绕相位求和，以得到解缠后的所述相位信号；其中z(m)为解缠后相位信号的第m个元素。
20.优选地，所述利用经验模态分解对所述相位信号进行信号分离，得到多个本征模态函数，包括：获取第一个本征模态函数、获取其它本征模态函数和设置停止条件；
21.所述获取第一个本征模态函数的具体步骤为：
22.计算所述相位信号z∈rm的局部极大值点坐标集i与局部极小值点坐标集j；
23.利用最邻近插值对局部极大值点{z(i)|i∈i}进行插值获取上包络信号a∈rm，并利用最邻近插值对局部极小值点{z(j)|j∈j}进行插值获取下包络信号b∈rm；
24.计算所述上包络信号与所述下包络信号的均值包络信号m1＝(a+b)/2；
25.令h1＝z-m1，判断h1是否为本征模态函数，若否，则令h1替代z，并跳转至步骤“计算所述相位信号z∈rm的局部极大值点坐标集i与局部极小值点坐标集j”，直到获得第一个本征模态函数；
26.所述获取其它本征模态函数的具体步骤为：
27.设rk＝r
k-1-hk，且令rk替代z，重复步骤“获取第一个本征模态函数”，直到获得第k+1个本征模态函数h
k+1
，其中k＝2,3,
…
；
28.所述设置停止条件的具体步骤为：
29.设置停止阈值t2＞0，计算判别式
30.若sd≤t2，则停止继续进行经验模态分解，此时记k＝k-1。
31.优选地，所述基于傅里叶分析的频谱分析方法，对所述本征模态函数分别进行运动参数估计，得到雷达目标运动特征，包括：
32.对所述本征模态函数依次进行傅里叶变换和取绝对值处理，得到频谱数据fk＝|fhk|；其中，f为离散傅里叶矩阵，所述离散傅里叶矩阵的元素为f
ij
＝exp(-2πj0(i-1)(j-1)/m)，其中，i,j∈{1,2,
…
,m}且j0为虚数单位，fk为所述频谱数据；hk为所述本征模态函数；
33.取所述频谱数据的最大值对应的频率值作为运动参数的估计，得到雷达目标运动特征wk＝argmaxfk,k＝1,2,
…
,k；其中wk为所述雷达目标运动特征。
34.一种雷达目标运动特征提取系统，包括：
35.预处理单元，用于对雷达一维距离像数据进行预处理，得到目标的相位信号；
36.信号分离单元，用于利用经验模态分解对所述相位信号进行信号分离，得到多个本征模态函数；
37.参数估计单元，用于基于傅里叶分析的频谱分析方法，对所述本征模态函数分别进行运动参数估计，得到雷达目标运动特征。
38.优选地，所述预处理单元包括：
39.求和模块，用于对所述雷达一维距离像x∈cm×n计算绝对值，再对所述绝对值进行逐列求和，得到求和参数其中，x为所述雷达一维距离像，y(n)为所述求和参数，c为复数域，m为一维距离像时间序列长度，n为一维距离像距离单元个数，m为一维距离像时间序列索引，n为一维距离像距离单元索引，n＝1,2,
…
,n，y∈rn；
40.提取模块，用于设置阈值设置阈值t1∈r，针对y∈rn，统计所述求和参数的指标集并根据所述指标集提取强散射单元的时序信号；其中γ为所述指标集，γ中的元素为雷达一维距离像x∈cm×n中强散射单元所在的列，k为一维距离像距离单元索引；
41.相位表示模块，用于令y＝arg(x)表示所述时序信号的相位数据；其中y＝arg(x)为提取复数相位的算子，且y∈rm；
42.第一计算模块，用于根据所述时序信号的相位数据计算相位梯度d(m)＝y(m)-y(m+1)；其中m＝1,2,
…
,m-1；
43.第二计算模块，用于根据所述时序信号的相位数据计算缠绕相位差：e(m)＝arctan(sin(y(m))/cos(y(m)))，其中m＝1,2,
…
,m-1；
44.初始化模块，用于令z(1)＝y(1)；z(1)为初始化数据；
45.第三计算模块，用于根据公式z(m)＝y(m)+e(m-1)对缠绕相位求和，以得到解缠后的所述相位信号；其中z(m)为解缠后相位信号的第m个元素。
46.优选地，所述信号分离单元包括：
47.第一获取模块，用于获取第一个本征模态函数；
48.第二获取模块，用于获取其它本征模态函数；
49.停止模块，用于设置停止条件；
50.所述第一获取模块包括：
51.极值计算子模块，用于计算所述相位信号z∈rm的局部极大值点坐标集i与局部极小值点坐标集j；
52.插值子模块，用于利用最邻近插值对局部极大值点{z(i)|i∈i}进行插值获取上包络信号a∈rm，并利用最邻近插值对局部极小值点{z(j)|j∈j}进行插值获取下包络信号b∈rm；
53.均值计算子模块，用于计算所述上包络信号与所述下包络信号的均值包络信号m1＝(a+b)/2；
54.第一获取子模块，用于令h1＝z-m1，判断h1是否为本征模态函数，若否，则令h1替代z，并跳转至步骤“计算所述相位信号z∈rm的局部极大值点坐标集i与局部极小值点坐标集j”，直到获得第一个本征模态函数；
55.所述第二获取模块包括：
56.第二获取子模块，用于设rk＝r
k-1-hk，且令rk替代z，重复步骤“获取第一个本征模态函数”，直到获得第k+1个本征模态函数h
k+1
，其中k＝2,3,
…
；
57.所述停止模块包括：
58.判别式计算子模块，用于设置停止阈值t2＞0，计算判别式
59.判断子模块，用于若sd≤t2时，停止继续进行经验模态分解，此时记k＝k-1。
60.优选地，所述参数估计单元包括：
61.傅里叶处理模块，用于对所述本征模态函数依次进行傅里叶变换和取绝对值处理，得到频谱数据fk＝|fhk|；其中，f为离散傅里叶矩阵，所述离散傅里叶矩阵的元素为f
ij
＝exp(-2πj0(i-1)(j-1)/m)，其中，i,j∈{1,2,
…
,m}且j0为虚数单位，fk为所述频谱数据；hk为所述本征模态函数；
62.估计模块，用于取所述频谱数据的最大值对应的频率值作为运动参数的估计，得到雷达目标运动特征wk＝argmaxfk,k＝1,2,
…
,k；其中wk为所述雷达目标运动特征。
63.根据本发明提供的具体实施例，本发明公开了以下技术效果：
64.本发明提供了一种雷达目标运动特征提取方法及系统，所述方法包括对雷达一维距离像数据进行预处理，得到目标的相位信号；利用经验模态分解对所述相位信号进行信号分离，得到多个本征模态函数；基于傅里叶分析的频谱分析方法，对所述本征模态函数分别进行运动参数估计，得到雷达目标运动特征。本发明能够适应雷达目标运动更加多样、高速运动、运动耦合等实际情况，提取的目标运动特征具有更高的精度。
附图说明
65.为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所
需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
66.图1为本发明提供的实施例中的方法流程图；
67.图2为本发明提供的实施例中的雷达目标运动特征提取流程图；
68.图3为本发明提供的实施例中的本征模态函数获取流程图；
69.图4为本发明提供的实施例中的雷达目标结构图；
70.图5为本发明提供的实施例中的相位解缠前后的第一幅相位信号图；
71.图6为本发明提供的实施例中的相位解缠前后的第二幅相位信号图；
72.图7为本发明提供的实施例中的经验模态分解结果。
具体实施方式
73.下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
74.在本文中提及“实施例”意味着，结合实施例描述的特定特征、结构或特性可以包含在本技术的至少一个实施例中。在说明书中的各个位置出现该短语并不一定均是指相同的实施例，也不是与其它实施例互斥的独立的或备选的实施例。本领域技术人员显式地和隐式地理解的是，本文所描述的实施例可以与其它实施例相结合。
75.本技术的说明书和权利要求书及所述附图中的术语“第一”、“第二”、“第三”和“第四”等是用于区别不同对象，而不是用于描述特定顺序。此外，术语“包括”和“具有”以及它们任何变形，意图在于覆盖不排他的包含。例如包含了一系列步骤、过程、方法等没有限定于已列出的步骤，而是可选地还包括没有列出的步骤，或可选地还包括对于这些过程、方法、产品或设备固有的其它步骤元。
76.本发明的目的是提供一种雷达目标运动特征提取方法及系统，能够适应雷达目标运动更加多样、高速运动、运动耦合等实际情况，提取的目标运动特征具有更高的精度。
77.为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
78.图1和图2分别为本发明提供的实施例中的方法流程图和雷达目标运动特征提取流程图，如图1和图2所示，本发明提供了一种雷达目标运动特征提取方法，包括：
79.步骤100：对雷达一维距离像数据进行预处理，得到目标的相位信号；
80.步骤200：利用经验模态分解对所述相位信号进行信号分离，得到多个本征模态函数；
81.步骤300：基于傅里叶分析的频谱分析方法，对所述本征模态函数分别进行运动参数估计，得到雷达目标运动特征。
82.优选地，所述步骤100包括：
83.对所述雷达一维距离像x∈cm×n计算绝对值，再对所述绝对值进行逐列求和，得到
求和参数其中，x为所述雷达一维距离像，y(n)为所述求和参数，c为复数域，m为一维距离像时间序列长度，n为一维距离像距离单元个数，m为一维距离像时间序列索引，n为一维距离像距离单元索引，n＝1,2,
…
,n，y∈rn；
84.设置阈值设置阈值t1∈r，针对y∈rn，统计所述求和参数的指标集并根据所述指标集提取强散射单元的时序信号；其中γ为所述指标集，γ中的元素为雷达一维距离像x∈cm×n中强散射单元所在的列，k为k为一维距离像距离单元索引；
85.令y＝arg(x)表示所述时序信号的相位数据；其中y＝arg(x)为提取复数相位的算子，且y∈rm；
86.根据所述时序信号的相位数据计算相位梯度d(m)＝y(m)-y(m+1)；其中m＝1,2,
…
,m-1；
87.根据所述时序信号的相位数据计算缠绕相位差：e(m)＝arctan(sin(y(m))/cos(y(m)))，其中m＝1,2,
…
,m-1；
88.令z(1)＝y(1)；z(1)为初始化数据；
89.根据公式z(m)＝y(m)+e(m-1)对缠绕相位求和，以得到解缠后的所述相位信号；其中z(m)为解缠后相位信号的第m个元素。
90.本实施例中第一个步骤为“数据预处理”步骤，设x∈cm×n为经过距离对准后的雷达一维距离像，数据预处理主要实现从雷达一维距离像数据中获取目标的相位信号，用于后续经验模态分解与运动特征提取。
91.其中具体包括：
92.s1.1：提取强散射单元
93.对雷达一维距离像x∈cm×n取绝对值，再逐列求和，即：
[0094][0095]
其中，n＝1,2,
…
,n，y∈rn。设置阈值t1∈r，针对y∈rn，统计其指标集：
[0096][0097]
则γ中的元素即为雷达一维距离像x∈cm×n中强散射单元所在的列。为了便于讨论，假设γ中仅有一个元素n0，此时令x(m)＝x(m,n0)表示强散射单元对应时序信号，其中m＝1,2,
…
,m。若γ为空集，则调整阈值t；若γ中元素数量大于1个，则分别针对每个元素重复进行仅有1个元素时的相同处理步骤即可。
[0098]
s1.2：相位解缠
[0099]
对于提取的强散射单元的时序信号x∈cm，令y＝arg(x)表示该信号的相位信号，
其中arg(
·
)为提取复数相位的算子，且y∈rm。由于直接提取相位其相位值会位于(-π,π]之间，需要利用相位解缠求解完整的相位信号，相位解缠的主要步骤如下：
[0100]
1)计算相位梯度：d(m)＝y(m)-y(m+1)，其中m＝1,2,
…
,m-1；
[0101]
2)计算缠绕相位差：e(m)＝arctan(sin(y(m))/cos(y(m)))，其中m＝1,2,
…
,m-1；
[0102]
3)初始化：令z(1)＝y(1)；
[0103]
4)对缠绕相位求和，得到解缠后的相位信号：z(m)＝y(m)+e(m-1)，其中m＝1,2,
…
,m-1；
[0104]
经过相位解缠后的相位信号z∈rm可用于后续信号分离与运动特征提取。
[0105]
优选地，所述步骤200包括：获取第一个本征模态函数、获取其它本征模态函数和设置停止条件；
[0106]
所述获取第一个本征模态函数的具体步骤为：
[0107]
计算所述相位信号z∈rm的局部极大值点坐标集i与局部极小值点坐标集j；
[0108]
利用最邻近插值对局部极大值点{z(i)|i∈i}进行插值获取上包络信号a∈rm，并利用最邻近插值对局部极小值点{z(j)|j∈j}进行插值获取下包络信号b∈rm；
[0109]
计算所述上包络信号与所述下包络信号的均值包络信号m1＝(a+b)/2；
[0110]
令h1＝z-m1，判断h1是否为本征模态函数，若否，则令h1替代z，并跳转至步骤“计算所述相位信号z∈rm的局部极大值点坐标集i与局部极小值点坐标集j”，直到获得第一个本征模态函数；
[0111]
所述获取其它本征模态函数的具体步骤为：
[0112]
设rk＝r
k-1-hk，且令rk替代z，重复步骤“获取第一个本征模态函数”，直到获得第k+1个本征模态函数h
k+1
，其中k＝2,3,
…
；
[0113]
所述设置停止条件的具体步骤为：
[0114]
设置停止阈值t2＞0，计算判别式
[0115]
若sd≤t2，则停止继续进行经验模态分解，此时记k＝k-1。
[0116]
本实施例中的第二个步骤为“经验模态分解”，对于相位信号z∈rm，首先利用经验模态分解实现不同运动形式对应信号的分离，具体步骤如下：
[0117]
s2.1：获取第一个本征模态函数
[0118]
获取第一个本征模态函数的流程图见图3，具体步骤如下：
[0119]
1)计算信号z∈rm的局部极大值点坐标集i与局部极小值点坐标集j；
[0120]
2)利用最邻近插值对局部极大值点{z(i)|i∈i}进行插值获取上包络信号a∈rm；同时利用最邻近插值对局部极小值点{z(j)|j∈j}进行插值获取下包络信号b∈rm；
[0121]
3)计算上包络信号与下包络信号的均值包络信号m1＝(a+b)/2；
[0122]
4)令h1＝z-m1；若h1为本征模态函数，则算法停止；否则令h1替代z，重复上述步骤1)-4)，直到获得第一个本征模态函数。
[0123]
第4)步中本征模态函数要满足的必要条件为：其极值点的个数p与过零点的个数q满足q-1≤p≤q+1。
[0124]
s2.2获取其它本征模态函数
[0125]
设信号z∈rm的第一个本征模态函数为h1，则令r1＝z-h1，且令r1替代z，重复上述步骤s2.1，直到获得第二个本征模态函数h2。
[0126]
令r2＝r
1-h2，且令r2替代z，重复上述步骤s2.1，直到获得第三个本征模态函数h3。
[0127]
以此类推，不妨设rk＝r
k-1-hk，且令rk替代z，重复上述步骤s2.1，直到获得第k+1个本征模态函数h
k+1
，其中k＝2,3,
…
。
[0128]
s2.3停止条件
[0129]
设置停止阈值t2＞0，计算判别式：
[0130][0131]
若sd≤t2，则停止继续进行经验模态分解，此时记k＝k-1，即将相位信号z∈rm共分解成了k个本征模态函数。
[0132]
优选地，所述步骤300包括：
[0133]
对所述本征模态函数依次进行傅里叶变换和取绝对值处理，得到频谱数据fk＝|fhk|；其中，f为离散傅里叶矩阵，所述离散傅里叶矩阵的元素为f
ij
＝exp(-2πj0(i-1)(j-1)/m)，其中，i,j∈{1,2,
…
,m}且j0为虚数单位，fk为所述频谱数据；hk为所述本征模态函数；
[0134]
取所述频谱数据的最大值对应的频率值作为运动参数的估计，得到雷达目标运动特征wk＝argmaxfk,k＝1,2,
…
,k；其中wk为所述雷达目标运动特征。
[0135]
具体的，本实施例中第三个步骤为运动特征提取，针对进行经验模态分解后的本征模态函数，分别进行运动参数估计，进而实现运动特征提取。在此主要利用基于傅里叶分析的频谱分析方法实现运动参数的估计。具体步骤如下：
[0136]
s3.1：频谱分析
[0137]
考虑本征模态函数hk∈rm，对其进行傅里叶变换，并取绝对值，即：
[0138]fk
＝|fhk|。
[0139]
其中，k＝1,2,
…
,k，f为离散傅里叶矩阵，其元素为f
ij
＝exp(-2πj0(i-1)(j-1)/m)，其中，i,j∈{1,2,
…
,m}且j0为虚数单位。
[0140]
s3.2：运动参数估计
[0141]
取fk的最大值对应的频率值作为运动参数的估计：
[0142]
wk＝argmaxfk,k＝1,2,
…
,k。
[0143]
此时，获得的雷达目标运动特征即为{w1,w2,
…
,wk}。
[0144]
本实施例还提供了另一具体实施方式，其中设雷达参数设置为：载频9.5ghz，带宽1ghz，脉冲重复频率300hz，其中目标为仿真的四个理想散射点，如图4所示，在雷达参数下仿真生成目标回波数据，并进行脉冲压缩与距离对齐，生成一维距离像x∈c
5000
×
1024
作为本方法的输入信号。具体步骤有：
[0145]
s1：数据预处理
[0146]
以经过距离对准后的雷达一维距离像x∈c
5000
×
1024
为输入，数据预处理主要实现从
雷达一维距离像数据中获取目标的相位信号，用于后续经验模态分解与运动特征提取。
[0147]
s1.1：提取强散射单元
[0148]
对雷达一维距离像x∈c
5000
×
1024
取绝对值，再逐列求和，即
[0149][0150]
其中，n＝1,2,
…
,1024，y∈r
1024
。设置阈值t1＝0.3，针对y∈r
1024
，统计其指标集：
[0151][0152]
则雷达一维距离像x∈c
5000
×
1024
中强散射单元所在的列为第503列与第518列。为了便于讨论，本例中仅考虑第503列强散射单元，第518列强散射单元采用类似的处理。此时令x(m)＝x(m,503)表示强散射单元对应时序信号，其中m＝1,2,
…
,5000。
[0153]
s1.2：相位解缠
[0154]
对于提取的强散射单元的时序信号x∈c
5000
，令y＝arg(x)表示该信号的相位信号，其中arg(
·
)为提取复数相位的算子，且y∈r
5000
。由于直接提取相位其相位值会位于(-π,π]之间，需要利用相位解缠求解完整的相位信号，相位解缠的主要步骤如下：
[0155]
1)计算相位梯度：d(m)＝y(m)-y(m+1)，其中m＝1,2,
…
,4999；
[0156]
2)计算缠绕相位差：e(m)＝arctan(sin(y(m))/cos(y(m)))，其中m＝1,2,
…
,4999；
[0157]
3)初始化：令z(1)＝y(1)；
[0158]
4)对缠绕相位求和，得到解缠后的相位信号：z(m)＝y(m)+e(m-1)，其中m＝1,2,
…
,4999；
[0159]
经过相位解缠后的相位信号z∈r
5000
如图5和图6所示，可用于后续信号分离与运动特征提取。
[0160]
s2：经验模态分解
[0161]
对于相位信号z∈r
5000
，首先利用经验模态分解实现不同运动形式对应信号的分离，具体步骤如下：
[0162]
s2.1：获取第一个本征模态函数
[0163]
按照图3，获取第一个本征模态函数；
[0164]
s2.2获取其它本征模态函数
[0165]
设信号z∈rm的第一个本征模态函数为h1，则令r1＝z-h1，且令r1替代z，重复上述步骤s2.1，直到获得第二个本征模态函数h2。
[0166]
以此类推，不妨设rk＝r
k-1-hk，且令rk替代z，重复上述步骤s2.1，直到获得第k+1个本征模态函数h
k+1
，其中k＝2,3,
…
。
[0167]
s2.3停止条件
[0168]
设置停止阈值t2＝0.2，计算判别式：
[0169][0170]
若sd≤t2，则停止继续进行经验模态分解，此时记k＝k-1＝2，即将相位信号z∈r
5000
共分解成了2个本征模态函数，结果见图7。
[0171]
s3：运动特征提取
[0172]
针对进行经验模态分解后的本征模态函数，分别进行运动参数估计，进而实现运动特征提取。在此主要利用基于傅里叶分析的频谱分析方法实现运动参数的估计。考虑上述两个本征模态函数，对其进行傅里叶变换，并取绝对值，进而取其最大值对应的频率值作为运动参数的估计，结果见下表，结果表明本方法能够提取目标自旋频率与进动频率两个运动特征，且运动参数的估计值与真实值非常接近，即精度较高。
[0173]
表1
[0174] 估计值真实值自旋频率9.75hz10hz进动频率0.49hz0.5hz
[0175]
本实施例还提供了一种雷达目标运动特征提取系统，包括：
[0176]
预处理单元，用于对雷达一维距离像数据进行预处理，得到目标的相位信号；
[0177]
信号分离单元，用于利用经验模态分解对所述相位信号进行信号分离，得到多个本征模态函数；
[0178]
参数估计单元，用于基于傅里叶分析的频谱分析方法，对所述本征模态函数分别进行运动参数估计，得到雷达目标运动特征。
[0179]
优选地，所述预处理单元包括：
[0180]
求和模块，用于对所述雷达一维距离像x∈cm×n计算绝对值，再对所述绝对值进行逐列求和，得到求和参数其中，x为所述雷达一维距离像，y(n)为所述求和参数，c为复数域，m为一维距离像时间序列长度，n为一维距离像距离单元个数，m为一维距离像时间序列索引，n为一维距离像距离单元索引，n＝1,2,
…
,n，y∈rn；
[0181]
提取模块，用于设置阈值设置阈值t1∈r，针对y∈rn，统计所述求和参数的指标集并根据所述指标集提取强散射单元的时序信号；其中γ为所述指标集，γ中的元素为雷达一维距离像x∈cm×n中强散射单元所在的列，k为一维距离像距离单元索引；
[0182]
相位表示模块，用于令y＝arg(x)表示所述时序信号的相位数据；其中y＝arg(x)为提取复数相位的算子，且y∈rm；
[0183]
第一计算模块，用于根据所述时序信号的相位数据计算相位梯度d(m)＝y(m)-y(m+1)；其中m＝1,2,
…
,m-1；
[0184]
第二计算模块，用于根据所述时序信号的相位数据计算缠绕相位差：e(m)＝arctan(sin(y(m))/cos(y(m)))，其中m＝1,2,
…
,m-1；
[0185]
初始化模块，用于令z(1)＝y(1)；z(1)为初始化数据；
[0186]
第三计算模块，用于根据公式z(m)＝y(m)+e(m-1)对缠绕相位求和，以得到解缠后的所述相位信号；其中z(m)为解缠后相位信号的第m个元素。
[0187]
优选地，所述信号分离单元包括：
[0188]
第一获取模块，用于获取第一个本征模态函数；
[0189]
第二获取模块，用于获取其它本征模态函数；
[0190]
停止模块，用于设置停止条件；
[0191]
所述第一获取模块包括：
[0192]
极值计算子模块，用于计算所述相位信号z∈rm的局部极大值点坐标集i与局部极小值点坐标集j；
[0193]
插值子模块，用于利用最邻近插值对局部极大值点{z(i)|i∈i}进行插值获取上包络信号a∈rm，并利用最邻近插值对局部极小值点{z(j)|j∈j}进行插值获取下包络信号b∈rm；
[0194]
均值计算子模块，用于计算所述上包络信号与所述下包络信号的均值包络信号m1＝(a+b)/2；
[0195]
第一获取子模块，用于令h1＝z-m1，判断h1是否为本征模态函数，若否，则令h1替代z，并跳转至步骤“计算所述相位信号z∈rm的局部极大值点坐标集i与局部极小值点坐标集j”，直到获得第一个本征模态函数；
[0196]
所述第二获取模块包括：
[0197]
第二获取子模块，用于设rk＝r
k-1-hk，且令rk替代z，重复步骤“获取第一个本征模态函数”，直到获得第k+1个本征模态函数h
k+1
，其中k＝2,3,
…
；
[0198]
所述停止模块包括：
[0199]
判别式计算子模块，用于设置停止阈值t2＞0，计算判别式
[0200]
判断子模块，用于若sd≤t2时，停止继续进行经验模态分解，此时记k＝k-1。
[0201]
优选地，所述参数估计单元包括：
[0202]
傅里叶处理模块，用于对所述本征模态函数依次进行傅里叶变换和取绝对值处理，得到频谱数据fk＝|fhk|；其中，f为离散傅里叶矩阵，所述离散傅里叶矩阵的元素为f
ij
＝exp(-2πj0(i-1)(j-1)/m)，其中，i,j∈{1,2,
…
,m}且j0为虚数单位，fk为所述频谱数据；hk为所述本征模态函数；
[0203]
估计模块，用于取所述频谱数据的最大值对应的频率值作为运动参数的估计，得到雷达目标运动特征wk＝argmaxfk,k＝1,2,
…
,k；其中wk为所述雷达目标运动特征。
[0204]
本发明的有益效果如下：
[0205]
(1)本发明能够提取目标自旋频率与进动频率两个运动特征，且运动参数的估计
值与真实值非常接近，即精度较高。
[0206]
(2)本发明能够适应雷达目标运动更加多样、高速运动、运动耦合等实际情况，提取的目标运动特征具有更高的精度。
[0207]
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
[0208]
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。