CN101299271A - 一种机动目标状态方程的多项式预测模型及跟踪方法 - Google Patents

Abstract

本发明属于机动目标跟踪技术领域，具体为一种机动目标状态方程的多项预测模型及跟踪算法。本发明根据匀变速运动的多项式形式，将变加速运动作为分段匀变速运动处理，提出了一种新的机动目标动态模型——多项式预测模型。该模型能够完全避免对描述未知运动的多项式系数的依赖，用本发明方法建立的状态方程始终能准确描述运动的动态。因此本发明提出的模型能够适用于各种运动目标状态方程，本发明提出的多项式预测模型的最优滤波算法适用于任何可以用多项式描述的机动目标状态估计问题，既不需要关于目标运动参数的先验信息，又不需要通过参数辨识等手段调整建立状态模型的参数，从而完全避免了由于状态模型不准确而带来的估计性能恶化问题。

Description

一种机动目标状态方程的多项式预测模型及跟踪方法

技术领域

本发明属于机动目标跟踪技术领域，具体涉及一种机动目标状态方程建模及跟踪的新方法。

背景技术

目标跟踪的关键在于能否有效地从观测中提取出有关目标状态的有用信息，而一个准确的目标动态模型则非常有利于这种信息的提取。总的说来，一个准确的模型甚至比大量的数据更加有效，在观测数据十分有限的情况下，模型的作用就更加重要^[1]。

几乎所有的机动目标跟踪算法都是基于模型的，即假设目标的运动及对目标的观测都可由已知的数学模型以足够高的精度表示。其中，最常用的状态空间模型一般用下式表示：

x_n+1＝f_n(x_n，u_n)+w_n            (1)

y_n＝h_n(x_n)+v_n

其中，x_n，y_n，u_n分别为目标在t_n时刻的状态、观测及控制输入，u_n的作用效果为目标的加速度变化。w_n和v_n分别为状态噪声和观测噪声。f_n和h_n是时变/时不变函数。

目标跟踪问题的一个主要难点就在于目标运动的不确定性。对跟踪者来说，被跟踪目标的准确动态模型通常是未知的。即使式(1)所示的状态空间模型可以作为目标状态的通用模型，但对不同的目标，跟踪者仍然很难确定目标实际控制变量u、函数f的具体形式和参数、噪声w的统计特性等^[1]。当目标的动态模型不准确时，估计算法往往无法发挥其效能，严重的情况下会出现发散。文献[2-4]详细分析了不准确的状态方程对卡尔曼滤波的影响。因此，目标动态模型的建立成了机动目标跟踪问题的首要任务。

在过去的几十年里，人们提出了目标运动的各种数学模型。文献[1]中，将各种文献中报道应用于机动目标跟踪的目标动态模型进行了总结，分析了这些模型的基本思想和假设，不同模型间内在的联系，并给出了这些模型各自的优缺点。总体上，目前对真实动态特性未知的机动目标建模有两种途径：其一是利用满足某种特性的随机过程对实际的非随机控制输入u进行逼近；其二是利用某些具有代表性的运动模型，通过调整设计参数来描述典型的目标轨迹。

在现有的目标动态模型中，有些模型形式简单而精确度不高，但应用简便灵活，如白噪声加速度模型；而有些更加精确的模型则形式复杂，非线性度大，需要有效的非线性滤波算法，计算量也很大，另外由于这种模型依赖于目标的状态，因此在目标状态估计不准确的情况下，将会导致跟踪性能更加恶化。

发明内容

本发明的目的在于提供一种形式简单并可避免对目标真实状态过分依赖的机动目标状态方程模型及跟踪的新方法。

本发明提出机动目标状态方程的多项式预测模型，形式简单，能够完全避免对描述未知运动的多项式系数的依赖，无论在什么时候描述未知运动的多项式系数发生了变化，该模型始终能准确描述该运动。本发明提出的机动目标状态方程的模型及滤波方法，适用于任何可以用多项式描述的机动目标跟踪问题，既不需要有关目标运动参数的先验信息，又不需要通过参数辨识等手段调整建立状态模型的参数。

假设信号x(n)是由下述方程描述的一个L阶多项式：

$x\left(n\right)=\underset{l=0}{\overset{L}{\mathrm{\Σ}}}p\left(l\right)n^l,---\left(2\right)$

其中，p(l)为多项式的系数，l＝0，...，L；

如果用信号x(n)的前K个时刻的值x(n-K+1)，...，x(n)来预测信号式(2)的将来值x(n+N)，即：

$x(n+N)=\underset{k=0}{\overset{K-1}{\mathrm{\Σ}}}h\left(k\right)x(n-k),---\left(3\right)$

显然，式(3)是一以h(k)(k＝0，...，K-1)为系数的FIR滤波器。

由式(2)和式(3)可知：

$\underset{l=0}{\overset{L}{\mathrm{\Σ}}}p\left(l\right){[n+N]}^l=\underset{k=0}{\overset{K-1}{\mathrm{\Σ}}}h\left(k\right)\underset{l=0}{\overset{L}{\mathrm{\Σ}}}p\left(l\right){[n-k]}^l,---\left(4\right)$

将式(4)展开为L+1个等式，有：

$p\left(l\right){[n+N]}^l=\underset{k=0}{\overset{K-1}{\mathrm{\Σ}}}h\left(k\right)p\left(l\right){[n-k]}^l,---\left(5\right)$

消去等式两边的p(l)，可得：

$\underset{k=0}{\overset{K-1}{\mathrm{\Σ}}}{k}^{l}h\left(k\right)={(-N)}^l---\left(6\right)$

式中，l＝0，...L。式(6)称为多项式预测滤波器系数h(k)的约束条件，但仅仅通过式(6)无法唯一确定h(k)，还需要其他的约束条件。

考虑到通常遇到的信号都是混有噪声的，因此为了使噪声通过该滤波器增益最小，通常希望滤波器的噪声增益：

$\mathrm{NG}=\frac{1}{\mathrm{\π}}\underset{0}{\overset{\mathrm{\π}}{\∫}}{\left|H\left(e^{\mathrm{j\ω}}\right)\right|}^2\mathrm{d\ω}---\left(7\right)$

最小。式(7)中，H(·)为滤波器的传递函数。对于式(3)描述的FIR滤波器，式(7)可写为：

$\mathrm{NG}=\underset{k=0}{\overset{K-1}{\mathrm{\Σ}}}{\left|h\left(k\right)\right|}^2---\left(8\right)$

利用拉格朗日数乘法，在式(6)的约束下，式(8)有最小的最优解h(k)^[5]，例如：

当N＝1，L＝1时，

$h\left(k\right)=\frac{4K-6k-4}{K(K-1)},---\left(9\right)$

当N＝1，L＝2时，

$h\left(k\right)=\frac{{9K}^2+(-27-36k)K+30k^2+42k+18}{K^3-{3K}^2+2K}---\left(10\right)$

其中，k＝0，1，...K-1。当N和L为其它值时，请详见文献[5]的计算方法及计算结果。

式(2)至(10)的推导告诉我们，当h(k)(k＝0，...，K-1)为式(8)的最小的最优解时(如由式(9)或(10)给出时)，多项式模型式(2)和FIR预测滤波器式(3)在数学上严格等效。h(k)(如式(9)和(10)同时也表明，式(3)的滤波器系数只与N，K，L有关，而与多项式模型式(2)的多项式系数p(l)(l＝0，...，L)无关。这意味着用式(3)的预测滤波器来建模符合多项式形式的信号时，除该多项式的阶以外，不需有关信号其它的任何先验信息。而一旦知道多项式的阶，式(3)的模型参数h(k)，即滤波器的系数可以由式(8)的最小最优解(如式(9)或(10))的预先计算获得。

对于作匀变速运动的目标，牛顿运动定律告诉我们目标的运动可由下述方程描述：

$J\left(t\right)=J\left(0\right)+v\left(0\right)\·t+\frac{1}{2}a\·t^2---\left(11\right)$

v(t)＝v(0)+a·t                (12)

其中，J(t)，v(t)为目标在t时刻的位移和速度，J(0)，v(0)为目标的初始位置和初始速度，a为目标的加速度。显然式(11)和(12)是用多项式描述的目标运动，简称它们为多项式信号。对于这样的多项式信号，前面的分析告诉我们，当式(11)和(12)被离散化后，在数学上它们可以用形如式(3)的预测滤波器严格等效。即：

$J(n+N)=\underset{k=0}{\overset{K-1}{\mathrm{\Σ}}}h\left(k\right)J(n-k)---\left(13\right)$

等效式(11)，而

$v(n+N)=\underset{k=0}{\overset{K-1}{\mathrm{\Σ}}}h\left(k\right)v(n-k)---\left(14\right)$

等效式(12)。其中：式(13)和(14)的系数h(k)(k＝0，...，K-1)可由式(9)或(10)预先计算获得。当式(2)给出机动目标位置和/或速度的观测时，即式(11)和(12)离散化的状态空间模型可以统一等效描述为：

X(n+1)＝F_ppfX(n)                (15a)

Y(n)＝H_ppfX(n)+W(n)            (15b)

其中，式(15b)是根据测量原理而确定的测量方程，X(n+1)＝[x(n+1)，x(n)，...，x(n-K+2)]^T，X(n)＝[x(n)，x(n-1)，...，x(n-K+1)]^T，观测向量为Y(n)＝[y(n)，y(n-1)，...，y(n-K+1)]^T，观测噪声向量为W(n)＝[w(n)，w(n-1)，...，w(n-K+1)]，方差E(W_jW_n ^T)＝diag(R(n)，...，R(n-K+1))，误差R(n)＝E(w(n)w^T(n))，y(n)＝cx(n)+w(n)，状态转移矩阵为 $F_{\mathrm{ppf}}=\left[\begin{array}{cccc}h\left(0\right)& h\left(1\right)& ...& h(K-1)\\ 1& 0& ...& 0\\ ...& ...& ...& ...\\ 0& ...& 1& 0\end{array}\right];$ h(0)，…h(K-1)由式(9)或(10)计算得出；状态观测矩阵为 $H_{\mathrm{ppf}}=\left[\begin{array}{cccc}c& 0& ...& 0\\ 0& c& ...& 0\\ ...& ...& ...& ...\\ 0& 0& 0& c\end{array}\right].$

C为与测量原理有关的测量函数，由实际测量对象确定。

需要强调的是，本发明中的机动目标模型式(15a、15b)没有像通常的状态模型那样，包含一个描述状态模型不确定性的噪声，这是因为我们描述目标运动的模型如式(9)和(10)所表明那样与式(11)和(12)的多项式系数无关，即式(15a)不会因为目标的初始位置J(0)，初始速度v(0)和目标的加速度a的不准确而引入不确定性。

由式(15a)可以看出，本发明提出的机动目标的多项式预测模型形式简单，并且将一个非线性的运动方程表示成了线性形式，因此如果观测方程也是线性的，那么将本发明提出的模型与经典卡尔曼滤波算法结合，即可得到目标的最优跟踪算法。可是，众所周知，在滤波器开始工作的初始时刻以及状态发生突变的时刻，卡尔曼滤波器需要经历一段过渡过程才能跟踪上目标的状态。而对于式(15)描述的机动目标多项式预测模型，虽然其状态模型(15a)可以描述任一常加速度情况下的运动规律，可是，如果机动目标在某一时刻其加速度的从一个常加速度突变到另一个常加速度，而包含加速度突变点的机动目标运动规律不符合多项式所能描述的运动规律。那么基于式(15)描述的机动目标的卡尔曼滤波跟踪需要经历一段过渡过程才能跟踪上目标的状态。针对这种状况，下面我们将讨论如何探测机动目标的加速度是否发生变化，如何缩短利用卡尔曼滤波器跟踪时的过渡过程，如何加速跟踪算法收敛到稳态值的方法。

假如在n+1时刻，机动目标在该时刻其加速度的从一个常加速度突变到另一个常加速度，而包含加速度突变点的机动目标运动规律不能用式(15a)精确描述，即n+1时刻的状态不能由前K个时刻的状态准确预测。那么我们将式(15a)可写为：

X(n+1)＝F_ppfX(n)+Q_ppf(n+1)            (16)

其中，Q_ppf(n+1)为附加状态噪声，即用它描述表示多项式预测的状态模型不能描述的预测误差。此时卡尔曼滤波的增益应该改写为：

$K(n+1)=\frac{P^{-}(n+1)H_{\mathrm{ppf}}^T}{(H_{\mathrm{ppf}}{P}^{-}(n+1)H_{\mathrm{ppf}}^T+R(n+1))}$

$=\frac{(F_{\mathrm{ppf}}P\left(n\right)F_{\mathrm{ppf}}^T+Q_{\mathrm{ppf}}(n+1))H_{\mathrm{ppf}}^T}{H_{\mathrm{ppf}}(F_{\mathrm{ppf}}P\left(n\right)F_{\mathrm{ppj}}^T+Q_{\mathrm{ppf}}(n+1))H_{\mathrm{ppj}}^T+R(n+1)}---\left(17\right)$

$=\frac{I}{H_{\mathrm{ppf}}+\frac{R(n+1)}{(F_{\mathrm{ppf}}P\left(n\right)F_{\mathrm{ppf}}^T+Q_{\mathrm{ppf}}(n+1))H_{\mathrm{ppf}}^T}}$

式中R(n+1)是观测噪声在时刻n+1的方差，P^-(n+1)是状态在n+1时刻的预测方差，P(n)是状态在n时刻的估计方差。

众所周知，当状态模型及其有关它的假设都是正确的时候，卡尔曼滤波器的新息为零均值的高斯白噪声^[6]，那么就可以通过新息序列的均值检验^[6]来判断多项式预测的状态模型是否可以准确描述目标的运动规律。卡尔曼滤波新息表达式可以写成：

$S(n+1)=Y(n+1)-H_{\mathrm{ppf}}{\hat{X}}^{-}(n+1)---\left(18\right)$

其中，Y(n+1)为测量向量，

为状态预测，而其新息方差则可以表示为：

P_s(n+1)＝H_ppf ^TP^-(n+1)H_ppf ^T+R(n+1)        (19)

根据文献[6]，S^T(n+1)P_s ^-1(n+1)S(n+1)为具有m个自由度的χ²变量，m为新息的维数。因此可以通过常用的χ²检验法来验证新息均值是否为零。即在显著性水平α之下，如果 $S^T(n+1)P_s^{-1}(n+1)S(n+1)>{\mathrm{\χ}}_{\mathrm{\α}}^2$ ，则认为新息均值不为零，即此时存在多项式预测的状态模型不能描述的加速度突变情况。

据式(17)，如果式(15)的模型是精确的，那么随着时间的增长，观测信息的不断增加，估计的精度将不断提高，因此状态估计值误差方差阵P(n)的范数将不断减小，因此新的观测数据对状态的修正作用逐渐地减弱。当式(15)的模型不准确时，例如加速度突变，尽管计算的误差方差阵的范数随着时间的增长而不断减小，但实际的估计值误差却可能不断的增大，此时，状态估计的预测不准确，而新的观测之中却含有较多的状态模型的信息，预测应设法增加这些观测值对状态估计的修正作用^[6]。由式(17)可以看出，Q_ppf(n+1)取值越大，新息加权值K(n+1)越大。因此可以在式(15)多项式预测的状态模型不能准确描述目标的运动规律时，增加观测值对状态估计的修正作用，以加速卡尔曼滤波器跟踪算法过渡过程结束。因此我们可以在显著性水平α，当 $S^T(n+1)P_s^{-1}(n+1)S(n+1)>{\mathrm{\χ}}_{\mathrm{\α}}^2$ 时。取Q_ppf(n+1)＝βR(n+1)，β为一个较大的正整数，使得 $\frac{R(n+1)}{(F_{\mathrm{ppf}}P\left(n\right)F_{\mathrm{ppf}}^T+Q_{\mathrm{ppf}}(n+1))H_{\mathrm{ppf}}^T}\≈0,$ K(n+1)≈I/H_ppf即可，其作用是在式(15)多项式预测的状态模型不能准确描述目标的运动规律时，用新息来加快卡尔曼滤波跟踪算法的收敛速度，而在同样的显著性水平α，当 $S^T(n+1)P_s^{-1}(n+1)S(n+1)\≤{\mathrm{\χ}}_{\mathrm{\α}}^2$ 时，即式(15)多项式预测的状态模型能正确描述目标的运动规律时，让Q_ppf(n+1)＝0。

当目标加速度两次突变的间隔小于算法的收敛时间，即一次突变发生后算法尚未收敛到稳态时，就发生了另一次突变，此时由于新息序列始终无法满足均值为零的条件，算法退化为状态模型不准确性情况下的卡尔曼滤波算法。根据文献[2]，此时如果状态噪声足够大，即在状态更新环节新息权值足够大，算法不会出现数学意义上的发散，但无法达到稳态最优。这种情况下，状态估计结果依赖于新息对状态预测结果的更新。也就是说，只要观测值是准确的，状态估计结果仍然可以跟踪上目标实际状态的变化，但估计结果的方差大于稳态值。

因此，基于多项式预测模型，能加速卡尔曼滤波器过渡过程的最优跟踪算法的计算步骤如下：

假设n＝0时的状态

和方差P₀已知，对n＝1，2，...

(1)根据文献[5]确定L和K，通常对一阶信号，取L＝1，K＝2，对二阶信号，取L＝2，K＝3，根据式(9)和(10)求出h(k)。

(2)预测：

${\hat{X}}^{-}(n+1)=F_{\mathrm{ppf}}\hat{X}\left(n\right)---\left(20\right)$

$P^{-}(n+1)=F_{\mathrm{ppf}}P\left(n\right)F_{\mathrm{ppf}}^T---\left(21\right)$

(3)新息检验

在显著性水平α之下，如果 $S^T(n+1)P_s^{-1}(n+1)S(n+1)>{\mathrm{\χ}}_{\mathrm{\α}}^2,$ 则认为新息均值不为零，Q_ppf(n+1)＝βR_ppf(n+1)，β为一个较大的正整数，使得K(n+1)≈I/H_ppf，否则，Q_ppf(n+1)＝0。

(4)更新

$P^{-}(n+1)=F_{\mathrm{ppf}}P\left(n\right)F_{\mathrm{ppf}}^T+Q_{\mathrm{ppf}}(n+1)---\left(22\right)$

$K(n+1)=P^{-}(n+1)H_{\mathrm{ppf}}^T(H_{\mathrm{ppf}}{P}^{-}(n+1)H_{\mathrm{ppf}}^T+R_{\mathrm{ppf}}(n+1))---\left(23\right)$

$\hat{X}(n+1)={\hat{X}}^{-}(n+1)+K(n+1)(Y(n+1)-H_{\mathrm{ppf}}{\hat{X}}^{-}(n+1))---\left(24\right)$

$P(n+1)=(I-K(n+1)H_{\mathrm{ppf}})P^{-}(n+1)$ (25)

$-(I-K(n+1)H_{\mathrm{ppf}})\left(F_{\mathrm{ppf}}P\left(n\right)F_{\mathrm{ppf}}^T\right)$

本发明根据匀变速运动的多项式形式，将变加速运动作为分段匀变速运动处理，提出了一种新的机动目标动态模型-多项式预测模型。该模型能够完全避免对描述未知运动的多项式系数的依赖，无论在什么时候描述未知运动的多项式系数发生了变化，用本发明方法建立的状态方程始终能准确描述运动的动态。因此本发明提出的模型能够适用于各种运动目标状态方程，包括低空、高空、超高空运动的，点迹目标及视频目标的状态方程；水下、水面、地面运动的点迹目标及视频目标的状态方程等。本发明提出的多项式预测模型的最优滤波算法适用于任何可以用多项式描述的机动目标状态估计问题，既不需要关于目标运动参数的先验信息，又不需要通过参数辨识等手段调整建立状态模型的参数，从而完全避免了由于状态模型不准确而带来的估计性能恶化问题。

附图说明

图1为机动目标跟踪仿真1，仿真中采用本文方法与交互多模型方法进行了比较，其中，(a)为x方向位置估计标准差，(b)为x方向速度估计标准差，(c)为y方向位置估计标准差，(d)为y方向速度估计标准差。

图2为机动目标跟踪100次蒙特卡洛仿真1。其中，(a)为100次蒙特卡洛仿真x方向位置估计平均标准差，(b)为100次蒙特卡洛仿真x方向速度估计平均标准差，(c)为100次蒙特卡洛仿真y方向位置估计平均标准差，(d)为100次蒙特卡洛仿真y方向速度估计平均标准差。

图3为机动目标跟踪仿真2，仿真中采用本文方法与交互多模型方法进行了比较，其中(a)为x方向位置估计标准差，(b)为x方向速度估计标准差，(c)为y方向位置估计标准差，(d)为y方向速度估计标准差；

图4为机动目标跟踪100次蒙特卡洛仿真2，其中，(a)为100次蒙特卡洛仿真x方向位置估计平均标准差，(b)为100次蒙特卡洛仿真x方向速度估计平均标准差，(c)为100次蒙特卡洛仿真y方向位置估计平均标准差，(d)为100次蒙特卡洛仿真y方向速度估计平均标准差。

具体实施方式

对于机动目标跟踪仿真1，假设目标在0-60秒，过程噪声标准差为0.2米/秒²，61-120秒过程噪声标准差为10米/秒²，120-180秒过程噪声标准差为0.2米/秒²，测量噪声标准差为100米。首先采用L＝2，K＝3的多项式模型对目标位置建模，采用L＝1，K＝2的多项式模型对目标速度建模；然后分别估计目标的位置和速度，即在一个循环内，先利用我们多项式模型和其相应的跟踪算法式(20)-(25)，根据测量值估计目标的位置；再把位置估计值的差分作为目标速度的测量值，同样利用多项式模型和其相应的跟踪算法估计目标速度。根据步骤3进行新息均值判断时，α＝0.01。当新息均值不为零时， $Q_{\mathrm{ppf}}={\mathrm{\β\σ}}_w^2{I}_{2K\×2K}.$ 其中，I_2K×2K为2K×2K的单位矩阵，位置估计中，β＝100，速度估计中，β＝100。从图1-2中可以看出，本发明算法的估计标准差与交互多模型算法相当。

对于机动目标跟踪仿真2，假设目标在0-180秒内，过程噪声标准差始终为10米/秒²，其余参数与仿真1中相同，从图3-4，可以看出，本发明算法的估计标准差小于交互多模型算法。

为了验证本发明提出的机动目标状态模型的合理性，以及多项式预测模型最优估计算法在机动目标跟踪中的实用性和有效性，对本发明算法和目前最为有效的机动目标跟踪算法-交互多模型算法^[7，8]进行了多种情况下的计算机仿真和比较，从图1-图4可以看出，本发明算法可以取得较好的估计效果。

参考文献

[1]X.Rong Li and Vesselin P.Jilkov，Survey of Maneuvering Target Tracking.Part I：DynamicModels.IEEE Trans.On Aerospace and electronic systems，Vol.39，No.4，Oct.2003.

[2]Fitzgerald RJ.Divergence of the kalman filter[J].IEEE Trans.on Automatic Control.1971，16(6)：736-747.

[3]邱恺，黄国荣，陈天如，杨亚莉.卡尔曼滤波过程的稳定性研究(J).系统工程与电子技术.2005(1)，Vol.27，No.1：33-35.

[4]Pearson J，Goodall R，EasthamM，et al.Investigation of kalman filter divergence usingrobust stability techniques[A].Decision and Control 1997[C].Proceedings of the 36th IEEEConference on，1997，5：4892-4893.

[5]P.Heinonen，Y Neuvo.FIR-median hybrid filters with predictive FIR substructures[J]，IEEETransactions on Acoustics，Speech，and Signal processing，1988，36(6)：892-899.

[6]张金槐，线性模型参数估计及其改进，长沙：国防科技大学出版社，1999。

[7]E.Mazor，A.Averbuch，Y.bar-Shalom et al；Interacing Multiple Model Methods in TargetTracking：A Survey.IEEE Trans.On AES，1998，Vol.34，No.1：103-123

[8]T.Kirubarajan，Y.Bar-Shalom，Kalman Filter Versus IMM Estimator：When Do We Need theLatter？IEEE Trans.On Aerospace and Electronic Systems，Vol.39，No.4，Oct.200

Claims (3)

1、一种机动目标状态方程的多项式预测模型，其特征在于假设信号x(n)是由下述方程描述的一个L阶多项式：

$x\left(n\right)=\underset{l=0}{\overset{L}{\mathrm{\Σ}}}p\left(l\right)n^1,---\left(2\right)$

其中，p(l)为多项式的系数，l＝0，...，L；

如果用信号x(n)的前K个时刻的值x(n-K+1)，...，x(n)来预测信号式(2)的将来值x(n+N)，即：

$x(n+N)=\underset{k=0}{\overset{K-1}{\mathrm{\Σ}}}h\left(k\right)x(n-k),---\left(3\right)$

对于作匀变速运动的目标，其运动由下述方程描述：

$J\left(t\right)=J\left(0\right)+v\left(0\right)\·t+\frac{1}{2}a\·t^2---\left(11\right)$

v(t)＝v(0)+a·t       (12)

其中，J(t)，v(t)为目标在t时刻的位移和速度，J(0)，v(0)为目标的初始位置和初始速度，a为目标的加速度；当式(11)和(12)被离散化后，在数学上它们用预测滤波器严格等效，即：

$J(n+N)=\underset{k=0}{\overset{K-1}{\mathrm{\Σ}}}h\left(k\right)J(n-k)---\left(13\right)$

等效式(11)，

$(n+N)=\underset{k=0}{\overset{K-1}{\mathrm{\Σ}}}h\left(k\right)v(n-k)---\left(14\right)$

等效式(12)，其中：式(13)和(14)的系数h(k)，k＝0，...，K-1，由式(9)或(10)预先计算获得；

当式(2)给出机动目标位置和/或速度的观测时，式(11)和(12)离散化的状态空间模型统一等效描述为：

X(n+1)＝F_ppfX(n)       (15a)

Y(n)＝H_ppfX(n)+W(n)    (15b)

其中，式(15b)是根据测量原理而确定的测量方程，X(n+1)＝[x(n+1)，x(n)，...，x(n-K+2)]^T，X(n)＝[x(n)，x(n-1)，...，x(n-K+1)]^T，观测向量为Y(n)＝[y(n)，y(n-1)，...，y(n-K+1)]^T，观测噪声向量为W(n)＝[w(n)，w(n-1)，...，w(n-K+1)]，方差E(W_jW_n ^T)＝diag(R(n)，...，R(n-K+1))，误差R(n)＝E(w(n)w^T(n))，y(n)＝cx(n)+w(n)，状态转移矩阵为 $F_{\mathrm{ppf}}=\left[\begin{array}{cccc}h\left(0\right)& h\left(1\right)& \·\·\·& h(K-1)\\ 1& 0& \·\·\·& 0\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ 0& \·\·\·& 1& 0\end{array}\right];$ h(0)，…h(K-1)由式(9)或(10)计算得出；状态观测矩阵为 $H_{\mathrm{ppf}}=\left[\begin{array}{cccc}c& 0& \·\·\·& 0\\ 0& c& \·\·\·& 0\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ 0& 0& 0& c\end{array}\right];$

这里C为与测量原理有关的测量函数，由实际测量对象确定；

当N＝1，L＝1时，

$h\left(k\right)=\frac{4K-6k-4}{K(K-1)},---\left(9\right)$

当N＝1，L＝2时，

$h\left(k\right)=\frac{9K^2+(-27-36k)K+30k^2+42k+18}{K^3-3K^2+2K}---\left(10\right)$

其中，k＝0，1，...K-1。

2、根据权利要求1所述多项式预测模型，其特征在于所述的(15a)进一步改写为

X(n+1)＝F_ppfX(n)+Q_ppf(n+1)    (16)

其中，Q_ppf(n+1)为附加状态噪声，即用它描述表示多项式预测的状态模型不能描述的预测误差。

3、一种基于如权利要求1所述的多项式预测模型的跟踪算法，其特征在于具体步骤如下：

假设n＝0时的状态

和方差P₀已知，对n＝1，2，...

(1)根据文献[5]确定L和K，对一阶信号，取L＝1，K＝2，对二阶信号，取L＝2，K＝3，根据式(9)和(10)求出h(k)。

(2)预测：

${\hat{X}}^{-}(n+1)=F_{\mathrm{ppf}}\hat{X}\left(n\right)---\left(20\right)$

$p^{-}(n+1)=F_{\mathrm{ppf}}P\left(n\right)F_{\mathrm{ppf}}^T---\left(21\right)$

(3)新息检验

在显著性水平α之下，如果 $S^T(n+1)P_s^{-1}(n+1)S(n+1)>{\mathrm{\χ}}_{\mathrm{\α}}^2$ ，则认为新息均值不为零，Q_ppf(n+1)＝βR_ppf(n+1)，β为一个较大的正整数，使得K(n+1)≈I/H_ppf，否则，Q_ppf(n+1)＝0；

(4)更新

$P^{-}(n+1)=F_{\mathrm{ppf}}P\left(n\right)F_{\mathrm{ppf}}^T+Q_{\mathrm{ppf}}(n+1)---\left(22\right)$

$K(n+1)=P^{-}(n+1)H_{\mathrm{ppf}}^T(H_{\mathrm{ppf}}{P}^{-}(n+1)H_{\mathrm{ppf}}^T+R_{\mathrm{ppf}}(n+1))---\left(23\right)$

$\hat{X}(n+1)=\widehat{X^{-}}(n+1)+K(n+1)(Y(n+1)-H_{\mathrm{ppf}}{\hat{X}}^{-}(n+1))---\left(24\right)$

$P(n+1)=(I-K(n+1)H_{\mathrm{ppf}})P^{-}(n+1)---\left(25\right).$

$=(I-K(n+1)H_{\mathrm{ppf}})\left(F_{\mathrm{ppf}}P\left(n\right)F_{\mathrm{ppf}}^T\right)$

Images