变分贝叶斯概率数据关联算法

恽鹏 吴盘龙 李星秀 何山

目标跟踪在军事和民用领域均具有广泛的应用[1−3],而在现实工程中,例如,在海上目标跟踪[4]和空中目标跟踪[5]等领域,观测设备在获取与目标相关的真实量测同时会获取与目标不相关的虚警数据.虚警数据也被称为杂波[6],由于无法在混有杂波信息的数据集中直接识别真实量测,基于真实量测设计的滤波器[7]无法对目标进行跟踪.此外密集的杂波信号会影响真实目标的回波信号,导致雷达探测能力下降进而影响目标的跟踪精度.因此如何提高杂波环境下目标跟踪精度是当前目标跟踪领域亟需解决的问题之一.

许多学者均对杂波环境下的目标跟踪问题展开研究,并取得了不错的研究成果.对于单目标跟踪系统,最近邻算法[8−10]以及概率数据关联算法(Probabilistic data association algorithm,PDA)[11−14]是较为经典的解决方法.最近邻算法通过构建距离指标函数选取与预测中心最近的数据作为真实量测,然而当杂波距离预测中心较近时,其往往会错误地将杂波选取为真实量测,导致跟踪精度下降.PDA算法本质上是一种全邻算法,其认为关联波门内的每一个数据均可作为真实量测,并通过新息函数为每一个数据赋予相应的权重.由于该方法相比于最近邻算法可降低杂波对状态更新的影响,因此具有更高的目标跟踪精度.PDA 算法的设计思想同样被应用于解决杂波环境下机动目标跟踪问题,相关的算法包括交互多模型概率数据关联算法[15](Interactive multi-model probabilistic data association,IMM-PDA)、联合交互多模型概率数据关联算法[16]以及联合交互多模型距离加权概率数据关联算法[17]等.

从贝叶斯统计原理的角度来看,杂波环境下目标真实状态后验概率密度函数应为真实量测作为观测信息情况下的状态后验概率密度函数.由于无法直接识别真实量测,因此目标真实状态后验概率密度函数无法获取并且每一个数据作为真实量测而获取的状态后验概率密度函数均可能为真实状态后验概率密度函数.此外,当检测概率不为1 时,所有数据与目标不相关时获取的状态后验概率密度函数也可能成为真实状态后验概率密度函数.代理概率密度函数表示某一变量可能的概率分布,所有数据关联形式下获取的状态后验概率密度函数可构成一个状态代理概率密度函数集合.基于该集合的统计特征以实现状态近似后验概率密度函数的获取是一个合理的选择.PDA 算法通过提取每一个数据−状态后验概率密度函数的均值与协方差并基于加权融合的方式以实现状态近似后验概率密度函数的获取.这类近似方法虽然简单且易于实现,但是在数据数目较多的情形下会存在计算效率低的缺陷.本文在变分贝叶斯框架[18]下设计了一种新的概率数据关联算法(Variational Bayesian based probabilistic data association algorithm,VB-PDA),该算法首先将关联事件视为随机变量并利用多项分布对两者进行建模,随后基于量测集、目标状态、关联事件的联合概率密度函数求取关联事件的后验概率密度函数,最后将关联事件的后验概率密度函数引入变分贝叶斯框架中以获取状态近似后验概率密度函数.VB-PDA算法本质上是一种在PDA 建模框架下通过权重Kullback-Leibler (KL)平均[19]准则来实现状态近似后验概率密度函数获取的方法.相比于PDA 算法,VB-PDA 算法在提高计算效率的同时,获取的近似后验概率密度函数也更能体现状态代理概率密度函数集的统计特征.相关仿真实验结果验证了算法的有效性.

本文内容安排如下: 第1 节对相关问题进行了描述;第2 节对PDA 算法进行了介绍以及对VBPDA 算法进行了推导;第3 节从权重KL 平均角度对VB-PDA 算法进行了分析;第4 节对算法计算复杂度进行了分析;第5 节通过相关仿真实验对提出的算法进行验证;第6 节对本文进行了总结.

1 问题描述

假设线性离散系统状态空间模型的状态方程和量测方程分别为:

由于无法直接区分真实量测与杂波,状态真实后验概率密度函数无法获取.如何获取合适的状态近似后验概率密度函数是本文的研究目的.

2 算法介绍

不同的杂波分布会导致算法具备不同的数学形式.假设目标的检测概率为pd,门概率为pg.系统视场 Θ 被视为关联门,杂波在 Θ 内服从均匀分布,该区域体积为Vk+1.杂波数目服从泊松分布,杂波密度为λk+1.Vk+1、pd和λk+1被认为是已知的.pg表示真实量测在视场内的概率.其可由如下公式获取:

式中,R 表示实数域.本文介绍的PDA 算法以及设计的VB-PDA 算法均是基于上述杂波环境给出的.

2.1 PDA 算法

2.2 VB-PDA 算法

在VB-PDA 算法中,量测集中任意一个数据为真实量测的概率相同被视为先验信息.此外,pd、pg以及杂波相关信息同样被视为已知的先验信息并被使用构建相应的概率密度函数.将关联事件视为随机变量,由于关联事件的数目是有限的,因此可用多项分布进行建模.此外,关联事件数目的确定意味着量测数目同样确定.因此基于已有的先验信息,关联事件的先验概率密度函数可被建模为:

可以发现,基于本文建模框架下推导的关联事件权重模型与PDA 算法相同.利用VB 算法[18]对状态近似后验密度函数qVB(xk+1) 进行求解,其求解公式如下:

基于推导过程可以发现,状态后验概率密度函数在贝叶斯框架下可直接被获取,并且该概率密度函数服从高斯混合分布.然而随着状态更新的过程不断进行,该后验概率密度函数的高斯分量会呈几何式增长,导致其计算复杂度过大.考虑到实时性的影响,采用VB 算法获取状态近似后验概率密度函数.VB-PDA 算法获取的状态近似后验概率密度函数服从高斯分布,这样的结果会为算法在后续的状态更新过程中提供便利.此外,VB-PDA 算法获取的状态协方差小于状态预测协方差以及PDA 算法获取的状态协方差,因此该算法可有效避免奇异现象的产生.

3 VB-PDA 算法讨论

本节从代理概率密度函数集的角度说明VBPDA 算法的优势.权重KL 平均被认为是最能体现代理概率密度函数集统计特征的概率密度函数[19].权重KL 平均(Ψ) 定义为:

式中,L (·) 表示下界函数,VB 算法是一种基于平均场理论实现参数近似后验概率密度函数获取的方法,qVB(ΨVB) 可被分解为:

根据式(25)、式(27)和式(39),式(37)的等价形式为:

根据式(36)和式(42)可知,VB-PDA 算法获取的参数集近似后验概率密度函数即为PDA 算法建模框架下的状态代理概率密度函数集的权重KL平均,因此该算法获取的状态近似后验概率密度函数比PDA 算法获取的状态近似后验概率密度函数更能体现代理概率密度函数集统计特征.

4 计算复杂度分析

算法的计算复杂度同样是影响算法跟踪性能的重要指标.本文将加减运算与乘除运算次数作为指标对4种算法的计算复杂度进行分析.表1 展示了PDA算法、距离加权概率数据关联算法[17](Distance weighted probabilistic data association,DWPDA)以及VB-PDA 算法在一步状态更新过程中所需的加减运算与乘除运算次数.可以看出,VBPDA 算法比PDA 算法以及DW-PDA 算法拥有更低的算法复杂度.当杂波数目较大多时,VB-PDA算法计算效率优势更加明显.

表1 一步状态更新过程中所需的加减运算与乘除运算次数Table 1 The number of addition and subtraction operations and multiplication and division operations required in the process of one-step state update

5 仿真实验

5.1 仿真场景 1

利用目标跟踪案例对算法的跟踪性能进行测试.假设目标在二维平面进行匀速直线运动,利用匀速模型(Constant velocity,CV)对该目标进行跟踪并且量测设备可以获取目标的位置信息.因此,目标的状态信息xk,状态转移矩阵Fk,量测矩阵Hk+1分别为:

式中,ak和bk表示目标位置信息,和表示目标速度信息,t=1 s 表示采样时间.系统过程噪声协方差和系统噪声量测协方差为:

目标真实初始状态信息x0为[100 m,100 m,1 m/s,1 m/s],初始状态协方差P0|0为diag{100 m2,100 m2,1 m2/s2,1 m2/s2},目标初始状态信息x0|0为x0与协方差为P0|0高斯白噪声的和[20],该高斯白噪声由Matlab 程序随机产生.因此在每一次蒙特卡洛仿真实验中,目标初始状态信息是不同的.

在仿真实验中,将均方根误差(Root mean square error,RMSE)和平均时间均方根误差(Time average root mean square error,TRMSE)作为指标来比较各算法的估计性能,目标位置和速度的RMSE和TRMSE 计算公式如下:

图2 仿真场景1 下 3种算法的速度 RMSEFig.2 The RMSEof velocity from three algorithms in scenario 1

假设杂波分布在 [ 0∼400 m,0∼400 m] 区域里服从均匀分布,因此Vk+1= 160 000 m2.杂波数目服从泊松分布[21].检测概率为0.95,杂波数目统计期望为20.图1～2 展示了3种算法的位置与速度RMSE,表2 展示了3种算法的位置与速度TRMSE.由表2 可以看出,VB-PDA 算法的位置估计精度在3种算法中最好,而DW-PDA 算法的速度估计精度在3种算法中最好.从TRMSE 数据可以发现,相比于PDA 算法,VB-PDA 算法在位置和速度估计精度上分别提高了2.18%和0.11%.DW-PDA算法是在权重计算的过程中对PDA 算法的改进,而VB-PDA 算法是从状态后验概率密度函数获取的过程中对PDA 算法的改进.由于两种算法的改进方式不同,出现DW-PDA 算法在速度估计精度优于VB-PDA 算法这一现象属于正常情形.表3展示了4种算法一次蒙特卡洛实验所需的计算时间,可以看出VB-PDA 算法的所需的计算时间最少,这也验证了第4 节分析的结果.从数据中可以发现,相比于PDA 算法,VB-PDA 算法在计算效率上提高了17.37%.

图1 场景 1 下 3种算法的位置 RMSEFig.1 The RMSE of position from three algorithms in scenario 1

表2 场景 1 下 3种算法的 TRMSETable 2 The TRMSE of three algorithms in scenario 1

表3 场景 1 下一次蒙特卡洛仿真实验所需的计算时间Table 3 The computational time at one Monte Carlo simulation experiment in scenario 1

5.2 仿真场景 2

考虑基于机动目标跟踪案例对提出算法的跟踪性能进行测试,利用交互多模型(Interactive multimodel,IMM)算法[22]与上述3种算法结合对目标进行跟踪.IMM 算法中采用两个带不同转弯速率匀速转弯模型(CT)对该目标进行跟踪并且量测设备可以获取目标的位置信息.因此,目标的状态信息xk,子模型状态转移矩阵Fk,j(j=1,2),量测矩阵Hk+1,j(j=1,2) 分别为:

其中,子模型的转弯速率分别为w1=−π/20 rad/s,w2=π/20 rad/s. 系统过程噪声协方差Qk,j和量测噪声协方差Rk+1,j设置为:

目标真实初始状态信息x0为 [100 m,100 m,3 m/s,3 m/s],目标在1～50 s 做转弯速率为π/100 rad/s的匀速转弯运动,在51～100 s 做转弯速率为−π/100 rad/s 的匀速转弯运动.假设2 个子模型初始权重均为0.5,马尔科夫转移矩阵为子模型初始状态协方差为diag{100 m2,100 m2,1 m2/s2,1m2/s2}.为验证算法的有效性,对IMM-PDA 算法、IMM-DW-PDA 算法和IMM-VB-PDA 算法进行测试.

假设杂波分布在 [ 0∼400 m,0∼400 m] 区域里服从均匀分布,因此Vk+1= 160 000 m2.杂波数目服从泊松分布.检测概率为0.9,杂波数目统计期望为10.图3～4 给出了3种算法获取的位置与速度RMSE.表4 给出了3种算法获取的位置与速度TRMSE.由表4 可以看出,在此类场景下,IMMVB-PDA 算法的估计精度最佳.从TRMSE 数据可以发现,相比于IMM-PDA 算法,IMM-VB-PDA算法在位置和速度估计精度上分别提高了5.26%和4.8%.

表4 仿真场景2 下3种算法的TRMSETable 4 The TRMSE of three algorithms in scenario 2

图3 仿真场景 2 下 3种算法的位置 RMSEFig.3 The RMSE of position from three algorithms in scenario 2

6 结束语

图4 仿真场景 2 下 3种算法的速度 RMSEFig.4 The RMSE of velocity from three algorithms in scenario 2

PDA 算法采用从状态代理概率密度函数集里提取均值信息与协方差信息的方式获取状态近似后验概率密度函数.该方法虽然便捷,但在数据量较大情形下存在计算效率低的缺陷.针对这一问题,本文提出了一种VB-PDA 算法,该算法本质上是在PDA 算法建模框架下提取相应的权重KL 平均以获取状态近似后验概率密度函数.相比PDA 算法,VB-PDA 算法在提高计算效率的同时,获取了近似程度更高的状态后验概率密度函数.相关仿真实验结果验证了VB-PDA 算法的有效性.由于本文提出的算法是针对单目标跟踪问题设计的,如何将类似的思想拓展到多目标跟踪领域是未来的研究目标之一.