CN102721951A - 一种高机动目标跟踪方法 - Google Patents

Abstract

本发明公开了一种高机动目标跟踪方法，主要解决现有技术中目标高机动带来的模型不匹配和跟踪精度低的问题。其实现过程是：(1)运动状态和模型参数的初始化；(2)建立改进jerk模型；(3)目标状态预测；(4)根据雷达的接收数据和状态预测值，求观测值的预测误差；(5)确定渐消因子；(6)确定滤波器增益；(7)对目标的状态进行更新，实现对目标的跟踪。本发明在保持一般机动目标跟踪精度的同时，提高了对目标状态突变时的跟踪性能，可用于高机动目标的跟踪。

Description

一种高机动目标跟踪方法

技术领域

本发明属于通信技术领域，更进一步涉及雷达技术领域中一种高机动目标跟踪方法。本发明可用于实现高机动目标的跟踪，提高对目标突发机动时的跟踪性能。

背景技术

随着现代航空航天技术的飞速发展，各种飞行器的机动性能大幅提高，机动形式不再局限于加速度阶跃、蛇行机动等形式，也可能会出现加速度变化率的阶跃机动等更复杂的机动。在此背景下，如何提高高机动目标的跟踪性能成为一个越来越重要的问题。

目标跟踪的问题实际上就是目标状态的跟踪滤波问题，即根据雷达已获得的目标量测数据对目标状态进行精确的估计。建立目标的运动模型和自适应滤波是机动目标跟踪的两个关键部分，高机动目标的精确跟踪在过去和现在都是一个难题，最根本的原因在于目标的高机动使得目标的运动模型与目标实际动力学模型不匹配，导致跟踪滤波器发散，跟踪性能严重下降。因此，改善目标模型与目标实际运动模型的匹配度和跟踪滤波器的性能是提高高机动目标的跟踪性能的主要解决办法。

乔向东和王宝树在论文“高度机动目标的‘当前’统计jerk模型”(系统工程与电子技术，2002年第24卷第6期)中提出了一种新的“当前”统计jerk模型及其跟踪滤波方法。该方法首先分析了jerk模型的局限性，提出了“当前”统计jerk模型，最后基于“当前”统计jerk模型应用卡尔曼滤波，实现对目标的跟踪。该方法存在的不足是，仅适用于实现对一般机动目标的跟踪，由于卡尔曼跟踪滤波对高机动目标存在自身鲁棒性和跟踪精度低的问题，对于高机动目标，此方法受到限制。

宋强、何友和杨俭在论文“基于强跟踪滤波器的jerk模型目标跟踪算法”(海军航空工程学院学报，2007年第22卷第3期)中提出了一种强跟踪滤波算法。该方法首先分析了jerk模型自适应跟踪算法的不足，提出了强跟踪滤波的思想。该方法存在的不足也是仅适用于实现对一般机动目标的跟踪，对于高机动目标，由于jerk模型中关于目标加速度变化率为零均值的假设是不符合实际的，这一缺点将使得目标的运动模型与目标实际动力学模型不匹配，导致跟踪滤波器发散，跟踪性能下降。

发明内容

本发明的目的在于克服上述现有技术的不足，提出了一种基于改进jerk模型的强跟踪滤波器的高机动目标跟踪方法。

实现本发明的目的的基本思路是：首先，基于目标实际运动中的加速度变化率是非零均值的时间相关随机过程，建立改进jerk模型；其次，根据建立的运动模型，对目标状态进行预测；再次，根据雷达接收数据的观测值，确定观测值的预测误差；然后，由观测值的预测误差确定渐消因子，并由渐消因子自适应调整跟踪滤波器增益；最后，对目标状态进行更新，完成对机动目标的跟踪。

实现本发明目的的具体步骤如下：

(1)运动状态和模型参数的初始化

1a)利用雷达观测的前三个观测值确定初始状态和初始协方差；

1b)渐消因子的初始值设为1；

1c)确定机动频率和目标加速度变化率的方差；

(2)建立改进jerk模型

2a)用下式对目标的运动特性建立运动状态方程：

$X(k+1)=F\left(k\right)X\left(k\right)+G\left(k\right)\stackrel{\‾}{j}+W\left(k\right)$

其中，X(k+1)为k+1时刻状态向量，k为采样时刻；F(k)为状态转移矩阵；X(k)为k时刻状态向量， $X\left(k\right)={\left[\begin{array}{cccc}x\left(k\right)& \stackrel{\·}{x\left(k\right)}& \stackrel{\·\·}{x\left(k\right)}& \stackrel{\·\·\·}{x\left(k\right)}\end{array}\right]}^T,$ []^T为向量转置操作，x(k)表示目标k时刻位置；

为位置的一阶求导，表示目标k时刻速度；

为位置的二阶求导，表示目标k时刻加速度；

为位置的三阶求导，表示目标k时刻加速度变化率；G(k)为控制矩阵；

为目标加速度变化率的均值；W(k)为过程噪声，其均值和协方差分别为0和Q(k)的高斯白噪声系列；

2b)用下式建立目标的观测方程：

Z(k+1)＝H(k+1)X(k+1)+V(k+1)

其中，Z(k+1)为目标在k+1时刻的雷达观测值，k为采样时刻；H(k+1)为观测矩阵；X(k+1)为k+1时刻状态向量；V(k+1)为观测噪声，其均值和协方差分别为0和R(k)的高斯白噪声系列，且与过程噪声W(k)相互独立；

(3)目标状态预测

3a)由建立的模型和初始值完成目标状态的一步预测：

$\hat{X}(k+1|k)=F\left(k\right)\hat{X}\left(k\right|k)+G\left(k\right)\stackrel{\‾}{j}$

其中，

表示k时刻预测目标在k+1时刻的状态，k为采样时刻，|表示条件操作符；F(k)为状态转移矩阵；

表示目标k时刻状态更新值；G(k)为控制矩阵；为目标加速度变化率的均值；

3b)按照下式完成目标状态协方差的一步预测：

P(k+1|k)＝λ(k+1)F(k)P(k|k)F^T(k)+Q(k)

其中，P(k+1|k)表示k时刻预测目标在k+1时刻的状态协方差，k为采样时刻；λ(k+1)为k+1时刻的渐消因子；F(k)为状态转移矩阵；P(k|k)表示目标k时刻状态协方差更新值；F^T(k)为状态转移矩阵的转置；Q(k)为过程噪声的协方差；

(4)根据雷达的接收数据和状态预测值，求观测值的预测误差

4a)由状态预测值按照下式确定目标在k时刻的观测值的预测：

$\hat{Z}(k+1|k)=H(k+1)\hat{X}(k+1|k)$

其中，

表示k时刻预测目标在k+1时刻的观测值，k为采样时刻；H(k+1)为观测矩阵；

表示k时刻预测目标在k+1时刻的状态；

4b)由雷达接收数据按照下式确定目标在k时刻的观测值的预测误差：

$\tilde{Z}(k+1)=Z(k+1)-\hat{Z}(k+1|k)$

其中，

表示目标观测值的预测误差，k为采样时刻；Z(k+1)为雷达接收数据在k+1时刻的目标观测值；

表示k时刻预测目标在k+1时刻的观测值；

(5)确定渐消因子

用预测误差确定渐消因子，其计算公式如下：

$\mathrm{\λ}(k+1)=\left\{\begin{array}{cc}\mathrm{\η}(k+1)& \mathrm{\η}(k+1)>1\\ 1& \mathrm{\η}(k+1)\≤1\end{array}\right.$

其中，λ(k+1)为渐消因子，k为采样时刻；η(k+1)为可由预测误差协方差矩阵求得的渐消因子实时变化值；

(6)确定滤波器增益

6a)按照下式确定预测误差的协方差矩阵：

S(k+1)＝H(k+1)P(k+1|k)H^T(k+1)+R(k+1)

其中，S(k+1)为预测误差的协方差矩阵，k为采样时刻；H(k+1)为观测矩阵；P(k+1|k)表示k时刻预测目标在k+1时刻的状态协方差；H^T(k+1)为观测矩阵的转置；R(k+1)为观测噪声的协方差矩阵；

6b)由预测误差的协方差矩阵按照下式确定滤波器增益：

K(k+1)＝P(k+1|k)H^T(k+1)S^-1(k+1)

其中，K(k+1)为滤波器增益，k为采样时刻；P(k+1|k)表示k时刻预测目标在k+1时刻的状态协方差；H^T(k+1)为观测矩阵的转置；S^-1(k+1)为预测误差的协方差矩阵的逆；

(7)对目标的状态进行更新

7a)按照下式完成目标状态的更新：

$\hat{X}(k+1|k+1)=\hat{X}(k+1|k)+K(k+1)\tilde{Z}(k+1)$

其中，

表示目标k+1时刻状态更新值，k为采样时刻；

表示k时刻预测目标在k+1时刻的状态；K(k+1)为滤波器增益；

表示目标观测值的预测误差；

7b)按照下式完成目标状态协方差的更新：

P(k+1|k+1)＝P(k+1|k)-K(k+1)S(k+1)K^T(k+1)

其中，P(k+1|k+1)表示目标k+1时刻状态协方差更新值，k为采样时刻；P(k+1|k)表示k时刻预测目标在k+1时刻的状态协方差；K(k+1)为滤波器增益；S(k+1)为预测误差的协方差矩阵；K^T(k+1)为滤波器增益矩阵的转置；

(8)判断状态更新是否完成，若更新没有完成，则执行步骤(3)；若更新完成，则结束。

本发明与现有技术相比具有以下优点：

第一，本发明充分结合目标实际运动中的加速度变化率是非零均值的时间相关随机过程，对jerk模型进行改进，克服了现有技术中jerk模型与高机动目标实际运动模型不匹配的缺点，使得本发明能实现高机动目标的跟踪。

第二，本发明在跟踪滤波器中引入渐消因子，使滤波器能够自适应调节增益，克服了现有技术中卡尔曼跟踪滤波对高机动目标存在自身鲁棒性和跟踪精度低的缺点，使得本发明可以提高高机动目标的跟踪性能。

第三，本发明中的渐消因子能够随着观测值的预测误差自适应调整，使跟踪滤波器能够保持对一般机动目标跟踪精度的同时，还能够使滤波器在目标发生突发机动时，自适应的调节增益，增大目标状态估计的补偿值，提高高机动目标的跟踪性能，因此，本发明的使用范围更广。

附图说明

图1为本发明的流程图；

图2为本发明的仿真图。

具体实施方式

下面结合附图1对本发明的步骤作进一步的详细说明。

步骤1，运动状态和模型参数的初始化

雷达观测的前三个时刻的量测值为Z(1)、Z(2)和Z(3)，将雷达观测值Z(1)、Z(2)和Z(3)，利用差分法通过下式确定初始状态

$\hat{X}\left(3\right|3)=\left[\begin{array}{c}Z\left(3\right)\\ [Z\left(3\right)-Z\left(2\right)]/T\\ ((Z\left(3\right)-Z\left(2\right))/T-(Z\left(2\right)-Z\left(1\right))/T)/T\end{array}\right]$

其中，

为初始状态；T为采样间隔。

初始协方差P(3|3)可由下式确定：

$P\left(3\right|3)=E\left[\hat{X}\left(3\right|3){\hat{X}}^T\left(3\right|3)\right]$

其中，P(3|3)表示初始协方差；表示初始状态；

表示初始状态的转置。

本发明实施例中将渐消因子的初始值设为1，并将机动频率的初始值设为α＝1/20。

步骤2，建立改进jerk模型

基于目标的加速度变化率是非零均值的时间相关随机过程建立目标运动状态方程：

$\stackrel{\·\·\·}{x\left(t\right)}=\stackrel{\‾}{j}+j\left(t\right)$

其中，

为目标加速度变化率，

为目标加速度变化率的均值，j(t)为零均值的指数相关的有色噪声模型，其相关函数为：

$R_j\left(\mathrm{\τ}\right)=E\left[j\left(t\right)j(t+\mathrm{\τ})\right]={\mathrm{\σ}}_j^2{e}^{-\mathrm{\α}\left|\mathrm{\τ}\right|}$

其中，R_j(τ)表示相关函数，

表示目标加速度变化率方差，α表示jerk机动频率，反应了目标的机动特性；对有色噪声j(t)做白化处理，得到：

$\stackrel{\·}{j\left(t\right)}=-\mathrm{\αj}\left(t\right)+w\left(t\right)$

其中，w(t)为均值为0、方差为

的高斯白噪声。由式和式

得到目标运动的连续状态方程：

$\stackrel{\·}{j\left(t\right)}=-\mathrm{\α}\stackrel{\·\·\·}{x\left(t\right)}+\mathrm{\α}\stackrel{\‾}{j}+w\left(t\right)$

以周期T采样，离散化后，得到目标的运动状态方程为：

$X(k+1)=F\left(k\right)X\left(k\right)+G\left(k\right)\stackrel{\‾}{j}+W\left(k\right)$

其中，X(k+1)为k+1时刻状态向量，k为采样时刻；F(k)为状态转移矩阵，表达式如下：

$F\left(k\right)=\left[\begin{array}{cccc}1& T& T^2/2& p1\\ 0& 1& T& q1\\ 0& 0& 1& r1\\ 0& 0& 0& s1\end{array}\right]$

p1、q1、r1和s1为上述状态转移矩阵F(k)的元素，表达式如下：

p1＝(2-2αT+α²T²-2e^-αT)/(2α³)

q1＝(αT-1+e^-αT)/α²

r1＝(1-e^-αT)/α

s1＝e^-αT

X(k)为k时刻状态向量， $X\left(k\right)={\left[\begin{array}{cccc}x\left(k\right)& \stackrel{\·}{x\left(k\right)}& \stackrel{\·\·}{x\left(k\right)}& \stackrel{\·\·\·}{x\left(k\right)}\end{array}\right]}^T,$ []^T为向量转置操作，x(k)表示目标k时刻位置；

为位置的一阶求导，表示目标k时刻速度；

为位置的二阶求导，表示目标k时刻加速度；

为位置的三阶求导，表示目标k时刻加速度变化率；G(k)为控制矩阵，表达式如下：

$G\left(k\right)=\left[\begin{array}{c}\frac{1}{2{\mathrm{\α}}^2}(2T-{\mathrm{\αT}}^2+\frac{{\mathrm{\α}}^2{T}^3}{3}-\frac{2(1-e^{-\mathrm{\αT}})}{\mathrm{\α}})\\ \frac{1}{\mathrm{\α}}(-T+\frac{{\mathrm{\αT}}^2}{2}+\frac{1-e^{-\mathrm{\αT}}}{\mathrm{\α}})\\ T-\frac{1-e^{-\mathrm{\αT}}}{\mathrm{\α}}\\ 1-e^{-\mathrm{\αT}}\end{array}\right]$

W(k)为过程噪声，其均值和协方差分别为0和Q(k)的高斯白噪声系列，Q(k)为对称矩阵，表达式如下：

$Q\left(k\right)=2{\mathrm{\α\σ}}_j^2\left[\begin{array}{cccc}{q}_{11}& q_{12}& q_{13}& q_{14}\\ q_{21}& q_{22}& q_{23}& q_{24}\\ q_{31}& q_{32}& q_{33}& q_{34}\\ q_{41}& q_{42}& q_{43}& q_{44}\end{array}\right]$

q_ij(i＝1～4；j＝1～4)为上述矩阵Q(k)的元素，表达式如下：

q₁₁＝1/2α⁷[α⁵T⁵/10-α⁴T⁴/2+4α³T³/3-2α²T²-3+4e^-αT+2α²T²e^-αT-e^-2αT]

q₁₂＝1/2α⁶[1-2αT+2α²T²-α³T³+α⁴T⁴/4+e^-2αT+2αT-2e^-αT-α²T²e^-αT]

q₁₃＝1/2α⁵[2αT-α²T²-α³T³/3-3-2e^-2αT+4e^-αT+α²T²e^-αT]

q₁₄＝1/2α⁵[1+e^-2αT-2e^-αT-α²T²e^-αT]

q₂₂＝1/2α⁵[1-e^-2αT+2αT+2α³T³/3-2α²T²-4αTe^-αT]

q₂₃＝1/2α⁴[e^-2αT+1-2e^-αT+2αTe^-αT-2αT+α²T²]

q₂₄＝1/2α³[1-e^-2αT-2αTe^-αT]

q₃₃＝1/2α³[4e^-αT-3-e^-2αT+2αT]

q₃₄＝1/2α²[e^-2αT+1-2e^-αT]

q₄₄＝1/2α[1-e^-2αT]

建立目标的观测方程：

Z(k+1)＝H(k+1)X(k+1)+V(k+1)

其中，Z(k+1)为目标在k+1时刻的雷达观测值，k为采样时刻；H(k+1)为观测矩阵；X(k+1)为k+1时刻状态向量；V(k+1)为观测噪声，其均值和协方差分别为0和R(k)的高斯白噪声系列，且与过程噪声W(k)相互独立。

步骤3，目标状态的预测

由建立的模型和初始值按照下式完成目标状态的一步预测：

$\hat{X}(k+1|k)=F\left(k\right)\hat{X}\left(k\right|k)+G\left(k\right)\stackrel{\‾}{j}$

其中，

表示k时刻预测目标在k+1时刻的状态，k为采样时刻；F(k)为状态转移矩阵；

表示目标k时刻状态更新值；G(k)为控制矩阵；为目标加速度变化率的均值。

按照下式完成目标状态协方差的一步预测：

P(k+1|k)＝λ(k+1)F(k)P(k|k)F^T(k)+Q(k)

其中，P(k+1|k)表示k时刻预测目标在k+1时刻的状态协方差，k为采样时刻；λ(k+1)为k+1时刻的渐消因子；F(k)为状态转移矩阵；P(k|k)表示目标k时刻状态协方差更新值；F^T(k)为状态转移矩阵的转置；Q(k)为过程噪声的协方差。

在协方差的一步预测中引入渐消因子λ(k+1)，使得预测协方差可以随着渐消因子的变化而自适应的调整，提高算法对目标机动的适应性。

步骤4，根据雷达的接收数据和状态预测值，求观测值的预测误差

由状态预测值按照下式确定目标在k时刻的观测值的预测：

$\hat{Z}(k+1|k)=H(k+1)\hat{X}(k+1|k)$

其中，

表示k时刻预测目标在k+1时刻的观测值，k为采样时刻；H(k+1)为观测矩阵；表示k时刻预测目标在k+1时刻的状态；

由雷达接收数据按照下式确定目标在k时刻的观测值的预测误差：

$\tilde{Z}(k+1)=Z(k+1)-\hat{Z}(k+1|k)$

其中，

表示目标观测值的预测误差，k为采样时刻；Z(k+1)为雷达接收数据在k+1时刻的目标观测值；表示k时刻预测目标在k+1时刻的观测值；

观测值的预测误差的反应了建立的模型与目标实际运动模型之间的匹配度，观测值的预测误差越小，说明建立的模型越好。

步骤5，确定渐消因子

用预测误差确定渐消因子，其计算公式如下：

$\mathrm{\λ}(k+1)=\left\{\begin{array}{cc}\mathrm{\η}(k+1)& \mathrm{\η}(k+1)>1\\ 1& \mathrm{\η}(k+1)\≤1\end{array}\right.$

其中，λ(k+1)为渐消因子，k为采样时刻；η(k+1)为渐消因子实时变化值，可由下式确定：

$\mathrm{\η}(k+1)=\frac{\mathrm{tr}\left[N(k+1)\right]}{\mathrm{tr}\left[M(k+1)\right]}$

其中，

N(k+1)＝V₀(k+1)-βR(k+1)-H(k+1)Q(k)H^T(k+1)

M(k+1)＝H(k+1)F(k)P(k|k)F^T(k)H^T(k+1)

$V_0(k+1)=E\left[\tilde{Z}(k+1){\tilde{Z}}^T(k+1)\right]$

$=\left\{\begin{array}{cc}\tilde{Z}\left(1\right){\tilde{Z}}^T\left(1\right)& k=0\\ \frac{{\mathrm{\ρV}}_0\left(k\right)+\tilde{Z}(k+1){\tilde{Z}}^T(k+1)}{1+\mathrm{\ρ}}& k\≥1\end{array}\right.$

其中，

为初始观测值预测误差，0＜ρ≤1是遗忘因子，本发明实施例中取ρ＝0.95，β是弱化因子，一般取大于1，本发明实施例中取β＝5。

在跟踪一般机动目标时，观测值的预测误差较小，使得渐消因子λ(k+1)为1，此时算法退化为基于jerk模型的卡尔曼自适应滤波算法，保持了对一般机动目标的跟踪精度。在目标发生突发机动时，强跟踪滤波器根据观测值的预测误差的增大而增大渐消因子，自适应的调节增益，增大了目标状态估计的补偿值，从而提高目标状态突变时的跟踪性能。

步骤6，确定滤波器增益

按照下式确定预测误差的协方差矩阵：

S(k+1)＝H(k+1)P(k+1|k)H^T(k+1)+R(k+1)

其中，S(k+1)为预测误差的协方差矩阵，k为采样时刻；H(k+1)为观测矩阵；P(k+1|k)表示k时刻预测目标在k+1时刻的状态协方差；H^T(k+1)为观测矩阵的转置；R(k+1)为观测噪声的协方差矩阵。

由预测误差的协方差矩阵按照下式确定滤波器增益：

K(k+1)＝P(k+1|k)H^T(k+1)S^-1(k+1)

其中，K(k+1)为滤波器增益，k为采样时刻；P(k+1|k)表示k时刻预测目标在k+1时刻的状态协方差；H^T(k+1)为观测矩阵的转置；S^-1(k+1)为预测误差的协方差矩阵的逆。

步骤7，对目标状态进行更新

按照下式完成目标状态的更新：

$\hat{X}(k+1|k+1)=\hat{X}(k+1|k)+K(k+1)\tilde{Z}(k+1)$

其中，

表示目标k+1时刻状态更新值，k为采样时刻；表示k时刻预测目标在k+1时刻的状态；K(k+1)为滤波器增益；

表示目标观测值的预测误差。

按照下式完成目标状态协方差的更新：

P(k+1|k+1)＝P(k+1|k)-K(k+1)S(k+1)K^T(k+1)

其中，P(k+1|k+1)表示目标k+1时刻状态协方差更新值，k为采样时刻；P(k+1|k)表示k时刻预测目标在k+1时刻的状态协方差；K(k+1)为滤波器增益；S(k+1)为预测误差的协方差矩阵；K^T(k+1)为滤波器增益矩阵的转置。

通过对目标的运动状态进行更新，实现对目标的每一步跟踪。

步骤8，判断状态更新是否完成，若更新没有完成，则执行步骤(3)；若更新完成，则结束。

下面结合附图2对本发明的效果做进一步说明。

1.仿真条件

设目标的真实航迹的起始状态为[60000，-420，0，0，1000，0，0，0]，对目标进行120s的连续观测，目标发生机动时刻及加速度的变化方式如下：

在38s时，a_x＝0，a_y＝-10，j_x＝0，j_y＝0；在61s时，a_x＝-5，a_y＝20，j_x＝0，j_y＝0；在71s时，a_x＝-10，a_y＝10，j_x＝0，j_y＝0；在91s时，a_x＝0，a_y＝0，j_x＝0，j_y＝0在100s时，a_x＝0，a_y＝0，j_x＝0.5，j_y＝0。其中，a_x，a_y，j_x和j_y分别表示目标的X方向的加速度，Y方向的加速度，X方向的加速度变化率和Y方向的加速度变化率。

在笛卡尔坐标系下，采取500次Monte-Carlo仿真实验，设雷达采样间隔T＝1s，雷达的测距精度为15m，观测方位角的精度为0.17°。

误差的计算公式为：

$\mathrm{RMSE}=\frac{1}{M}\underset{j=1}{\overset{M}{\mathrm{\Σ}}}|X_i^j\left(k\right)-{\hat{X}}_i^j\left(k\right|k)|$

其中，M为Monte-Carlo仿真次数，

表示第j次仿真中第k时刻状态向量X(k)的第i个分量，

表示第j次仿真中第k时刻状态估计量

的第i个分量。

2.仿真内容

采用本发明方法和jerk模型自适应跟踪方法分别对目标的X方向位置、Y方向位置、X方向速度和Y方向速度进行跟踪估计，并比较跟踪效果。

3.仿真分析

图2(a)是本发明和jerk模型自适应跟踪方法对目标的X方向位置的跟踪估计效果比较图。图中实线表示采用本发明对目标的X方向位置进行跟踪估计结果曲线，虚线表示采用jerk模型自适应跟踪方法对目标的X方向位置进行跟踪估计结果曲线。从图中可以看出，整体趋势上本发明的稳定性和跟踪性能都要优于jerk模型的自适应跟踪方法，特别是在100s后，目标在X方向作变加速度运动，机动性更强时，此时本发明的跟踪性能明显优于jerk模型的自适应跟踪方法。

图2(b)是本发明和jerk模型自适应跟踪方法对目标的Y方向位置的跟踪估计效果比较图。图中实线表示采用本发明对目标的Y方向位置进行跟踪估计结果曲线，虚线表示采用jerk模型自适应跟踪方法对目标的Y方向位置进行跟踪估计结果曲线。从图中可以看出，整体趋势上本发明的跟踪误差小于jerk模型的自适应跟踪方法。

图2(c)是本发明和jerk模型自适应跟踪方法对目标的X方向速度的跟踪估计效果比较图。图中实线表示采用本发明对目标的X方向速度进行跟踪估计结果曲线，虚线表示采用jerk模型自适应跟踪方法对目标的X方向速度进行跟踪估计结果曲线。从图中可以看出，整体趋势上本发明的跟踪误差小于jerk模型的自适应跟踪方法，特别是在100s后，目标在X方向作变加速度运动，机动性更强时，此时本发明的跟踪性能明显优于jerk模型的自适应跟踪方法。

图2(d)是本发明和jerk模型自适应跟踪方法对目标的Y方向速度的跟踪估计效果比较图。图中实线表示采用本发明对目标的Y方向速度进行跟踪估计结果曲线，虚线表示采用jerk模型自适应跟踪方法对目标的Y方向速度进行跟踪估计结果曲线。从图中可以看出，整体趋势上本发明的跟踪误差小于jerk模型的自适应跟踪方法。

Claims (4)

1.一种高机动目标跟踪方法，包括如下步骤：

(1)运动状态和模型参数的初始化

1a)利用雷达观测的前三个观测值确定初始状态和初始协方差；

1b)渐消因子的初始值设为1；

1c)确定机动频率和目标加速度变化率的方差；

(2)建立改进jerk模型

2a)用下式对目标的运动特性建立运动状态方程：

$X(k+1)=F\left(k\right)X\left(k\right)+G\left(k\right)\stackrel{\‾}{j}+W\left(k\right)$

其中，X(k+1)为k+1时刻状态向量，k为采样时刻；F(k)为状态转移矩阵；X(k)为k时刻状态向量， $X\left(k\right)={\left[\begin{array}{cccc}x\left(k\right)& \stackrel{\·}{x\left(k\right)}& \stackrel{\·\·}{x\left(k\right)}& \stackrel{\·\·\·}{x\left(k\right)}\end{array}\right]}^T,$ []^T为向量转置操作，x(k)表示目标k时刻位置；

为位置的一阶求导，表示目标k时刻速度；

为位置的二阶求导，表示目标k时刻加速度；

为位置的三阶求导，表示目标k时刻加速度变化率；G(k)为控制矩阵；

为目标加速度变化率的均值；W(k)为过程噪声，其均值和协方差分别为0和Q(k)的高斯白噪声系列；

2b)用下式建立目标的观测方程：

Z(k+1)＝H(k+1)X(k+1)+V(k+1)

其中，Z(k+1)为目标在k+1时刻的雷达观测值，k为采样时刻；H(k+1)为观测矩阵；X(k+1)为k+1时刻状态向量；V(k+1)为观测噪声，其均值和协方差分别为0和R(k)的高斯白噪声系列，且与过程噪声W(k)相互独立；

(3)目标状态预测

3a)由建立的模型和初始值完成目标状态的一步预测：

$\hat{X}(k+1|k)=F\left(k\right)\hat{X}\left(k\right|k)+G\left(k\right)\stackrel{\‾}{j}$

其中，表示k时刻预测目标在k+1时刻的状态，k为采样时刻，|表示条件操作符；F(k)为状态转移矩阵；

表示目标k时刻状态更新值；G(k)为控制矩阵；为目标加速度变化率的均值；

3b)按照下式完成目标状态协方差的一步预测：

P(k+1|k)＝λ(k+1)F(k)P(k|k)F^T(k)+Q(k)

其中，P(k+1|k)表示k时刻预测目标在k+1时刻的状态协方差，k为采样时刻；λ(k+1)为k+1时刻的渐消因子；F(k)为状态转移矩阵；P(k|k)表示目标k时刻状态协方差更新值；F^T(k)为状态转移矩阵的转置；Q(k)为过程噪声的协方差；

(4)根据雷达的接收数据和状态预测值，求观测值的预测误差

4a)由状态预测值按照下式确定目标在k时刻的观测值的预测：

$\hat{Z}(k+1|k)=H(k+1)\hat{X}(k+1|k)$

其中，

表示k时刻预测目标在k+1时刻的观测值，k为采样时刻；H(k+1)为观测矩阵；

表示k时刻预测目标在k+1时刻的状态；

4b)由雷达接收数据按照下式确定目标在k时刻的观测值的预测误差：

$\tilde{Z}(k+1)=Z(k+1)-\hat{Z}(k+1|k)$

其中，

表示目标观测值的预测误差，k为采样时刻；Z(k+1)为雷达接收数据在k+1时刻的目标观测值；

表示k时刻预测目标在k+1时刻的观测值；

(5)确定渐消因子

用预测误差确定渐消因子，其计算公式如下：

$\mathrm{\λ}(k+1)=\left\{\begin{array}{cc}\mathrm{\η}(k+1)& \mathrm{\η}(k+1)>1\\ 1& \mathrm{\η}(k+1)\≤1\end{array}\right.$

其中，λ(k+1)为渐消因子，k为采样时刻；η(k+1)为可由预测误差协方差矩阵求得的渐消因子实时变化值；

(6)确定滤波器增益

6a)按照下式确定预测误差的协方差矩阵：

S(k+1)＝H(k+1)P(k+1|k)H^T(k+1)+R(k+1)

其中，S(k+1)为预测误差的协方差矩阵，k为采样时刻；H(k+1)为观测矩阵；P(k+1|k)表示k时刻预测目标在k+1时刻的状态协方差；H^T(k+1)为观测矩阵的转置；R(k+1)为观测噪声的协方差矩阵；

6b)由预测误差的协方差矩阵按照下式确定滤波器增益：

K(k+1)＝P(k+1|k)H^T(k+1)S^-1(k+1)

其中，K(k+1)为滤波器增益，k为采样时刻；P(k+1|k)表示k时刻预测目标在k+1时刻的状态协方差；H^T(k+1)为观测矩阵的转置；S^-1(k+1)为预测误差的协方差矩阵的逆；

(7)对目标的状态进行更新

7a)按照下式完成目标状态的更新：

$\hat{X}(k+1|k+1)=\hat{X}(k+1|k)+K(k+1)\tilde{Z}(k+1)$

其中，表示目标k+1时刻状态更新值，k为采样时刻；

表示k时刻预测目标在k+1时刻的状态；K(k+1)为滤波器增益；

表示目标观测值的预测误差；

7b)按照下式完成目标状态协方差的更新：

P(k+1|k+1)＝P(k+1|k)-K(k+1)S(k+1)K^T(k+1)

其中，P(k+1|k+1)表示目标k+1时刻状态协方差更新值，k为采样时刻；P(k+1|k)表示k时刻预测目标在k+1时刻的状态协方差；K(k+1)为滤波器增益；S(k+1)为预测误差的协方差矩阵；K^T(k+1)为滤波器增益矩阵的转置；

(8)判断状态更新是否完成，若更新没有完成，则执行步骤(3)；若更新完成，则结束。

2.根据权利要求1所述的一种高机动目标跟踪方法，其特征在于，步骤1a)中所述初始状态是利用雷达的前三个量测值运用差分法确定的。

3.根据权利要求1所述的一种高机动目标跟踪方法，其特征在于，步骤1a)中所述初始协方差的计算方法如下：

$P\left(3\right|3)=E\left[\hat{X}\left(3\right|3){\hat{X}}^T\left(3\right|3)\right]$

其中，P(3|3)表示初始协方差；表示初始状态；

表示初始状态的转置。

4.根据权利要求1所述的一种高机动目标跟踪方法，其特征在于，步骤1c)中所述机动频率的取值为0～1之间的任意值。

Images