CN109597042B - 基于奇异谱分析的目标进动频率估计方法 - Google Patents

Abstract

本发明属于雷达信号处理技术领域，公开了基于奇异谱分析的目标进动频率估计方法。该方法利用奇异谱分析对目标的原始RCS序列进行分析：首先获取目标的原始RCS序列，构造轨迹矩阵；再利用PCA对轨迹矩阵进行分解，得到主分量和主分量对应的空间转换矩阵；然后利用主分量和空间转换矩阵，选取重构参数，对轨迹矩阵进行重构，提取重构的RCS序列；最后利用重构RCS序列，对目标的进动频率进行估计。本发明能在估计目标进动频率前，去除RCS序列中的噪声和非平稳分量，进而对目标进动频率进行估计，能够去除目标进动频率估计过程中的主要误差源，提升目标进动频率估计时的抗噪性能，提升算法稳定性。

Description

基于奇异谱分析的目标进动频率估计方法

技术领域

本发明涉及雷达信号处理技术领域，尤其涉及基于奇异谱分析的目标进动频率估计方法。

背景技术

雷达散射截面积(英文全称：Radar Cross Section,英文缩写：RCS)反映了目标对雷达发射电磁波的散射能力，与入射电磁波频率、角度、目标尺寸、结构和表面涂层有密切关系，因而目标RCS特性广泛应用于目标分类与识别之中。目标RCS特性主要包括频率特性和角度特性：目标RCS会对不同频率的入射电磁波呈现不同的散射特征，根据目标尺寸和电磁波长的相对大小，在实际的观测场景中，运动目标相对雷达的视角不断变化，目标RCS值也随之变化，构成目标RCS序列，由此可以提取出目标的一些运动特征，如进动频率，进动角等。

对于包含锥旋运动的空间目标而言，进动是目标运动的主要组成部分，而目标进动频率主要与目标的质量分布有关，是目标分类和识别中非常重要的特征量。然而实际情况中，由于目标运动情况复杂，且锥旋运动的空间目标尺寸一般较小，目标RCS本身很小，使得雷达观测的目标RCS序列信噪比较低，进而导致基于低信噪比的目标RCS序列估计目标的进动频率困难较大。

针对目标进动频率估计问题，国内外学者开展了一系列的研究。自相关函数(英文全称:Auto Correlation Function，英文缩写:AUTOC)法可以用于噪声环境下的语音信息特征提取。AUTOC法直接使用目标RCS观测序列的自相关函数进行周期分析，可以较为直观地反映目标RCS观测序列的周期性，但是由于噪声和目标复杂运动带来的非平稳分量的影响，使用AUTOC法难以准确估计目标的频率特征，例如提取目标进动频率时容易出现整数倍的误差。平均幅度差函数(英文全称:Average Magnitude Difference Function,英文缩写:AMDF)法是提取语音特征中重复性和周期性的高效方法。AMDF法从另一个角度提取目标RCS观测序列的周期特征，但同样由于噪声和非平稳分量的影响，该方法存在目标进动频率估计结果容易出现整数倍误差的问题，此外，该方法计算得到的AMDF函数往往存在大量的虚假值，对目标RCS序列周期特征提取造成很大的影响。在AUTOC法和AMDF法的基础上提出的循环自相关函数-循环平均幅度差函数法(英文全称:Circular Auto Correlation-Circular Average Magnitude Difference Function,英文缩写:CAUTOC-CAMDF)，该方法在前两种方法的改进基础上进行结合,可以一定程度上避免虚假值的影响，增强抗噪性能，但不能完全去除噪声和非平稳分量带来的估计误差。此外，一种通用的周期特征提取方法是基于变换域数据的傅里叶变换方法，该方法利用傅里叶变换(英文全称:FourierTransform，英文简称：FT)的方法对进动频率进行估计，需要较长的观测时间以保证估计精度，此时，长观测时间和低采样率条件下目标运动的非平稳分量对傅里叶变换估计的准确性影响均会变大，降低目标进动频率的估计精度。

综上所述，目标RCS序列中的噪声分量和非平稳分量，在现有的进动频率估计方法使用中会引起严重的估计误差，而非平稳分量在长观测时间或低采样率条件下也会引起估计误差，二者均是造成现有方法无法从目标RCS序列中准确估计进动频率参数的主要原因。

发明内容

为了克服现有技术的不足，本发明提供一种基于奇异谱分析的目标进动频率估计方法，可以在估计目标进动频率前，去除RCS序列中的噪声和非平稳分量，进而能够避免在目标进动频率估计过程中引入误差，提升目标进动频率估计时的抗噪性能，提升算法稳定性，为雷达目标特征提取和识别的准确性提供良好的保障。

为达到上述目的，本发明的实施例采用如下技术方案：

步骤1，获取目标的原始雷达散射截面积RCS序列，利用目标的原始RCS序列，构造轨迹矩阵R。

步骤2，利用主成分分析法PCA对轨迹矩阵R进行分解，得到轨迹矩阵R的主分量Y及主分量Y的空间转换矩阵A。

步骤3，利用主分量Y和空间转换矩阵A，选取重构参数q，对轨迹矩阵R进行重构，得到轨迹矩阵R的重构矩阵

并从重构矩阵

中提取重构RCS序列。

步骤4，利用重构RCS序列，对目标进动频率进行估计。

基于本发明提供的基于奇异谱分析的目标进动频率估计方法，利用奇异谱分析对目标的原始RCS序列进行分析：首先获取目标的原始RCS序列，构造轨迹矩阵；再利用PCA对轨迹矩阵进行分解，得到主分量和主分量对应的空间转换矩阵；然后利用主分量和空间转换矩阵，选取合适的重构参数，对轨迹矩阵进行重构，提取重构的RCS序列；最后利用重构RCS序列，对目标的进动频率进行估计。

本发明方法使用PCA对轨迹矩阵进行分解，得到对角化矩阵和主分量，对角化矩阵是主分量的协方差矩阵，对角化矩阵的对角线元素对应着主分量的行向量的方差，该方差等价于原始RCS序列经奇异谱分析后各个分量能量的大小。根据特征值的大小关系，可知分量能量大小是依次递减的，其中，代表序列周期变化的分量能量最大，而其他分量，如噪声分量和非平稳运动分量能量较小。也就是说，本发明方法利用所得到的对角化矩阵和主分量对轨迹矩阵进行重构，进而从重构矩阵中提取RCS序列，在重构的过程中能够去除原始RCS序列中存在的噪声分量和非平稳分量，这样一来，在利用新提取的RCS序列估计目标进动频率时，即可去除目标进动频率估计过程中的主要误差源，提升目标进动频率估计时的抗噪性能，并提升算法稳定性，为雷达目标特征提取和识别的准确性提供良好的保障。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例提供的基于奇异谱分析的目标进动频率估计方法流程示意图；

图2为奇异谱分析去除噪声和非平稳分量效果图；

图3为本发明实施例提供的方法与传统方法的采样信噪比-估计误差对比图；其中，图3(a)为本发明方法和FT的采样信噪比-估计误差对比图，图3(b)为本发明方法和CAUTOC-CAMDF的采样信噪比-估计误差对比图；

图4为本发明实施例提供的方法与传统方法的RCS序列采样时间-估计误差对比图；其中图4(a)为本发明方法和FT的RCS序列采样时间-估计误差对比图，图4(b)为本发明方法和CAUTOC-CAMDF的RCS序列采样频率-估计误差对比图；

图5为本发明实施例提供的方法与传统方法的RCS序列采样频率-估计误差对比图；其中，其中图5(a)为本发明方法和FT的RCS序列采样频率-估计误差对比图，图5(b)为本发明方法和CAUTOC-CAMDF的RCS序列采样频率-估计误差对比图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

图1为本发明实施例提供的基于奇异谱分析的目标进动频率估计方法流程示意图。

参加图1，本发明实施例提供的基于奇异谱分析的目标进动频率估计方法包括以下步骤：

步骤1，获取目标的原始雷达散射截面积RCS序列，利用目标的原始RCS序列，构造轨迹矩阵R。

具体的，步骤1包括以下步骤：

(1.1)通过雷达采集数据，获取目标的原始RCS序列。

(1.2)确定构造轨迹矩阵的参数：轨迹矩阵R的目标的原始RCS序列长度N，信号滑动窗口的长度L，信号样本数M。

其中，L需保证信号滑动窗口至少包含两个采样周期的序列信息，信号样本数M＝N-L+1，若N为偶数，则L＝N/2；若N为奇数，则L＝(N+1)/2。

(1.3)令目标的原始RCS序列中的第n个单元为RCS(n)，n∈{0,1,...,N-1}，构造轨迹矩阵R：

其中，

m∈{0,1,…,M-1}；上标T表示矩阵的转置。

步骤2，利用主成分分析法PCA对轨迹矩阵R进行分解，得到轨迹矩阵R的主分量Y及主分量Y的空间转换矩阵A。

具体的，步骤2包括以下步骤：

(2.1)分别计算

中所有元素的平均值

得到轨迹矩阵R的平均向量

其中，

(2.2)计算轨迹矩阵R的自相关矩阵C_R:

其中，E_1×(N-M+1)为1×(N-M+1)阶的矩阵，矩阵中的元素全为1，自相关矩阵C_R为M×M阶的矩阵。

(2.3)计算自相关矩阵C_R的对角化矩阵C_Y和空间转换矩矩阵A：

C_Y＝AC_RA^T＝diag[λ₀,λ₁,…,λ_M-1]；

其中，diag[λ₀,λ₁,…,λ_M-1]为对角线元素为λ₀,λ₁,…λ_M-1的矩阵，λ₀≥λ₁≥…≥λ_M-1，λ₀,λ₁,…λ_M-1为自相关矩阵C_R的所有特征值，空间转换矩阵A的行是λ₀,λ₁,…λ_M-1对应的特征向量；对角化矩阵C_Y为M×M阶的矩阵，自相关矩阵为M×M阶的矩阵。

(2.4)计算主分量Y，

主分量Y为M×(N-M+1)阶的矩阵。

需要说明的是，对角化矩阵C_Y是主分量Y的协方差矩阵，主分量

主分量的行向量

对应着对角化矩阵C_Y的方差λ₀,λ₁,…λ_M-1，由于已经进行去均值处理，该方差等价于序列中奇异谱分析后各个分量能量的大小。根据特征值的大小关系，可知分量能量大小是依次递减的，其中，代表序列周期变化的分量能量最大，而其他分量，如噪声分量和非平稳运动分量能量较小。

步骤3，利用主分量Y和空间转换矩阵A，选取重构参数q，对轨迹矩阵R进行重构，得到轨迹矩阵R的重构矩阵

并从重构矩阵

中提取重构RCS序列。

具体的，步骤3包括以下步骤：

(3.1)选取重构参数q＝1，将M×M阶矩阵A截断为M×1阶矩阵A_q，计算所述重构矩阵

其中，所述重构矩阵

为一维向量，

(3.2)从所述重构矩阵

中提取重构RCS序列：令所述重构RCS序列第l个单元为RCS(l)，将所述重构矩阵

中的元素作为所述重构RCS序列中对应位置处的单元，令RCS(l)＝r_l，l∈{1,2,...,N-M+1}。

需要说明的是，在对轨迹矩阵R进行重构的过程中，使用有损失的重构去除能量较低的噪声分量和非平稳分量，保留能量最大的周期分量。因此，选取重构参数q＝1对轨迹矩阵R进行重构，进而从重构矩阵中提取重构RCS序列，能够去除目标RCS序列中的噪声分量和非平稳分量。

示例性的，参见图2，图2(a)为原始的RCS序列，可以看出其中包含了噪声和非平稳分量，从中提取频率信息十分困难；图2(b)为利用本发明实施例方法对图2(a)中的序列去除噪声和非平稳分量得到的结果，由图2(b)可以看出，去除噪声分量和非平稳分量的效果显著。

步骤4，利用重构RCS序列，对目标进动频率进行估计。

具体来说，可利用传统周期估计方法，如CAUTOC-CAMDF和傅里叶变换法，对目标进动频率进行估计。

即，本发明实施例的一种具体的实施方式中，步骤4具体可以包括以下步骤：

(4.1a)构造重构RCS序列的循环自相关函数：

其中，mod(l+k,N-M+1)表示将l+k按照模N-M+1取余。

(4.2a)构造重构后RCS序列的循环平均幅度差函数：

(4.3a)构造目标函数：

(4.4a)确定所述目标函数取最大值时对应的k的取值，假设为k₀，，k₀∈{0,1,...,N-M+1}。则重构RCS序列的周期值为：

T＝k₀/f_s，

其中，f_s为目标的原始RCS序列的采样频率。

(4.5a)根据所述重构RCS序列的周期值，计算得到所述目标进动频率：

f＝1/T。

本发明实施例的另一种具体的实施方式中，步骤4具体可以包括以下步骤：

(4.1b)构造重构RCS序列的离散傅里叶变换式：

其中，j为虚数单位，j²＝-1。

(4.2b)确定所述目标函数取最大值时对应的k的取值，假设为k₀，

则所述重构目标进动频率估计值表示为：

其中，f_s为所述目标的原始RCS序列的采样频率。

本发明方法使用PCA对轨迹矩阵进行分解，得到对角化矩阵和主分量，对角化矩阵是主分量的协方差矩阵，对角化矩阵的对角线元素对应着主分量的行向量的方差，该方差等价于原始RCS序列经奇异谱分析后各个分量能量的大小。根据特征值的大小关系，可知分量能量大小是依次递减的，其中，代表序列周期变化的分量能量最大，而其他分量，如噪声分量和非平稳运动分量能量较小。也就是说，本发明方法利用所得到的对角化矩阵和主分量对轨迹矩阵进行重构，进而从重构矩阵中提取RCS序列，在重构的过程中能够去除原始RCS序列中存在的噪声分量和非平稳分量，这样一来，在利用新提取的RCS序列估计目标进动频率时，即可去除目标进动频率估计过程中的主要误差源，提升目标进动频率估计时的抗噪性能，并提升算法稳定性，为雷达目标特征提取和识别的准确性提供良好的保障。

进一步的，以下通过仿真实验对本发明实施例方法上述有益效果进行验证：

实验一：

1、仿真条件

设定目标的原始RCS序列采样频率f_s＝10Hz，目标进动频率f＝1.5Hz，目标的原始RCS序列采样时长t＝30s。

2、仿真实验内容及结果分析

基于目标的原始RCS序列，采用传统的FT和传统的CAUTOC-CAMDF及本发明方法进行目标进动频率估计，并对比两种方法的估计误差，绘制本发明方法和FT的信噪比-估计误差对比图以及本发明方法和SSA-FT的信噪比-估计误差对比图。

仿真结果如图3所示，其中图3(a)为为本发明方法和FT的信噪比估计误差对比图，图3(b)为本发明方法和CAUTOC-CAMDF的信噪比估计误差对比图；观察图3(a)和图3(b)，可以得到，当采样序列的信噪比较低时，对于同样大小的信噪比，本发明提供的方法得到的目标进动频率的估计误差较小，具有更好的抗噪性能。

实验二

1、仿真条件

设定目标的原始RCS序列采样频率f_s＝10Hz，目标进动频率f＝1.5Hz，目标的原始RCS序列信噪比SNR＝5dB.

2、仿真实验内容及结果分析

基于目标的原始RCS序列，采用传统的FT和传统的CAUTOC-CAMDF及本发明方法进行目标进动频率估计，并对比两种方法的估计误差，绘制本发明方法和FT的采样时间-估计误差对比图以及本发明方法和CAUTOC-CAMDF的采样时间-估计误差对比图。

仿真结果如图4所示，其中图4(a)为本发明方法和FT的RCS序列采样时间-估计误差对比图，图4(b)为本发明方法和CAUTOC-CAMDF的RCS序列采样频率-估计误差对比图；观察图4(a)和图4(b)，可以得到，仿真结果表明，对于相同的RCS序列采样时间，本发明提供的方法估计的目标进动频率的估计误差较小，当采样时间变化时，并且对不同的RCS序列采样时间估计更加稳定，大幅减小了传统方法的不稳定性。

实验三

1、仿真条件

设定目标的原始RCS序列采样时长t＝30s，目标进动频率f＝1.5Hz，目标的原始RCS序列信噪比SNR＝5dB。

2、仿真实验内容及结果分析

基于目标的原始RCS序列，采用传统的FT和传统的CAUTOC-CAMDF及本发明方法进行目标进动频率估计，并对比两种方法的估计误差，绘制本发明方法和FT的采样频率-估计误差对比图以及本发明方法和CAUTOC-CAMDF的采样频率-估计误差对比图。

仿真结果如图5所示，其中图5(a)为本发明方法和FT的RCS序列采样频率-估计误差对比图，图5(b)为本发明方法和CAUTOC-CAMDF的RCS序列采样频率-估计误差对比图；观察图5(a)和图5(b)，可以得到，对于相同的RCS序列采样频率，本发明提供的方法估计的目标进动频率的估计误差较小，即对于采样频率的变化不敏感，具有较好的稳定性。

本领域普通技术人员可以理解：实现上述方法实施例的全部或部分步骤可以通过程序指令相关的硬件来完成，前述的程序可以存储于一计算机可读取存储介质中，该程序在执行时，执行包括上述方法实施例的步骤；而前述的存储介质包括：ROM、RAM、磁碟或者光盘等各种可以存储程序代码的介质。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。

Claims (6)

1.基于奇异谱分析的目标进动频率估计方法，其特征在于，包括以下步骤：

步骤1，获取目标的原始雷达散射截面积RCS序列，利用所述目标的原始RCS序列，构造轨迹矩阵R；

步骤2，利用主成分分析法PCA对所述轨迹矩阵R进行分解，得到所述轨迹矩阵R的主分量Y及所述主分量Y的空间转换矩阵A；

步骤3，利用所述主分量Y和所述空间转换矩阵A，选取重构参数q，对所述轨迹矩阵R进行重构，得到所述轨迹矩阵R的重构矩阵

并从所述重构矩阵

中提取重构RCS序列；

步骤4，利用所述重构RCS序列，对目标进动频率进行估计。

2.根据权利要求1所述的方法，其特征在于，所述步骤1具体包括：

(1.1)通过雷达采集数据，获取所述目标的原始RCS序列；

(1.2)确定构造轨迹矩阵的参数：所述轨迹矩阵R的所述目标的原始RCS序列长度N，信号滑动窗口的长度L，信号样本数M；

其中，L需保证信号滑动窗口至少包含两个采样周期的序列信息，信号样本数M＝N-L+1，若N为偶数，则L＝N/2；若N为奇数，则L＝(N+1)/2；

(1.3)令所述目标的原始RCS序列中的第n个单元为RCS(n)，n∈{0,1,…,N-1}，构造所述轨迹矩阵R：

其中，

m∈{0,1,…,M-1}；上标T表示矩阵的转置。

3.根据权利要求1所述的方法，其特征在于，所述步骤2具体包括：

(2.1)分别计算

中所有元素的平均值

得到所述轨迹矩阵R的平均向量

其中，

(2.2)计算所述轨迹矩阵R的自相关矩阵C_R:

其中，E_1×(N-M+1)为1×(N-M+1)阶的矩阵，矩阵中的元素全为1，所述自相关矩阵C_R为M×M阶的矩阵；

(2.3)计算所述自相关矩阵C_R的对角化矩阵C_Y和所述空间转换矩矩阵A：

C_Y＝AC_RA^T＝diag[λ₀,λ₁,…,λ_M-1]；

其中，diag[λ₀,λ₁,…,λ_M-1]表示对角线元素为λ₀,λ₁,…λ_M-1的矩阵，λ₀≥λ₁≥…≥λ_M-1，λ₀,λ₁,…λ_M-1为所述自相关矩阵C_R的所有特征值，所述空间转换矩阵A的行是λ₀,λ₁,…λ_M-1对应的特征向量；所述对角化矩阵C_Y为M×M阶的矩阵，所述自相关矩阵为M×M阶的矩阵；

(2.4)计算所述主分量Y，

所述主分量Y为M×(N-M+1)阶的矩阵。

4.根据权利要求1所述的方法，其特征在于，所述步骤3具体包括：

(3.1)选取重构参数q＝1，将M×M阶矩阵A截断为M×1阶矩阵A_q，计算所述重构矩阵

其中，所述重构矩阵

为一维向量，

(3.2)从所述重构矩阵

中提取重构RCS序列：令所述重构RCS序列第l个单元为RCS(l)，将所述重构矩阵

中的元素作为所述重构RCS序列中对应位置处的单元，令RCS(l)＝r_l，l∈{1,2,...,N-M+1}。

5.如权利要求1述的方法，其特征在于，步骤4包括以下子步骤：

(4.1a)构造所述重构RCS序列的循环自相关函数：

其中，mod(l+k,N-M+1)表示将l+k按照模N-M+1取余；

(4.2a)构造所述重构后RCS序列的循环平均幅度差函数：

(4.3a)构造目标函数：

(4.4a)确定所述目标函数取最大值时对应的k的取值，假设为k₀，k₀∈{0,1,...,N-M+1}，则所述重构RCS序列的周期值为：

T＝k₀/f_s，

其中，f_s为所述目标的原始RCS序列的采样频率；

(4.5a)根据所述重构RCS序列的周期值，计算得到所述目标进动频率：

f＝1/T。

6.根据权利要求1所述的方法，其特征在于，所述步骤4包括以下子步骤：

(4.1b)构造所述重构RCS序列的离散傅里叶变换式：

其中，j为虚数单位，j²＝-1，；

(4.2b)确定所述目标函数取最大值时对应的k的取值，假设为k₀，

则所述重构目标进动频率估计值表示为：

其中，f_s为所述目标的原始RCS序列的采样频率。

Images