一种基于变维滤波算法的机动目标跟踪方法与流程

本发明属于信号处理技术，具体涉及变维滤波算法。

背景技术：

机动目标跟踪属于信号处理领域的一项技术，对国防和经济建设具有重要作用，一直以来很多学者致力于研究目标跟踪算法。但由于现代飞行器的机动性不断增强，对目标进行可靠并准确的跟踪已经成为一大难题。

目前机动目标跟踪算法大致可以分为两类，一类是无须机动检测的自适应跟踪算法，另一类是具有机动检测的跟踪算法。所谓机动检测，实质是利用目标的残差对目标状态进行判别和调整。目前常用的机动检测算法包括可调白噪声模型、变维滤波算法、输出估计算法等。其中，变维滤波算法是在1982年由bar-shalom和birmiwal在variabledimensionfilterformaneuveringtargettracking中提出。变维滤波算法采用“开关”型切换，根据检测手段确定目标的机动情况并在非机动模型和机动模型之间进行切换。原本的变维滤波算法的优点在于不依赖目标机动的先验假设，但滤波器在非机动状态和机动状态之间的切换往往导致目标状态估计产生较大误差，从而增大了跟踪误差。

技术实现要素：

本发明的目的在于提出一种基于变维滤波算法的机动目标跟踪方法，提高对机动目标的跟踪性能。

实现本发明的基本思路是：利用衰减记忆平均值来检测机动的发生，一旦机动假设被接受，则利用当前时刻的量测值和机动发生时刻的状态估计初始化机动时刻的状态估计，同时更新该时刻的状态协方差矩阵，以使机动模型下的滤波快速收敛。

本发明技术方案为一种基于变维滤波算法的机动目标跟踪方法,该方法包括：

步骤1：建立目标运动模型

1a)分别建立非机动模型和机动模型的运动状态方程：

xk＝fxk-1+wk-1

其中，xk和xk-1分别表示非机动模型在k和k-1时刻的状态向量，xk和表示目标k时刻位置和速度，和分别表示机动模型在k和k-1时刻的状态向量；和分别表示目标k时刻位置、速度和加速度，f和f^m分别为非机动和机动时的状态转移矩阵；wk-1和为k-1时刻离散时间白噪声序列，均值为零，方差分别为qδkj和q^mδkj，q和q^m表示非机动模型和机动模型的噪声协方差矩阵，δkj表示二维冲激函数，当k＝j的时候为1，否则为0，k，j随机产生，上标m表示处于机动模式情况下；

1b)用下式建立目标非机动模型和机动模型的量测方程：

zk＝hxk+vk

其中，zk表示k时刻雷达的量测值；h和h^m分别为非机动模型和机动模型的量测矩阵；vk为量测噪声，其均值为0，方差为rk，并与过程噪声相互独立，rk为量测噪声的协方差矩阵；

步骤2：非机动模型的状态预测

2a)由步骤1建立的非机动模型和上一时刻状态估计完成目标状态的一步预测：

其中，表示k-1时刻预测目标在k时刻的状态；表示目标k-1时刻状态更新值；

2b)目标状态协方差矩阵的一步预测由下式决定：

pk|k-1＝fpk-1|k-1f^t+q

其中，pk|k-1表示k-1时刻预测k时刻的状态协方差；pk-1|k-1表示k-1时刻状态协方差更新值；f^t为非机动状态转移矩阵的转置；q表示噪声协方差矩阵；

2c)由状态预测值按下式确定目标在k-1时刻量测值的预测：

其中，表示k-1时刻预测目标在k时刻的量测值；表示k-1时刻预测目标在k时刻的状态；

2d)利用下式确定滤波新息(残差)：

其中，vk为k时刻的滤波新息，zk为k时刻量测值；

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

sk＝hpk|k-1h^t+rk

其中，sk为预测误差的协方差矩阵；pk|k-1表示k-1时刻预测目标在k时刻的状态协方差；矩阵h^t为量测矩阵的转置；rk为量测噪声的协方差矩阵；

步骤3：检测目标机动是否发生

3a)按照下式计算滤波新息的衰减记忆平均值：

μ＝1-1/s

其中，ρk和ρk-1分别表示k和k-1时刻的滤波新息的衰减记忆平均值，为归一化滤波新息的平方,且服从自由度为1/(1-μ)的卡方分布；μ为折扣因子，s为滑窗长度，v'k表示vk的转置；

3b)按照下式检测机动发生：

p{ρk≤ρmax}＝1-α

其中，ρmax为机动假设的阈值，α为所设定的显著性水平；如果ρk超过下式所设定的阈值ρmax，则接受机动发生的假设，在阈值点上估计器从非机动模型转换为机动模型，进入步骤5；否则，拒绝接受机动发生的假设，进入步骤4；

步骤4：非机动模型的状态更新

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

其中，kk为k时刻的滤波器增益；pk|k-1表示k-1时刻预测目标在k时刻的状态协方差；sk为预测误差的协方差矩阵；h^t为量测矩阵的转置；

4b)按照下式完成对目标状态的更新：

其中，表示目标k时刻状态更新值；表示k-1时刻预测目标在k时刻的状态；zk表示目标量测值；表示k-1时刻预测目标在k时刻的量测值；

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

其中，pk|k表示目标k时刻状态协方差更新值；为滤波器增益矩阵的转置；

4d)重复步骤2～3；

步骤5：机动模式的初始化

5a)当k时刻检测到机动，认为k-s时刻目标开始机动，因此，增加加速度维度，对机动模型在k-s时刻的状态估计进行初始化，s为设定的固定数值，表达式如下：

其中，表示机动模型在k-s时刻的状态更新值，表示非机动模型在k-s时刻对位置的估计，表示非机动模型在k-s时刻对速度的估计，zk表示就一维状态下k时刻的量测值，t表示雷达采样间隔；

5b)按照下式完成对k-s时刻状态协方差矩阵的初始化：

其中，表示k-s时刻机动模型的状态协方差矩阵，表示中第i行第j列的元素，pk-s|k-s表示k-s时刻非机动模型的状态协方差矩阵，pk-s|k-s(i,j)表示pk-s|k-s中第i行第j列的元素；

步骤6：机动模型的状态预测

6a)由建立的机动模型和上一时刻状态更新值完成目标状态的一步预测：

其中，表示k-1时刻预测目标在k时刻的状态；表示目标k-1时刻状态更新值；f^m表示机动时的状态转移矩阵；

6b)目标状态协方差矩阵的一步预测由下式决定：

其中，表示k-1时刻预测k时刻的状态协方差；表示k-1时刻状态协方差更新值；q^m表示噪声协方差矩阵；

6c)由状态预测值按下式确定目标在k-1时刻量测值的预测：

其中，表示k-1时刻预测目标在k时刻的量测值；表示k-1时刻预测目标在k时刻的状态；h^m表示目标机动时的量测矩阵；

6d)利用下式确定滤波新息(残差)：

其中，为k时刻的滤波新息，zk为k时刻量测值；

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

其中，为预测误差的协方差矩阵；表示k-1时刻预测目标在k时刻的状态协方差；(h^m)^t表示量测矩阵的转置；rk为量测噪声的协方差矩阵；

步骤7：检测机动是否结束

7a)当出现加速度突然下降到0的情况，即机动突然结束，当机动模型的新息超过95％的置信区域时，切换到非机动模型，进入步骤2；

7b)按下式计算非机动检测的统计量：

其中，δk表示k时刻加速度估计显著性检验的统计量，表示k时刻加速度分量的估计，表示机动模型的协方差矩阵与加速度分量相对应的块，表示在长度为p的滑窗上加速度估计显著性检验统计量之和；

7c)当落在所设定的阈值以下时，认为加速度是不显著的，切换到非机动模型，进入步骤2；否则，认为目标保持机动状态，进入步骤8；

步骤8：机动模式的状态更新

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

其中，为k时刻的滤波器增益；表示k-1时刻预测目标在k时刻的状态协方差；(h^m)^t表示量测矩阵的转置；为预测误差的协方差矩阵；

8b)按照下式完成对目标状态的更新：

其中，表示目标k时刻状态更新值；表示k-1时刻预测目标在k时刻的状态；zk表示目标量测值；表示k-1时刻预测目标在k时刻的量测值；

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

其中，表示目标k时刻状态协方差更新值；为滤波器增益矩阵的转置；

8d)重复步骤6～7。

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

第一、本发明与变维滤波算法相比，当k时刻检测到机动，认为非机动模型下k-s时刻的状态估计是合理的，并利用k-s时刻的状态估计以及k时刻的量测值计算k-s时刻的加速度，从而提高了加速度估计的准确性，减小了模型切换时的误差。

第二、本发明与变维滤波算法相比，利用k-s时刻加速度估计对机动模型的状态协方差矩阵中相关项进行调整，减小了加速度方差，从而使滤波器能够快速准确地调整滤波状态，加快了滤波收敛的速度。

附图说明

图1为本发明的流程示意图；

图2为本发明与变维滤波算法的仿真误差对比图。

具体实施方式

下面结合图1，本发明实现步骤如下：

步骤1，建立目标运动模型

以周期t采样，对目标分别建立非机动模型和机动模型的运动状态方程：

xk＝fxk-1+wk-1

其中，xk和xk-1分别表示非机动模型在k和k-1时刻的状态向量，xk和表示目标k时刻位置和速度，和分别表示机动模型在k和k-1时刻的状态向量，和分别表示目标k时刻位置、速度和加速度，f和f^m分别为非机动和机动时的状态转移矩阵，表达式如下：

wk-1和为离散时间白噪声序列，均值为零，方差分别为为e[wkwj]＝qkδkj和

本发明实例中，给定则过程噪声方差矩阵分别为：

非机动模型和机动模型的量测方程为：

zk＝hxk+vk

其中，zk表示k时刻雷达的量测值；h和h^m分别为非机动模型和机动模型的量测矩阵，表达式如下：

h＝[10]

h^m＝[100]

vk为量测噪声，其均值为0，方差为rk，并与过程噪声相互独立。

步骤2，非机动模型的状态预测

由建立的非机动模型和上一时刻状态更新值完成目标状态的一步预测：

其中，表示k-1时刻预测目标在k时刻的状态；表示目标k-1时刻状态更新值；f表示状态转移矩阵。

目标状态协方差矩阵的一步预测由下式决定

pk|k-1＝fpk-1|k-1f^t+q

其中，pk|k-1表示k-1时刻预测k时刻的状态协方差；pk-1|k-1表示k-1时刻状态协方差更新值，f^t表示状态转移矩阵的转置，q表示噪声方差矩阵。

由状态预测值按下式确定目标在k-1时刻量测值的预测：

其中，表示k-1时刻预测目标在k时刻的量测值；表示k-1时刻预测目标在k时刻的状态，h为量测矩阵。

利用下式确定滤波新息(残差)：

其中，vk为k时刻的滤波新息，zk为k时刻量测值。

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

sk＝hpk|k-1h^t+rk

其中，sk为预测误差的协方差矩阵；pk|k-1表示k-1时刻预测目标在k时刻的状态协方差；h^t为量测矩阵的转置；rk为量测噪声的协方差矩阵。

步骤3，检测目标机动是否发生

按照下式计算滤波新息的衰减记忆平均值：

μ＝1-1/s

其中，ρk和ρk-1分别表示k和k-1时刻的滤波新息的衰减记忆平均值，服从自由度为1/(1-μ)的卡方分布；μ为折扣因子，s为滑窗长度，本发明实例中s＝5。

按照下式检测机动发生：

p{ρk≤ρmax}＝1-α

其中，ρmax为机动假设的阈值，α为所设定的显著性水平；本发明实例中α取0.005，根据卡方分布表，自由度为5时对应的阈值ρmax＝16.75。

如果ρk超过下式所设定的阈值ρmax，则接受机动发生的假设，在阈值点上估计器从非机动模型转换为机动模型，进入步骤5；否则，拒绝接受机动发生的假设，进入步骤4；

步骤4，非机动模型的状态更新

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

其中，kk为k时刻的滤波器增益；pk|k-1表示k-1时刻预测目标在k时刻的状态协方差；h^t为量测矩阵的转置；[·]^-1表示矩阵的逆。

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

其中，表示目标k时刻状态更新值；表示k-1时刻预测目标在k时刻的状态；表示k-1时刻预测目标在k时刻的量测值；zk表示目标量测值，kk滤波器增益矩阵。

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

其中，pk|k表示目标k时刻状态协方差更新值；sk为预测误差的协方差矩阵；为滤波器增益矩阵的转置。

重复步骤2～3。

步骤5，机动状态的初始化

当k时刻检测到机动，认为k-s时刻目标开始机动，因此，增加加速度维度，对机动模型在k-s时刻的状态估计进行初始化，表达式如下：

其中，表示机动模型在k-s时刻的状态更新值，表示非机动模型在k-s时刻对位置的估计，表示非机动模型在k-s时刻对速度的估计。

按照下式完成对k-s时刻状态协方差矩阵的初始化：

其中，表示k-s时刻机动模型的状态协方差矩阵，表示中第i行第j列的元素，pk-s|k-s表示k-s时刻非机动模型的状态协方差矩阵，pk-s|k-s(i,j)表示pk-s|k-s中第i行第j列的元素。

步骤6，机动模型的状态预测

由建立的机动模型和上一时刻状态更新值完成目标状态的一步预测：

其中，表示k-1时刻预测目标在k时刻的状态；表示目标k-1时刻状态更新值；f^m表示目标状态转移矩阵。

目标状态协方差矩阵的一步预测由下式决定：

其中，表示k-1时刻预测k时刻的状态协方差；表示k-1时刻状态协方差更新值；(f^m)^t表示目标状态转移矩阵的转置；q^m表示噪声协方差矩阵。

由状态预测值按下式确定目标在k-1时刻量测值的预测：

其中，表示k-1时刻预测目标在k时刻的量测值；表示k-1时刻预测目标在k时刻的状态；h^m表示目标量测矩阵。

利用下式确定滤波新息(残差)：

其中，为k时刻的滤波新息，zk为k时刻量测值。

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

其中，为预测误差的协方差矩阵；表示k-1时刻预测目标在k时刻的状态协方差；(h^m)^t表示量测矩阵的转置。

步骤7，检测机动是否结束

当出现加速度突然下降到0的情况(即机动突然结束)，即机动模型的新息超过95％的置信区域时，切换到非机动模型，进入步骤2。

按下式计算非机动检测的统计量：

其中，δk表示k时刻加速度估计显著性检验的统计量，表示k时刻加速度分量的估计，表示机动模型的协方差矩阵与加速度分量相对应的块，表示在长度为p的滑窗上加速度估计显著性检验统计量之和，本发明实例中p取2。

当落在所设定的阈值以下时，认为加速度是不显著的，切换到非机动模型，进入步骤2；否则，认为目标保持机动状态，进入步骤8；

本发明实例中，在显著性水平为0.005时，根据卡方分布表，自由度为2对应的阈值

步骤8，机动模式的状态更新

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

其中，为k时刻的滤波器增益；表示k-1时刻预测目标在k时刻的状态协方差；(h^m)^t表示量测矩阵的转置；为预测误差的协方差矩阵。

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

其中，表示目标k时刻状态更新值；表示k-1时刻预测目标在k时刻的状态；zk表示目标量测值；表示k-1时刻预测目标在k时刻的量测值。

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

其中，表示目标k时刻状态协方差更新值；为滤波器增益矩阵的转置。

重复步骤6～7。

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

1.仿真条件

设目标的真实航迹的起始状态为[10000m,100m/s,0m/s²]，对目标进行180s的采样观测，目标的具体运动如下：

在1～30s内，目标做匀速直线运动，在31～60s内加速度由30m/s²减小到0，在61～90s时，目标做匀加速直线运动，在91～120s内加速度由0增加到30m/s²，在121～180s内做匀速直线运动。

在笛卡尔坐标系下，采用1000次蒙特卡洛实验，设雷达采样间隔t＝1s，雷达的量测位置方差r＝2500m²。

仿真的评价指标为均方根误差，即rmse，其计算公式为：

其中，和分别表示第i次运行k时刻状态的第j分量的真值和估计值。

2.仿真内容

采用本发明方法和变维滤波算法分别对目标的位置、速度进行跟踪估计，并比较跟踪效果。

3.仿真分析

图2(a)是本发明和变维滤波算法对一维目标进行跟踪时位置均方根误差，实线表示本发明的跟踪误差曲线，虚线表示变维滤波算法的跟踪误差曲线。从图中可以看出，当目标的机动性较弱时，本发明与变维滤波算法的位置误差水平相近；当目标机动性较强时，本发明能够有效地降低位置误差(相对降低约40％～50％)，并且滤波收敛速度要快于变维滤波算法。

图2(b)是本发明和变维滤波算法对一维目标进行跟踪时速度均方根误差，从图中可以看出，当目标的机动性较弱时，本发明与变维滤波算法的速度误差水平几乎一致；当目标机动性较强时，本发明对速度跟踪性能保持较好，而变维滤波算法几乎丧失对速度的跟踪能力。因此，本发明极大改善了变维滤波算法在目标机动时对速度的跟踪性能。