CN103308896A - 一种适于非引擎机动目标的高精度跟踪方法 - Google Patents

Abstract

本发明涉及一种适用于非引擎机动目标的高精度跟踪方法，其首先根据目标实际运动中加速度满足的非零均值时间相关随机过程特性，建立含有系统自适应参数的目标运动模型；其次，根据建立的目标运动模型，对目标运动特性进行预测；然后，利用目标状态预测值和雷达观测数据值计算目标当前状态估计值，最后，根据目标加速度估计值修正系统自适应参数，利用修正的系统自适应参数更新目标的运动模型，利用更新的目标运动模型进行下一次的预测与估计；本方法对具有非引擎机动随机特征的目标进行在线估计，提高了机动目标模型描述实际运动特性的准确性，提高了非引擎机动目标的跟踪精度。

Description

一种适于非引擎机动目标的高精度跟踪方法

技术领域

本发明涉及目标跟踪技术领域，尤其涉及雷达对非引擎机动目标的高精度跟踪方法。

背景技术

目标跟踪技术包括目标建模及利用测量数据进行估计两部分。传统的机动目标模型，如常速度模型、常加速度模型、Singer模型、Jerk模型及“当前统计模型”都是根据先验知识对目标运动特性进行建模。以下专利都是采用上述模型中的一种：基于粒子滤波和航迹管理的微弱目标检测前跟踪方法（201210106282）、一种高机动目标跟踪方法（201210138397）、一种基于多尺度维分解的目标跟踪方法（201110361072.5）、一种采用前向散射雷达的三维目标跟踪方法（201110247074）、卡尔曼滤波与经验模态分解有机结合的机动目标跟踪方法（201010235062）、一种基于分布式处理的无线传感器网络目标跟踪方法（200910037862），其中一种高机动目标跟踪方法（201210138397）使用的是改进的Jerk模型，需要已知目标的加速度机动频率及方差。由于驱动引擎本身的动力学约束，引擎驱动目标加速度机动频率及其随机特性是可测的，因此，尽管引擎目标的速度较快但加速度的机动特性可以较为准确的描述。但人、动物等非引擎驱动目标，加速度机动频率及其随机特性十分复杂而且种类繁多，目前尚没有每种运动方式的加速度机动频率及其随机特性的准确描述数据，这导致了这些方法在跟踪非引擎机动目标无法取得良好的跟踪效果。

基于嵌入式平台的自适应滤波目标跟踪定位方法（201110404800）公开了一种基于嵌入式平台的自适应滤波目标跟踪定位方法，其中涉及到根据三维空间运动的三阶拉格朗日级数展开地理位置观测数据进行建模，将运动模型建模模块输出的误差值进行残差数据处理，利用噪声空间估计关键参数，对增益矩阵进行参数优化，但其存在观测数据具有不确定性的缺点。

综上所述，目前关于目标的建模与跟踪的方法没有根据非引擎目标运动特征实时建模、在跟踪该类目标时无法取得良好的跟踪效果。

发明内容

本发明所要解决的技术问题是针对现有技术的不足，提供一种适用于非引擎机动目标的高精度跟踪方法，该方法能够利用测量数据对目标随机加速度运动特性进行实时建模并进行跟踪。

本发明解决上述技术问题的技术方案如下：

实现本发明的目的的基本思路是：首先，根据目标实际运动中加速度满足的非零均值时间相关随机过程特性，建立含有系统自适应参数的目标运动模型；其次，根据建立的目标运动模型，对目标运动特性进行预测；然后，利用目标状态预测值和雷达数据观测值对目标状态进行更新，最后，根据目标加速度估计值修正系统自适应参数，利用修正的系统自适应参数更新目标的运动模型，利用更新的目标运动模型进行下一次的预测与估计。

具体实现步骤如下：

步骤1：目标运动状态和系统自适应参数初始化；

步骤2：建立具有系统自适应参数的运动模型；

步骤3：根据建立的具有系统自适应参数的运动模型对目标状态进行预测，得到目标的状态预测值；

步骤4:根据目标的状态预测值和雷达数据观测值对目标状态进行更新，得到目标状态估计值；

步骤5：根据目标状态估计值计算目标加速度均值及目标加速度估计值；

步骤6：根据目标加速度估计值对系统自适应参数进行修正；

步骤7：利用步骤5中求得的目标加速度均值和步骤6中求得的修正后的系统自适应参数更新步骤2中的目标运动模型；

步骤8：重复步骤2至7，直至所有测量数据全部执行完毕，则结束。

在上述技术方案的基础上，本发明还可以做如下改进。

进一步，所述步骤1中目标运动状态和系统自适应参数初始化包括：

1.1设置状态初值

1.2系统自适应参数初值α=α₀和

其中α为机动频率，为加速度方差，

1.3自相关函数初值r₀(0)和r₀(1)的初值取为r₀(0)=0,r₀(1)=0，

1.4系统加速度分量初值

进一步，所述步骤2中建立具有系统自适应参数的运动模型的具体步骤包括：

2.1利用下式描述目标的运动特征

$x(k+1)=\mathrm{\Φ}(k+1,k)x\left(k\right)+U\left(k\right)\stackrel{\‾}{a}\left(k\right)+w\left(k\right)$

其中

为3维状态列向量，

分别是位移、速度和加速度，x(k+1)为k+1时刻目标的状态向量，k为采样时刻；Φ(k+1,k)为状态转移矩阵；x(k)为k时刻目标的状态向量；U(k)为控制矩阵；

为0时刻开始至k时刻目标的加速度均值；w(k)为过程噪声，其均值为0，方差为Q(k)；所述Φ(k+1,k)、U(k)及Q(k)中含有目标机动频率α和目标机动加速度方差

随着系统自适应参数的变化而变化。

2.2利用下式建立目标测量方程

y(k)=H(k)x(k)+v(k)

其中k为采样时刻，y(k)为目标在k时刻的雷达观测值，H(k)为测量矩阵，x(k)为k时刻目标的状态向量；v(k)为高斯测量白噪声，其方差为R，且与过程噪声w(k)相互独立。

进一步，所述步骤3中根据建立的具有系统自适应参数的运动模型对目标状态进行预测，得到目标的状态预测值的具体步骤如下：

3.1根据建立的具有系统自适应参数的运动模型和初始值完成目标状态的一步预测，预测方程式如下：

$\hat{x}\left(k\right|k-1)=\mathrm{\Φ}(k,k-1)\hat{x}(k-1|k-1)+U(k-1)\stackrel{\‾}{a}(k-1)$

其中

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

表示目标k-1时刻目标的状态估计值；U(k-1)为控制矩阵；

为从0时刻开始至k-1的加速度均值;

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

P(k|k-1)=Φ(k,k-1)P(k-1|k-1)Φ^T(k,k-1)+Q(k-1)

其中P(k|k-1)表示k-1时刻预测目标在k时刻的状态协方差，k为采样时刻，|表示条件操作符；P(k-1|k-1)表示k-1时刻目标的状态协方差的估计值；Φ(k,k-1)为状态转移矩阵；Q(k-1)为过程噪声协方差。

进一步，所述步骤4中根据目标的状态预测值和雷达数据观测值对目标状态进行更新，得到目标状态估计值的具体步骤如下：

4.1根据目标状态协方差预测值、测量矩阵及测量噪声方差按照下式计算滤波器增益

K(k)=P(k|k-1)H^T(k)[H(k)P(k|k-1)H^T(k)+R]^T

其中，K(k)为滤波器增益，k为采样时刻，P(k|k-1)表示k-1时刻预测目标在k时刻的状态协方差，H(k)为k时刻的测量矩阵，R为高斯测量白噪声的方差，H^T(k)为k时刻的测量矩阵的转置；

4.2利用目标状态预测值和雷达数据观测值计算目标当前状态估计值 $\hat{x}\left(k\right|k)=\hat{x}\left(k\right|k-1)+K\left(k\right)[y\left(k\right)-H\left(k\right)\hat{x}\left(k\right|k-1)]$

其中，

表示k时刻目标的状态估计值，

表示k-1时刻时预测目标在k时刻的状态，k为采样时刻，K(k)为k时刻滤波器增益，y(k)为雷达接收数据在k时刻的目标观测值，H(k)为k时刻的测量矩阵；

4.3按照下式计算目标状态协方差估计值，

P(k|k)=[I-K(k)H(k)]P(k|k-1)

其中I是3维单位矩阵，P(k|k)表示k时刻的目标状态协方差的估计值，k为采样时刻，K(k)为k时刻滤波器增益，H(k)为k时刻的测量矩阵，P(k|k-1)表示k-1时刻预测目标在k时刻的状态协方差。

进一步，所述步骤5中根据目标状态估计值计算目标加速度均值及目标加速度估计值的具体步骤如下：

步骤5.1：根据目标状态估计值计算目标加速度均值

$\stackrel{\‾}{a}\left(k\right)=\frac{1}{k}\underset{k=0}{\overset{k-1}{\mathrm{\Σ}}}\widehat{\stackrel{\·\·}{x}}\left(k\right|k)$

其中

为0至k时刻加速度均值，

为k时刻目标的状态估计值

的第三行值，k为采样时刻；

步骤5.2：按照下式获取系统k-1时刻和k时刻的加速度估计值

$\hat{a}\left(k\right),$

$\hat{a}(k-1)=\widehat{\stackrel{\·\·}{x}}(k-1|k-1)$

$\hat{a}\left(k\right)=\widehat{\stackrel{\·\·}{x}}\left(k\right|k)$

其中

为k-1时刻状态估计

的第三行值，

为k时刻状态估计的第三行值。

进一步，所述步骤6中根据加速度估计值对系统自适应参数进行修正的具体步骤如下：

根据采样时刻k值的大小，选择修正系统自适应参数α和的方法，若k小于等于4进入步骤6.1，若k大于4进入步骤6.2，

6.1当采样时刻k小于等于4时，按下式计算系统自适应参数α和

α=α₀,其中α₀为系统自适应参数α的初值，

如果 $\hat{a}\left(k\right)>0$ 则取 ${\mathrm{\σ}}_a^2=\frac{4-\mathrm{\π}}{\mathrm{\π}}{[a_M-\hat{a}\left(k\right)]}^2,$

如果 $\hat{a}\left(k\right)<0$ 则取 ${\mathrm{\&sigma;}}_a^2=\frac{4-\mathrm{\&pi;}}{\mathrm{\&pi;}}{[\hat{a}\left(k\right)-a_{-M}]}^2,$

如果则

取(0,10]之间的任意数，

其中，为k时刻目标加速度估计值，π为圆周率，取为3.14，a_M为正的常数，取为3，a_-M为与a_M绝对值相等的负常数，取为-3；

6.2当采样时刻k大于4时，按下式计算系统自适应参数α和

$r_k\left(1\right)=r_{k-1}\left(1\right)+\frac{1}{k}[30*\hat{a}\left(k\right)\hat{a}(k-1)-r_{k-1}\left(1\right)]$

$r_k\left(0\right)=r_{k-1}\left(0\right)+\frac{1}{k}[30*\hat{a}\left(k\right)\hat{a}\left(k\right)-r_{k-1}\left(0\right)]$

$\mathrm{\&alpha;}=-\frac{\ln r_k\left(1\right)-\ln r_k\left(0\right)}{T}$

T ${\mathrm{\&delta;}}_a^2=\frac{r_k\left(0\right)-\mathrm{\&alpha;}{r}_k\left(1\right)}{1-{\left(\frac{r_k\left(1\right)}{r_k\left(0\right)}\right)}^2}$

其中r_k(1)为k时刻目标加速度向前一步相关函数，r_k-1(1)为k-1时刻目标加速度向前一步相关函数，和

分别为k-1时刻和k时刻目标加速度估计值；r_k(0)为k时刻目标加速度自相关函数，r_k-1(0)为k-1时刻目标加速度自关函数；α为一个系统自适应参数，T为采样间隔；

为另一个系统自适应参数。

进一步，所述状态初始值x₀为3维全0列向量，维数为系统模型中状态向量的维数。

进一步，所述系统自适应参数的初始值α₀和

取任意正数。

本发明的有益效果是：

1.针对非引擎机动目标的加速度统计特性，建立含有系统自适应参数的目标运动模型，利用加速度估计值在线地实时修正系统自适应参数，保证了目标运动模型的实时准确性，实现高精度跟踪；

2.将具有自适应参数的目标运动模型与卡尔曼滤波器相结合，利用卡尔曼滤波器实时修正加速度估计值，利用加速度估计值在线地实时修正系统自适应参数，进而实时更新目标运动模型，保证了目标运动模型与实际目标运动的实时逼近，避免由于使用不准确观测数据进行建模可能导致滤波发散的缺点，克服了应用现有机动目标模型进行滤波时无法准确描述非引擎目标运动模型的缺点，适用于目标跟踪领域中的非引擎机动目标跟踪。

附图说明

图1为本发明所述一种适于非引擎机动目标的高精度跟踪方法流程图；

具体实施方式

以下结合附图对本发明的原理和特征进行描述，所举实例只用于解释本发明，并非用于限定本发明的范围。

如图1所示，一种适于非引擎机动目标的高精度跟踪方法，包括如下步骤：

步骤1：目标运动状态和系统自适应参数初始化

1.1设置状态初值

为3维全0列向量，维数为系统模型中状态向量的维数，

1.2系统自适应参数初值α=α₀和

取任意正数，本实施例中α₀取值

取值3，

1.3自相关函数初值r₀(0)和r₀(1)的初值取为

$r_0\left(1\right)={\stackrel{\&CenterDot;\&CenterDot;}{x}}_0=0,$

1.4系统加速度分量初值

本实施例中

取0；

步骤2：建立具有系统自适应参数的运动模型

2.1利用下式描述目标的运动特征

设目标加速度为非零均值的时间相关随机过程其中为加速度均值，a(t)为零均值指数相关有色噪声模型，其相关函数为

$R_a\left(\mathrm{\&tau;}\right)=E\left[a\left(t\right)a(t+\mathrm{\&tau;})\right]={\mathrm{\&delta;}}_a^2{e}^{-\mathrm{\&alpha;}\left|\mathrm{\&tau;}\right|}$

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

表示加速度方差；α为机动频率，反应目标的机动随机特性；

对有色噪声a(t)做白化处理，得到：

$\stackrel{\&CenterDot;}{a}\left(t\right)=-\mathrm{\&alpha;a}\left(t\right)+w\left(t\right)$

其中，w(t)为零均值白噪声，方差为

由 $\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)=\stackrel{\&OverBar;}{a}+a\left(t\right)$ 和 $\stackrel{\&CenterDot;}{a}\left(t\right)=-\mathrm{\&alpha;a}\left(t\right)+w\left(t\right)$ 得到目标运动的连续状态方程

$\stackrel{\&CenterDot;}{x}\left(t\right)=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 0& 0& -\mathrm{\&alpha;}\end{array}\right]x\left(t\right)+\left[\begin{array}{c}0\\ 0\\ \mathrm{\&alpha;}\end{array}\right]\stackrel{\&OverBar;}{a}\left(t\right)+\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]w\left(t\right)$

以周期T采样，离散化后系统目标运动满足以下方程

$x(k+1)=\mathrm{\&Phi;}(k+1,k)x\left(k\right)+U\left(k\right)\stackrel{\&OverBar;}{a}\left(k\right)+w\left(k\right)$

其中为3维状态列向量，

分别是位移、速度和加速度，x(k+1)为k+1时刻目标的状态向量，k为采样时刻；Φ(k+1,k)为状态转移矩阵；x(k)为k时刻目标的状态向量；U(k)为控制矩阵；

为0时刻开始至k时刻目标的加速度均值；w(k)为过程噪声，其均值为0，方差为Q(k)；所述Φ(k+1,k)、U(k)及Q(k)中含有目标机动频率α和目标机动加速度方差随着系统自适应参数的变化而变化；

其中状态转移矩阵Φ(k+1,k)的表达式如下：

$\mathrm{\&Phi;}(k+1,k)=\left[\begin{array}{ccc}1& T& \frac{\mathrm{\&alpha;T}-1+e^{-\mathrm{\&alpha;T}}}{{\mathrm{\&alpha;}}^2}\\ 0& 1& \frac{1-e^{-\mathrm{\&alpha;T}}}{\mathrm{\&alpha;}}\\ 0& 0& e^{-\mathrm{\&alpha;T}}\end{array}\right]$

控制矩阵U(k)的表达式如下：

$U\left(k\right)=\left[\begin{array}{c}\frac{1}{\mathrm{\&alpha;}}(-T+\frac{\mathrm{\&alpha;}{T}^2}{2}+\frac{1-e^{-\mathrm{\&alpha;T}}}{\mathrm{\&alpha;}})\\ T-\frac{1-e^{-\mathrm{\&alpha;T}}}{\mathrm{\&alpha;}}\\ 1-e^{-\mathrm{\&alpha;T}}\end{array}\right]$

过程噪声w(k)的方差Q(k)的表达式如下：

$Q\left(k\right)=E\left[w\left(k\right)w^T\left(k\right)\right]=2\mathrm{\&alpha;}{\mathrm{\&delta;}}_a^2\left[\begin{array}{ccc}{q}_{11}& q_{12}& q_{13}\\ q_{12}& q_{22}& q_{23}\\ q_{13}& q_{23}& q_{33}\end{array}\right]$

其中， $\left.\begin{array}{c}{q}_{11}=\frac{1}{2{\mathrm{\&alpha;}}^5}[1-e^{-2\mathrm{\&alpha;T}}+2\mathrm{\&alpha;T}+\frac{2{\mathrm{\&alpha;}}^3{T}^3}{3}-2{\mathrm{\&alpha;}}^2{T}^2-4\mathrm{\&alpha;T}{e}^{-\mathrm{\&alpha;T}}]\\ q_{12}=\frac{1}{2{\mathrm{\&alpha;}}^4}[e^{-2\mathrm{\&alpha;T}}+1-2e^{-\mathrm{\&alpha;T}}+2\mathrm{\&alpha;T}{e}^{-\mathrm{\&alpha;T}}-2\mathrm{\&alpha;T}+{\mathrm{\&alpha;}}^2{T}^2]\\ q_{13}=\frac{1}{2{\mathrm{\&alpha;}}^3}[1-e^{-2\mathrm{\&alpha;T}}-2\mathrm{\&alpha;}{\mathrm{Te}}^{-\mathrm{\&alpha;T}}]\\ q_{22}=\frac{1}{2{\mathrm{\&alpha;}}^3}[4e^{-\mathrm{\&alpha;T}}-3-e^{-2\mathrm{\&alpha;T}}+2\mathrm{\&alpha;T}]\\ q_{23}=\frac{1}{2{\mathrm{\&alpha;}}^2}[e^{-2\mathrm{\&alpha;T}}+1-2\mathrm{\&alpha;T}]\\ q_{33}=\frac{1}{2\mathrm{\&alpha;}}[1-e^{-2\mathrm{\&alpha;T}}]\end{array}\right\}$

E[·]表示求取均值；E[w(k)w^T(k)]表示求取零均值噪声w(k)的方差；

2.2利用下式建立目标测量方程

y(k)=H(k)x(k)+v(k)

其中k为采样时刻，y(k)为目标在k时刻的雷达观测值，H(k)为测量矩阵，x(k)为k时刻目标的状态向量；v(k)为高斯测量白噪声，其方差为R，且与过程噪声w(k)相互独立；

步骤3：根据建立的具有系统自适应参数的运动模型对目标状态进行预测，

3.1根据建立的具有系统自适应参数的运动模型和初始值完成目标状态的一步预测，预测方程式如下：

$\hat{x}\left(k\right|k-1)=\mathrm{\&Phi;}(k,k-1)\hat{x}(k-1|k-1)+U(k-1)\stackrel{\&OverBar;}{a}(k-1)$

其中

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

表示目标k-1时刻目标的状态估计值；U(k-1)为控制矩阵；

为从0时刻开始至k-1的加速度均值;

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

P(k|k-1)=Φ(k,k-1)P(k-1|k-1)Φ^T(k,k-1)+Q(k-1)

其中P(k|k-1)表示k-1时刻预测目标在k时刻的状态协方差，k为采样时刻，|表示条件操作符；P(k-1|k-1)表示k-1时刻目标的状态协方差的估计值；Φ(k,k-1)为状态转移矩阵；Q(k-1)为过程噪声协方差;

步骤4:根据目标的状态预测值、雷达数据观测值和目标的状态协方差预测值对目标状态进行更新

4.1根据目标状态协方差预测值、测量矩阵及测量噪声方差按照下式计算滤波器增益

K(k)=P(k|k-1)H^T(k)[H(k)P(k|k-1)H^T(k)+R]^T

其中，K(k)为滤波器增益，k为采样时刻，P(k|k-1)表示k-1时刻预测目标在k时刻的状态协方差，H(k)为k时刻的测量矩阵，R为高斯测量白噪声的方差，H^T(k)为k时刻的测量矩阵的转置；

4.2利用目标状态预测值和雷达数据观测值计算目标当前状态估计值

$\hat{x}\left(k\right|k)=\hat{x}\left(k\right|k-1)+K\left(k\right)[y\left(k\right)-H\left(k\right)\hat{x}\left(k\right|k-1)]$

其中，

表示k时刻目标的状态估计值，

表示k-1时刻时预测目标在k时刻的状态，k为采样时刻，K(k)为k时刻滤波器增益，y(k)为雷达接收数据在k时刻的目标观测值，H(k)为k时刻的测量矩阵；

4.3按照下式计算目标状态协方差的估计值

P(k|k)=[I-K(k)H(k)]P(k|k-1)

其中I是3维单位矩阵，P(k|k)表示k时刻的目标状态协方差的估计值，k为采样时刻，K(k)为k时刻滤波器增益，H(k)为k时刻的测量矩阵，P(k|k-1)表示k-1时刻预测目标在k时刻的状态协方差；

步骤5：根据目标状态估计值计算目标加速度均值及目标加速度估计值，

步骤5.1：利用下式计算目标加速度均值

$\stackrel{\&OverBar;}{a}\left(k\right)=\frac{1}{k}\underset{k=0}{\overset{k-1}{\mathrm{\&Sigma;}}}\widehat{\stackrel{\&CenterDot;\&CenterDot;}{x}}\left(k\right|k)$

其中为0至k时刻加速度均值，

为k时刻目标的状态估计值

的第三行值，k为采样时刻；

步骤5.2：按照下式获取系统k-1时刻和k时刻的加速度估计值

$\hat{a}(k-1)=\widehat{\stackrel{\&CenterDot;\&CenterDot;}{x}}(k-1|k-1)$

$\hat{a}\left(k\right)=\widehat{\stackrel{\&CenterDot;\&CenterDot;}{x}}\left(k\right|k)$

其中

为k-1时刻状态估计

的第三行值，

为k时刻状态估计的第三行值；

步骤6：根据目标加速度估计值对系统自适应参数进行修正，

根据采样时刻k值的大小，选择修正系统自适应参数α和

的方法，若k小于等于4进入步骤6.1，若k大于4进入步骤6.2，

6.1当采样时刻k小于等于4时，按下式计算系统自适应参数α和

α=α₀,其中α₀为系统自适应参数α的初值，

如果 $\hat{a}\left(k\right)>0$ 则取 ${\mathrm{\&sigma;}}_a^2=\frac{4-\mathrm{\&pi;}}{\mathrm{\&pi;}}{[a_M-\hat{a}\left(k\right)]}^2,$

如果 $\hat{a}\left(k\right)<0$ 则取 ${\mathrm{\&sigma;}}_a^2=\frac{4-\mathrm{\&pi;}}{\mathrm{\&pi;}}{[\hat{a}\left(k\right)-a_{-M}]}^2,$

如果

则

取(0,10]之间的任意数，

其中，

为k时刻目标加速度估计值，π为圆周率，取为3.14，a_M为正的常数，取为3，a_-M为与a_M绝对值相等的负常数，取为-3；

6.2当采样时刻k大于4时，按下式计算系统自适应参数α和

$r_k\left(1\right)=r_{k-1}\left(1\right)+\frac{1}{k}[30*\hat{a}\left(k\right)\hat{a}(k-1)-r_{k-1}\left(1\right)]$

$r_k\left(0\right)=r_{k-1}\left(0\right)+\frac{1}{k}[30*\hat{a}\left(k\right)\hat{a}\left(k\right)-r_{k-1}\left(0\right)]$

其中r_k(1)为k时刻目标加速度向前一步相关函数，r_k-1(1)为k-1时刻目标加速度向前一步相关函数，

和

分别为k-1时刻和k时刻目标加速度估计值；r_k(0)为k时刻目标加速度自相关函数，r_k-1(0)为k-1时刻目标加速度自关函数；

根据系统方程 $x(k+1)=\mathrm{\&Phi;}(k+1,k)x\left(k\right)+U\left(k\right)\stackrel{\&OverBar;}{a}\left(k\right)+w\left(k\right)$ 得到目标运动加速度方程满足如下一阶马尔科夫随机序列：

$\hat{a}(k+1)=\mathrm{\&beta;}\hat{a}\left(k\right)+w^a\left(k\right)$

其中

为k+1时刻的加速度，

为k时刻的加速度，β为离散后加速度随机序列的机动频率，w^a(k)为零均值白噪声离散序列，方差为

其中

为零均值白噪声w(t)的方差，β与α的关系为β=e^-αT；

一阶马尔科夫时间加速度序列满足以下参数关系：

$\mathrm{\&beta;}=\frac{r_k\left(1\right)}{r_k\left(0\right)}$ ${\mathrm{\&delta;}}_{\mathrm{aw}}^2=r_k\left(0\right)-\mathrm{\&alpha;}{r}_k\left(1\right)$

其中，r_k(1)为k时刻的加速度向前一步相关函数，r_k(0)为k时刻的加速度自相关函数，α与β分别为加速度及其离散化后加速度序列的机动频率，自适应参数可按照下式计算得到：

$\mathrm{\&alpha;}=-\frac{\ln r_k\left(1\right)-\ln r_k\left(0\right)}{T}$ ${\mathrm{\&delta;}}_a^2=\frac{r_k\left(0\right)-\mathrm{\&alpha;}{r}_k\left(1\right)}{1-{\left(\frac{r_k\left(1\right)}{r_k\left(0\right)}\right)}^2}$

其中，r_k(1)为k时刻的加速度向前一步相关函数，r_k-1(0)为k时刻的加速度自相关函数，ln为取以e为底的对数计算；α为一个系统自适应参数，T为采样间隔；

为另一个系统自适应参数；

步骤7：利用步骤5中求得的目标加速度均值和步骤6中求得的修正后的系统自适应参数更新步骤2中的目标运动模型；

步骤8：重复步骤2至7，直至所有测量数据全部执行完毕，则结束。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (9)

1.一种适于非引擎机动目标的高精度跟踪方法，其特征在于，包括如下步骤：

步骤1：目标运动状态和系统自适应参数初始化；

步骤2：建立具有系统自适应参数的运动模型；

步骤3：根据建立的具有系统自适应参数的运动模型对目标状态进行预测，得到目标的状态预测值；

步骤4:根据目标的状态预测值和雷达数据观测值对目标状态进行更新，得到目标状态估计值；

步骤5：根据目标状态估计值计算目标加速度均值及目标加速度估计值；

步骤6：根据目标加速度估计值对系统自适应参数进行修正；

步骤7：利用步骤5中求得的目标加速度均值和步骤6中求得的修正后的系统自适应参数更新步骤2中的目标运动模型；

步骤8：重复步骤2至7，直至所有测量数据全部执行完毕，则结束。

2.根据权利要求1所述一种适于非引擎机动目标的高精度跟踪方法，其特征在于，所述步骤1中目标运动状态和系统自适应参数初始化包括：

1.1设置状态初值

1.2系统自适应参数初值α=α₀和

其中α为机动频率，

为加速度方差，

1.3自相关函数初值r₀(0)和r₀(1)的初值取为r₀(0)=0,r₀(1)=0，

1.4系统加速度分量初值

3.根据权利要求2所述一种适于非引擎机动目标的高精度跟踪方法，其特征在于，所述步骤2中建立具有系统自适应参数的运动模型的具体步骤包括：

2.1利用下式描述目标的运动特征

$x(k+1)=\mathrm{\&Phi;}(k+1,k)x\left(k\right)+U\left(k\right)\stackrel{\&OverBar;}{a}\left(k\right)+w\left(k\right)$

其中

为3维状态列向量，

分别是位移、速度和加速度，x(k+1)为k-1时刻目标的状态向量，k为采样时刻；Φ(k+1,k)为状态转移矩阵；x(k)为k时刻目标的状态向量；U(k)为控制矩阵；

为0时刻开始至k时刻目标的加速度均值；w(k)为过程噪声，其均值为0，方差为Q(k)；所述Φ(k+1,k)、U(k)及Q(k)中含有目标机动频率α和目标机动加速度方差

随着系统自适应参数的变化而变化；

2.2利用下式建立目标测量方程

y(k)=H(k)x(k)+v(k)

其中k为采样时刻，y(k)为目标在k时刻的雷达观测值，H(k)为测量矩阵，x(k)为k时刻目标的状态向量；v(k)为高斯测量白噪声，其方差为R，且与过程噪声w(k)相互独立。

4.根据权利要求3所述一种适于非引擎机动目标的高精度跟踪方法，其特征在于，所述步骤3中根据建立的具有系统自适应参数的运动模型对目标状态进行预测，得到目标的状态预测值的具体步骤如下：

3.1根据建立的具有系统自适应参数的运动模型和初始值完成目标状态的一步预测，预测方程式如下：

$\hat{x}\left(k\right|k-1)=\mathrm{\&Phi;}(k,k-1)\hat{x}(k-1|k-1)+U(k-1)\stackrel{\&OverBar;}{a}(k-1)$

其中

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

表示目标k-1时刻目标的状态估计值；U(k-1)为控制矩阵；

为从0时刻开始至k-1的加速度均值;

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

P(k|k-1)=Φ(k,k-1)P(k-1|k-1)Φ^T(k,k-1)+Q(k-1)

其中P(k|k-1)表示k-1时刻预测目标在k时刻的状态协方差，k为采样时刻，|表示条件操作符；P(k-1|k-1)表示k-1时刻目标的状态协方差的估计值；Φ(k,k-1)为状态转移矩阵；Q(k-1)为过程噪声协方差。

5.根据权利要求4所述一种适于非引擎机动目标的高精度跟踪方法，其特征在于，所述步骤4中根据目标的状态预测值和雷达数据观测值对目标状态进行更新，得到目标状态估计值的具体步骤如下：

4.1根据目标状态协方差预测值、测量矩阵及测量噪声方差按照下式计算滤波器增益

K(k)=P(k|k-1)H^T(k)[H(k)P(k|k-1)H^T(k)+R]^T

其中，K(k)为滤波器增益，k为采样时刻，P(k|k-1)表示k-1时刻预测目标在k时刻的状态协方差，H(k)为k时刻的测量矩阵，R为高斯测量白噪声的方差，H^T(k)为k时刻的测量矩阵的转置；

4.2利用目标状态预测值和雷达数据观测值计算目标当前状态估计值 $\hat{x}\left(k\right|k)=\hat{x}\left(k\right|k-1)+K\left(k\right)[y\left(k\right)-H\left(k\right)\hat{x}\left(k\right|k-1)]$

其中，

表示k时刻目标的状态估计值，

表示k-1时刻时预测目标在k时刻的状态，k为采样时刻，K(k)为k时刻滤波器增益，y(k)为雷达接收数据在k时刻的目标观测值，H(k)为k时刻的测量矩阵；

4.3按照下式计算目标状态协方差估计值，

P(k|k)=[I-K(k)H(k)]P(k|k-1)

其中I是3维单位矩阵，P(k|k)表示k时刻的目标状态协方差的估计值，k为采样时刻，K(k)为k时刻滤波器增益，H(k)为k时刻的测量矩阵，P(k|k-1)表示k-1时刻预测目标在k时刻的状态协方差。

6.根据权利要求5所述一种适于非引擎机动目标的高精度跟踪方法，其特征在于，所述步骤5中根据目标状态估计值计算目标加速度均值及目标加速度估计值的具体步骤如下：

步骤5.1：根据目标状态估计值计算目标加速度均值

$\stackrel{\&OverBar;}{a}\left(k\right)=\frac{1}{k}\underset{k=0}{\overset{k-1}{\mathrm{\&Sigma;}}}\widehat{\stackrel{\&CenterDot;\&CenterDot;}{x}}\left(k\right|k)$

其中

为0至k时刻加速度均值，为k时刻目标的状态估计值的第三行值，k为采样时刻；

步骤5.2：按照下式获取系统k-1时刻和k时刻的加速度估计值

$\hat{a}\left(k\right),$

$\hat{a}(k-1)=\widehat{\stackrel{\&CenterDot;\&CenterDot;}{x}}(k-1|k-1)$

$\hat{a}\left(k\right)=\widehat{\stackrel{\&CenterDot;\&CenterDot;}{x}}\left(k\right|k)$

其中为k-1时刻状态估计

的第三行值，

为k时刻状态估计

的第三行值。

7.根据权利要求6所述一种适于非引擎机动目标的高精度跟踪方法，其特征在于，所述步骤6中根据目标加速度估计值对系统自适应参数进行修正的具体步骤如下：

根据采样时刻k值的大小，选择修正系统自适应参数α和

的方法，若k小于等于4进入步骤6.1，若k大于4进入步骤6.2，

6.1当采样时刻k小于等于4时，按下式计算系统自适应参数α和

α=α₀,其中α₀为系统自适应参数α的初值，

如果 $\hat{a}\left(k\right)>0$ 则取 ${\mathrm{\&sigma;}}_a^2=\frac{4-\mathrm{\&pi;}}{\mathrm{\&pi;}}{[a_M-\hat{a}\left(k\right)]}^2,$

如果 $\hat{a}\left(k\right)<0$ 则取 ${\mathrm{\&sigma;}}_a^2=\frac{4-\mathrm{\&pi;}}{\mathrm{\&pi;}}{[\hat{a}\left(k\right)-a_{-M}]}^2,$

如果

则

取(0,10]之间的任意数，

其中，

为k时刻目标加速度估计值，π为圆周率，取为3.14，a_M为正的常数，取为3，a_-M为与a_M绝对值相等的负常数，取为-3；

6.2当采样时刻k大于4时，按下式计算系统自适应参数α和

$r_k\left(1\right)=r_{k-1}\left(1\right)+\frac{1}{k}[30*\hat{a}\left(k\right)\hat{a}(k-1)-r_{k-1}\left(1\right)]$

$r_k\left(0\right)=r_{k-1}\left(0\right)+\frac{1}{k}[30*\hat{a}\left(k\right)\hat{a}\left(k\right)-r_{k-1}\left(0\right)]$

$\mathrm{\&alpha;}=-\frac{\ln r_k\left(1\right)-\ln r_k\left(0\right)}{T}$ ${\mathrm{\&delta;}}_a^2=\frac{r_k\left(0\right)-\mathrm{\&alpha;}{r}_k\left(1\right)}{1-{\left(\frac{r_k\left(1\right)}{r_k\left(0\right)}\right)}^2}$

其中r_k(1)为k时刻目标加速度向前一步相关函数，r_k-1(1)为k-1时刻目标加速度向前一步相关函数，

和

分别为k-1时刻和k时刻目标加速度估计值；r_k(0)为k时刻目标加速度自相关函数，r_k-1(0)为k-1时刻目标加速度自关函数；α为一个系统自适应参数，T为采样间隔；

为另一个系统自适应参数。

8.根据权利要求2所述一种适于非引擎机动目标的高精度跟踪方法，其特征在于，所述状态初始值x₀为3维全0列向量，维数为系统模型中状态向量的维数。

9.根据权利要求2所述一种适于非引擎机动目标的高精度跟踪方法，其特征在于，所述系统自适应参数的初始值α₀和

取任意正数。

Images