CN104007430A - 基于瞬时调频率估计的进动目标的微多普勒提取方法 - Google Patents

Abstract

本发明公开了一种基于瞬时调频率估计的进动目标的微多普勒提取方法,其步骤为：步骤1，建立空间进动锥体目标的等效散射点模型；雷达接收锥体目标的回波信号，利用等效散射点模型从回波信号中获取回波序列；步骤2，用松弛解线调频的方法估计回波序列的调频率；步骤3，通过随机抽样一致性算法处理回波序列的调频率，得到P条调频率曲线；步骤4，对P条调频率曲线分别做积分运算，得到P条瞬时微多普勒频率的估计曲线。本发明实现在低信噪比的情况下也能得到比较高的估计精度。本发明可用于分析空间进动锥体目标的微多普勒频率。

Description

基于瞬时调频率估计的进动目标的微多普勒提取方法

技术领域

本发明属于雷达技术领域，涉及瞬时多普勒频率的估计方法，尤其涉及一种基于瞬时调频率估计的进动目标的微多普勒提取方法，用于估计空间进动目标的微多普勒频率。

背景技术

表面光滑的目标在大气层外飞行时，为了保持姿态的稳定性，需要做自旋运动，当受到横向干扰时，在做自旋运动的同时还会绕空间另一定向轴做锥旋运动，这种运动形式被称为进动。美国海军实验室的V.C.Chen将进动定义为微运动的一种，并将微运动产生的雷达回波中的多普勒调制命名为微多普勒调制。由于目标的微运动可以反映其结构、尺寸、微动频率等重要物理特性，因此微多普勒调制对空间目标识别与参数估计有着重要的意义，而正确的从雷达回波中提取出微多普勒调制则是利用微动信息进行目标识别与参数估计的前提。

现有的瞬时多普勒频率估计方法大致可以分为非参数化方法和参数化方法两类。非参数化方法一般首先计算信号的时频分布，然后通过跟踪时频分布图中的峰值来得到频率随时间的变化。例如邵长宇等提出的基于多目标跟踪的空间锥体目标微多普勒频率提取方法，该方法首先应用经典的短时傅里叶变换得到回波的时频分布图，然后通过多目标跟踪技术跟踪时频分布图上的时频曲线，从而达到提取目标各等效散射中心微多普勒频率的目的，但这种方法中用到的时频分析的方式存在着时间分辨率与频率分辨率相矛盾的问题，难以得到较高的瞬时频率估计精度；瞬时频率估计的另一类方法是参数化方法，与非参数化方法的区别在于需要利用信号的先验信息建立一个参数化的模型来描述信号，例如韩勋等提出的基于时变自回归模型的瞬时频率估计方法，该方法首先用p阶时变自回归模型来表示目标的回波信号，求得模型的时变参数，然后通过求信号功率谱密度极点的相位来估计瞬时频率，该方法虽然在高信噪比情况下可以得到较高的估计精度，但是在低信噪比情况下估计精度会有所下降，另外还存在模型阶数选择的问题。

发明内容

针对上述现有技术的不足，本发明提出一种基于瞬时调频率估计的进动目标的微多普勒提取方法，实现在低信噪比的情况下也能得到比较高的估计精度。

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

一种基于瞬时调频率估计的进动目标的微多普勒提取方法，其特征在于，包括以下步骤：

步骤1，建立空间进动锥体目标的等效散射点模型；根据等效散射点模型得到各等效散射点瞬时微多普勒频率的导数公式；

步骤2，雷达接收锥体目标的回波信号，再从回波信号中获取回波序列S_r(t_m)；

步骤3，对回波序列S_r(t_m)做短时傅里叶变换，得到回波序列的时频图，通过回波序列的时频图得到可视等效散射点的个数G以及正弦曲线的周期数目M；

步骤4，将回波序列分为4LM段的观测时间的回波段，L为正整数；

步骤5，将每一观测时间的回波段近似为T个线性调频信号分量的叠加，设定T等于可视等效散射点的数目G，用松弛解线调频的方法得到每个线性调频信号分量的调频率估计值，总共得到4LM*G个线性调频信号分量的调频率估计值；*表示乘号；

步骤6，根据各等效散射点瞬时微多普勒频率的导数公式建立等效散射点瞬时微多普勒频率导数模型，然后用随机抽样一致性算法处理4LM*G个线性调频信号分量的调频率估计值，得到G条调频率曲线；

步骤7，对G条调频率曲线分别做积分运算，得到G条瞬时微多普勒频率的估计曲线。

上述技术方案的特点和进一步改进在于：

(1)步骤1包括以下子步骤：

1a)所述空间进动锥体目标为无尾翼的光滑锥体目标，其等效散射点为点P₁、点P₂和点P₃，设定点P₁、点P₂和点P₃为无尾翼的光滑锥体目标的等效散射点，点P₁和点P₂为锥体目标的底部边缘上的两点，点P₃为锥顶上的点，锥体目标高度为H，旋转中心距底面距离为h，底面半径为r，雷达位于雷达坐标系(U,V,W)的原点Q，锥体目标位于雷达远场区域，方位角为α，俯仰角为β，(X,Y,Z)表示与雷达坐标系平行的参考坐标系，锥体目标坐标系(x,y,z)，初始时刻锥体目标坐标系与雷达坐标系平行，两坐标系原点之间的距离为R₀，锥体目标坐标系的原点为锥体目标的质心，锥体目标坐标系的z轴为锥体目标的对称轴，锥旋轴OC在yOz平面内，建立如下式所示的新坐标系(X_n,Y_n,Z_n)：

Z_n＝M(t)·z₀

X_n＝rLOS×Z_n

Y_n＝Z_n×X_n

其中，M(t)表示锥体目标的微动矩阵，z₀表示初始时刻锥体目标对称轴的单位方向矢量，即z₀＝(0,0,1)^T，上标T是转置符号，rLOS表示雷达视线在初始时刻锥体目标坐标系中的单位方向向量，即rLOS＝(cosα·sinβ,sinα·sinβ,-cosβ)^T，设定α＝90°，则rLOS＝(0,sinβ,-cosβ)^T，符号×表示叉乘；

1b)设定锥体目标只存在进动或平动已经被完全补偿，即锥体目标除绕其对称轴Oz以角速度ω_s做自旋运动外，还绕锥旋轴OC以角速度ω_c做锥旋运动，轴OC和轴Oz之间的夹角为θ，称为进动角；

微动矩阵M(t)公式具体形式如下：

$M\left(t\right)=\left[\begin{array}{ccc}\cos \left({\mathrm{\ω}}_{c}t\right)& -\cos \mathrm{\θ}\sin \left({\mathrm{\ω}}_{c}t\right)& \sin \mathrm{\θ}\sin \left({\mathrm{\ω}}_{c}t\right)\\ \cos \mathrm{\θ}\sin \left({\mathrm{\ω}}_{c}t\right)& 1-{\cos }^2\mathrm{\θ}[1-\cos \left({\mathrm{\ω}}_{c}t\right)]& \sin \mathrm{\θ}\cos \mathrm{\θ}[1-\cos \left({\mathrm{\ω}}_{c}t\right)]\\ -\sin \mathrm{\θ}\sin \left({\mathrm{\ω}}_{c}t\right)& \sin \mathrm{\θ}\cos \mathrm{\θ}[1-\cos \left({\mathrm{\ω}}_{c}t\right)]& 1-{\sin }^2\mathrm{\θ}[1-\cos \left({\mathrm{\ω}}_{c}t\right)]\end{array}\right]$

在锥体目标飞行过程中，由于存在遮挡效应，使得一些等效散射点无法被雷达波照射到，如附图3所示，在等效散射点P₁、P₂和P₃的新坐标系(X_n,Y_n,Z_n)中，γ为锥体目标的半锥角，rLOS表示雷达视线，表示姿态角，也就是锥体目标对称轴与雷达视线的夹角，等效散射点的遮挡效应由姿态角与锥体目标的半锥角γ共同决定；

设定姿态角在[0,π]范围内变化，各个等效散射点的遮挡情况如下表3所示：

表3

上表中Y表示遮挡，N表示不遮挡；从上表中看到，点P₁总是不被遮挡，点P₂在时被遮挡，而点P₃在时被遮挡；

1c)等效散射点P₁、P₂和P₃距离雷达的瞬时径向距离公式为：

$\left.\begin{array}{c}{R}_1\left(t\right)=R_0-r\·\sqrt{1-{\cos }^2\mathrm{\θ}{\cos }^2(\mathrm{\β}+\mathrm{\θ})}+h\·\cos \mathrm{\θ}\cos (\mathrm{\β}+\mathrm{\θ})\\ +(h\·\sin \mathrm{\θ}\sin (\mathrm{\β}+\mathrm{\θ})+r\·\frac{\sin \mathrm{\θ}\cos \mathrm{\θ}\sin (\mathrm{\β}+\mathrm{\θ})\cos (\mathrm{\β}+\mathrm{\θ})}{\sqrt{1-{\cos }^2\mathrm{\θ}{\cos }^2(\mathrm{\β}+\mathrm{\θ})}})\cos \left({\mathrm{\ω}}_{c}t\right))\\ R_2\left(t\right)=R_0+r\·\sqrt{1-{\cos }^2\mathrm{\θ}{\cos }^2(\mathrm{\β}+\mathrm{\θ})}+h\·\cos \mathrm{\θ}\cos (\mathrm{\β}+\mathrm{\θ})\\ +(h\·\sin \mathrm{\θ}\sin (\mathrm{\β}+\mathrm{\θ})-r\·\frac{\sin \mathrm{\θ}\cos \mathrm{\θ}\sin (\mathrm{\β}+\mathrm{\θ})\cos (\mathrm{\β}+\mathrm{\θ})}{\sqrt{1-{\cos }^2\mathrm{\θ}{\cos }^2(\mathrm{\β}+\mathrm{\θ})}})\cos \left({\mathrm{\ω}}_{c}t\right))\\ R_3\left(t\right)=R_0-(H-h)\cos \mathrm{\θ}\cos (\mathrm{\β}+\mathrm{\θ})-(H-h)\sin \mathrm{\θ}\sin (\mathrm{\β}+\mathrm{\θ})\cos \left({\mathrm{\ω}}_{c}t\right)\end{array}\right\}---\left(1\right)$

其中，R₁(t)表示等效散射点P₁在任意时刻t距离雷达的瞬时径向距离，R₂(t)表示等效散射点P₂在任意时刻t距离雷达的瞬时径向距离，R₃(t)表示等效散射点P₃在任意时刻t距离雷达的瞬时径向距离，R₀表示初始时刻锥体目标与雷达之间的距离，H表示锥体目标高度，h表示旋转中心距底面的距离，r表示底面圆半径，θ表示进动角，β表示俯仰角，ω_c表示锥旋角速度；

等效散射点P₁、P₂和P₃的瞬时微多普勒频率公式为：

$\left.\begin{array}{c}{\mathrm{fd}}_1\left(t\right)=\frac{2}{\mathrm{\λ}}{\mathrm{\ω}}_c(h\sin \mathrm{\θ}\sin (\mathrm{\β}+\mathrm{\θ})+r\frac{\sin \mathrm{\θ}\cos \mathrm{\θ}\sin (\mathrm{\β}+\mathrm{\θ})\cos (\mathrm{\β}+\mathrm{\θ})}{\sqrt{1-{\cos }^2\mathrm{\θ}{\cos }^2(\mathrm{\β}+\mathrm{\θ})}})\sin \left({\mathrm{\ω}}_{c}t\right)\\ {\mathrm{fd}}_2\left(t\right)=\frac{2}{\mathrm{\λ}}{\mathrm{\ω}}_c(h\sin \mathrm{\θ}\sin (\mathrm{\β}+\mathrm{\θ})-r\frac{\sin \mathrm{\θ}\cos \mathrm{\θ}\sin (\mathrm{\β}+\mathrm{\θ})\cos (\mathrm{\β}+\mathrm{\θ})}{\sqrt{1-{\cos }^2\mathrm{\θ}{\cos }^2(\mathrm{\β}+\mathrm{\θ})}})\sin \left({\mathrm{\ω}}_{c}t\right)\\ {\mathrm{fd}}_3\left(t\right)=-\frac{2}{\mathrm{\λ}}(H-h){\mathrm{\ω}}_c\sin \mathrm{\θ}\sin (\mathrm{\β}+\mathrm{\θ})\sin \left({\mathrm{\ω}}_{c}t\right)\end{array}\right\}---\left(2\right)$

其中，fd₁(t)表示等效散射点P₁在任意时刻t的瞬时微多普勒频率，fd₂(t)表示等效散射点P₂在任意时刻t的瞬时微多普勒频率，fd₃(t)表示等效散射点P₃在任意时刻t的瞬时微多普勒频率，H表示锥体目标高度，h表示旋转中心距底面的距离，r表示底面圆半径，θ表示进动角，β表示俯仰角，ω_c表示锥旋角速度，λ表示雷达发射信号的波长；

等效散射点P₁、P₂和P₃的瞬时微多普勒频率的导数公式为：

$\left.\begin{array}{c}{\mathrm{fd}}_1^{\′}\left(t\right)=\frac{2}{\mathrm{\λ}}{{\mathrm{\ω}}_c}^2(h\sin \mathrm{\θ}\sin (\mathrm{\β}+\mathrm{\θ})+r\frac{\sin \mathrm{\θ}\cos \mathrm{\θ}\sin (\mathrm{\β}+\mathrm{\θ})\cos (\mathrm{\β}+\mathrm{\θ})}{\sqrt{1-{\cos }^2\mathrm{\θ}{\cos }^2(\mathrm{\β}+\mathrm{\θ})}})\cos \left({\mathrm{\ω}}_{c}t\right)\\ {\mathrm{fd}}_2^{\′}\left(t\right)=\frac{2}{\mathrm{\λ}}{{\mathrm{\ω}}_c}^2(h\sin \mathrm{\θ}\sin (\mathrm{\β}+\mathrm{\θ})-r\frac{\sin \mathrm{\θ}\cos \mathrm{\θ}\sin (\mathrm{\β}+\mathrm{\θ})\cos (\mathrm{\β}+\mathrm{\θ})}{\sqrt{1-{\cos }^2\mathrm{\θ}{\cos }^2(\mathrm{\β}+\mathrm{\θ})}})\cos \left({\mathrm{\ω}}_{c}t\right)\\ {\mathrm{fd}}_3^{\′}\left(t\right)=-\frac{2}{\mathrm{\λ}}(H-h){{\mathrm{\ω}}_c}^2\sin \mathrm{\θ}\sin (\mathrm{\β}+\mathrm{\θ})\cos \left({\mathrm{\ω}}_{c}t\right)\end{array}\right\}---\left(3\right)$

其中，fd₁′(t)表示等效散射点P₁在任意时刻t的瞬时微多普勒频率的导数，fd₂′(t)表示等效散射点P₂在任意时刻t的瞬时微多普勒频率的导数，fd₃′(t)表示等效散射点P₃在任意时刻t的瞬时微多普勒频率的导数，H表示锥体目标高度，h表示旋转中心距底面的距离，r表示底面圆半径，θ表示进动角，β表示俯仰角，ω_c表示锥旋角速度，λ表示雷达发射信号的波长。

(2)步骤2具体包括以下子步骤：

2a)雷达系统发射线性调频信号，并接收到空间进动锥体目标的回波信号；

2a1)雷达发射线性调频信号，线性调频信号公式如下：

$S_t(\hat{t},t_m)=\mathrm{rect}\left(\frac{\hat{t}}{t_p}\right)\exp (j2{\mathrm{\πf}}_{c}t+\mathrm{j\πv}{\hat{t}}^2)$

其中， $\mathrm{rect}\left(u\right)=\left\{\begin{array}{cc}1& \left|u\right|\≤\frac{1}{2}\\ 0& \left|u\right|>\frac{1}{2}\end{array}\right.,\hat{t}=t-m{T}_r$ 表示快时间，t表示全时间，T_r表示脉冲重复时间，t_m＝m T_r表示慢时间，m为整数，f_c表示载频，ν表示调频率，t_p表示脉冲宽度；

2a2)设定锥体目标在不同遮挡情况下，锥体目标上的可视等效散射点个数为G，经过一定的时延，根据线性调频信号公式得到接收到的回波信号公式如下：

$S_r(\hat{t},t_m)=\underset{g=1}{\overset{G}{\mathrm{\Σ}}}{\mathrm{\σ}}_g\mathrm{rect}\left(\frac{\hat{t}-\frac{{2R}_g\left(t_m\right)}{C}}{t_p}\right)*\exp (j2\mathrm{\π}{f}_c(t-\frac{{2R}_g\left(t_m\right)}{C})+\mathrm{j\πv}{(\hat{t}-\frac{{2R}_g\left(t_m\right)}{C})}^2)$

其中，R_g(t_m)表示第m个周期第g个可视等效散射点距离雷达的瞬时径向距离，可视等效散射点也就是未被遮挡的等效散射点，σ_g表示第g个可视等效散射点的后向散射系数，C表示光速，g大于等于1并且小于等于G，G表示可视等效散射点的个数，最大值是3；

2a3)将各等效散射点距离雷达的瞬时径向距离公式(1)代入步骤2a2)中的接收到的回波信号公式，得到空间进动锥体目标的回波信号；

2b)将回波信号去载频，则去载频后的回波信号的表达式如下：

$S_r(\hat{t},t_m)=\underset{g=1}{\overset{G}{\mathrm{\Σ}}}{\mathrm{\σ}}_g\mathrm{rect}\left(\frac{\hat{t}-\frac{{2R}_g\left(t_m\right)}{C}}{t_p}\right)*\exp (-j2\mathrm{\π}{f}_c\frac{{2R}_g\left(t_m\right)}{C})+\mathrm{j\πv}{(\hat{t}-\frac{{2R}_g\left(t_m\right)}{C})}^2)$

其中， $\mathrm{rect}\left(u\right)=\left\{\begin{array}{cc}1& \left|u\right|\≤\frac{1}{2}\\ 0& \left|u\right|>\frac{1}{2}\end{array}\right.,$ R_g(t_m)表示第m个周期第g个可视等效散射点距离雷达的瞬时径向距离，σ_g表示第g个可视等效散射点的后向散射系数，C表示光速，f_c表示载频，ν表示调频率，t_p表示脉冲宽度，表示快时间，t_m＝m T_r表示慢时间，m为整数，T_r表示脉冲重复时间；

2c)令其中，R_g(t_m)表示第m个周期第g个可视等效散射点距离雷达的瞬时径向距离，C表示光速，τ_g(t_m)表示第m个周期第g个可视等效散射点的时延，代入去载频后的回波信号的表达式之后，去载频后的回波信号的表达式化简为：

$S_r(\hat{t},t_m)=\underset{g=1}{\overset{G}{\mathrm{\Σ}}}{\mathrm{\σ}}_g\mathrm{rect}\left(\frac{\hat{t}-{\mathrm{\τ}}_g\left(t_m\right)}{t_p}\right)\exp (-j2\mathrm{\π}{f}_c{\mathrm{\τ}}_g\left(t_m\right)+\mathrm{j\πv}{(\hat{t}-{\mathrm{\τ}}_g\left(t_m\right))}^2)$

2d)对去载频后的回波信号进行脉压处理，得到回波信号进行脉压之后的表达式：

$S_r(\hat{t},t_m)=\underset{g=1}{\overset{G}{\mathrm{\Σ}}}{\mathrm{\σ}}_g\exp (-j2\mathrm{\π}{f}_c{\mathrm{\τ}}_g\left(t_m\right))\sin c\left(\mathrm{\πB}(\hat{t}-{\mathrm{\τ}}_g\left(t_m\right))\right)$

其中，σ_g表示第g个可视等效散射点的后向散射系数，f_c表示载频，τ_g(t_m)表示第m个周期第g个可视等效散射点的时延，B表示信号带宽，表示快时间；

2e)取回波信号进行脉压之后的表达式中每一慢时间内的最大值，得到回波序列S_r(t_m)，回波序列S_r(t_m)的表达式为：

$S_r\left(t_m\right)=\underset{g=1}{\overset{G}{\mathrm{\Σ}}}{\mathrm{\σ}}_g\exp (-j2\mathrm{\π}{f}_c{\mathrm{\τ}}_g\left(t_m\right))=\underset{g=1}{\overset{G}{\mathrm{\Σ}}}{\mathrm{\σ}}_g\exp (-j4\mathrm{\π}{f}_c\frac{R_g\left(t_m\right))}{C})$

其中，σ_g表示第g个可视等效散射点的后向散射系数，R_g(t_m)表示第m个周期第g个可视等效散射点距离雷达的瞬时径向距离，τ_g(t_m)表示第m个周期第g个可视等效散射点的时延，C表示光速，f_c表示载频。

(3)步骤5具体包括以下子步骤；

5a)将每一观测时间的回波段近似为T个线性调频信号分量的叠加，T个线性调频信号分量的叠加之后，叠加后每一观测时间的回波段信号模型表示为：

$S\left(t_m\right)=\underset{q=1}{\overset{T}{\mathrm{\Σ}}}{\mathrm{\δ}}_q\exp \left(j2\mathrm{\π}(f_{0_q}{t}_m+\frac{1}{2}{\mathrm{\μ}}_q{t_m}^2)\right)$

其中，δ_q表示第q个线性调频信号分量的复幅度，表示第q个线性调频信号分量的初始频率，μ_q表示第q个线性调频信号分量的调频率，T表示线性调频信号分量个数，设定T等于可视等效散射点的数目G，t_m＝m T_r表示慢时间，m为整数，T_r表示脉冲重复时间，S(t_m)表示叠加之后每一观测时间的回波段；

5b)先设定线性调频信号分量个数T＝1，即单分量情况，设单分量线性调频信号为：

$s\left(t_m\right)=\mathrm{\δexp}\left(j2\mathrm{\π}(f_0{t}_m+\frac{1}{2}{{\mathrm{\μt}}_m}^2)\right)---\left(4\right)$

其中，δ表示单分量线性调频信号的复幅度，f₀表示单分量线性调频信号的初始频率，μ表示单分量线性调频信号的调频率，t_m＝m T_r表示慢时间，m为整数，T_r表示脉冲重复时间；

用解线调频的方法得到单分量线性调频信号的复幅度δ的估计值初始频率f₀的估计值调频率μ的估计值

5c)对于每一观测时间的回波段的T个线性调频信号分量，用松弛解线调频算法估计每个线性调频信号分量的复幅度、初始频率和调频率，得到每个线性调频信号分量的复幅度的估计值初始频率的估计值和调频率的估计值

(4)步骤6具体包括以下子步骤：

6a)将各等效散射点瞬时微多普勒频率的导数公式(3)简化后得到等效散射点瞬时微多普勒频率导数模型表达式a·cos(2π·bx)＝y；

6b)用随机抽样一致性算法处理每一观测时间的回波段的调频率的估计值，把估计出来的调频率区分为局内点和局外点，然后由局内点确定一条调频率曲线；

6c)根据可视等效散射点的个数G，确定G条调频率曲线；

若G＝1，根据6b)得到一条调频率曲线；

若G＝2，则可视等效散射点的个数为两个，由6b)中得到的局外点估计出模型a·cos(2π·bx)＝y中的另一组参数a和b，得到另一条调频率曲线；

若G＝3，则可视等效散射点的个数为三个，由6b)中的局内点确定一条调频率曲线；对6b)中区分的局外点再用一次随机抽样一致性算法，把6b)中区分的局外点再区分为两类不同的点，由这两类不同的点估计出模型a·cos(2π·bx)＝y中的两组不同的参数a和b，进而确定另外两条不同的调频率曲线。

(5)子步骤6b)具体包括以下子步骤：

6b1)随机选取4LM*G个线性调频信号分量的调频率估计值中的两个点，根据这两个点以及建立的等效散射点瞬时微多普勒频率导数模型a·cos(2π·bx)＝y，估计模型中的未知参数a和b，选择的两个点被设定为局内点，没有选择的点为局外点；

6b2)通过等效散射点瞬时微多普勒频率导数模型测试4LM*G个线性调频信号分量的调频率估计值中的其他点；也就是将除6b1)随机选择的两个点局内点之外的局外点代入等效散射点瞬时微多普勒频率导数模型，如果等式a·cos(2π·bx)＝y+Δy中的错误Δy小于检测阈值ζ，则将该点添加到局内点集合；

6b3)用所有局内点去重新估计模型的参数a和b，再统计局内点的个数d，以及所有局内点与模型之间的错误之和E；

6b4)重复6b1)‐6b3)，得到局内点的个数记为d'，局内点与模型之间的错误之和E'；

6b5)如果迭代之后局内点的个数d'大于迭代之前局内点的个数d，并且迭代之后的错误之和E'小于迭代之前的错误之和E，则更新局内点的个数d和局内点与模型之间的错误E，并且更新模型的参数a和b以及局内点的集合；否则，转至6b6)；

6b6)以预定的次数重复执行6b4)‐6b5)；

6b7)运行预定的次数后，输出该模型的参数a和b以及局内点的集合，由参数a和b确定一条调频率曲线。

与现有技术相比，本发明具有突出的实质性特点和显著的进步。本发明与现有方法相比，具有以下优点：

第一，对噪声不敏感，在低信噪比情况下也能得到比较高的估计精度。

第二，虽然本发明也用到了短时傅里叶变换的方法，但并不是在时频分析的基础上进行的，进行短时傅里叶变换只是为了确定可视等效散射点的个数，因此不存在非参数化方法中时间分辨率与频率分辨率相矛盾的问题。

附图说明

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

图1是本发明的流程图；

图2是空间进动锥体目标的等效散射点模型图；

图3是遮挡情况示意图；

图4是回波序列的时频图；

图5是瞬时微多普勒频率的估计曲线图；

图6是瞬时多普勒频率的估计曲线与理论曲线的对比图。

具体实施方式

参照图1，说明本发明的一种基于瞬时调频率估计的进动目标的微多普勒提取方法，本发明用于估计空间进动锥体目标的微多普勒频率，其具体步骤如下：

步骤1，建立空间进动锥体目标的等效散射点模型；根据等效散射点模型得到各等效散射点瞬时微多普勒频率的导数公式。

对于无尾翼的光滑锥体目标，在等效散射点模型中有三个散射点起主要作用，分别是锥顶和底部边缘上的两点，底部边缘上的两点为入射平面与底部边缘的交点，所谓入射平面就是锥体目标对称轴和雷达视线所确定的平面。

如附图2所示，建立空间进动锥体目标的等效散射点模型。

1a)设定点P₁、点P₂和点P₃为无尾翼的光滑锥体目标的等效散射点，点P₁和点P₂为锥体目标的底部边缘上的两点，点P₃为锥顶上的点，锥体目标高度为H，旋转中心距底面距离为h，底面半径为r，雷达位于雷达坐标系(U,V,W)的原点Q，锥体目标位于雷达远场区域，方位角为α，俯仰角为β，(X,Y,Z)表示与雷达坐标系平行的参考坐标系，锥体目标坐标系(x,y,z)，初始时刻锥体目标坐标系与雷达坐标系平行，两坐标系原点之间的距离为R₀，锥体目标坐标系的原点为锥体目标的质心，锥体目标坐标系的z轴为锥体目标的对称轴，锥旋轴OC在yOz平面内，建立如下式所示的新坐标系(X_n,Y_n,Z_n)：

Z_n＝M(t)·z₀

X_n＝rLOS×Z_n

Y_n＝Z_n×X_n

其中，M(t)表示锥体目标的微动矩阵，z₀表示初始时刻锥体目标对称轴的单位方向矢量，即z₀＝(0,0,1)^T，上标T是转置符号，rLOS表示雷达视线在初始时刻锥体目标坐标系中的单位方向向量，即rLOS＝(cosα·sinβ,sinα·sinβ,-cosβ)^T，设定α＝90°，则rLOS＝(0,sinβ,-cosβ)^T，符号×表示叉乘。

1b)设定锥体目标只存在进动或平动已经被完全补偿，即锥体目标除绕其对称轴Oz以角速度ω_s做自旋运动外，还绕锥旋轴OC以角速度ω_c做锥旋运动，轴OC和轴Oz之间的夹角为θ，称为进动角。

微动矩阵M(t)公式具体形式如下：

$M\left(t\right)=\left[\begin{array}{ccc}\cos \left({\mathrm{\ω}}_{c}t\right)& -\cos \mathrm{\θ}\sin \left({\mathrm{\ω}}_{c}t\right)& \sin \mathrm{\θ}\sin \left({\mathrm{\ω}}_{c}t\right)\\ \cos \mathrm{\θ}\sin \left({\mathrm{\ω}}_{c}t\right)& 1-{\cos }^2\mathrm{\θ}[1-\cos \left({\mathrm{\ω}}_{c}t\right)]& \sin \mathrm{\θ}\cos \mathrm{\θ}[1-\cos \left({\mathrm{\ω}}_{c}t\right)]\\ -\sin \mathrm{\θ}\sin \left({\mathrm{\ω}}_{c}t\right)& \sin \mathrm{\θ}\cos \mathrm{\θ}[1-\cos \left({\mathrm{\ω}}_{c}t\right)]& 1-{\sin }^2\mathrm{\θ}[1-\cos \left({\mathrm{\ω}}_{c}t\right)]\end{array}\right]$

在锥体目标飞行过程中，由于存在遮挡效应，使得一些等效散射点无法被雷达波照射到，如附图3所示，在等效散射点P₁、P₂和P₃的新坐标系(X_n,Y_n,Z_n)中，γ为锥体目标的半锥角，rLOS表示雷达视线，表示姿态角，也就是锥体目标对称轴与雷达视线的夹角，等效散射点的遮挡效应由姿态角与锥体目标的半锥角γ共同决定。

设定姿态角在[0,π]范围内变化，各个等效散射点的遮挡情况如下表3所示：

表3

上表中Y表示遮挡，N表示不遮挡。从上表中看到，点P₁总是不被遮挡，点P₂在时被遮挡，而点P₃在时被遮挡。

1c)等效散射点P₁、P₂和P₃距离雷达的瞬时径向距离公式为：

$\left.\begin{array}{c}{R}_1\left(t\right)=R_0-r\·\sqrt{1-{\cos }^2\mathrm{\θ}{\cos }^2(\mathrm{\β}+\mathrm{\θ})}+h\·\cos \mathrm{\θ}\cos (\mathrm{\β}+\mathrm{\θ})\\ +(h\·\sin \mathrm{\θ}\sin (\mathrm{\β}+\mathrm{\θ})+r\·\frac{\sin \mathrm{\θ}\cos \mathrm{\θ}\sin (\mathrm{\β}+\mathrm{\θ})\cos (\mathrm{\β}+\mathrm{\θ})}{\sqrt{1-{\cos }^2\mathrm{\θ}{\cos }^2(\mathrm{\β}+\mathrm{\θ})}})\cos \left({\mathrm{\ω}}_{c}t\right))\\ R_2\left(t\right)=R_0+r\·\sqrt{1-{\cos }^2\mathrm{\θ}{\cos }^2(\mathrm{\β}+\mathrm{\θ})}+h\·\cos \mathrm{\θ}\cos (\mathrm{\β}+\mathrm{\θ})\\ +(h\·\sin \mathrm{\θ}\sin (\mathrm{\β}+\mathrm{\θ})-r\·\frac{\sin \mathrm{\θ}\cos \mathrm{\θ}\sin (\mathrm{\β}+\mathrm{\θ})\cos (\mathrm{\β}+\mathrm{\θ})}{\sqrt{1-{\cos }^2\mathrm{\θ}{\cos }^2(\mathrm{\β}+\mathrm{\θ})}})\cos \left({\mathrm{\ω}}_{c}t\right))\\ R_3\left(t\right)=R_0-(H-h)\cos \mathrm{\θ}\cos (\mathrm{\β}+\mathrm{\θ})-(H-h)\sin \mathrm{\θ}\sin (\mathrm{\β}+\mathrm{\θ})\cos \left({\mathrm{\ω}}_{c}t\right)\end{array}\right\}---\left(1\right)$

其中，R₁(t)表示等效散射点P₁在任意时刻t距离雷达的瞬时径向距离，R₂(t)表示等效散射点P₂在任意时刻t距离雷达的瞬时径向距离，R₃(t)表示等效散射点P₃在任意时刻t距离雷达的瞬时径向距离，R₀表示初始时刻锥体目标与雷达之间的距离，H表示锥体目标高度，h表示旋转中心距底面的距离，r表示底面圆半径，θ表示进动角，β表示俯仰角，ω_c表示锥旋角速度。

等效散射点P₁、P₂和P₃的瞬时微多普勒频率公式为：

$\left.\begin{array}{c}{\mathrm{fd}}_1\left(t\right)=\frac{2}{\mathrm{\λ}}{\mathrm{\ω}}_c(h\sin \mathrm{\θ}\sin (\mathrm{\β}+\mathrm{\θ})+r\frac{\sin \mathrm{\θ}\cos \mathrm{\θ}\sin (\mathrm{\β}+\mathrm{\θ})\cos (\mathrm{\β}+\mathrm{\θ})}{\sqrt{1-{\cos }^2\mathrm{\θ}{\cos }^2(\mathrm{\β}+\mathrm{\θ})}})\sin \left({\mathrm{\ω}}_{c}t\right)\\ {\mathrm{fd}}_2\left(t\right)=\frac{2}{\mathrm{\λ}}{\mathrm{\ω}}_c(h\sin \mathrm{\θ}\sin (\mathrm{\β}+\mathrm{\θ})-r\frac{\sin \mathrm{\θ}\cos \mathrm{\θ}\sin (\mathrm{\β}+\mathrm{\θ})\cos (\mathrm{\β}+\mathrm{\θ})}{\sqrt{1-{\cos }^2\mathrm{\θ}{\cos }^2(\mathrm{\β}+\mathrm{\θ})}})\sin \left({\mathrm{\ω}}_{c}t\right)\\ {\mathrm{fd}}_3\left(t\right)=-\frac{2}{\mathrm{\λ}}(H-h){\mathrm{\ω}}_c\sin \mathrm{\θ}\sin (\mathrm{\β}+\mathrm{\θ})\sin \left({\mathrm{\ω}}_{c}t\right)\end{array}\right\}---\left(2\right)$

其中，fd₁(t)表示等效散射点P₁在任意时刻t的瞬时微多普勒频率，fd₂(t)表示等效散射点P₂在任意时刻t的瞬时微多普勒频率，fd₃(t)表示等效散射点P₃在任意时刻t的瞬时微多普勒频率，H表示锥体目标高度，h表示旋转中心距底面的距离，r表示底面圆半径，θ表示进动角，β表示俯仰角，ω_c表示锥旋角速度，λ表示雷达发射信号的波长。

等效散射点P₁、P₂和P₃的瞬时微多普勒频率的导数公式为：

$\left.\begin{array}{c}{\mathrm{fd}}_1^{\′}\left(t\right)=\frac{2}{\mathrm{\λ}}{{\mathrm{\ω}}_c}^2(h\sin \mathrm{\θ}\sin (\mathrm{\β}+\mathrm{\θ})+r\frac{\sin \mathrm{\θ}\cos \mathrm{\θ}\sin (\mathrm{\β}+\mathrm{\θ})\cos (\mathrm{\β}+\mathrm{\θ})}{\sqrt{1-{\cos }^2\mathrm{\θ}{\cos }^2(\mathrm{\β}+\mathrm{\θ})}})\cos \left({\mathrm{\ω}}_{c}t\right)\\ {\mathrm{fd}}_2^{\′}\left(t\right)=\frac{2}{\mathrm{\λ}}{{\mathrm{\ω}}_c}^2(h\sin \mathrm{\θ}\sin (\mathrm{\β}+\mathrm{\θ})-r\frac{\sin \mathrm{\θ}\cos \mathrm{\θ}\sin (\mathrm{\β}+\mathrm{\θ})\cos (\mathrm{\β}+\mathrm{\θ})}{\sqrt{1-{\cos }^2\mathrm{\θ}{\cos }^2(\mathrm{\β}+\mathrm{\θ})}})\cos \left({\mathrm{\ω}}_{c}t\right)\\ {\mathrm{fd}}_3^{\′}\left(t\right)=-\frac{2}{\mathrm{\λ}}(H-h){{\mathrm{\ω}}_c}^2\sin \mathrm{\θ}\sin (\mathrm{\β}+\mathrm{\θ})\cos \left({\mathrm{\ω}}_{c}t\right)\end{array}\right\}---\left(3\right)$

其中，fd₁′(t)表示等效散射点P₁在任意时刻t的瞬时微多普勒频率的导数，fd₂′(t)表示等效散射点P₂在任意时刻t的瞬时微多普勒频率的导数，fd₃′(t)表示等效散射点P₃在任意时刻t的瞬时微多普勒频率的导数，H表示锥体目标高度，h表示旋转中心距底面的距离，r表示底面圆半径，θ表示进动角，β表示俯仰角，ω_c表示锥旋角速度，λ表示雷达发射信号的波长。

步骤2，雷达接收锥体目标的回波信号，再从回波信号中获取回波序列S_r(t_m)。

2a)雷达系统发射线性调频信号，经过一定的时延，接收到空间进动锥体目标的回波信号。

2a1)雷达发射线性调频信号，线性调频信号公式如下：

$S_t(\hat{t},t_m)=\mathrm{rect}\left(\frac{\hat{t}}{t_p}\right)\exp (j2{\mathrm{\πf}}_{c}t+\mathrm{j\πv}{\hat{t}}^2)$

其中， $\mathrm{rect}\left(u\right)=\left\{\begin{array}{cc}1& \left|u\right|\≤\frac{1}{2}\\ 0& \left|u\right|>\frac{1}{2}\end{array}\right.,\hat{t}=t-m{T}_r$ 表示快时间，t表示全时间，T_r表示脉冲重复时间，t_m＝m T_r表示慢时间，m为整数，f_c表示载频，ν表示调频率，t_p表示脉冲宽度。

2a2)设定锥体目标在不同遮挡情况下，锥体目标上的可视等效散射点个数为G，经过一定的时延，根据线性调频信号公式得到接收到的回波信号公式如下：

$S_r(\hat{t},t_m)=\underset{g=1}{\overset{G}{\mathrm{\Σ}}}{\mathrm{\σ}}_g\mathrm{rect}\left(\frac{\hat{t}-\frac{{2R}_g\left(t_m\right)}{C}}{t_p}\right)*\exp (j2\mathrm{\π}{f}_c(t-\frac{{2R}_g\left(t_m\right)}{C})+\mathrm{j\πv}{(\hat{t}-\frac{{2R}_g\left(t_m\right)}{C})}^2)$

其中，R_g(t_m)表示第m个周期第g个可视等效散射点距离雷达的瞬时径向距离，可视等效散射点也就是未被遮挡的等效散射点，σ_g表示第g个可视等效散射点的后向散射系数，C表示光速，g大于等于1并且小于等于G，G表示可视等效散射点的个数，最大值是3。

2a3)将各等效散射点距离雷达的瞬时径向距离公式(1)代入步骤2a2)中的接收到的回波信号公式，得到空间进动锥体目标的回波信号。

2b)将回波信号去载频，则去载频后的回波信号的表达式如下：

$S_r(\hat{t},t_m)=\underset{g=1}{\overset{G}{\mathrm{\Σ}}}{\mathrm{\σ}}_g\mathrm{rect}\left(\frac{\hat{t}-\frac{{2R}_g\left(t_m\right)}{C}}{t_p}\right)*\exp (-j2\mathrm{\π}{f}_c\frac{{2R}_g\left(t_m\right)}{C})+\mathrm{j\πv}{(\hat{t}-\frac{{2R}_g\left(t_m\right)}{C})}^2)$

其中， $\mathrm{rect}\left(u\right)=\left\{\begin{array}{cc}1& \left|u\right|\≤\frac{1}{2}\\ 0& \left|u\right|>\frac{1}{2}\end{array}\right.,$ R_g(t_m)表示第m个周期第g个可视等效散射点距离雷达的瞬时径向距离，σ_g表示第g个可视等效散射点的后向散射系数，C表示光速，f_c表示载频，ν表示调频率，t_p表示脉冲宽度，表示快时间，t_m＝m T_r表示慢时间，m为整数，T_r表示脉冲重复时间。

2c)令其中，R_g(t_m)表示第m个周期第g个可视等效散射点距离雷达的瞬时径向距离，C表示光速，τ_g(t_m)表示第m个周期第g个可视等效散射点的时延，代入去载频后的回波信号的表达式之后，去载频后的回波信号的表达式化简为：

$S_r(\hat{t},t_m)=\underset{g=1}{\overset{G}{\mathrm{\Σ}}}{\mathrm{\σ}}_g\mathrm{rect}\left(\frac{\hat{t}-{\mathrm{\τ}}_g\left(t_m\right)}{t_p}\right)\exp (-j2\mathrm{\π}{f}_c{\mathrm{\τ}}_g\left(t_m\right)+\mathrm{j\πv}{(\hat{t}-{\mathrm{\τ}}_g\left(t_m\right))}^2)$

2d)对去载频后的回波信号进行脉压处理，得到回波信号进行脉压之后的表达式：

$S_r(\hat{t},t_m)=\underset{g=1}{\overset{G}{\mathrm{\Σ}}}{\mathrm{\σ}}_g\exp (-j2\mathrm{\π}{f}_c{\mathrm{\τ}}_g\left(t_m\right))\sin c\left(\mathrm{\πB}(\hat{t}-{\mathrm{\τ}}_g\left(t_m\right))\right)$

其中，σ_g表示第g个可视等效散射点的后向散射系数，f_c表示载频，τ_g(t_m)表示第m个周期第g个可视等效散射点的时延，B表示信号带宽，表示快时间。

2e)取回波信号进行脉压之后的表达式中每一慢时间内的最大值，得到回波序列S_r(t_m)，回波序列S_r(t_m)的表达式为：

$S_r\left(t_m\right)=\underset{g=1}{\overset{G}{\mathrm{\Σ}}}{\mathrm{\σ}}_g\exp (-j2\mathrm{\π}{f}_c{\mathrm{\τ}}_g\left(t_m\right))=\underset{g=1}{\overset{G}{\mathrm{\Σ}}}{\mathrm{\σ}}_g\exp (-j4\mathrm{\π}{f}_c\frac{R_g\left(t_m\right))}{C})$

其中，σ_g表示第g个可视等效散射点的后向散射系数，R_g(t_m)表示第m个周期第g个可视等效散射点距离雷达的瞬时径向距离，τ_g(t_m)表示第m个周期第g个可视等效散射点的时延，C表示光速，f_c表示载频。

步骤3，对回波序列S_r(t_m)做短时傅里叶变换，得到回波序列的时频图，通过回波序列的时频图得到可视等效散射点的个数G以及正弦曲线的周期数目M。

具体的，也就是将时频图中显示的完整的正弦曲线的数目作为可视等效散射点的个数G，时频图中显示的正弦曲线的周期数目记为M；

虽然本发明也用到了短时傅里叶变换的方法，但并不是在时频分析的基础上进行的，进行短时傅里叶变换只是为了确定可视等效散射点的个数，因此不存在非参数化方法中时间分辨率与频率分辨率相矛盾的问题。

步骤4，将回波序列分为4LM段的观测时间的回波段，L为正整数。

步骤5，将每一观测时间的回波段近似为T个线性调频信号分量的叠加，设定T等于可视等效散射点的数目G，用松弛解线调频的方法得到每个线性调频信号分量的调频率估计值，总共得到4LM*G个线性调频信号分量的调频率估计值；*表示乘号。

5a)将每一观测时间的回波段近似为T个线性调频信号分量的叠加，T个线性调频信号分量的叠加之后，叠加后每一观测时间的回波段信号模型表示为：

$S\left(t_m\right)=\underset{q=1}{\overset{T}{\mathrm{\Σ}}}{\mathrm{\δ}}_q\exp \left(j2\mathrm{\π}(f_{0_q}{t}_m+\frac{1}{2}{\mathrm{\μ}}_q{t_m}^2)\right)$

其中，δ_q表示第q个线性调频信号分量的复幅度，表示第q个线性调频信号分量的初始频率，μ_q表示第q个线性调频信号分量的调频率，T表示线性调频信号分量个数，设定T等于可视等效散射点的数目G，t_m＝m T_r表示慢时间，m为整数，T_r表示脉冲重复时间，S(t_m)表示叠加之后每一观测时间的回波段。

5b)先设定线性调频信号分量个数T＝1，即单分量情况，设单分量线性调频信号为：

$s\left(t_m\right)=\mathrm{\δexp}\left(j2\mathrm{\π}(f_0{t}_m+\frac{1}{2}{{\mathrm{\μt}}_m}^2)\right)---\left(4\right)$

其中，δ表示单分量线性调频信号的复幅度，f₀表示单分量线性调频信号的初始频率，μ表示单分量线性调频信号的调频率，t_m＝m T_r表示慢时间，m为整数，T_r表示脉冲重复时间。

用解线调频的方法得到单分量线性调频信号的复幅度δ的估计值初始频率f₀的估计值调频率μ的估计值

解线调频的方法具体见文献(邢孟道，保铮，冯大政.基于调幅-线性调频信号参数估计的机动目标成像方法.现代雷达,2000(6):44～49)，具体通过以下5b1)、5b2)和5b3)实现：

5b1)设定f_η(t_m)＝s(t_m)·exp(-jπηt_m ²)；η表示调频率变量，f_η(t_m)是调频率变量的函数，s(t_m)表示单分量线性调频信号，t_m表示慢时间；

5b2)改变调频率的变量η并对每次f_η(t_m)作傅里叶变换并取模，得到关于初始频率f₀和调频率的变量η的二维分布图；

5b3)从二维分布图的峰值点位置得到单分量线性调频信号的初始频率f₀的估计值和调频率μ的估计值将f_η(t_m)作傅里叶变换并取模之后的最大值所对应的复常数作为复幅度δ的估计值

5c)对于每一观测时间的回波段的T个线性调频信号分量，用松弛解线调频算法估计每个线性调频信号分量的复幅度、初始频率和调频率，得到每个线性调频信号分量的复幅度的估计值初始频率的估计值和调频率的估计值

松弛解线调频的方法是参考解线调频RELAX算法(具体见孙长印，保铮.雷达成像的近似二维模型及其超分辨算法.电子学报,1999,27(12):84～87)；具体的采用松弛解线调频算法估计每个线性调频信号分量的复幅度、初始频率和调频率的步骤如下：

5c1)设定线性调频信号分量个数K＝1，通过解线调频的方法得到设定的第一个线性调频信号分量的复幅度的估计值初始频率的估计值和调频率的估计值具体来说，也就是通过步骤5b1)、5b2)和5b3)实现。

设定将回波序列分为4LM段的观测时间的回波段后，每一回波段为S^*(t_m)，构建循环变量S_k(t_m)，表达式为下式：

$S_k\left(t_m\right)=S^{*}\left(t_m\right)-\underset{l=1,l\≠k}{\overset{K}{\mathrm{\Σ}}}{\widehat{\mathrm{\δ}}}_l\exp \left(j2\mathrm{\π}({\hat{f}}_{0_l}{t}_m+\frac{1}{2}{\widehat{\mathrm{\μ}}}_l{t_m}^2)\right)---\left(5\right)$

其中，表示设定的第l个线性调频信号分量的复幅度的估计值，表示设定的第l个线性调频信号分量的初始频率的估计值，表示设定的第l个线性调频信号分量的调频率的估计值，l＝1,2,...,K，k是K的变量，K表示设定的线性调频信号分量个数，K＝1,2,...,T；

需要说明的是，本发明中l表示在进行松弛解线调频的方法时，该松弛解线调频的方法中设定的线性调频信号分量的变量，K表示松弛解线调频的方法所设定的线性调频信号分量个数。

5c2)设定线性调频信号分量个数K＝2，利用5c1)求得的复幅度的估计值初始频率的估计值和调频率的估计值代入公式(5)，得到循环变量S₂(t_m)；

将循环变量S₂(t_m)近似为单分量线性调频信号s(t_m)，单分量线性调频信号s(t_m)见公式(4)，通过解线调频的方法得到设定的第二个线性调频信号分量复幅度的估计值初始频率的估计值和调频率的估计值

利用求得的设定的第二个线性调频信号分量复幅度的估计值初始频率的估计值和调频率的估计值代入公式(5)计算循环变量S₁(t_m)，将循环变量S₁(t_m)近似为单分量线性调频信号s(t_m)，进而重新确定设定的第一个线性调频信号分量复幅度的估计值初始频率的估计值和调频率的估计值

重复5c2)，直至满足收敛判据；收敛判据为比较相邻两次迭代前后残余信号能量的变化量，如果这个变化量小于某个预设的值ξ＝10^-3，则认为该步骤已收敛；否则，重复5c2)。

应当注意，在用解线调频的方法时，调频率的步长应充分小，(也就是调频率的变量的改变值应充分小)且利用FFT来进行傅里叶变换时应在信号后面填补充分长的零，这在很大程度上影响着估计的准确度。为减少计算量，先利用较大的调频频率步长求出一个粗略的估计，然后再利用较小的步长在估计值附近进行优化。

5c3)设定线性调频信号分量个数K＝3，利用5c2)求得的设定的第一个线性调频信号分量复幅度的估计值初始频率的估计值和调频率的估计值以及设定的第二个线性调频信号分量复幅度的估计值初始频率的估计值和调频率的估计值代入公式(5)计算循环变量S₃(t_m)；

将循环变量S₃(t_m)近似为单分量线性调频信号s(t_m)，单分量线性调频信号s(t_m)见公式(4)；通过解线调频的方法得到设定的第三个线性调频信号分量复幅度的估计值初始频率的估计值和调频率的估计值

利用设定的第三个线性调频信号分量复幅度的估计值初始频率的估计值和调频率的估计值以及设定的第一个线性调频信号分量复幅度的估计值初始频率的估计值和调频率的估计值代入公式(5)计算循环变量S₂(t_m)，将循环变量S₂(t_m)近似为单分量线性调频信号s(t_m)，进而重新确定设定的第二个线性调频信号分量复幅度的估计值初始频率的估计值和调频率的估计值

利用设定的第三个线性调频信号分量复幅度的估计值初始频率的估计值和调频率的估计值以及设定的第二个线性调频信号分量复幅度的估计值初始频率的估计值和调频率的估计值代入公式(5)计算循环变量S₁(t_m)，将循环变量S₁(t_m)近似为单分量线性调频信号s(t_m)，重新确定出设定的第一个线性调频信号分量复幅度的估计值初始频率的估计值和调频率的估计值

重复第5c3)，直至满足收敛判据。收敛判据为比较相邻两次迭代前后残余信号能量的变化量，如果这个变化量小于预设的值ξ＝10^-3，则认为该步骤已收敛；否则，重复5c3)；。

继续以上步骤，直到设定的线性调频信号分量个数K等于线性调频信号分量个数T结束。

在本发明中，把每一观测时间段中间时刻作为横坐标，估计出来的调频率作为纵坐标，这样得到若干个离散的点，同一时刻有G个点，同一时刻的这G个点表示的是这一时刻可视等效散射点瞬时微多普勒频率的导数。

步骤6，根据各等效散射点瞬时微多普勒频率的导数公式建立等效散射点瞬时微多普勒频率导数模型，然后用随机抽样一致性算法处理4LM*G个线性调频信号分量的调频率估计值，得到G条调频率曲线。

随机抽样一致性算法参考沈乐章的文章《随机抽样一致性算法RANSAC源程序和教程》。

6a)将各等效散射点瞬时微多普勒频率的导数公式(3)简化后得到等效散射点瞬时微多普勒频率导数模型表达式a·cos(2π·bx)＝y。

6b)用随机抽样一致性算法处理每一观测时间的回波段的调频率的估计值，把估计出来的调频率区分为局内点和局外点，然后由局内点确定一条调频率曲线。

具体的步骤6b)包括如下：

6b1)随机选取4LM*G个线性调频信号分量的调频率估计值中的两个点，根据这两个点以及建立的等效散射点瞬时微多普勒频率导数模型a·cos(2π·bx)＝y，估计模型中的未知参数a和b，选择的两个点被设定为局内点，没有选择的点为局外点；

6b2)通过等效散射点瞬时微多普勒频率导数模型测试4LM*G个线性调频信号分量的调频率估计值中的其他点；也就是将除6b1)随机选择的两个点局内点之外的局外点代入等效散射点瞬时微多普勒频率导数模型，如果等式a·cos(2π·bx)＝y+Δy中的错误Δy小于检测阈值ζ，则将该点添加到局内点集合；

6b3)用所有局内点去重新估计模型的参数a和b，再统计局内点的个数d，以及所有局内点与模型之间的错误之和E；

6b4)重复6b1)-6b3)，得到局内点的个数记为d'，局内点与模型之间的错误之和E'；

6b5)如果迭代之后局内点的个数d'大于迭代之前局内点的个数d，并且迭代之后的错误之和E'小于迭代之前的错误之和E，则更新局内点的个数d和局内点与模型之间的错误E，并且更新模型的参数a和b以及局内点的集合；否则，转至6b6)；

6b6)以预定的次数重复执行6b4)-6b5)；

6b7)运行预定的次数后，输出该模型的参数a和b以及局内点的集合，由参数a和b确定一条调频率曲线；

6c)根据可视等效散射点的个数G，确定G条调频率曲线。

若G＝1，根据6b)得到一条调频率曲线；

若G＝2，则可视等效散射点的个数为两个，由6b)中得到的局外点估计出模型a·cos(2π·bx)＝y中的另一组参数a和b，得到另一条调频率曲线；

若G＝3，则可视等效散射点的个数为三个，由6b)中的局内点确定一条调频率曲线；对6b)中区分的局外点再用一次随机抽样一致性算法，把6b)中区分的局外点再区分为两类不同的点，由这两类不同的点估计出模型a·cos(2π·bx)＝y中的两组不同的参数a和b，进而确定另外两条不同的调频率曲线。

步骤7，对G条调频率曲线分别做积分运算，得到G条瞬时微多普勒频率的估计曲线。

下面结合仿真实验对本发明的效果做进一步说明。

1、实验条件

实验数据为MATLAB仿真回波数据，锥体目标高为H＝0.96m，旋转中心距底面距离为h＝0.32m，底面半径为r＝0.25m，雷达方位角为α＝90°，俯仰角为β＝60°，锥体目标做进动，自旋频率f_ss＝4Hz，锥旋频率为f_cc＝2Hz，进动角为θ＝10°。雷达发射线性调频信号，载频为f_c＝10GHz，带宽为B＝1MHz，脉冲宽度为t_p＝10μs，雷达脉冲重复频率为prf＝1KHz，驻留时间为TT＝2s。

2、实验内容

2.1)对回波序列做短时傅里叶变换，得到回波序列的时频图如附图4所示；横坐标是时间，纵坐标是瞬时多普勒频率。图中显示两条完整的正弦曲线，表明可视等效散射点的个数为两个，在驻留时间内有4个周期。

2.2)将回波序列分成32段观测时间的回波段，将每一观测时间的回波段近似为两个线性调频信号分量的叠加，用松弛解线调频的方法估计每个线性调频信号分量的调频率，并用随机抽样一致性方法处理估计的调频率，得到两条调频率的曲线，对这两条调频率曲线分别做积分运算，得到瞬时微多普勒频率的估计曲线，结果如附图5所示，横坐标是时间，纵坐标是瞬时多普勒频率。

2.3)其他条件不变，加入高斯白噪声，取不同的信噪比进行实验。

3、实验结果分析

从附图4看出，有两条完整的正弦曲线，表明可视等效散射点的个数为两个，由于设定的俯仰角为β＝60°，进动角为θ＝10°，锥体目标高度为H＝0.96m，底面半径为r＝0.25m，所以姿态角的变化范围为[40°,60°]，而半锥角γ＝14°，满足对照步骤1b)中的表3，点P₂在驻留时间内被遮挡，图中显示的两条正弦曲线分别表示的是等效散射点P₃和P₁的瞬时微多普勒频率随时间的变化。

根据等效散射点P₃和P₁的瞬时微多普勒频率公式(2)得到瞬时微多普勒频率的理论曲线。

将瞬时微多普勒频率的估计曲线以及理论曲线画在同一幅图中，如附图6所示，横坐标是时间，纵坐标是瞬时多普勒频率。图中实线表示的是等效散射点P₃的瞬时多普勒频率的估计曲线，虚线表示的是等效散射点P₁的瞬时多普勒频率的估计曲线，点划线表示的是等效散射点P₃的瞬时多普勒频率的理论曲线，点线表示的是等效散射点P₁的瞬时多普勒频率的理论曲线。由图中看出，瞬时多普勒频率的估计曲线跟理论曲线基本重合。

为了更加具体地描述估计结果，现引入下式来描述估计的正确率：

$A=(1-\frac{\underset{k}{\mathrm{\Σ}}|{\mathrm{IF}}_r\left(k\right)-{\mathrm{IF}}_e\left(k\right)|}{\underset{k}{\mathrm{\Σ}}\left|{\mathrm{IF}}_r\left(k\right)\right|})*100\%$

其中，IF_e(k)表示瞬时多普勒频率的估计值，IF_r(k)表示瞬时多普勒频率的理论值。

现取附图6中横坐标每隔0.001s所对应的值分别作为P₃和P₁瞬时多普勒频率的估计值以及瞬时多普勒频率的理论值。然后将P₃的瞬时多普勒频率的估计值和理论值代入正确率表达式，得到P₃的估计正确率为96.62％，将P₁的瞬时多普勒频率的估计值和理论值代入正确率表达式，得到P₁的估计正确率为93.68％。

不同信噪比情况下，两个等效散射点P₁和P₃的瞬时多普勒频率的估计正确率见下表1：

表1

现有技术中提出的基于时变自回归模型的瞬时频率估计方法得出的估计正确率见下表2：

表2

SNR(dB)	20	15	10
				P3(％)	93.9	92.5	91.8
P1(％)	89.7	88..5	86.1

从表1和表2中可以看出，在低信噪比情况下，采用本方法估计的正确率也比较高。

Claims (6)

1.一种基于瞬时调频率估计的进动目标的微多普勒提取方法，其特征在于，包括以下步骤： 

步骤1，建立空间进动锥体目标的等效散射点模型；根据等效散射点模型得到各等效散射点瞬时微多普勒频率的导数公式； 

步骤2，雷达接收锥体目标的回波信号，再从回波信号中获取回波序列S_r(t_m)； 

步骤3，对回波序列S_r(t_m)做短时傅里叶变换，得到回波序列的时频图，通过回波序列的时频图得到可视等效散射点的个数G以及正弦曲线的周期数目M； 

步骤4，将回波序列分为4LM段的观测时间的回波段，L为正整数； 

步骤5，将每一观测时间的回波段近似为T个线性调频信号分量的叠加，设定T等于可视等效散射点的数目G，用松弛解线调频的方法得到每个线性调频信号分量的调频率估计值，总共得到4LM*G个线性调频信号分量的调频率估计值；*表示乘号； 

步骤6，根据各等效散射点瞬时微多普勒频率的导数公式建立等效散射点瞬时微多普勒频率导数模型，然后用随机抽样一致性算法处理4LM*G个线性调频信号分量的调频率估计值，得到G条调频率曲线； 

步骤7，对G条调频率曲线分别做积分运算，得到G条瞬时微多普勒频率的估计曲线。 

2.根据权利要求1所述的基于瞬时调频率估计的进动目标的微多普勒提取方法，其特征在于，步骤1具体包括以下子步骤： 

1a)所述空间进动锥体目标为无尾翼的光滑锥体目标，其等效散射点为点P₁、点P₂和点P₃，点P₁、点P₂和点P₃为无尾翼的光滑锥体目标的等效散射点，点P₁和点P₂为锥体目标的底部边缘上的两点，点P₃为锥顶上的点，锥体目标高度为H，旋转中心距底面距离为h，底面半径为r，雷达位于雷达坐标系(U,V,W)的原点Q，锥体目标位于雷达远场区域，方位角为α，俯仰角为β，(X,Y,Z)表示与雷达坐标系平行的参考坐标系，锥体目标坐标系(x,y,z)，初始时刻锥体目标坐标系与雷达坐标系平行，两坐标系原点之间的距离为R₀，锥体目标坐标系的原点为锥体目标的质心，锥体目标坐标系的z轴为锥体目标的对称轴，锥旋轴OC在yOz平面内，建立如下式所示的新坐标系(X_n,Y_n,Z_n)： 

Z_n＝M(t)·z₀

X_n＝rLOS×Z_n

Y_n＝Z_n×X_n

其中，M(t)表示锥体目标的微动矩阵，z₀表示初始时刻锥体目标对称轴的单位方向矢量，即z₀＝(0,0,1)^T，上标T是转置符号，rLOS表示雷达视线在初始时刻锥体目标坐标系中的单位方向向量，即rLOS＝(cosα·sinβ,sinα·sinβ,-cosβ)^T，设定α＝90°，则rLOS＝(0,sinβ,-cosβ)^T，符号×表示叉乘； 

1b)设定锥体目标只存在进动或平动已经被完全补偿，即锥体目标除绕其对称轴Oz以角速度ω_s做自旋运动外，还绕锥旋轴OC以角速度ω_c做锥旋运动，轴OC和轴Oz之间的夹角为θ，称为进动角； 

微动矩阵M(t)公式具体形式如下： 

在锥体目标飞行过程中，由于存在遮挡效应，使得一些等效散射点无法被雷达波照射到，如附图3所示，在等效散射点P₁、P₂和P₃的新坐标系(X_n,Y_n,Z_n)中，γ为锥体目标的半锥角，rLOS表示雷达视线，表示姿态角，也就是锥体目标对称轴与雷达视线的夹角，等效散射点的遮挡效应由姿态角与锥体目标的半锥角γ共同决定； 

设定姿态角在[0,π]范围内变化，各个等效散射点的遮挡情况如下表3所示： 

表3 

上表中Y表示遮挡，N表示不遮挡；从上表中看到，点P₁总是不被遮挡，点P₂在 时被遮挡，而点P₃在时被遮挡； 

1c)等效散射点P₁、P₂和P₃距离雷达的瞬时径向距离公式为： 

其中，R₁(t)表示等效散射点P₁在任意时刻t距离雷达的瞬时径向距离，R₂(t)表示等效散射点P₂在任意时刻t距离雷达的瞬时径向距离，R₃(t)表示等效散射点P₃在任意时刻t距离雷达的瞬时径向距离，R₀表示初始时刻锥体目标与雷达之间的距离，H表示锥体目标高度，h表示旋转中心距底面的距离，r表示底面圆半径，θ表示进动角，β表示俯仰角，ω_c表示锥旋角速度； 

等效散射点P₁、P₂和P₃的瞬时微多普勒频率公式为： 

其中，fd₁(t)表示等效散射点P₁在任意时刻t的瞬时微多普勒频率，fd₂(t)表示等效散射点P₂在任意时刻t的瞬时微多普勒频率，fd₃(t)表示等效散射点P₃在任意时刻t的瞬时微多普勒频率，H表示锥体目标高度，h表示旋转中心距底面的距离，r表示底面圆半径，θ表示进动角，β表示俯仰角，ω_c表示锥旋角速度，λ表示雷达发射信号的波长； 

等效散射点P₁、P₂和P₃的瞬时微多普勒频率的导数公式为： 

其中，fd₁′(t)表示等效散射点P₁在任意时刻t的瞬时微多普勒频率的导数，fd₂′(t)表示等效散射点P₂在任意时刻t的瞬时微多普勒频率的导数，fd₃′(t)表示等效散射点P₃在任意时刻t的瞬时微多普勒频率的导数，H表示锥体目标高度，h表示旋转中心距底面的距离，r表示底面圆半径，θ表示进动角，β表示俯仰角，ω_c表示锥旋角速度，λ表示雷达发射信号的波长。 

3.根据权利要求2所述的基于瞬时调频率估计的进动目标的微多普勒提取方法，其特征在于，步骤2具体包括以下子步骤： 

2a)雷达系统发射线性调频信号，并接收到空间进动锥体目标的回波信号； 

2a1)雷达发射线性调频信号，线性调频信号公式如下： 

其中，表示快时间，t表示全时间，T_r表示脉冲重复时间，t_m＝m T_r表示慢时间，m为整数，f_c表示载频，ν表示调频率，t_p表示脉冲宽度； 

2a2)设定锥体目标在不同遮挡情况下，锥体目标上的可视等效散射点个数为G，经过一定的时延，根据线性调频信号公式得到接收到的回波信号公式如下： 

其中，R_g(t_m)表示第m个周期第g个可视等效散射点距离雷达的瞬时径向距离，可视等效散射点也就是未被遮挡的等效散射点，σ_g表示第g个可视等效散射点的后向散射系 数，C表示光速，g大于等于1并且小于等于G，G表示可视等效散射点的个数，最大值是3； 

2a3)将各等效散射点距离雷达的瞬时径向距离公式(1)代入步骤2a2)中的接收到的回波信号公式，得到空间进动锥体目标的回波信号； 

2b)将回波信号去载频，则去载频后的回波信号的表达式如下： 

其中，R_g(t_m)表示第m个周期第g个可视等效散射点距离雷达的瞬时径向距离，σ_g表示第g个可视等效散射点的后向散射系数，C表示光速，f_c表示载频，ν表示调频率，t_p表示脉冲宽度，表示快时间，t_m＝m T_r表示慢时间，m为整数，T_r表示脉冲重复时间； 

2c)令其中，R_g(t_m)表示第m个周期第g个可视等效散射点距离雷达的瞬时径向距离，C表示光速，τ_g(t_m)表示第m个周期第g个可视等效散射点的时延，代入去载频后的回波信号的表达式之后，去载频后的回波信号的表达式化简为： 

2d)对去载频后的回波信号进行脉压处理，得到回波信号进行脉压之后的表达式： 

其中，σ_g表示第g个可视等效散射点的后向散射系数，f_c表示载频，τ_g(t_m)表示第m个周期第g个可视等效散射点的时延，B表示信号带宽，表示快时间； 

2e)取回波信号进行脉压之后的表达式中每一慢时间内的最大值，得到回波序列S_r(t_m)，回波序列S_r(t_m)的表达式为： 

其中，σ_g表示第g个可视等效散射点的后向散射系数，R_g(t_m)表示第m个周期第g个可视等效散射点距离雷达的瞬时径向距离，τ_g(t_m)表示第m个周期第g个可视等效散射点的时延，C表示光速，f_c表示载频。 

4.根据权利要求1所述的基于瞬时调频率估计的进动目标的微多普勒提取方法，其特征在于，步骤5具体包括以下子步骤； 

5a)将每一观测时间的回波段近似为T个线性调频信号分量的叠加，T个线性调频信号分量的叠加之后，叠加后每一观测时间的回波段信号模型表示为： 

其中，δ_q表示第q个线性调频信号分量的复幅度，表示第q个线性调频信号分量的初始频率，μ_q表示第q个线性调频信号分量的调频率，T表示线性调频信号分量个数，设定T等于可视等效散射点的数目G，t_m＝m T_r表示慢时间，m为整数，T_r表示脉冲重复时间，S(t_m)表示叠加之后每一观测时间的回波段； 

5b)先设定线性调频信号分量个数T＝1，即单分量情况，设单分量线性调频信号为： 

其中，δ表示单分量线性调频信号的复幅度，f₀表示单分量线性调频信号的初始频率，μ表示单分量线性调频信号的调频率，t_m＝m T_r表示慢时间，m为整数，T_r表示脉冲重复时间； 

用解线调频的方法得到单分量线性调频信号的复幅度δ的估计值初始频率f₀的估计值调频率μ的估计值

5c)对于每一观测时间的回波段的T个线性调频信号分量，用松弛解线调频算法估计每个线性调频信号分量的复幅度、初始频率和调频率，得到每个线性调频信号分量的复幅度的估计值初始频率的估计值和调频率的估计值。

5.根据权利要求2所述的基于瞬时调频率估计的进动目标的微多普勒提取方法，其特 征在于，步骤6具体包括以下子步骤： 

6a)将各等效散射点瞬时微多普勒频率的导数公式(3)简化后得到等效散射点瞬时微多普勒频率导数模型表达式a·cos(2π·bx)＝y； 

6b)用随机抽样一致性算法处理每一观测时间的回波段的调频率的估计值，把估计出来的调频率区分为局内点和局外点，然后由局内点确定一条调频率曲线； 

6c)根据可视等效散射点的个数G，确定G条调频率曲线； 

若G＝1，根据6b)得到一条调频率曲线； 

若G＝2，则可视等效散射点的个数为两个，由6b)中得到的局外点估计出模型a·cos(2π·bx)＝y中的另一组参数a和b，得到另一条调频率曲线； 

若G＝3，则可视等效散射点的个数为三个，由6b)中的局内点确定一条调频率曲线；对6b)中区分的局外点再用一次随机抽样一致性算法，把6b)中区分的局外点再区分为两类不同的点，由这两类不同的点估计出模型a·cos(2π·bx)＝y中的两组不同的参数a和b，进而确定另外两条不同的调频率曲线。 

6.根据权利要求5所述的基于瞬时调频率估计的进动目标的微多普勒提取方法，其特征在于，子步骤6b)具体包括以下子步骤： 

6b1)随机选取4LM*G个线性调频信号分量的调频率估计值中的两个点，根据这两个点以及建立的等效散射点瞬时微多普勒频率导数模型a·cos(2π·bx)＝y，估计模型中的未知参数a和b，选择的两个点被设定为局内点，没有选择的点为局外点； 

6b2)通过等效散射点瞬时微多普勒频率导数模型测试4LM*G个线性调频信号分量的调频率估计值中的其他点；也就是将除6b1)随机选择的两个点局内点之外的局外点代入等效散射点瞬时微多普勒频率导数模型，如果等式a·cos(2π·bx)＝y+Δy中的错误Δy小于检测阈值ζ，则将该点添加到局内点集合； 

6b3)用所有局内点去重新估计模型的参数a和b，再统计局内点的个数d，以及所有局内点与模型之间的错误之和E； 

6b4)重复6b1)-6b3)，得到局内点的个数记为d'，局内点与模型之间的错误之和E'； 

6b5)如果迭代之后局内点的个数d'大于迭代之前局内点的个数d，并且迭代之后的错误之和E'小于迭代之前的错误之和E，则更新局内点的个数d和局内点与模型之间的错 误E，并且更新模型的参数a和b以及局内点的集合；否则，转至6b6)； 

6b6)以预定的次数重复执行6b4)-6b5)； 

6b7)运行预定的次数后，输出该模型的参数a和b以及局内点的集合，由参数a和b确定一条调频率曲线。