CN104360336B - 一种自适应提取雷达目标微动周期的新方法 - Google Patents

Abstract

该发明公开了一种自适应提取雷达目标微动周期的新方法，本发明属于仅利用雷达接收目标回波数据对所探测运动目标的微动周期进行估计的技术范畴，涉及利用雷达目标回波数据的时频图的轮廓提取和曲线拟合方法对目标的微动周期进行较精确估计的技术。该方法利用雷达接收的目标回波数据先通过短时傅里叶变换得到信号的时频分布，然后对得到的时频分布曲线进行轮廓提取，接着通过自适应设置门限的方法提取时频曲线变化的轮廓，之后采用多项式拟合的方法对得到的粗轮廓进行细节平滑，最后通过对平滑后的波形图检测谱峰并给出目标微动周期的估计。一种低成本、自适应、稳健的高精度雷达目标微动周期估计新方法。

Description

一种自适应提取雷达目标微动周期的新方法

技术领域

本发明属于仅利用雷达接收目标回波数据对所探测运动目标的微动周期进行估计的技术范畴，涉及利用雷达目标回波数据的时频图的轮廓提取和曲线拟合方法对目标的微动周期进行较精确估计的技术。

背景技术

雷达目标识别问题在雷达研究和应用领域中具有重要的意义。雷达发射电磁波作用到目标物体上，再反射回来被雷达所接收，其回波信号将携带目标的本身的特征信息。当目标物体与雷达间存在相对运动时，回波的载波频率就会发生偏移，这被称为多普勒效应(Doppler effect)。如果雷达目标除了主体的平移运动外还存在整体或部分的震动或转动等微运动，那么可能会产生一个额外的频率调制，使得回波信号在多普勒频移的基础上产生边带，这就是微多普勒效应(micro-Doppler effect)。微多普勒效应使得我们能够通过研究回波信号的时频分布信息分析目标物体的动力学特性，并且提供了一种精细分析微动目标的新方法。

文献Victor C.Chen，Fayin Lin，Micro-Doppler effect in radar:phenomenon,model,and simulation study,IEEETransactions on aerospace and electronicsystems，Vol.42,No.1，pp.1-21,2006.对雷达目标进行各种微动对回波的调制情况进行了理论推导和分析。研究结果表面，雷达目标做不同的微动运动会对雷达回波的相位进行不同调制，该理论为利用雷达回波信号分析目标的微动特性提供了理论依据。上述文献仅给出了弹头、人等雷达目标的微多普勒运动模型和时频分析结果，验证了短时傅里叶变换在微多普勒分析中的有效性，但并没有给出基于该时频图的微多普勒周期的有效提取方法。针对微动周期提取问题，国内的国防科技大学、电子科技大学、南京14所、西安电子科技大学等高校也进行了初步研究如文献：金光虎，朱玉鹏，高勋章，黎湘，信号处理[J]，第25卷，第5期，pp.772-776,2009.，并取得了一定的效果，但大部分算法当时频曲线较为模糊、且存在毛刺时，算法的微动周期提取效果较差。

针对上述问题，本发明创造性的给出了利用雷达目标回波数据的时频图轮廓的自适应提取和曲线拟合方法对目标的微动周期进行估计的新方法，通过仿真和实测数据验证，验证了本发明所提算法的有效性和可行性。对较为模糊的时频曲线，本发明提出的微动周期提取方法具有较好的稳健性。

发明内容

在对雷达目标回波进行时频变换的基础上，对得到的时频图运用局部特征提取和曲线拟合方法实现雷达回波信号的周期估计，本发明的目的是提供一种自适应、低复杂度、稳健的雷达目标微动周期提取的新方法。

为了实现本发明的目的，将采用以下技术方案：一种适用于微动目标的雷达回波信号周期估计方法，该方法利用雷达接收的目标回波数据先通过短时傅里叶变换(STFT)得到信号的时频分布，然后对得到的时频分布曲线进行轮廓(局部细节)提取，接着通过自适应设置门限的方法提取时频曲线变化的轮廓，之后采用多项式拟合的方法对得到的粗轮廓进行细节平滑，最后通过对平滑后的波形图检测谱峰并给出目标微动周期的估计。因而该方法包括以下步骤：

步骤1，回波数据采集：雷达位于坐标原点，以一定的PRF发射窄脉冲，采集雷达目标回波信号数据；提取每次脉冲作用下目标回波信号一维距离像幅值的最大值和该幅值最大值出现的位置，形成一矢量数据x(n),n＝1,2,…,N，其中N表示目标回波的个数；

步骤2，时频分析：对该矢量数据进行短时傅里叶变换，获得维数为m×n的时频变换矩阵X(m为频率轴变量，n为时间轴变量)；记录下时间轴的单位时间间距Δt。

步骤3，时频曲线的提取，获取包含微动信息的局部时频图：对t＝1,2,…,n，保存时频变换矩阵X中每一列中数据模的最大值P(t)及该最大值在该列的位置序号Loc(t)；以位置序号Loc(t)中的最大值和最小值为上限和下限对原始时频图进行裁剪，只保留时频曲线出现的局部时频图，得到包含时频曲线信息的时频图X₁，其维数为m₁×n,m₁＜＜m；

步骤4，对时频图X₁进行去毛刺和细化：对X₁的每一列数据计算其均值μ_t和标准差σ_t，并对每一列数据预设门限为p_t＝μ_t+λσ_t，其中λ∈[0,1/2]之间的数值，每一列内的各值与门限进行比较后，保存所有大于该门限值的位置序号矢量其中F_t的第j个元素f_j∈N⁺为位置序号的值，是正整数；j＝1,…,M_t，M_t为第t列数据中大于门限p_t的元素个数，并取所有序号的均值形成新的位置向量Y＝[y₁,…,y_t]^T，该矢量中第t个元素y_t的计算公式如下：其中，floor(·)表示向下取整操作；同时，找出时频图X₁中所有元素模的最大值，对t＝1,2,…,n，把时频图X₁中位置(t,y_t)的像素值用最大值替换原值，其他位置全部用1填补，形成处理后的时频图X₂，该时频图X₂为时频曲线的大致轮廓；

步骤5，多项式平滑：利用多项式平滑算法对位置矢量Y进行平滑，得到平滑后矢量Y₁；

步骤6，搜索平滑后矢量Y₁的谱峰；

步骤7，估计周期：微动周期T的估计为两个相邻谱峰之间的间隔L乘以时频图时间轴相邻点的间距Δt，即当时频曲线只存一个谱峰及一个谱谷时，周期的估计为相邻谱峰和谱谷之间时间差的2倍。

所述步骤2中的短时傅里叶变换，是一种线性的时频变换，它的表达形式为：

$X(n,\mathrm{\ω})=\underset{m=-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\Σ}}x\left(m\right)w(n-m)e^{-j\mathrm{\ω}m}$

其中，x(n)为短时傅里叶变换前的已知信号，w(n)为短时傅里叶变换的窗函数，可以根据对估计的均值和方差的不同需要选择Bartlett、Kaiser、Hamming、Blackman等窗函数,X(n,ω)为信号短时傅里叶变换后的时频分布图。

所述步骤4中的门限选择的原则是：方差σ_i越大，则λ选择趋近于0的值，反之，λ则选择趋近与1/2的数值。

所述步骤5中多项式拟合，表达式为

p(x)＝p₁xⁿ+p₂x^n-1+…+p_nx+p_n+1

其中n为拟合的多项式阶数，p₁到p_n+1为拟合的多项式系数。

所述步骤5中多项式拟合表达式的阶数n＝10。

本发明利用短时傅里叶变换(STFT)处理雷达接收的目标回波数据，得到雷达回波信号的功率谱，再通过局部特征提取和多项式曲线拟合对功率谱曲线进行细化和平滑，最后通过谱峰搜索的方法求出周期，解决了微多普勒效应下雷达回波信号周期不易求出的问题，对回波周期估计有着较高精确度。在微多普勒条件下雷达回波信号周期的求解对于分析目标的微动特性有着重要意义，本发明方法解决了雷达回波信号周期求解的问题，具有实际的应用价值。

附图说明

图1为雷达和目标微动的几何示意图。

图2为本发明的流程图。

图3为短时傅立叶变换的时频图。

图4为对时频图曲线进行预处理后的结果。

图5为去毛刺和细化后的微动特性曲线。

图6为多项式平滑后的微动特性曲线。

具体实施方式

下面结合附图对本发明详细说明。

如图1所示(U,V,W)是雷达坐标系，雷达位于坐标系的原点Q。假设雷达目标为圆柱体且位于雷达坐标系中的O点。另外，假设有一个随其运动的目标坐标系(x,y,z)，相对雷达坐标系存在平移和转动。为了描述目标的转动，引入一个参考坐标系(X,Y,Z)，它与目标坐标系有着相同的原点和同样的平移运动，但是没有相对于雷达坐标系的转动，如图1所示。

假设目标是刚体，其相对雷达的平移速度为并且在参考坐标系中的旋转角速度为可表示为于是，对于目标物体上的一个参考点P，将从P移动到P'。这个过程可以看作P先平移到P”，再从P”旋转到P'，最终位置P'可以表示为

目标的多普勒频率为

$f_D=\frac{2f}{c}{\[\stackrel{\→}{V}+\stackrel{\→}{\mathrm{\ω}}\×\stackrel{\→}{r}\]}_{radial}---\left(2\right)$

其中，f是雷达发射信号的载波频率，c是光速，下标_radial代表沿雷达径向。上式第一项是由目标平移运动产生的多普勒频移，第二项则是目标旋转运动的微多普勒效应，表示为：

$f_{micro-Doppler}=\frac{2f}{c}{\[\stackrel{\→}{\mathrm{\ω}}\×\stackrel{\→}{r}\]}_{radial}---\left(3\right)$

由此可见，如果能够求出雷达接收信号的多普勒频移信息，就可以精确分析目标的微动特征。周期提取算法的流程图如图2所示.。

下面以雷达具体回波数据为例进行分析：

已知目标做进动运动，微动周期为0.65s。雷达工作频率为10GHz，PRF为1ms，在每个观测周期连续发射60个脉冲，总共观测27组回波数据，总回波点数为1620点的复数。对1620点数据进行短时傅立叶变换。窗函数选择160点的Bartlett窗(约为总数据长度的1/10)，短时傅立叶变换的原始时频图如图3所示；经过预处理后的时频图如图4所示；利用步骤四进行去毛刺和细化后的曲线如图5所示，显然该图无法提取周期结果。平滑多项式拟合表达式的阶数n选择为10阶，阶数越高，得到的曲线的去毛刺效果越好，但会带来计算量的增加，不利于实验结果表面10阶已经足够保证去毛刺效果。进一步利用多项式平滑后的结果如图6所示，从该图很容易看出周期性。在图6的基础上经过谱峰搜索，得到相邻谱峰之间的间隔L＝700-54＝646，再乘以时间间隔Δt＝0.001s，最终的微动周期估计结果为和真实微动周期的误差为0.004s。上述实测数据结果验证了本发明算法的有效性。

Claims (5)

1.一种自适应提取雷达目标微动周期的新方法，该方法包括：

步骤1，回波数据采集：雷达位于坐标原点，以一定的PRF发射窄脉冲，采集雷达目标回波信号数据；提取每次脉冲作用下目标回波信号一维距离像幅值的最大值和该幅值的最大值出现的位置，形成一矢量数据x(n),n＝1,2,…,N，其中N表示目标回波的个数；

步骤2，时频分析：对该矢量数据进行短时傅里叶变换，获得维数为m×n的时频变换矩阵X，其中m为频率轴变量，n为时间轴变量；记录下时间轴的单位时间间距△t；

步骤3，时频曲线的提取，获取包含微动信息的局部时频图：对t＝1,2,…,n，保存时频变换矩阵X中每一列中数据模的最大值P(t)及该最大值在该列的位置序号Loc(t)；以位置序号Loc(t)中的最大值和最小值为上限和下限对原始时频图进行裁剪，只保留时频曲线出现的局部时频图，得到包含时频曲线信息的紧凑时频图X₁，其维数为m₁×n,m₁＜＜m；

步骤4，对紧凑时频图X₁进行去毛刺和细化：对紧凑时频图X₁的每一列数据计算其均值μ_t和标准差σ_t，并对每一列数据预设门限为p_t＝μ_t+λσ_t，其中λ∈[0,1/2]之间的数值，每一列内的各值与门限进行比较后，保存所有大于该门限值的位置序号矢量其中F_t的第j个元素f_j∈N⁺为位置序号的值，是正整数；j＝1,…,M_t，M_t为第t列数据中大于门限p_t的元素个数，并取所有序号的均值形成新的位置向量Y＝[y₁,…,y_t]^T，该矢量中第t个元素y_t的计算公式如下：其中，floor(·)表示向下取整操作；同时，找出紧凑时频图X₁中所有元素模的最大值，对t＝1,2,…,n，把紧凑时频图X₁中位置(t,y_t)的像素值用最大值替换原值，其他位置全部用1填补，形成处理后的时频图X₂，该时频图X₂为时频曲线的大致轮廓；

步骤5，多项式平滑：利用多项式平滑算法对位置矢量Y进行平滑，得到平滑后矢量Y₁；

步骤6，搜索平滑后矢量Y₁的谱峰；

步骤7，估计周期：微动周期T的估计为两个相邻谱峰之间的间隔L乘以时频图时间轴相邻点的间距△t，即当时频曲线只存一个谱峰及一个谱谷时，周期的估计为相邻谱峰和谱谷之间时间差的2倍。

2.如权利要求1所述的一种自适应提取雷达目标微动周期的新方法，其特征在于所述步骤2中的短时傅里叶变换，是一种线性的时频变换，它的表达形式为：

$X(n,\mathrm{\ω})=\underset{m=-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\Σ}}x\left(m\right)w(n-m)e^{-j\mathrm{\ω}m}$

其中，x(n)为短时傅里叶变换前的已知信号，w(n)为短时傅里叶变换的窗函数，可以根据对估计的均值和方差的不同需要选择Bartlett、Kaiser、Hamming、Blackman等窗函数,X(n,ω)为信号短时傅里叶变换后的时频分布图。

3.如权利要求1所述的一种自适应提取雷达目标微动周期的新方法，其特征在于所述步骤4中的门限选择的原则是：方差σ_i越大，则λ选择趋近于0的值，反之，λ则选择趋近于1/2的数值。

4.如权利要求1所述的一种自适应提取雷达目标微动周期的新方法，其特征在于所述步骤5中多项式拟合，表达式为

p(x)＝p₁xⁿ+p₂x^n-1+…+p_nx+p_n+1

其中n为拟合的多项式阶数，p₁到p_n+1为拟合的多项式系数。

5.如权利要求4所述的一种自适应提取雷达目标微动周期的新方法，其特征在于所述步骤5中多项式拟合表达式的阶数n＝10。