面向非线性MTT的多模型泊松多伯努利混合滤波算法

陈嵩杰,李 波,张 露

(辽宁工业大学 电子与信息工程学院,辽宁 锦州 121001)

0 引 言

多目标跟踪[1](Multi Target Tracking,MTT)是一个重要的研究方向,适用于航空航天、空中交管和自动驾驶等领域[2,3].考虑MTT共轭性所反映的目标先验与后验分布概率密度具有相同形式的特点[4],基于多伯努利共轭先验理论的带标签广义标签多伯努利(Generalized Label Multi-Bernoulli,GLMB)滤波[5]和不带标签的泊松多伯努利混合(Poisson Multi-Bernoulli Mixture Filter,PMBM)滤波[6]已有研究.其中,PMBM滤波可为目标概率密度提供精确闭式解,计算代价和精度较GLMB滤波器更优.文献[7]提出了一种多模型多伯努利滤波算法,给出了该算法在线性高斯模型的实现;文献[8]提出了一种多伯努利混合滤波的高斯实现;文献[9]提出了一种基于多模型泊松多伯努利混合滤波的高斯逆威沙特实现,利用高斯逆威沙特分布估计目标质心,并引入强跟踪滤波修正状态模型协方差矩阵.

针对非线性模型的MTT问题,现有文献主要采用序贯蒙特卡洛(Sequential Monte Carlo,SMC)和高斯混合(Gaussian Mixture,GM)两种方式实现.其中,SMC实现[10,11]存在粒子退化问题,滤波计算繁琐,且跟踪速度慢;GM实现[12]的计算复杂程度相对较低,适用于线性高斯模型的MTT.基于此,文献[13]提出了一种非线性模型的扩展卡尔曼滤波(Extended Kalman Filter,EKF)的PMBM算法,但EKF在强非线性系统中存在较低的鲁棒性,其雅克比矩阵计算也较为困难;文献[14]提出了一种非线性模型的(Cubature Kalman Filter-Multi-Bernoulli,CKF-MB)和(Square-rooted CKF-MB,SCKF-MB)滤波算法,但在低检测概率时,目标势估计存在一定偏差;文献[15]提出了一种用于非线性模型δ-GLMB滤波的积分卡尔曼高斯混合实现,但有效性不高.针对非线性高杂波环境的MTT,文献[16]提出了一种均方根容积卡尔曼滤波与δ-GLMB结合的GM实现,但仍存在平衡滤波运算时间和跟踪精度难以平衡的问题.

针对MTT的非线性与低检测概率问题,本文在非线性模型下融合SCKF与高斯混合多模型泊松多伯努利混合滤波(GM Multiple Model PMBM,GM-MM-PMBM),提出一种基于SCKF的GM-MM-PMBM滤波算法.首先,推导出MM-PMBM预测和更新的GM实现.然后,采用误差协方差的平方根方式进行递推,应用等权值容积点集计算多目标密度函数的均值与协方差矩阵.最后,结合SCKF对GM-MM-PMBM滤波高斯项进行预测和更新,采用最优子模式[17,18](Optimal Sub-Pattern Assignment,OSPA)距离评价该算法的性能.

1 相关工作

1.1 随机有限集

在随机有限集(Random Finite Set,RFS)理论框架中[19],假设单目标状态x∈nx,多目标状态X∈F(nx),X为单目标状态向量的集合,F(nx)为nx空间中所有有限子集集合,目标状态可由量测集合Z∈F(nz)中的量测观察.

假设k时刻的状态集为X={x1,…,xn};测量集为Z={Z1,…,Zn}∪κ(z),κ(z)为杂波量测集合且满足κ(z)=λc(z),λ为泊松率,c(z)为杂波空间概率分布函数,Zi为目标i产生的量测集合,κ(z)与Z1…Zn相互独立.

1.2 PMBM共轭先验

PMBM滤波可由泊松和多伯努利混合两部分表示[20].其中,泊松部分用于表示未检目标;多伯努利混合部分则表示已检目标[21].假设X=Xu+∪Xd为目标集合,Xu和Xd分别为未检目标和已检目标,则共轭先验可定义为:

(1)

(2)

(3)

式中,wi,a为权重,+∪表示不相交并集;fp(·)为泊松部分概率密度函数;μ(·)为泊松强度;fmbm(·)为多伯努利混合部分概率密度函数;A为全局假设索引集合,a为所有全局假设索引,n为所有伯努利分量的数目.

于是,第a个全局假设中第i个伯努利量为:

(4)

式中,ri,a(·)和fi,a(·)分别为第a个全局假设中第i个的目标存在概率和伯努利分量.

1.3 标准PMBM算法

1.3.1 预测过程

由式(1)～式(3)可知,PMBM滤波器的预测过程包含泊松部分和多伯努利混合部分.假设k-1时刻泊松强度为μk-1(xk-1),则k时刻泊松强度预测密度函数为:

μk|k-1(xk)=λb(xk)+〈μk-1(xk-1),pSf(xk|xk-1)〉

(5)

式中,λb(·)为新生强度,pS为目标存活概率,f(xk|xk-1)为状态转移函数,[·,·]表示函数的内积.

同时,多伯努利部分预测参数可表示为:

(6)

1.3.2 更新过程

更新过程包含3个部分:未检目标的泊松强度更新、首次检测目标的泊松强度更新、已检目标的伯努利更新[22].

未检目标的泊松强度更新为:

(7)

首次检测目标的伯努利更新为:

(8)

(9)

已检目标(漏检目标和量测目标)的伯努利更新为:

(10)

(11)

可以看出,泊松强度更新和伯努利更新过程是分别进行的,首次检测目标的伯努利更新和已检目标的伯努利更新差异是前者中的伯努利分量含泊松分量.

2 多伯努利混合滤波非线性实现

由于标准PMBM算法在递归步骤中存在难以求解RFS积分的问题[23],因此基于高斯线性模型推导出该算法的GM实现,将积分形式传递的伯努利参数转化为高斯项加权求和形式.

为表述简洁,下文中的各函数主要以变量ξ定义,且假设检测概率pd(·)和存活概率pS(·)与模型相关,满足pd(xk,ξ)=pd(ξ),pS(xk,ξ)=pS(ξ).于是,PMBM算法的GM实现的预测和更新递推如下.

2.1 GM-MM-PMBM算法

2.1.1 预测过程

基于高斯分布N(·),新生目标泊松强度的GM形式可记为:

(12)

假设k-1时刻泊松强度的GM形式为:

(13)

则k时刻预测强度的GM形式为:

(14)

式中,Fk(ξ)为模型ξ的状态转移矩阵,Qk(ξ)为模型ξ的过程噪声协方差矩阵,fk|k-1(ξ|ξ′)为模型ξ′～模型ξ的转移概率.

假设k-1时刻伯努利密度函数为:

(15)

则预测过程伯努利分量的权重保持不变,其存在概率满足:

(16)

于是,多伯努利部分预测密度函数为:

(17)

2.1.2 更新过程

未检目标的泊松强度更新:假设泊松预测强度GM形式由式(19)给出.

(18)

则更新后的后验强度为:

μk|k(xk,ξ)=(1-pd)μk|k-1(xk,ξ)

(19)

首次检测到目标的泊松强度更新:

(20)

(21)

式中,相关参数定义如下:

已检目标的伯努利更新:对于漏检目标定义下述公式.

(22)

(23)

(24)

对于量测目标则有:

(25)

(26)

(27)

式中,相关参数定义如下:

2.2 基于SCKF的GM-MM-PMBM算法

考虑运动模型ξ,假设目标状态方程和量测方程分别为:

(28)

式中,xk和zk分别为状态向量和量测向量;fk(·)和hk(·)分别为状态函数和量测函数;uk-1是均值为0、协方差矩阵为Qk-1的过程噪声;vk是均值为0、协方差矩阵为Rk的量测噪声.

假设k-1时刻误差协方差的平方根满足下三角阵Γk-1|k-1=Tria(A),Tria(·)为矩阵的三角化运算,于是基于SCKF的非线性滤波过程如下:

2.2.1 预测过程

由Cholesky分解计算预测高斯项的协方差矩阵:

(29)

利用三阶球面-径向容积准则生成容积点χq,k-1及权重wq=1/2n,q=1,2,…,2n,这里系统状态维数设为2n.于是,计算非线性状态函数的传递容积点:

(30)

计算状态预测均值:

(31)

计算预测误差协方差矩阵的平方根:

(32)

2.2.2 更新过程

产生容积点χq,+并计算非线性量测函数传递容积点:

(33)

计算量测预测均值:

(34)

计算新息协方差平方根:

(35)

计算互协方差矩阵:

(36)

计算卡尔曼增益:

(37)

计算状态均值更新:

(38)

计算协方差矩阵的平方根更新:

(39)

3 仿真验证与结果分析

3.1 仿真情形及参数设置

目标运动模型ξ设为匀速(Constant Velocity,CV)运动和匀速转弯(Constant Turn,CT)运动,其状态转移方程分别为:

式中,wk-1是均值为0、协方差矩阵为Qk-1的过程噪声.

目标非线性量测模型方程为:

观测矩阵和量测噪声协方差矩阵分别为:

式中,标准差σε=10,I2为二维单位方阵.此外,目标存活概率为pS=0.99,由泊松分布对杂波建模.模型1定义为CV运动,设过程噪声为5m/s2;模型2定义为CT运动,设过程噪声为5 m/s2,逆时针转动角度为3°/s;模型3定义为CT运动,设过程噪声为5m/s2,顺时针转动角度为3°/s.于是,定义CV和CT运动的线性状态转移矩阵如下:

假设目标新生模型为泊松形式,新生目标泊松强度GM为:

所有存活目标在前20s进行CV运动,在第21～ 40s内按逆时针方向由CV运动转换为CT运动,在第60～80s内按顺时针方向转换为CT运动,模型间转移概率矩阵定义为:

图1给出了检测区域和3个目标的非线性运动轨迹.

图1 目标真实轨迹Fig.1 Target real trajectory

3.2 实验结果分析

为比较本文SCKF-GM-MM-PMBM算法与EKF-PMBM算法[13]、SCKF-GM-MM-MB算法[14]、SMC-MM-PMBM算法的跟踪性能,本实验进行100次SMC仿真.

3.2.1 固定杂波密度及高检测概率情形下

假设检测概率为pd=0.96,杂波密度为λC=0.6×10-3/m2,目标势估计及误差估计仿真结果如图2所示.

图2 高检测概率下势估计Fig.2 Cardinality estimation under high detection probability

可以看出,随着时间的变化,各算法都能较好完成势估计.EKF-PMBM算法、SCKF-GM-MM-MB算法和本文算法都能接近目标真实数目,而经典SMC-MM-PMBM算法在粒子数目较少时,会出现低于目标真实值的估计情形;并在第15s和第30s时,SMC-MM-PMBM算法估计时延较长且对目标势估计有明显偏差,EKF-PMBM算法和SCKF-GM-MM-MB算法的目标数目估计时延较低,而本文算法的估计时延和消耗较最低,能有效地估计目标势变化.

本文采用OSPA距离作为MTT性能的评价指标,该距离越大表示算法的总体精度越差,设定参误差调节参数c为100,阶次参数p为1[24].各滤波算法的OSPA距离对比如图3所示.可以看出,各算法具有相似的OSPA误差,在第5s和第18s时,SCKF-GM-MM-MB算法给出的OSPA距离略优于SCKF-GM-MM-PMBM算法.若目标在第20s出现,在第60s和第80s时目标消失,本文算法优于SCKF-GM-MM-MB算法,其原因为PMBM算法中泊松部分保留了漏检分量和刚消失的目标分量,能及时检测新生目标并延迟对消亡目标的响应.EKF-PMBM算法和SMC-MM-PMBM算法的OSPA距离相比于另两种算法偏大,其原因为强非线性下EKF需先进行线性处理,忽略了高阶项,出现发散问题,导致跟踪目标丢失;而SMC的采样粒子较多易发生粒子退化现象,若当采样粒子数目较少时,目标跟踪性能和运算速度降低.

图3 OSPA距离对比Fig.3 OSPA distance comparison

3.2.2 固定杂波密度及低检测概率情形下

假设目标检测概率为pd=0.6,仿真结果如图4所示.

图4 低检测概率下势估计Fig.4 Cardinality estimation under low detection probability

对比高检测概率环境,可以看出SCKF-GM-MM-PMBM算法的势估计与SCKF-GM-MM-MB算法、EKF-PMBM算法和SMC-MM-PMBM算法有着明显的差异,SCKF-GM-MM-PMBM算法具有更好的势估计误差.

各算法的OSPA距离对比如图5所示.可以看出,当第62s目标消亡时,本文算法的OSPA距离误差较SCKF-GM-MM-MB算法大,其原因为较低的Pd使本文算法在更新过程中的泊松分量权重增大,无法及时删除刚刚消失的目标状态.

图5 OSPA距离对比Fig.5 OSPA distance comparison

表1综合比较了OSPA距离误差和单步平均运行时间.可以看出,本文算法在低检测概率时下具有更好的跟踪性能.

表1 不同检测概率的OSPA距离和单步平均运行时间Table 1 OSPA distance and single-step average running time under different detection probabilities

3.2.3 不同杂波密度情形下

为进一步验证所提算法的跟踪性能,给出了各算法在不同杂波密度λC=(0.6,1.2,1.8,2.4,3.0)×10-3/m2时的势估计误差与OSPA距离对比.

图6和图7分别为势估计误差对比图与OSPA距离对比图.可以看出,当杂波密度增大时,各算法的势估计误差与OSPA距离误差增大,EKF-PMBM算法仅次于SCKF-GM-MM-MB算法,其原因为弱非线性环境下EKF-PMBM算法可优先考虑,在强非线性情形下EKF算法与线性假设原则冲突,导致滤波发散并产生估计误差,但本文算法在不同杂波密度情况下,其势估计误差和OSPA距离误差都要优于其他3种算法.

图6 势估计误差对比Fig.6 Cardinality estimation and error comparison

图7 OSPA距离对比Fig.7 OSPA distance comparison

表2对比了当检测概率pd=0.96时不同杂波密度环境下的OSPA距离与单步平均运行时间.可以看出,当同一检测概率时,各算法的OSPA距离和运行时间都随杂波密度的增加而增大.当不同杂波时,本文算法具有较高跟踪精度,运行时间低于其他算法.当杂波密度较高时,由于EKF-PMBM算法进行雅克比矩阵计算,导致运行时间增加;而SMC-MM-PMBM算法所需的采样粒子数目不断增加,导致运行时间增加.

表2 不同杂波密度下的OSPA距离和单步平均运行时间Table 2 OSPA distance and single-step average running time under different clutter densities

综上,本文算法在不同检测概率情形下能以较高的精度估计目标的势和运动状态.尤其是在低检测概率情形下,其跟踪精度和运行时间优于其他3种算法;在高检测概率不同杂波密度情形下,各算法的OSPA距离与运行时间具有一定差距,本文算法的OSPA距离优于其他3种算法.

4 总 结

本文针对MTT存在的非线性与低检测概率问题,提出了MM-PMBM算法的GM实现,在低检测概率环境下为目标概率密度提供精确闭式解,其跟踪精度和运行时间上相比MM-MB算法有所提高;在非线性系统模型下,融合SCKF滤波框架,避免了计算协方差矩阵的平方根,其数值稳定性和滤波性能得以保证;在不同杂波密度情形下,该算法可准确估计目标的势与运动状态,提高了MTT的跟踪精度.由不同情形下的仿真实验结果可知,本文算法在非线性情形中OSPA距离误差与运行时间均最小.

针对MTT的运行时间,接下来将深入研究本文算法用于解决机动目标、扩展目标和群目标等比较复杂的问题,进一步节省运行时间.