一种基于高阶粒子滤波的目标微多普勒曲线提取方法与流程

本发明属于雷达目标识别技术领域，涉及一种基于高阶粒子滤波的目标微多普勒曲线提取方法，适用于针对正弦形式调制的微多普勒曲线提取。

背景技术：

微多普勒曲线通常被认为能够独一无二地表征目标的运动特性，具有为雷达目标的辨识与分类提供可靠信息的能力；然而，如何有效地提取目标的微多普勒曲线在实际应用中仍是一项挑战。

通常，提取微多普勒曲线的方法大致可分为两类，第一类是根据微多普勒信号模型建立子信号字典，将目标的雷达回波与建立的字典进行匹配，从而估计得到相应的模型参数；这类方法需要目标微多普勒信号模型的先验信息，当模型参数较多时，需要构建高维度的子信号字典，从而导致运算量较大；以自旋目标为例，即使自旋目标平动得到了完全的补偿，子信号字典也是三维的，因为其包含了三个未知的模型参数；第二类方法是利用时频分析方法来获得时间-瞬时多普勒图，并根据时间-瞬时多普勒图中的瞬时多普勒曲线来估计目标的模型参数，常用的时频分析方法有短时傅里叶变换(shorttimefouriertransform，stft)、wigner-ville分布、s方法、时变自回归模型(time-varyingautoregressivemodel，tvar)等；所述常用的时频分析方法的计算复杂度通常相对较低，但是需要清晰的时间-瞬时多普勒图，且要求在单个距离单元内的散射中心个数尽可能地少；前者通常很难得到保证，后者可以通过加大雷达发射信号的带宽实现。

radon变换是目前低信噪比环境中从时频面上提取正弦曲线的一种普遍使用的方法，其需要假设目标平动被完全补偿；然而，这一假设通常很难达到，尤其对高频雷达而言，微小的补偿误差也可能导致非常大的微多普勒调制；因此，有必要通过三阶或者更高阶的速度来描述目标的平动，但这会导致radon变换的计算复杂度急剧增加；实际上，在低信噪比条件下从时频面上提取微多普勒曲线的问题与微弱目标的检测前跟踪(trackingbeforedetection，tbd)问题一致，对目标回波幅度的起伏非常稳健，只是其主要针对一阶马尔可夫过程，无法直接扩展到微多普勒曲线的提取中。

技术实现要素：

针对上述现有技术存在的不足，本发明的目的在于提出一种基于高阶粒子滤波的目标微多普勒曲线提取方法，该种基于高阶粒子滤波的目标微多普勒曲线提取方法能够用于正弦形式调制的微多普勒曲线提取，并且能够有效估计正弦微多普勒曲线的参数。

本发明的实现思路是：将正弦调制的瞬时多普勒曲线看作为一个高阶马尔可夫过程，建立目标的动态方程；将目标雷达回波的短时傅里叶变换结果看作为观测，建立目标的观测模型；利用核平滑方法估计目标的静态模型参数，再采用辅助粒子滤波方法估计目标的状态参数，从而提取得到目标的微多普勒曲线。

为达到上述技术目的，本发明采用以下技术方案予以实现。

一种基于高阶粒子滤波的目标微多普勒曲线提取方法，其特征在于，包括以下步骤：

步骤1，确定雷达，所述雷达检测范围内存在目标，将目标相对于雷达的运动分解为目标平动和目标微动，进而得到m次目标雷达回波信号的短时傅里叶变换结果；其中，m为获取目标雷达回波信号的总次数，m为大于或等于1的正整数；

步骤2，使用高阶粒子滤波法建立第k时刻目标瞬时多普勒频率的状态模型；其中，0≤k≤k-1，k表示m次目标雷达回波信号的短时傅里叶变换结果的离散时间长度，k为大于1的正整数；

其中所述m次目标雷达回波信号的短时傅里叶变换结果包含k×n个离散短时傅里叶变换单元，n表示m次目标雷达回波信号的短时傅里叶变换结果包含的离散频率点总个数，k表示m次目标雷达回波信号的短时傅里叶变换结果的离散时间长度；

步骤3，根据第k时刻目标瞬时多普勒频率的状态模型，确定第k时刻n个离散短时傅里叶变换单元的目标观测值构成的观测向量，进而得到第k时刻目标的观测模型；其中，n≤m；

初始化：令k表示第k时刻，k＝0,1,2,…,k-1，k的初始值为0，k表示m次目标雷达回波信号的短时傅里叶变换结果的离散时间长度；

步骤4，根据第k时刻目标的观测模型，计算得到第k+1时刻目标平动的瞬时多普勒频率对应的状态变量估计值和第k+1时刻目标微动的瞬时多普勒频率对应的状态变量估计值

步骤5，根据第k+1时刻目标平动的瞬时多普勒频率对应的状态变量估计值和第k+1时刻目标微动的瞬时多普勒频率对应的状态变量估计值计算得到第k+1时刻目标的瞬时多普勒频率估计值

步骤6，令k的值分别取0至k-1，重复步骤4至步骤5，进而分别得到第1时刻目标的瞬时多普勒频率估计值至第k时刻目标的瞬时多普勒频率估计值然后根据第1时刻目标的瞬时多普勒频率估计值至第k时刻目标的瞬时多普勒频率估计值绘制得到曲线，得到的所述曲线为提取的目标微多普勒曲线。

本发明的有益效果：将正弦调制的微多普勒曲线提取的过程看作是状态估计问题，将目标的瞬时多普勒频率看作为状态建立动态方程，目标雷达回波的短时傅里叶变换结果看作为观测建立观测方程，利用辅助粒子滤波结合核平滑方法估计得到目标的状态参数，从而实现对正弦微多普勒曲线的提取。

附图说明

下面结合附图说明和具体实施方式对本发明作进一步详细说明。

图1为本发明的一种基于高阶粒子滤波的目标微多普勒曲线提取方法流程图；

图2为电磁仿真的目标模型示意图；

图3为电磁仿真目标的时间-距离像示意图；

图4a为无噪声情况下短时傅里叶变换方法对光滑自旋目标的微多普勒分析结果图；

图4b为5db信噪比下短时傅里叶变换方法对光滑自旋目标的微多普勒分析结果图；

图4c为0db信噪比下短时傅里叶变换方法对光滑自旋目标的微多普勒分析结果图；

图4d为-5db信噪比下短时傅里叶变换方法对光滑自旋目标的微多普勒分析结果图；

图5a为无噪声情况下的瞬时多普勒曲线提取的结果图；

图5b为本发明对5db信噪比下的瞬时多普勒曲线提取的结果图；

图5c为本发明对0db信噪比下的瞬时多普勒曲线提取的结果图；

图5d为本发明对-5db信噪比下的瞬时多普勒曲线提取的结果图。

具体实施方式

参照图1，为本发明的一种基于高阶粒子滤波的目标微多普勒曲线提取方法流程图；其中所述基于高阶粒子滤波的目标微多普勒曲线提取方法，包括以下步骤：

步骤1，建立正弦微多普勒调制的雷达回波模型及其离散短时傅里叶变换的信号形式。

确定雷达，所述雷达检测范围内存在目标，将目标相对于雷达的运动分解为目标平动和目标微动，目标平动指目标沿着雷达视线方向所做的径向运动，目标微动指目标或目标的组成部分除质心平动以外的振动、转动和加速度等微小运动，目标微动过程中目标上的各点相对于雷达都会产生不同的运动形态，这也就对应了不同的多普勒频移量，记为目标微动的瞬时多普勒频率；然后分别设定目标平动的初始速度为v0，目标平动的二阶速度为a，目标平动的三阶速度为b；设定目标微动的频率为ω，目标微动的幅度为e、目标微动的初始相位为φ0；将雷达发射信号经目标反射回到雷达的回波信号，记为目标雷达回波信号，分别设定目标雷达回波信号的幅度为ρ，设定目标雷达回波信号的初始相位为ψ0，设定雷达发射信号的波长为λ，目标雷达回波信号的脉冲重复周期为tr；进而计算得到第m次目标雷达回波信号x(m)，其表达式为：

其中，令m为目标雷达回波信号的序号，m＝0,…,m-1，第m次对应mtr时刻，m为获取目标雷达回波信号的总次数，m为大于或等于1的正整数；u(m)为第m次目标雷达回波信号x(m)的噪声，且u(m)服从均值为0、方差为的复高斯分布；exp表示指数函数，j表示虚数单位，ρ表示目标雷达回波信号的幅度，v0表示目标平动的初始速度，tr表示目标雷达回波信号的脉冲重复周期。

采用短时傅立叶变换方法对第0次目标雷达回波信号x(0)至第m-1次目标雷达回波信号x(m-1)进行短时傅立叶变换，进而得到m次目标雷达回波信号的短时傅里叶变换结果。

分别设定k为m次目标雷达回波信号的短时傅里叶变换结果的离散时间变量，l为m次目标雷达回波信号的短时傅里叶变换结果的离散频率变量，则m次目标雷达回波信号的短时傅里叶变换结果的第(k,l)个离散短时傅里叶变换单元的信号形式为z(k,l)：

其中，x[(n-nol)k+n]表示第(n-nol)k+n次目标雷达回波信号，n＝0,…,n-1，w(n)表示m次目标雷达回波信号的短时傅里叶变换结果中第n个离散频率点处的窗函数，n表示m次目标雷达回波信号的短时傅里叶变换结果包含的离散频率点总个数，nol表示窗函数w(n)滑动时与m次目标雷达回波信号的重叠点数，n≤m，k表示m次目标雷达回波信号的短时傅里叶变换结果的离散时间长度，表示向下取整，n、m、k分别为大于或等于1的正整数。

步骤2，将目标的瞬时多普勒频率看作为状态，使用高阶粒子滤波法建立第k时刻目标瞬时多普勒频率的状态模型。

其中，目标的瞬时多普勒频率分为两部分来对待-平动部分引起的瞬时多普勒频率与微动部分引起的瞬时多普勒频率，对其分别建立动态方程，从而得到第k时刻目标瞬时多普勒频率的状态模型。

根据目标的雷达回波信号(1)，可得第m次目标雷达回波信号x(m)的瞬时多普勒频率γm可表示为：

其中，ω＝ωtr，v0表示目标平动的初始速度，a表示目标平动的二阶速度，tr表示目标雷达回波信号的脉冲重复周期，λ表示雷达发射信号的波长，b表示目标平动的三阶速度，e表示目标微动的幅度，ω表示目标微动的频率。

为了采用粒子滤波来提取目标的瞬时多普勒频率曲线，将第0次目标雷达回波信号x(0)的瞬时多普勒频率γ0至第m-1次目标雷达回波信号x(m-1)的瞬时多普勒频率γm-1看作为状态，考虑本发明目标观测值为m次目标雷达回波信号的短时傅里叶变换结果，并将目标观测值的总观测时间长度记为k'，k'与k取值相等；设定短时傅里叶变换的步长为△t，则k+1时刻目标平动的状态方程为：

其中，dk+1为第k+1时刻目标平动的初始速度v0引起的瞬时多普勒频率，dk为第k时刻目标平动的初始速度v0引起的瞬时多普勒频率，ηk+1为第k+1时刻目标平动的二阶速度a引起的瞬时多普勒频率，ηk为第k时刻目标平动的二阶速度a引起的瞬时多普勒频率，为第k+1时刻目标平动的三阶速度b引起的瞬时多普勒频率，为第k时刻目标平动的三阶速度b引起的瞬时多普勒频率；为第k+1时刻目标平动的三阶速度b引起的瞬时多普勒频率的过程噪声，且服从均值为0、方差为的高斯分布。

令第k+1时刻目标平动的状态向量为βk+1，上标t表示转置；目标平动的状态转移矩阵为ψ，第k+1时刻目标平动的过程噪声向量为ξk+1，为第k+1时刻目标平动的三阶速度b引起的瞬时多普勒频率的过程噪声；进而得到第k+1时刻目标平动的状态方程为：

βk+1＝ψβk+ξk+1(6)

令第k+1时刻目标微动的状态变量为gk+1，则第k+1时刻目标微动的状态方程为：

其中，tk+1＝tk+△t，gk表示第k时刻目标微动的状态变量，gk＝a'cos(ωtk+φ0)，tk+1表示第k+1时刻，tk表示第k时刻，△t表示设定的时间间隔，φ0表示目标微动的初始相位；ω＝ωtr，tr表示目标雷达回波信号的脉冲重复周期，λ表示雷达发射信号的波长，e表示目标微动的幅度，ω表示目标微动的频率，sin表示正弦函数，cos表示余弦函数；sgn(·)表示符号函数，满足：

其中，由于目标观测值在时间维度上是离散化的，所以式(7)中的微分部分无法直接计算，需要采用一阶或高阶导数近似计算获得。

采用p阶导数近似后，1≤p≤p，p表示设定的最高阶数，p为大于或等于1的正整数，本实施例中设定的最高阶数p取值为4；进而计算得到p阶导数近似后第k+1时刻目标微动的状态方程：

其中，ω＝ωtr，tr表示目标雷达回波信号的脉冲重复周期，ω表示目标微动的频率，△t表示短时傅里叶变换的步长，d^p(gk)表示第k时刻目标微动的状态变量gk的第p阶导数，sgn表示符号函数；分别将第k时刻目标微动的状态变量gk的第1阶导数记为d⁽¹⁾(gk)，将第k时刻目标微动的状态变量gk的第2阶导数记为d⁽²⁾(gk)，将第k时刻目标微动的状态变量gk的第3阶导数记为d⁽³⁾(gk-1)，将第k时刻目标微动的状态变量gk的第4阶导数记为d⁽⁴⁾(gk-1)，其表达式分别为：

d⁽¹⁾(gk)＝(gk-gk-1)(10a)

根据式(9)可以看出，gk与gk-p到gk-1这p个时刻的历史状态相关，所以目标微动的瞬时多普勒频率是一个高阶动态过程；需要指出的是，第k时刻目标微动的状态变量gk的第p阶导数仅当观测时刻k≥p时存在。

结合式(5)给出的第k时刻目标平动的状态方程，并考虑实际应用中目标微动的频率与幅度均存在随机扰动，因此，建立第k时刻目标瞬时多普勒频率的状态模型为：

其中，gk+1,p表示p阶导数近似后第k+1时刻目标微动的状态变量，gk,p表示p阶导数近似后第k时刻目标微动的状态变量，βk+1表示第k+1时刻目标平动的状态变量，βk表示第k时刻目标平动的状态变量，ξk+1表示第k+1时刻目标平动的过程噪声向量，υk+1表示p阶导数近似后第k+1时刻目标微动的状态变量gk+1,p的过程噪声，且υk+1服从均值为0、方差为的高斯分布；d^p(gk)表示第k时刻目标微动的状态变量gk的第p阶导数，sgn表示符号函数，ψ表示目标平动的状态转移矩阵。

此时第k时刻目标的瞬时多普勒频率γk可表示为：

γk＝tβk+gk(12)

其中，t表示设定的第k时刻的状态向量，且t＝[100]。

步骤3，根据第k时刻目标瞬时多普勒频率的状态模型，将m次目标雷达回波信号的短时傅里叶变换结果看作为观测，建立第k时刻目标的观测模型。

目标观测值为m次目标雷达回波信号的短时傅里叶变换结果，所述m次目标雷达回波信号的短时傅里叶变换结果包含k×n个离散短时傅里叶变换单元，n表示m次目标雷达回波信号的短时傅里叶变换结果包含的离散频率点总个数，k表示m次目标雷达回波信号的短时傅里叶变换结果的离散时间长度。

令zk为第k时刻n个离散短时傅里叶变换单元的目标观测值构成的观测向量，zk＝[z(k,0),…,z(k,l),…,z(k,n-1)]^t，z(k,l)表示m次目标雷达回波信号的短时傅里叶变换结果的第(k,l)个离散短时傅里叶变换单元的信号形式，0≤k≤k-1,0≤l≤n-1，上标t表示转置。

给定第k时刻目标瞬时多普勒频率的状态变量为{gk,p,βk}，设定包含在p阶导数近似后第k时刻目标微动的状态变量gk,p中的静态模型参数为θk，θk＝[a,ω]^t，令第(k,l)个离散短时傅里叶变换单元的点散布函数为h(k,l；gk,p,θk,βk)，令hk为第k时刻k×n个离散短时傅里叶变换单元的点散布函数向量。

此外，第(k,l)个离散短时傅里叶变换单元的点散布函数h(k,l；gk,p,θk,βk)可能存在未知的相移exp{jφk}，φk表示第(k,l)个离散短时傅里叶变换单元的相移相位，φk服从[0,2π)上的均匀分布；考虑第(k,l)个离散短时傅里叶变换单元的加性噪声nk的存在，nk服从均值为0、方差为的复高斯分布，in表示n维单位矩阵；第k时刻n个离散短时傅里叶变换单元的目标观测值构成的观测向量zk满足：

zk＝exp{jφk}hk+nk(13)

观测方程采用目标的联合似然比函数(即存在目标时的似然函数与不存在目标时的似然函数之比)来描述；假设第(k,l)个离散短时傅里叶变换单元的点散步函数的相移相位φ已知，则第k时刻n个离散短时傅里叶变换单元的目标观测值构成的观测向量zk的目标条件联合似然比函数为l(zk|gk,p,θk,βk,φk)：

其中，φk表示第(k,l)个离散短时傅里叶变换单元的相移相位；将l(zk|gk,p,θk,βk,φ)在第(k,l)个离散短时傅里叶变换单元的相移相位φk的分布范围[0,2π)上进行积分，进而得到第k时刻n个离散短时傅里叶变换单元的目标观测值构成的观测向量zk的联合似然函数比l(zk|gk,p,θk,βk)，其表达式为：

其中，上标h表示共轭转置，i0(·)表示修正bessel函数，θk表示包含在p阶导数近似后第k时刻目标微动的状态变量gk,p中的静态模型参数，表示p阶导数近似后第k+1时刻目标微动的状态变量gk+1,p的过程噪声υk+1的方差，dφk表示φk的微分。

所述第k时刻n个离散短时傅里叶变换单元的目标观测值构成的观测向量zk的联合似然函数比l(zk|gk,p,θk,βk)为第k时刻目标的观测模型。

初始化：令k表示第k时刻，k＝0,1,2,…,k-1，k的初始值为0。

步骤4，根据第k时刻目标的观测模型，利用辅助粒子滤波算法辅助粒子滤波方法估计第k时刻目标雷达回波信号的状态变量和静态模型参数{gk,p,θk,βk}，得到第k+1时刻目标平动的瞬时多普勒频率对应的状态变量估计值和第k+1时刻目标微动的瞬时多普勒频率对应的状态变量估计值

具体地，根据第k时刻目标的观测模型，利用辅助粒子滤波算法结合核平滑方法估计第k时刻目标雷达回波信号的状态变量和静态模型参数{gk,p,θk,βk}，θk表示包含在p阶导数近似后第k时刻目标微动的状态变量gk,p中的静态模型参数，θk＝[a,ω]^t，ω＝ωtr。

设定第k时刻粒子集合为所述第k时刻粒子集合为中包含ns个粒子，ns表示粒子总个数，ns为大于0的正整数；表示第k时刻第i个粒子中目标微动的瞬时多普勒频率对应的状态变量，表示第k时刻第i个粒子中目标平动的瞬时多普勒频率对应的状态变量，表示第k时刻第i个粒子中的目标静态模型参数；其中，第0时刻粒子集合为表示第0时刻第i个粒子中目标微动的瞬时多普勒频率对应的状态变量，表示第0时刻第i个粒子中目标平动的瞬时多普勒频率对应的状态变量，表示第0时刻第i个粒子中的目标静态模型参数。

(1)分别计算第k时刻第i个粒子对应的先验估计量即第k+1时刻第i个粒子中目标微动的瞬时多普勒频率对应的状态变量的先验估计量为第k+1时刻第i个粒子中目标平动的瞬时多普勒频率对应的状态变量的先验估计量为第k时刻第i个粒子中的目标静态模型参数的先验估计量为其表达式分别为：

其中，～表示从给定分布中产生随机数，h表示设定的平滑因子，0<h<1，表示包含在p阶导数近似后第k时刻目标微动的状态变量gk,p中的静态模型参数θk的均值估计量，表示第k时刻第i个粒子的第一级权重；当k取值为0时，

(2)计算第k+1时刻第i个粒子的第一级权重为：

其中，∝表示比例关系，表示对应第k+1时刻n个离散短时傅里叶变换单元的目标观测值构成的观测向量zk+1的联合似然函数比，表示第k时刻第i个粒子对应的先验估计量。

(3)对第k+1时刻第i个粒子的第一级权重作归一化处理，得到归一化处理后第k+1时刻第i个粒子的第一级权重

(4)令i分别取1至ns，重复执行(3)，进而分别得到归一化处理后第k+1时刻第1个粒子的第一级权重至归一化处理后第k+1时刻第ns个粒子的第一级权重并记为归一化处理后第k+1时刻ns个粒子的第一级权重

为了改善粒子滤波器的退化问题，根据归一化处理后第k时刻ns个粒子的第一级权重并利用对称重采样算法对第k时刻粒子集合进行重采样，以减少小权值粒子的数目，提高有效粒子数目，进而得到第k时刻新的粒子集合所述第k时刻新的粒子集合中包含ns个粒子，ns表示粒子总个数，ns为大于0的正整数；表示第k时刻第i^j个粒子中目标微动的瞬时多普勒频率对应的状态变量，表示第k时刻第i^j个粒子中目标平动的瞬时多普勒频率对应的状态变量，表示第k时刻第i^j个粒子中的目标静态模型参数，i^j表示第j次从1…ns个粒子中重复抽取后得到的粒子标号，i^j∈{1,…,ns}，∈表示属于。

(5)确定均值为方差为h²vk+1的高斯分布表示第k时刻第i^j个粒子中的目标静态模型参数的先验估计量，h表示设定的平滑因子，0<h<1；并从均值为方差为h²vk+1的高斯分布中随机产生第k+1时刻第j个粒子中的静态模型参数

其中，θk+1表示包含在p阶导数近似后第k+1时刻目标微动的状态变量gk+1,p中的静态模型参数，vk表示第k时刻ns个粒子的方差估计量，

(6)从概率分布中随机产生第k+1时刻第j个粒子的目标微动状态变量以及从概率分布中随机产生第k+1时刻第j个粒子的目标平动状态变量

其中，表示给定及后第k+1时刻目标微动的瞬时多普勒频率对应的状态变量的后验概率分布，gk+1,p表示p阶导数近似后第k+1时刻目标微动的状态变量，表示给定后第k+1时刻目标平动的瞬时多普勒频率对应的状态变量的后验概率分布，表示第k时刻第i^j个粒子中目标平动的瞬时多普勒频率对应的状态变量。

(7)估计第k+1时刻第j个粒子的第二级权重j＝1,…,ns：

其中，表示给定第k+1时刻新的粒子集合对应第k+1时刻n个离散短时傅里叶变换单元的目标观测值构成的观测向量zk+1的联合似然函数比，表示对应第k+1时刻n个离散短时傅里叶变换单元的目标观测值构成的观测向量zk+1的后验概率，表示第k+1时刻新的粒子集合中第i^j个粒子对应的先验估计量，表示第k+1时刻第i^j个粒子中目标微动的瞬时多普勒频率对应的状态变量，表示第k+1时刻第i^j个粒子中目标平动的瞬时多普勒频率对应的状态变量，表示第k+1时刻第i^j个粒子中的目标静态模型参数，zk+1表示第k+1时刻n个离散短时傅里叶变换单元的目标观测值构成的观测向量。

(8)对第k+1时刻第j个粒子的第二级权重作归一化处理，得到归一化处理后第k+1时刻第j个粒子的第二级权重

(9)利用第k+1时刻新的粒子集合的加权平均值作为第k+1时刻ns个粒子的状态变量估计值和第k+1时刻ns个粒子的模型静态参数估计值其表达式分别为：

其中，表示第k+1时刻目标微动的瞬时多普勒频率对应的状态变量估计值，表示第k+1时刻目标平动的瞬时多普勒频率对应的状态变量估计值，表示包含在p阶导数近似后第k+1时刻目标微动的状态变量gk+1,p中的静态模型参数估计值，表示第k+1时刻第j个粒子中目标微动的瞬时多普勒频率对应的状态变量，表示第k+1时刻第j个粒子中目标平动的瞬时多普勒频率对应的状态变量，表示第k+1时刻第j个粒子中的目标静态模型参数，θk+1表示包含在p阶导数近似后第k+1时刻目标微动的状态变量gk+1,p中的静态模型参数。

步骤5，根据第k+1时刻目标平动的瞬时多普勒频率对应的状态变量估计值和第k+1时刻目标微动的瞬时多普勒频率对应的状态变量估计值计算得到第k+1时刻目标的瞬时多普勒频率估计值

具体地，将状态变量和模型静态参数的估计值代入公式(12)，计算得到第k+1时刻目标的瞬时多普勒频率估计值为其表达式为：

步骤6，令k的值分别取0至k-1，重复步骤4至步骤5，进而分别得到第1时刻目标的瞬时多普勒频率估计值至第k时刻目标的瞬时多普勒频率估计值然后根据第1时刻目标的瞬时多普勒频率估计值至第k时刻目标的瞬时多普勒频率估计值绘制得到曲线，得到的所述曲线为提取的目标微多普勒曲线。

通过以下电磁仿真实验对本发明效果作进一步验证说明。

(1)电磁仿真实验参数设置

本实验采用专业3维电磁仿真软件cststudiosuite2015产生雷达观测样本。目标采用简单的光滑圆锥体模型，如图2所示，材质为理想良导体；其中，圆锥体的底面半径为0.1m，高度为0.3m。

仿真的雷达参数设置如下：载频为10.01ghz，频率范围为9.5ghz～10.5ghz，脉冲重复频率为360hz，驻留时间为0.5s；目标的自旋角频率为2πrad/s，平动的残余速度为v0＝1.5m/s,二阶残余速度为a＝0.01m/s²，三阶残余速度为b＝-0.09m/s³；信噪比定义为单次回波的平均能量与噪声功率之比；微多普勒分析通过stft方法实现，其采用副瓣衰减系数为100db、长度为128点的chebyshev窗函数；辅助粒子滤波的阶数为5，产生的粒子数为500。

(2)电磁仿真实验内容

1)检验旋锥体目标锥顶与锥底对应的散射中心在距离上是可分离的，如图3所示，目标的距离-时间像通过产生一组宽带数据来获得；图3中，横坐标为时间，单位为(s)，纵坐标为距离，坐标为(m)。

2)检验stft方法在不同信噪比下的微多普勒分析性能如图4a、图4b、图4c、图4d所示，其中，图4a、图4b、图4c、图4d分别是stft方法对无噪声、5db、0db、及-5db信噪比下光滑自旋目标的微多普勒分析结果图；图4a、图4b、图4c、图4d中的横坐标都为时间，单位为(s)，纵坐标都为多普勒频率，单位为(hz)。

3)检验本发明在不同信噪比下的瞬时多普勒曲线提取性能如图5a、图5b、图5c、图5d所示，其中，图5a、图5b、图5c、图5d分别是本发明在无噪声、5db、0db、及-5db信噪比下的瞬时多普勒曲线提取结果图；图5a、图5b、图5c、图5d中的横坐标都为时间，单位为(s)，纵坐标都为多普勒频率，单位为(hz)。

(3)电磁仿真实验结果分析：

1)图3展示了根据电磁仿真数据获得的距离-时间像随时间的变化情况，从图3中可清晰地看到两条平行的直线，表明自旋锥体目标锥顶与锥底对应的散射中心在距离上是可分离的。

2)从图4a、图4b、图4c、图4d中可看出，随着信噪比的降低，stft方法对电磁仿真数据的微多普勒分析性能逐步下降，时频图中的微多普勒曲线的畸变程度逐渐加大；当信噪比为-5db时，时频图中的微多普勒曲线无法清晰显现，几乎被噪声淹没。

3)从图5a、图5b、图5c、图5d中可看出，本发明方法能够有效提取噪声环境中的微多普勒曲线；当无噪声或信噪比为5db时，提取得到的微多普勒曲线结果非常接近真值，当信噪比为0db或-5db时，提取得到的微多普勒曲线存在部分失真，但结果基本满意。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。