CN111693954B - 空间旋转群目标微多普勒分离方法 - Google Patents

Abstract

提出一种空间旋转群目标微多普勒分离方法，具体包括下列步骤：分析旋转目标微动模型和回波模型，构造不同长度下回波的相关矩阵，对该矩阵做奇异值分解，求解奇异值比与回波分段长度的数学关系；利用回波长度‑奇异值比值之间的关系，估计各目标对应的旋转周期；根据第二步估计得到的子目标回波周期，采用奇异值分解的方法分离各子目标对应的回波；利用时频重排的方法，获取各子目标时频图，估计各散射中心旋转中心，对时频图进行旋转处理，求得各散射中心时频图；利用维特比Viterbi算法对各子散射中心的时频图像进行提取，提取其时间‑微多普勒曲线。该方法可以简便、高效、准确的估计出群目标中各目标微动周期并提取各散射中心微多普勒信息，为后续的微动特征提取及微动目标识别提供技术支撑。

Description

空间旋转群目标微多普勒分离方法

技术领域

本发明属于信号处理技术领域，具体涉及一种对空间旋转群目标微多普勒分离算法。

背景技术

美国科学家V.C.Chen于2000年首次提出微动的概念，并对基于微动特性的目标特征提取与识别展开深入的研究。微动是指目标主体之外的其他部件进行的与目标主体不同的运动，典型的微动包括人行走过程中手臂的运动、直升机飞行过程中旋翼的运动、坦克行进时履带的运动等。微动也会产生相应的多普勒，被称为微多普勒。微多普勒作为微动目标的独有特征，能够有效反应目标特性，因而得到广泛的研究。

V.C.Chen将典型的微动方式归结为自旋、锥旋、摆动以及进动等。其中自旋、锥旋、摆动等目标的微多普勒均表现为正弦形式，因此将这些目标可以统一归类为旋转目标。根据散射中心理论，在光学区，目标回波可以等效为多个散射中心的叠加。因此，对目标的特性分析往往要通过各散射中心的微动特性进行分析。同时，由于空间中往往存在着多个不同旋转周期的微动目标，多目标、多散射中心的微多普勒信息往往在时频图上相互交叠，严重影响单个散射中心微多普勒的提取。因此，研究空间旋转群目标微多普勒分离算法具有重要的意义。

自微多普勒的概念被应用到目标特征提取方面以来，陆续产生了许多微多普勒分离方法。Thayaparan T在《Micro-Doppler radar signatures for intelligent targetreconition》(Defence Research and Development Canada，2004)中提出用峰值法来对散射中心的微多普勒进行提取，然而由于频率分辨率的限制和噪声的影响，该算法提取的微多普勒精度相对较差。LI P在《Separation of Micro-Doppler Signals Based on TimeFrequency Filter and Viterbi Algorithm》(Signal，Image and Video Processing，2011，7(3)：1-13)提出基于Viterbi和时频滤波的算法对振动散射中心微多普勒进行提取，然而该算法对于多散射中心微多普勒混叠时提取效果相对较差。赵盟盟在《一种用于空间群目标分辨的滑动窗轨迹跟踪算法》(宇航学报，2015，36(10)：1187-1194)对回波时频图采用形态学处理的方法抑制旁瓣，然后提出一种滑动窗轨迹跟踪的算法分离出相互缠绕的微多普勒，然而当对时频图进行形态学处理时，会造成部分散射中心时频图的断裂，从而影响后续的微多普勒提取。王义哲在《弹道中段微多普勒分离与提取仿真研究》(系统仿真学报，2017，29(6)：1201-1209)中提取采用基于完备总体经验模态分解和改进自适应Viterbi算相结合的群目标微多普勒提取算法，该算法受限完备总体经验模态分解容易产生模态混叠和端点效应，容易导致散射中心的微多普勒提取存在误差。

发明内容

针对现有技术存在的问题，本发明提供一种空间旋转群目标微多普勒分离方法，具体包括下列步骤：

第一步：分析旋转目标微动模型和回波模型，构造不同长度下回波的相关矩阵，对该矩阵做奇异值分解，求解奇异值比与回波分段长度的数学关系

采用窄带雷达发射单载频信号对目标进行观测；若雷达发射信号为s_i(t)

式中，t表示时间，T_c为脉冲宽度，f_c为雷达发射电磁波的载频；

电磁波遇到目标后，发生反射现象；若目标和雷达之间的距离变化为R(t)，则雷达接收到的目标回波与发射回波之间产生时延为2R(t)/c，相应的雷达回波表示为：

式中，c表示光速；

获得目标原始回波之后，需要对回波进行零中频处理消除载波的影响；同时，为了利用回波的相位信息，采用I/Q双通道正交解调处理，得到回波信号的复数表示：

式中，σ为散射中心的散射系数，j表示虚数单位；

对于空间旋转目标而言，由于雷达对其观测的条件一般满足远场条件，其回波可以等效为目标上几个有效的散射中心回波之和；若雷达观测范围内存在多个目标，每个目标上存在多个散射中心，则雷达回波进一步表示为

式中，L表示回波中的包含目标总数，l表示目标序号，M_l表示第l个目标对应的散射中心总数，i表示散射中心的序号，σ_li表示第l个目标上第i个散射中心的散射系数，R_li(t)表示第l个目标上第i个散射中心的微距离变化；

这里只对空间旋转散射中心的微动模型进行分析；

通过对回波进行采样处理，得N个采样时刻的回波的离散化表示s_r(n)为

式中，n表示采样点的序号，t_s表示采样间隔；

定义一个线性递增的整数变量N_t，N_t∈[1 N]，N_t∈Z，Z表示整数集合，将长度为N的回波s_r(n)划分N₀段长度为N_t的回波，并将这些回波排列为一个回波矩阵C；

式中，S_r，1，Nt、s_r，2，Nt、

分别表示经过分段处理后长度为N_t的第1段回波、第2段回波以及第N₀段回波，N₀的计算方法为

式中，函数ceil(·)表示靠左取整；从上式中看出，回波矩阵C的行数为N₀，列数为N_t；

对回波矩阵C进行奇异值分解SVD，

式中，U为左奇异矩阵，S为右奇异矩阵，∑为该矩阵的奇异值矩阵；∑中对角线元素为C的奇异值，即为C^HC对应的特征值的

次方，()^H表示对矩阵进行共轭转置矩阵处理；矩阵∑中非对角线上的元素均为0，矩阵S中每一列为矩阵C^HC的特征值对应的特征向量；分析可知，当N_t＝kN_s，k＝1，2，…，k表示正整数，N_s表示一个周期的回波所对应的长度；也就是当分段长度N_t为子目标周期的整数倍时，回波矩阵C中各列之间只是相差一个系数因子，此时奇异值矩阵中除了第一个奇异值较大外，其他N₀-1个奇异值会比较小；而当回波长度不满足子目标周期的整数倍时，回波矩阵C的奇异值普遍较小；因此，通过奇异值分布的差异，能够实现对目标旋转周期的有效估计；

如果将奇异值矩阵∑表示为

∑＝diag(σ₁，σ₂，...σ_N0) (9)

式中，diag(·)表示对角矩阵，即该矩阵只有对角线上的元素非零，其他位置的元素均为零；σ₁，σ₂，...σ_N0表示该矩阵中的N₀个元素，也就是回波矩阵C的N₀个奇异值；

定义奇异值比值κ＝σ₁/σ₃，求得不同回波长度N_t下对应κ值，能够求得回波长度-奇异值比值之间的关系；

第二步：利用回波长度-奇异值比值之间的关系，估计各目标对应的旋转周期；

利用先验知识对N_t进行约束；考虑到微动目标的微动频率ω一般满足约束条件ω＝[2π 10π]rad/s，因此将N_t的搜索限定在该范围内，也就是

在对奇异值比值进行分析的过程中，采用峰值搜索的方法确定各子目标的微动周期；同时需要注意的是，如果某个位置是目标周期的倍频分量，该处也将以极大值的状态出现在奇异值比值序列当中；因此在对微动频率进行估计的过程中，需要考虑倍频分量的影响；

第三步：根据第二步估计得到的子目标回波周期，采用奇异值分解的方法分离各子目标对应的回波；

根据第二步估计得到的L个子目标的旋转周期分别为T₁、T₂、…、T_L，根据式(32)计算这L个旋转周期所对应的回波长度分别为N₁、N₂、…、N_L；

通过式(32)计算出各个周期对应的回波长度之后，分别对N_t取值N₁、N₂、…、N_L，根据式(27)构造对应的相关矩阵C₁、C₂、…C_L，采用式(29)分别对C₁、C₂、…C_L进行奇异值分解，得到左、右奇异矩阵和奇异值矩阵；选择右奇异矩阵中第一个列向量作为子目标的回波，即可实现子目标回波的有效分离；

将式(26)中的回波模型进行进一步简化处理；定义g_li为第l个目标上第i个散射中心一个周期内的回波对应的基波分量，

式中，N_l表示第l个目标一个周期对应的回波长度；则一个周期内，该散射中心的回波表示为

s_li(n)＝σ_lig_li n＝1，2，3…，N_l (12)

以第1个目标为待分离对象，对子目标回波分离原理做具体分析如下；

令N_t＝N₁，根据式(27)构造回波矩阵C₁，

对于C₁进行奇异值分解，也就是求C^H ₁C₁对应的特征值及特征向量；

式中，M₁表示第1个目标对应的散射中心个数，I表示单位矩阵；

从式(36)中看出，C^H ₁C₁最大特征值对应的特征向量为a₁表示为

式中，q_1i表示第1个目标上第i个散射中心基波分量的幅度系数；比较式(37)和式(34)看出，特征向量a₁与第1个目标回波为仅在幅度上存在一定差异，其中基波分量相同，所以其相位中包含的微动信息相同；因此，通过上述方法能够从回波中分离出第一个目标的回波；

以目标1为例，通过上述分析，可以有效分离目标1对应的回波；根据估计得到的回波长度N₁、N₂、…、N_L，对L个目标均执行式(35)-(37)的操作，直至分离出所有子目标的回波；

第四步：利用时频重排的方法，获取各子目标时频图，估计各散射中心旋转中心，对时频图进行旋转处理，求得各散射中心时频图

采用时频重排的方法对各子目标回波进行时频分析，获取子目标的时频图；对时频图上零频位置进行局部峰值搜索，将这些局部峰值位置作为各散射中心对应旋转中心；以一个周期为变换空间，将时频图绕旋转中心进行180°旋转变换，得到各散射中心对应的时频图；

第五步：利用维特比Viterbi算法对各子散射中心的时频图像进行提取，提取其时间-微多普勒曲线；

Viterbi算法是基于目标时频信息的两个特点进行微多普勒曲线提取的；

特点1：散射中心真实微多普勒与其对应时频图像中位置上的幅度尽可能大；

特点2：相邻时刻之间散射中心的微多普勒变化尽可能小，不会出现突变或者其他剧烈变化的情况；

根据上述两个特点，对Viterbi算法的原理说明如下；

假设存在两个采样时刻n₁和n₂，采样时刻n∈[n₁ n₂]，k(n)表示路径，K表示所有路径的集合，则Viterbi算法的目标函数对应的最小路径

满足式(39)对应的约束条件；

式中，g(k(n)，k(n+1))表示对应于|k(n+1)-k(n)|的非减函数，|·|表示求绝对值，h(TF(n，k(n))为一个非增函数，TF(n，k(n))表示某个采样时刻为n、路径为k(n)的点在时频图中的幅值；

首先将n时刻的时频图对应频谱按由大到小降序排列

式中，m对应时频图中频率单元的序号，f_i表示频率值，i＝1，2，…M₀，M₀表示n时频图中频率单元的个数；

定义h(·)满足式(41)中的表达式；

h(TF(n，f_i))＝i-1 (18)

该函数表明，对于n时刻的频谱，峰值越大，则该点越有可能是其真实微多普勒的位置；

g(k(n)，k(n+1))对应于第二个特点，即要求相邻时刻内的微多普勒的变化较小，满足式(42)中的表达式；

式中，Δ表示微多普勒分辨率，由信号的重复频率决定，u为代价倍数；从上式看出，随着两个时刻微多普勒差异的增大，其代价函数越来越大；

采用Viterbi算法对第四步中得到的各个散射中心时频图处理，根据提取到的最优路径，计算散射中心对应的微多普勒；至此，群目标中各散射中心的微多普勒得到有效分离。

在本发明的一个实施例中，第二步具体包括以下步骤：

Step1：考虑到奇异值比值序列的分布特性，采用峰值搜索方法对奇异值比值序列进行峰值搜索，确定第一个峰值所对应的N_t，得到第一个目标对应的回波长度N₁；考虑到回波长度和周期之间满足：T₁＝N₁t_s，利用第一个峰值对应的N₁估计出周期T₁；

Step2：定义wid为窗口长度，以上一步中N₁及其倍频分量对应位置为中心，将这些位置所在的窗口内所有奇异值比值置零；然后在剩余的奇异值比值中继续进行搜索，得到相应的峰值，并重复上式对周期进行估计；

Step3：根据比较门限，若此时搜索得到的峰值与第一峰值的比例小于门限，则停止搜索；否则，重复上述步骤，直到满足条件为止；

定义门限

ε＝κ_{max 1}/κ_{max i} (20)

式中，κ_{max 1}为奇异值比值得最大值，κ_{max i}为第i次周期估计对应的峰值，当ε＜0.2，认为没有包含其他周期分量；同时，对每次搜索对应的倍频分量处的奇异值比进行分析，如果倍频处不是峰值，那么将该处对应的周期从原先的周期估计中删除。

在本发明的另一个实施例中，第四步具体包括以下步骤：

Step1：根据第三步分离得到的子目标回波，采用时频重排的方法获取各子目标时频图；

时频重排是在短时傅里叶变换的基础上对信号的一种时频分析方法；对于信号x(t)，经过时频重排后，时频图上任意一点(t′，f′)对应的幅值R_X(t′，f′；h)为所有重排至此位置的时频图值之和；时频重排的具体原理如式(38)所示；

式中，t′表示时频重排后的时间单元、f′为重排后的频率单元、t表示短时傅里叶变换后的时间单元、f为短时傅里叶变换后的频率单元、

表示时间重排算子，

表示频率重排算子，δ(·)表示为冲激函数，X(t，f；h)表示短时傅里叶变换的结果，h表示时频分析中采用的窗函数；

Step2：沿第l子目标时频图零频位置，搜索局部峰值点；根据搜索结果，得到M_l个局部峰值点，分别对应M_l旋转中心；

Step3：选择其中一个旋转中心，截取一个周期内的回波时频图，将截取的时频图绕该旋转中心旋转180°，即可得到该散射中心的时频信息；考虑到旋转目标的微多普勒表现为正弦形式，通过对时频图进行旋转处理，该散射中心的微多普勒曲线能够实现倍数级的强化；相应的，非该旋转中心对应散射中心的时频图像则在时频图中被弱化；与此同时，为了得到更为清晰的时频信息，对旋转变换后得到的时频图可以采用恒虚警降噪的方法进行进一步处理；

Step4：存储上一步得到的时频信息，然后从原时频图中减去该散射中心对应的时频信息；

Step5：重复Step2-Step4，直到完成所有散射中心的提取；

上述提取方法需要对目标的观测时间至少大于一个周期。

本发明为了解决空间旋转群目标各散射中心微多普勒交叉混叠、难以提取的问题，提出一种空间旋转群目标散射中心微多普勒分离方法。通过对构造旋转群目标微动回波矩阵，并对回波矩阵进行奇异值分解，估计各目标微动周期，利用周期差异分离各子目标回波，采用时频重排的方法获取各子目标时频图，对时频图进行旋转变换，提取各散射中心时频图，利用Viterbi算法提取各散射中心微多普勒，最终分离出各散射中心的微多普勒。

本发明的方法利用不同长度下回波矩阵奇异值的分布特性，实现对群目标中各个子目标微动周期的估计，然后利用微动周期对应长度下回波矩阵最大奇异值对应的特征向量，分离出每个目标的微动信息；根据旋转散射中心微多普勒的正弦特性，搜索时频图旋转中心，实现对单个散射中心时频图时频信息的增强处理和对噪声的抑制；最后利用Viterbi算法对各散射中心的微多普勒进行提取。本发明可以实现旋转群目标多散射中心微多普勒的有效分离和提取，为后续的目标参数估计和识别打下基础。

附图说明

图1为本发明空间旋转群目标微多普勒分离的示意性流程图；

图2为空间群目标回波时频图；

图3为奇异值比值随信号分段长度变化示意图；

图4为对特征值对应特征向量进行时频分析后得到的各个子目标时频图；

图5为对时频图进行旋转处理后提取到的各散射中心时频图；

图6为最终提取到得到的散射中心微多普勒。

具体实施方式

下面结合实施例、附图对本发明作进一步描述。

本发明的方法是：第一步：采用回波长度递增的方法构建回波矩阵，对回波矩阵进行奇异值分解，求取奇异值比值；第二步：通过对奇异值比值序列的峰值位置进行搜索，分析目标周期及其倍频分量之间的关系，估计各子目标对应的微动周期；第三步：根据估计得到的微动周期，重新构建回波矩阵并进行奇异值分解，分离出各个子目标的回波；第四步：利用各散射中心时频图对应的旋转中心，对单个周期内散射中心的时频图进行旋转变换，提取出各散射中心对应的时频图；第五步，利用Viterbi算法对各散射中心时频图进行处理，提取各散射中心对应的微多普勒。

具体步骤如下：

第一步：分析旋转目标微动模型和回波模型，构造不同长度下回波的相关矩阵，对该矩阵做奇异值分解，求解奇异值比与回波分段长度的数学关系；

本发明中采用窄带雷达发射单载频信号对目标进行观测。若雷达发射信号为s_i(t)

式中，t表示时间，T_c为脉冲宽度，f_c为雷达发射电磁波的载频。

电磁波遇到目标后，发生反射现象。若目标和雷达之间的距离变化为R(t)，则雷达接收到的目标回波与发射回波之间产生时延为2R(t)/c，相应的雷达回波可以表示为：

式中，c表示光速。

获得目标原始回波之后，需要对回波进行零中频处理消除载波的影响。同时，为了利用回波的相位信息，通常采用I/Q双通道正交解调处理，得到回波信号的复数表示：

式中，σ为散射中心的散射系数，j表示虚数单位。

对于空间旋转目标而言，由于雷达对其观测的条件一般满足远场条件，其回波可以等效为目标上几个有效的散射中心回波之和。若雷达观测范围内存在多个目标，每个目标上存在多个散射中心，则雷达回波可以进一步表示为

式中，L表示回波中的包含目标总数，l表示目标序号，M_l表示第l个目标对应的散射中心总数，i表示散射中心的序号，σ_li表示第l个目标上第i个散射中心的散射系数，R_li(t)表示第l个目标上第i个散射中心的微距离变化。

本发明中主要针对空间旋转目标，因此这里只对空间旋转散射中心的微动模型进行分析。需要说明的是空间旋转目标通常在进行微动的同时还进行着轨道运动，轨道运动对回波时频信息的影响可以通过多种方式补偿，这里不考虑其影响。

通过对回波进行采样处理，得N个采样时刻的回波的离散化表示s_r(n)为

式中，n表示采样点的序号，t_s表示采样间隔。

定义一个线性递增的整数变量N_t，N_t∈[1 N]，N_t∈Z，Z表示整数集合，将长度为N的回波s_r(n)划分N₀段长度为N_t的回波，并将这些回波排列为一个回波矩阵C。

式中，s_r，1，Nt、s_r，2，Nt、

分别表示经过分段处理后长度为N_t的第1段回波、第2段回波以及第N₀段回波，N₀的计算方法为

式中，函数ceil(·)表示靠左取整。从上式中可以看出，回波矩阵C的行数为N₀，列数为N_t。

对回波矩阵C进行奇异值分解(Singular value decomposition，SVD)，

式中，U为左奇异矩阵，S为右奇异矩阵，∑为该矩阵的奇异值矩阵。∑中对角线元素为C的奇异值，即为C^HC对应的特征值的

次方，()^H表示对矩阵进行共轭转置矩阵处理。矩阵∑中非对角线上的元素均为0，矩阵S中每一列为矩阵C^HC的特征值对应的特征向量。分析可知，当N_t＝kN_s，k＝1，2，…，k表示正整数，N_s表示一个周期的回波所对应的长度。即当分段长度N_t为子目标周期的整数倍时，回波矩阵C中各列之间只是相差一个系数因子，此时奇异值矩阵中除了第一个奇异值较大外，其他N₀-1个奇异值会比较小。而当回波长度不满足子目标周期的整数倍时，回波矩阵C的奇异值普遍较小。因此，通过奇异值分布的差异，可以实现对目标旋转周期的有效估计。

如果将奇异值矩阵∑表示为

∑＝diag(σ₁，σ₂，...σ_N0) (30)

式中，diag(·)表示对角矩阵，即该矩阵只有对角线上的元素非零，其他位置的元素均为零。σ₁，σ₂，...σ_N0表示该矩阵中的N₀个元素，即回波矩阵C的N₀个奇异值。

定义奇异值比值κ＝σ₁/σ₃，求得不同回波长度N_t下对应κ值，即可求得回波长度-奇异值比值之间的关系。

第二步：利用回波长度-奇异值比值之间的关系，估计各目标对应的旋转周期。

在第一步中，假设回波对应的长度N，如果使得N_t从0到N依次递增，同时又进行奇异值分解，运算量相对较大。因此，可以利用先验知识对N_t进行约束。考虑到微动目标的微动频率ω一般满足约束条件ω＝[2π 10π]rad/s，因此将N_t的搜索限定在该范围内，即

在对奇异值比值进行分析的过程中，主要采用峰值搜索的方法确定各子目标的微动周期。同时需要注意的是，如果某个位置是目标周期的倍频分量，该处也将以极大值的状态出现在奇异值比值序列当中。因此在对微动频率进行估计的过程中，需要考虑倍频分量的影响。

综上所述，将第二步目标旋转周期估计的过程总结为以下三个步骤：

Step1：考虑到奇异值比值序列的分布特性，采用峰值搜索方法对奇异值比值序列进行峰值搜索(张希会.多项式相位信号的检测与参数估计[D].电子科技大学，2012.)，确定第一个峰值所对应的N_t，得到第一个目标对应的回波长度N₁。考虑到回波长度和周期之间满足：T₁＝N₁t_s，利用第一个峰值对应的N₁估计出周期T₁；

Step2：定义wid为窗口长度，以上一步中N₁及其倍频分量对应位置为中心，将这些位置所在的窗口内所有奇异值比值置零。然后在剩余的奇异值比值中继续进行搜索，得到相应的峰值，并重复上式对周期进行估计；

Step3：根据比较门限，若此时搜索得到的峰值与第一峰值的比例小于门限，则停止搜索；否则，重复上述步骤，直到满足条件为止。

定义门限

ε＝κ_{max 1}/κ_{max i} (31)

式中，κ_{max 1}为奇异值比值得最大值，κ_{max i}为第i次周期估计对应的峰值，一般认为当ε＜0.2，则认为没有包含其他周期分量。同时，对每次搜索对应的倍频分量处的奇异值比进行分析，如果倍频处不是峰值，那么将该处对应的周期从原先的周期估计中删除。

第三步：根据第二步估计得到的子目标回波周期，采用奇异值分解的方法分离各子目标对应的回波。

根据第二步估计得到的L个子目标的旋转周期分别为T₁、T₂、…、T_L，根据式(32)计算这L个旋转周期所对应的回波长度分别为N₁、N₂、…、N_L。

通过式(32)计算出各个周期对应的回波长度之后，分别对N_t取值N₁、N₂、…、N_L，根据式(27)构造对应的相关矩阵C₁、C₂、…C_L，采用式(29)分别对C₁、C₂、…C_L进行奇异值分解，得到左、右奇异矩阵和奇异值矩阵。选择右奇异矩阵中第一个列向量作为子目标的回波，即可实现子目标回波的有效分离。

为了说明上文中右奇异矩阵中第一个列向量即为该周期下子目标的回波，将式(26)中的回波模型进行进一步简化处理。定义g_li为第l个目标上第i个散射中心一个周期内的回波对应的基波分量，

式中，N_l表示第l个目标一个周期对应的回波长度。则一个周期内，该散射中心的回波可以表示为

s_li(n)＝σ_lig_li n＝1，2，3…，N_l (34)

以第1个目标为待分离对象，对子目标回波分离原理做具体分析如下。

令N_t＝N₁，根据式(27)构造回波矩阵C₁，

对于C₁进行奇异值分解，即求C^H ₁C₁对应的特征值及特征向量。

式中，M₁表示第1个目标对应的散射中心个数，I表示单位矩阵。

从式(36)中可以看出，C^H ₁C₁最大特征值对应的特征向量为a₁可以表示为

式中，q_1i表示第1个目标上第i个散射中心基波分量的幅度系数。比较式(37)和式(34)可以看出，特征向量a₁与第1个目标回波为仅在幅度上存在一定差异，其中基波分量相同，所以其相位中包含的微动信息相同。因此，通过上述方法能够从回波中分离出第一个目标的回波。

以目标1为例，通过上述分析，可以有效分离目标1对应的回波。根据估计得到的回波长度N₁、N₂、…、N_L，对L个目标均执行式(35)-(37)的操作，直至分离出所有子目标的回波。

第四步：利用时频重排的方法，获取各子目标时频图，估计各散射中心旋转中心，对时频图进行旋转处理，求得各散射中心时频图

采用时频重排的方法对各子目标回波进行时频分析(李楠，曲长文，平殿发，苏峰.基于时频重排与WHT的多分量LFM信号识别[J].航天电子对抗，2009，25(06)：33-36.)，获取子目标的时频图。对时频图上零频位置进行局部峰值搜索(张希会.多项式相位信号的检测与参数估计[D].电子科技大学，2012.)，将这些局部峰值位置作为各散射中心对应旋转中心。以一个周期为变换空间，将时频图绕旋转中心进行180°旋转变换(程云鹏，张凯院，徐仲.矩阵论[M].西北工业大学出版社，2006.)，得到各散射中心对应的时频图。

综上所述，将本发明中第四步总结为以下几个步骤：

Step1：根据第三步分离得到的子目标回波，采用时频重排的方法获取各子目标时频图(李楠，曲长文，平殿发，苏峰.基于时频重排与WHT的多分量LFM信号识别[J].航天电子对抗，2009，25(06)：33-36.)；

时频重排是在短时傅里叶变换的基础上对信号的一种时频分析方法。对于信号x(t)，经过时频重排后，时频图上任意一点(t′，f′)对应的幅值R_X(t′，f′；h)为所有重排至此位置的时频图值之和。时频重排的具体原理如式(38)所示；

式中，t′表示时频重排后的时间单元、f′为重排后的频率单元、t表示短时傅里叶变换后的时间单元、f为短时傅里叶变换后的频率单元、

表示时间重排算子，

表示频率重排算子，δ(·)表示为冲激函数，X(t，f；h)表示短时傅里叶变换的结果，h表示时频分析中采用的窗函数。

Step2：沿第l子目标时频图零频位置，搜索局部峰值点(张希会.多项式相位信号的检测与参数估计[D].电子科技大学，2012.)。根据搜索结果，得到M_l个局部峰值点，分别对应M_l旋转中心；

Step3：选择其中一个旋转中心，截取一个周期内的回波时频图，将截取的时频图绕该旋转中心旋转180°，即可得到该散射中心的时频信息。考虑到旋转目标的微多普勒表现为正弦形式，通过对时频图进行旋转处理(程云鹏，张凯院，徐仲.矩阵论[M].西北工业大学出版社，2006.)，该散射中心的微多普勒曲线可以实现倍数级的强化。相应的，非该旋转中心对应散射中心的时频图像则在时频图中被弱化。与此同时，为了得到更为清晰的时频信息，对旋转变换后得到的时频图可以采用恒虚警降噪的方法进行进一步处理(吴顺君，梅晓春.雷达信号处理和数据处理[M].电子工业出版社，2008.)；

Step4：存储上一步得到的时频信息，然后从原时频图中减去该散射中心对应的时频信息；

Step5：重复Step2-Step4，直到完成所有散射中心的提取。

通过上述步骤，各子散射中心的时频图可以得到较好的提取。需要说明的是，上述提取方法需要对目标的观测时间至少大于一个周期。考虑到实际观测情况，一般满足这一观测条件，因此在这里不再进行说明。

第五步：利用维特比(Viterbi)算法对各子散射中心的时频图像进行提取，提取其时间-微多普勒曲线。

Viterbi算法是一种经典的散射中心时频提取方法，它原理简单，计算量适中，可以有效提取散射中心时频信息。尽管Viterbi也存在一些缺点，比如当散射中心时频图交叠严重时，会产生关联错误。考虑到本发明中该算法处理的时频图均只包含单个散射中心的时频信息，因此该算法的劣势不会影响本发明中的提取效果。

Viterbi算法是基于目标时频信息的两个特点进行微多普勒曲线提取的。

特点1：散射中心真实微多普勒与其对应时频图像中位置上的幅度尽可能大；

特点2：相邻时刻之间散射中心的微多普勒变化尽可能小，不会出现突变或者其他剧烈变化的情况。

根据上述两个特点，对Viterbi算法的原理说明如下。

假设存在两个采样时刻n₁和n₂，采样时刻n∈[n₁ n₂]，k(n)表示路径，K表示所有路径的集合，则Viterbi算法的目标函数对应的最小路径

满足式(39)对应的约束条件。

式中，g(k(n)，k(n+1))表示对应于|k(n+1)-k(n)|的非减函数，|·|表示求绝对值，h(TF(n，k(n))为一个非增函数，TF(n，k(n))表示某个采样时刻为n、路径为k(n)的点在时频图中的幅值。

首先将n时刻的时频图对应频谱按由大到小降序排列

式中，m对应时频图中频率单元的序号，f_i表示频率值，i＝1，2，…M₀，M₀表示n时频图中频率单元的个数。

在本发明中，定义h(·)满足式(41)中的表达式。

h(TF(n，f_i))＝i-1 (41)

该函数表明，对于n时刻的频谱，峰值越大，则该点越有可能是其真实微多普勒的位置。

g(k(n)，k(n+1))对应于第二个特点，即要求相邻时刻内的微多普勒的变化较小，满足式(42)中的表达式。

式中，Δ表示微多普勒分辨率，一般由信号的重复频率决定，u为代价倍数。从上式可以看出，随着两个时刻微多普勒差异的增大，其代价函数越来越大。

采用Viterbi算法对第四步中得到的各个散射中心时频图处理，根据提取到的最优路径，计算散射中心对应的微多普勒。至此，群目标中各散射中心的微多普勒得到有效分离。

实例：含旋转部件目标微动信号分离

仿真参数设定：假设空间存在两个旋转目标，目标1中包含两个散射中心，旋转半径分别为r₁₁＝0.6m，r₁₂＝0.2m；在z轴上的坐标分别为z₁₁＝0.8m，z₁₂＝0.6m；初始相位分别为

目标旋转轴与雷达视线夹角为α₁＝45°。目标2中包含两个散射中心，旋转半径分别为r₂₁＝0.3m，r₂₂＝0.7m；在z轴上的坐标分别为z₂₁＝0.6m，z₂₂＝-0.2m；初始相位分别为

旋转角频率分别为ω₁＝8πrad/s，ω₂＝5πrad/s；目标旋转轴与雷达视线夹角为α₁＝85°。四个散射中心的系数之比为σ₁₁∶σ₁₂∶σ₂₁∶σ₂₂＝1∶1.3∶1∶1.3。雷达载频f_c＝10GHz，脉冲重复频率PRF＝2KHz，观测时间为4s，信噪比为10dB。

采用上述设定的参数进行仿真实验，图2为两个目标回波的时频图；图3为目标回波的分段点数和奇异值比值之间的关系；图4分离后的两个子目标时频图；图5各个子散射中心对应的时频图；图6为采用Viterbi算法分离得到的各个散射中心微多普勒。比较图2和图6可以看出，利用本发明可以较好分离群目标中各个散射中心的微多普勒，有利于下一步目标微动参数的估计。

Claims (3)

1.一种空间旋转群目标微多普勒分离方法，其特征在于，具体包括下列步骤：

第一步：分析旋转目标微动模型和回波模型，构造不同长度下回波的相关矩阵，对该矩阵做奇异值分解，求解奇异值比与回波分段长度的数学关系

采用窄带雷达发射单载频信号对目标进行观测；若雷达发射信号为s_i(t)

式中，t表示时间，T_c为脉冲宽度，f_c为雷达发射电磁波的载频；

电磁波遇到目标后，发生反射现象；若目标和雷达之间的距离变化为R(t)，则雷达接收到的目标回波与发射回波之间产生时延为2R(t)/c，相应的雷达回波表示为：

式中，c表示光速；

获得目标原始回波之后，需要对回波进行零中频处理消除载波的影响；同时，为了利用回波的相位信息，采用I/Q双通道正交解调处理，得到回波信号的复数表示：

式中，σ为散射中心的散射系数，j表示虚数单位；

对于空间旋转目标而言，由于雷达对其观测的条件一般满足远场条件，其回波可以等效为目标上几个有效的散射中心回波之和；若雷达观测范围内存在多个目标，每个目标上存在多个散射中心，则雷达回波进一步表示为

式中，L表示回波中的包含目标总数，l表示目标序号，M_l表示第l个目标对应的散射中心总数，i表示散射中心的序号，σ_li表示第l个目标上第i个散射中心的散射系数，R_li(t)表示第l个目标上第i个散射中心的微距离变化；

这里只对空间旋转散射中心的微动模型进行分析；

通过对回波进行采样处理，得N个采样时刻的回波的离散化表示s_r(n)为

式中，n表示采样点的序号，t_s表示采样间隔；

定义一个线性递增的整数变量N_t，N_t∈[1 N]，N_t∈Z，Z表示整数集合，将长度为N的回波s_r(n)划分N₀段长度为N_t的回波，并将这些回波排列为一个回波矩阵C；

式中，s_r，1，Nt、s_r，2，Nt、

分别表示经过分段处理后长度为N_t的第1段回波、第2段回波以及第N₀段回波，N₀的计算方法为

式中，函数ceil(·)表示靠左取整；从上式中看出，回波矩阵C的行数为N₀，列数为N_t；

对回波矩阵C进行奇异值分解SVD，

式中，U为左奇异矩阵，S为右奇异矩阵，∑为该矩阵的奇异值矩阵；∑中对角线元素为C的奇异值，即为C^HC对应的特征值的

次方，()^H表示对矩阵进行共轭转置矩阵处理；矩阵∑中非对角线上的元素均为0，矩阵S中每一列为矩阵C^HC的特征值对应的特征向量；分析可知，当N_t＝kN_s，k＝1，2，…，k表示正整数，N_s表示一个周期的回波所对应的长度；也就是当分段长度N_t为子目标周期的整数倍时，回波矩阵C中各列之间只是相差一个系数因子，此时奇异值矩阵中除了第一个奇异值较大外，其他N₀-1个奇异值会比较小；而当回波长度不满足子目标周期的整数倍时，回波矩阵C的奇异值普遍较小；因此，通过奇异值分布的差异，能够实现对目标旋转周期的有效估计；

如果将奇异值矩阵∑表示为

∑＝diag(σ₁，σ₂，...σ_N0) (51)

式中，diag(·)表示对角矩阵，即该矩阵只有对角线上的元素非零，其他位置的元素均为零；σ₁，σ₂，...σ_N0表示该矩阵中的N₀个元素，也就是回波矩阵C的N₀个奇异值；

定义奇异值比值κ＝σ₁/σ₃，求得不同回波长度N_t下对应κ值，能够求得回波长度-奇异值比值之间的关系；

第二步：利用回波长度-奇异值比值之间的关系，估计各目标对应的旋转周期；

利用先验知识对N_t进行约束；考虑到微动目标的微动频率ω一般满足约束条件ω＝[2π10π]rad/s，因此将N_t的搜索限定在该范围内，也就是

在对奇异值比值进行分析的过程中，采用峰值搜索的方法确定各子目标的微动周期；同时需要注意的是，如果某个位置是目标周期的倍频分量，该处也将以极大值的状态出现在奇异值比值序列当中；因此在对微动频率进行估计的过程中，需要考虑倍频分量的影响；

第三步：根据第二步估计得到的子目标回波周期，采用奇异值分解的方法分离各子目标对应的回波；

根据第二步估计得到的L个子目标的旋转周期分别为T₁、T₂、…、Y_L，根据式(32)计算这L个旋转周期所对应的回波长度分别为N₁、N₂、…、N_L；

通过式(32)计算出各个周期对应的回波长度之后，分别对N_t取值N₁、N₂、…、N_L，根据式(27)构造对应的相关矩阵C₁、C₂、…C_L，采用式(29)分别对C₁、C₂、…C_L进行奇异值分解，得到左、右奇异矩阵和奇异值矩阵；选择右奇异矩阵中第一个列向量作为子目标的回波，即可实现子目标回波的有效分离；

将式(26)中的回波模型进行进一步简化处理；定义g_li为第l个目标上第i个散射中心一个周期内的回波对应的基波分量，

式中，N_l表示第l个目标一个周期对应的回波长度；则一个周期内，该散射中心的回波表示为

s_li(n)＝σ_lig_li n＝1，2，3…，N_l (54)

以第1个目标为待分离对象，对子目标回波分离原理做具体分析如下；

令N_t＝N₁，根据式(27)构造回波矩阵C₁，

对于C₁进行奇异值分解，也就是求C^H ₁C₁对应的特征值及特征向量；

式中，M₁表示第1个目标对应的散射中心个数，I表示单位矩阵；

从式(36)中看出，C^H ₁C₁最大特征值对应的特征向量为a₁表示为

式中，q_1i表示第1个目标上第i个散射中心基波分量的幅度系数；比较式(37)和式(34)看出，特征向量a₁与第1个目标回波为仅在幅度上存在一定差异，其中基波分量相同，所以其相位中包含的微动信息相同；因此，通过上述方法能够从回波中分离出第一个目标的回波；

以目标1为例，通过上述分析，可以有效分离目标1对应的回波；根据估计得到的回波长度N₁、N₂、…、N_L，对L个目标均执行式(35)-(37)的操作，直至分离出所有子目标的回波；

第四步：利用时频重排的方法，获取各子目标时频图，估计各散射中心旋转中心，对时频图进行旋转处理，求得各散射中心时频图

采用时频重排的方法对各子目标回波进行时频分析，获取子目标的时频图；对时频图上零频位置进行局部峰值搜索，将这些局部峰值位置作为各散射中心对应旋转中心；以一个周期为变换空间，将时频图绕旋转中心进行180°旋转变换，得到各散射中心对应的时频图；

第五步：利用维特比Viterbi算法对各子散射中心的时频图像进行提取，提取其时间-微多普勒曲线；

Viterbi算法是基于目标时频信息的两个特点进行微多普勒曲线提取的；

特点1：散射中心真实微多普勒与其对应时频图像中位置上的幅度尽可能大；

特点2：相邻时刻之间散射中心的微多普勒变化尽可能小，不会出现突变或者其他剧烈变化的情况；

根据上述两个特点，对Viterbi算法的原理说明如下；

假设存在两个采样时刻n₁和n₂，采样时刻n∈[n₁ n₂]，k(n)表示路径，K表示所有路径的集合，则Viterbi算法的目标函数对应的最小路径

满足式(39)对应的约束条件；

式中，g(k(n)，k(n+1))表示对应于|k(n+1)-k(n)|的非减函数，|·|表示求绝对值，h(TF(n，k(n))为一个非增函数，TF(n，k(n))表示某个采样时刻为n、路径为k(n)的点在时频图中的幅值；

首先将n时刻的时频图对应频谱按由大到小降序排列

式中，m对应时频图中频率单元的序号，f_i表示频率值，i＝1，2，…M₀，M₀表示n时频图中频率单元的个数；

定义h(·)满足式(41)中的表达式；

h(TF(n，f_i))＝i-1 (60)

该函数表明，对于n时刻的频谱，峰值越大，则该点越有可能是其真实微多普勒的位置；

g(k(n)，k(n+1))对应于第二个特点，即要求相邻时刻内的微多普勒的变化较小，满足式(42)中的表达式；

式中，Δ表示微多普勒分辨率，由信号的重复频率决定，u为代价倍数；从上式看出，随着两个时刻微多普勒差异的增大，其代价函数越来越大；

采用Viterbi算法对第四步中得到的各个散射中心时频图处理，根据提取到的最优路径，计算散射中心对应的微多普勒；至此，群目标中各散射中心的微多普勒得到有效分离。

2.如权利要求1所述的空间旋转群目标微多普勒分离方法，其特征在于，第二步具体包括以下步骤：

Step1：考虑到奇异值比值序列的分布特性，采用峰值搜索方法对奇异值比值序列进行峰值搜索，确定第一个峰值所对应的N_t，得到第一个目标对应的回波长度N₁；考虑到回波长度和周期之间满足：T₁＝N₁t_s，利用第一个峰值对应的N₁估计出周期T₁；

Step2：定义wid为窗口长度，以上一步中N₁及其倍频分量对应位置为中心，将这些位置所在的窗口内所有奇异值比值置零；然后在剩余的奇异值比值中继续进行搜索，得到相应的峰值，并重复上式对周期进行估计；

Step3：根据比较门限，若此时搜索得到的峰值与第一峰值的比例小于门限，则停止搜索；否则，重复上述步骤，直到满足条件为止；

定义门限

ε＝κ_max1/κ_maxi (62)

式中，κ_max1为奇异值比值得最大值，κ_maxi为第i次周期估计对应的峰值，当ε＜0.2，认为没有包含其他周期分量；同时，对每次搜索对应的倍频分量处的奇异值比进行分析，如果倍频处不是峰值，那么将该处对应的周期从原先的周期估计中删除。

3.如权利要求1所述的空间旋转群目标微多普勒分离方法，其特征在于，第四步具体包括以下步骤：

Step1：根据第三步分离得到的子目标回波，采用时频重排的方法获取各子目标时频图；

时频重排是在短时傅里叶变换的基础上对信号的一种时频分析方法；对于信号x(t)，经过重排后，时频图上任意一点(t′，f′)对应的幅值R_X(t′，f′；h)为所有重排至此位置的时频图值之和；时频重排的具体原理如式(38)所示；

式中，t′表示时频重排后的时间单元、f′为重排后的频率单元、t表示短时傅里叶变换后的时间单元、f为短时傅里叶变换后的频率单元、

表示时间重排算子，

表示频率重排算子，δ(·)表示为冲激函数，X(t，f；h)表示短时傅里叶变换的结果，h表示时频分析中采用的窗函数；

Step2：沿第l子目标时频图零频位置，搜索局部峰值点；根据搜索结果，得到M_l个局部峰值点，分别对应M_l旋转中心；

Step3：选择其中一个旋转中心，截取一个周期内的回波时频图，将截取的时频图绕该旋转中心旋转180°，即可得到该散射中心的时频信息；考虑到旋转目标的微多普勒表现为正弦形式，通过对时频图进行旋转处理，该散射中心的微多普勒曲线能够实现倍数级的强化；相应的，非该旋转中心对应散射中心的时频图像则在时频图中被弱化；与此同时，为了得到更为清晰的时频信息，对旋转变换后得到的时频图可以采用恒虚警降噪的方法进行进一步处理；

Step4：存储上一步得到的时频信息，然后从原时频图中减去该散射中心对应的时频信息；

Step5：重复Step2-Step4，直到完成所有散射中心的提取；

上述提取方法需要对目标的观测时间至少大于一个周期。

Images