侯利明,连峰,王伟
(1.西安交通大学智能网络与网络安全教育部重点实验室,710049,西安;2.陕西省组合与智能导航重点实验室,710068,西安;3.中国电子科技集团公司第二十研究所导航重点实验室,710068,西安)
近20年来,随机有限集(RFS)理论[1-2]被广泛地应用于多目标跟踪领域。与传统的多目标跟踪算法[3-5]相比,基于RFS理论的多目标跟踪算法不需要明确量测与目标之间复杂繁琐的对应关系,尤其适用于多目标跟踪问题。目前,基于RFS的多目标跟踪滤波器主要有概率假设密度(PHD)滤波器[6]、势概率假设密度(CPHD)滤波器[7]、多目标多伯努利(MeMBer)滤波器[1]和δ-广义标记多伯努利(δ-GLMB)滤波器[8-9]。特别地,δ-GLMB滤波器不仅在跟踪精度方面优于其他RFS滤波器,且能自动生成目标航迹,是首个理论可证的具有精确闭式解的多目标贝叶斯跟踪滤波器。由于δ-GLMB滤波器的递推过程是严格封闭的,因此它对目标个数及其状态的估计精度高于其他同类滤波器,但同样面临计算量大的问题。该滤波器的计算复杂度理论上与目标的个数呈线性关系,与量测的个数呈3次方关系[10],但在其实现过程中仍需要花费较大的计算代价用于假设排序问题。文献[8-9]给出了该滤波器在线性高斯模型下的高斯混合(GM)实现和非线性模型下的序贯蒙特卡洛(SMC)实现方法,其中,由于SMC-δ-GLMB滤波算法受粒子数和粒子退化等问题的影响,尽管可以通过增加粒子数以获得较高的跟踪精度,但大量的粒子同样意味着需要消耗更大的计算量,难以实现δ-GLMB滤波器对目标的快速跟踪。在实际的多目标跟踪应用中,非线性模型则更为普遍,如雷达和声呐等量测均为非线性[11],因此研究非线性模型下δ-GLMB滤波器的多目标跟踪问题具有非常重要的意义。
GM-δ-GLMB滤波器在线性高斯模型下具有闭式解,并且可以将GM-δ-GLMB滤波器与扩展卡尔曼(EK)[11]和无味卡尔曼(UK)[12]相结合用于解决非线性模型下的多目标跟踪问题。由于EK算法只是简单地将非线性系统模型进行局部线性化处理,存在严重的模型描述误差,当模型的非线性程度较高时,该滤波算法可能会导致被跟踪目标丢失,且在对非线性函数进行Taylor级数展开时,需要求解Jacobians矩阵和Hessians矩阵(二阶EKF),而这2种矩阵的计算容易导致数值不稳定,且计算量较大,更有甚者,它们根本就不存在[13]。UK算法则是通过选择一组确定的采样点来近似非线性函数的均值和方差,精度可达到二阶,但其采样点及权系数一般由状态维数和凭惯例或经验而设定的参数共同来决定,很难达到较高精度的跟踪效果。关于EK和UK滤波算法的上述限制条件,直接阻碍了δ-GLMB滤波器向非线性多目标跟踪应用领域的推广。
本文将积分卡尔曼(QK)[14-15]非线性滤波算法与GM-δ-GLMB滤波器相结合,提出了具有更加严谨的数学理论基础且适用性更强的QK-GM-δ-GLMB滤波器。为了验证该算法的有效性,通过在不同的杂波密度和检测概率条件下的仿真实验,将所提算法的滤波精度和实时性能分别与EK-GM-δ-GLMB、UK-GM-δ-GLMB和SMC-δ-GLMB 3种滤波器做了较详细的对比。实验结果表明:QK-GM-δ-GLMB滤波器的跟踪精度是4种滤波器中最优的,尽管其时间消耗有所增加,但与另外两种δ-GLMB滤波器的GM实现算法均保持在同一量级,这在许多实际应用中是可以接受的。
假设每个目标均服从如下的非线性运动和量测模型,即
xk=f(xk-1)+wk-1
(1)
zk=h(xk)+vk
(2)
式中:k为观测时刻;xk为k时刻目标的状态向量;zk为k时刻传感器的观测向量;f(·)和h(·)分别为已知的非线性状态转移函数和非线性观测函数;wk-1和vk分别表示相互独立的过程噪声和量测噪声,并假设它们服从均值为零、协方差矩阵分别为Qk-1和Rk的高斯分布。
积分卡尔曼滤波算法的主要思想是基于Gauss -Hermite数值积分规则求取一组带权重的积分点,利用这些积分点求取密度函数的均值和协方差估计[14-15]。假设向量x~N(x;0,Inx),则nx维的积分规则可由下式表示
(3)
通过构造一个具有零对角元素的对称三对角矩阵J来求取其特征值,其中,εl表示矩阵J的第l个特征值,且该矩阵其他元素可由下式求得
(4)
根据上述所求取的一组带权重积分点,可以推导出基于QK的δ-GLMB滤波器的高斯混合实现过程,具体预测步骤如下。假设k时刻的多目标后验概率密度具有如下的δ-GLMB RFS形式[9]
(5)
式中:ξ为直到k时刻的所有关联映射对应于标签集合I的关联历史,ξ=(θ1,…,θk),θk表示k时刻航迹标签到传感器量测的关联映射;X表示带标签的多伯努利状态集合;数对(I,ξ)表示当前时刻所有航迹的标签集合I∈F(L)与直到当前时刻的标签-量测关联历史ξ∈Ξ之间的一种对应关系;w(I,ξ)表示该假设为真的概率;F(L)表示由标签空间L的子集所构成的集合;离散空间Ξ表示航迹标签到量测的关联过程。直到k时刻的所有关联历史ξ条件下标签为l对应目标航迹的广义标记多伯努利分量的空间分布密度函数由J(ξ)(l)个高斯项组成
(6)
由δ-GLMB分布的共轭先验性[8]可知,k+1时刻的多目标预测分布密度也具有δ-GLMB形式,即
(7)
(8)
式中:pB(x,l)表示新生目标的状态空间密度;HL(l)与HB(l)分别表示存活目标和新生目标的示性函数。标签为l∈L的存活航迹的空间分布密度具有如下的高斯混合形式
(9)
(10)
(11)
(12)
为简单起见,假设目标航迹存活概率独立于标签l和状态x,即pS(x,l)=pS,航迹消亡概率为qS(x,l)=1-pS(x,l)=1-pS。
设新生目标概率密度可表示为如下的LMB RFS形式
(13)
(14)
(15)
(16)
式中:I+⊆L+=L∪B表示预测目标的标签空间,而新生目标航迹的假设wB(J)则与具体的应用场景相关。
用求得的积分点bl和相应权重wl来近似式(9)中第j个高斯项的均值和协方差,具体实现步骤如下。
步骤1对式(6)中组成标签为l的航迹中的第j个高斯项的协方差矩阵进行Cholesky分解
(17)
步骤2计算组成标签l的第j个高斯项的积分点
(18)
步骤3非线性状态方程传递积分点
(19)
步骤4计算预测状态均值
(20)
步骤5计算状态预测协方差矩阵
(21)
(22)
则更新后的后验密度仍然为δ-GLMB形式
δI+(L(X))[p(ξ,θ)(·|Z)]X}
(23)
式中Θ表示航迹标签到量测关联θ的映射空间,即L+→{0,1,…,|Z|},0表示漏检。此时的后验标签集合与预测航迹标签集合相等,仍为I+。更新假设的概率为
(24)
令k+1时刻标签l的航迹关联映射为θ(l),若θ(l)>0,表示该航迹被传感器收到的量测集合zθ(l)更新;若θ(l)=0,则表示没有量测与之关联,该航迹漏检。总的量测更新广义似然函数可表示为
(25)
假设检测概率独立于标签l和状态x,即pD(x,l)=pD,漏检概率qD(x,l)=1-pD(x,l)=1-pD。
对于更新航迹,标签为l的航迹关联映射满足θ(l)>0,式(25)的广义似然函数可具体表示为
(26)
将式(22)与式(26)相乘并进行归一化,即可得航迹l更新后的δ-GLMB空间分布密度
(27)
(28)
归一化常数
(29)
(30)
式中:q(ξ,j)(z(θ(l)),l)表示目标(x,l)生成量测z(θ(l))的似然概率;κ(z(θ(l)))表示量测过程中的杂波密度函数。
对于漏检航迹,因没有量测与之关联,故θ(l)=0,此时的广义似然函数可重新表示为
(31)
(32)
将式(22)与式(31)相乘并除以式(32),可得标签l的航迹漏检时的空间分布密度,此时的假设更新概率、状态均值和状态协方差矩阵分别对应于该航迹的预测空间密度式(22)中的各参数,即
(33)
(34)
(35)
步骤1对式(22)中组成航迹l的第j个高斯项的协方差矩阵进行Cholesky分解
(36)
步骤2计算标签l的第j个高斯项的积分点
(37)
步骤3非线性量测方程传递积分点
(38)
步骤4计算预测量测均值
(39)
步骤5计算量测预测协方差矩阵
S(ξ,j)=
(40)
步骤6计算互协方差矩阵
(41)
步骤7计算卡尔曼增益
K(ξ,j)=G(ξ,j)[S(ξ,j)]-1
(42)
步骤8均值更新
(43)
步骤9协方差矩阵更新
(44)
在二维平面内,设观测区域范围为[0°,180°]×[0 m,2 000 m],观测时长为50 s,传感器采样间隔T=1 s,最多可能同时有8个目标随机出现在该区域的任何位置,且目标出现和消失的时刻都是不确定的。目标真实运动轨迹如图1所示。
图1 极坐标下多目标的真实运动轨迹
目标的运动模型为CT模型[4],且观测数据存在角度和距离的观测噪声,k时刻目标的状态向量xk=[px,k,vx,k,py,k,vy,k,ωk]T由位置分量(px,k,py,k)、速度分量(vx,k,vy,k)和转弯速率ωk组成。CT模型的状态转移方程可表示为
xk=F(ωk-1)xk-1+Gkwk-1
(45)
F(ωk-1)=
(46)
(47)
含方位角和径向距离噪声的非线性量测方程为
(48)
目标的存活概率pS,k=0.99,检测概率pD,k=0.98。杂波模型服从强度为κk=λcVu(z)的泊松RFS,其中λc表示平均每帧杂波的个数,V为传感器监控区域的面积,u(·)为该区域内的均匀概率分布。
采用剪切和合并技术限制高斯项数目的增长,设航迹的截断概率阈值T′=10-5,高斯项权重的阈值T″=10-5,合并阈值U=4,航迹最大值Mmax=300,每条航迹对应高斯项个数的最大值Jmax=100,其具体剪切和合并过程与文献[16]所述一致。
实验通过100次蒙特卡洛(MC)仿真验证了所提算法的有效性,采用最优子模式分配(OSPA)距离DOSPA[17]来评价几种多目标跟踪算法的性能。本实验中,将OSPA表达式的误差调节参数c和阶次参数p分别设置为100、2。与已有的EK-GM-δ-GLMB、UK-GM-δ-GLMB和SMC-δ-GLMB滤波算法关于平均OSPA距离和平均单步运行时间在不同检测概率和不同的杂波强度条件下做了比较,其中,在SMC实现过程中,每个假设航迹被分配的最大粒子数为10 000,最小粒子数为2 000。
仿真实验运行环境软件和硬件参数为Matlab R2017b,Windows 7 64 bit,Intel Core i5-4570 CPU 3.20 GHz,RAM 4.00 GB。
图2给出了4种算法在检测概率pD为0.98、杂波密度λc为3.2×10-3(rad·m)-1的实验场景中经100次MC仿真实验后DOSPA随时间变化曲线的对比。从图2可以看出,4种算法的DOSPA具有相同的变化趋势,说明在存在目标消失和新生目标出现的场景中,它们均可以准确跟踪目标,但在跟踪精度上有所差别。
图2 100次MC仿真实验后DOSPA的对比
由图2可知,各滤波器跟踪精度依QK-GM-δ-GLMB、EK-GM-δ-GLMB、UK-GM-δ-GLMB、SMC-δ-GLMB顺序依次降低,具体原因分析如下。
(1)在δ-GLMB滤波器的更新过程中涉及到量测信息与航迹标签之间的关联问题,而这种关联问题存在较大的不确定性,随着量测信息的不断增加,量测-标签之间的关联假设呈指数增长,造成大量粒子被用来逼近这些小权重的关联假设,使得大量的计算资源被用于粒子的重采样过程,从而导致跟踪精度与计算时间效率都较低。SMC-δ-GLMB算法的跟踪精度在很大程度上取决于粒子数的大小,随着粒子数的增加,SMC-δ-GLMB算法的实现精度也会随之升高,但同时意味着需要更长的计算时间,不易实现快速跟踪。
(2)UK-GM-δ-GLMB算法的滤波参数值一般是凭经验和惯例而得到的,针对复杂和高状态维数的多目标跟踪问题往往难以凑效。
(3)在本次实验中,EK-GM-δ-GLMB算法的跟踪精度仅次于QK-GM-δ-GLMB算法,这是因为针对弱非线性问题,EK算法是首要考虑的方法,当非线性程度较高时,很容易出现滤波发散问题,从而导致被跟踪目标的丢失。
(4)本文所提的QK-GM-δ-GLMB滤波器具有较高的跟踪精度,因为其所有滤波参数均是基于Gauss-Hermite积分规则得到,所以得到的参数更加准确,从而可实现较高精度的目标跟踪。
为进一步详细比较4种算法的目标跟踪精度和实时性,在表1和表2中分别列出了不同场景中各滤波器经100次MC仿真实验后的DOSPA和单步运行时间的比较结果。
表1给出了当杂波密度λc为0.8×10-3(rad·m)-1时,在检测概率pD为0.80、0.90、0.98的条件下4种滤波器的DOSPA和单步运行时间的实验数据。由表1可以看出,在相同的杂波密度场景下,各滤波器的DOSPA和单步运行时间都随检测概率的增加而减小,其中,在高检测概率场景中,QK-GM-δ-GLMB算法的跟踪精度最高,而在较低检测概率场景中QK-GM-δ-GLMB与EK-GM-δ-GLMB算法同时具有最高跟踪精度,且实时性优于SMC-δ-GLMB算法。
表2给出了在检测概率pD为0.98的条件下,当杂波密度λc分别为0.8、1.6、3.2×10-3(rad·m)-1时4种滤波器的DOSPA和平均单步运行时间的实验数据。由表2可以看出,在同一检测概率条件下,各滤波器的DOSPA和单步运行时间都随着杂波密度的增强而增加,其中,在不同的杂波强度下,QK-GM-δ-GLMB算法具有最高的跟踪精度,同样其实时性表现处于UK-GM-δ-GLMB算法和SMC-δ-GLMB算法之间。当杂波密度较高时,SMC-δ-GLMB算法的计算代价是巨大的。
表1 不同检测概率条件下的DOSPA和平均单步运行时间
表2 不同杂波密度下的DOSPA和平均单步运行时间
综上所述,本文所提的QK-GM-δ-GLMB滤波算法在不同的场景中都能够以较高精度对目标进行准确跟踪,而其对应的时间消耗在实际应用中处于可接受的范围内,其适用场景也得到进一步扩展。
在非线性高斯多目标运动场景下,给出了QK-GM-δ-GLMB滤波算法的递推过程,QK-GM-δ-GLMB滤波算法具有良好的滤波性能,是一种适应性强、估计精度高的跟踪算法,其时间开销在实际的多目标跟踪应用中是完全可接受的,所以QK-GM-δ-GLMB滤波算法具有一定的工程应用价值。
由于QK-GM-δ-GLMB滤波算法在非线性场景中具有较高的跟踪精度,下一步将致力于将其应用于机动目标、扩展目标和群目标等较复杂的多目标跟踪问题。