昆明理工大学信息工程与自动化学院, 云南 昆明 650500
收稿/录用/修回: 2018-07-04/2018-09-26/2018-12-11
基金项目: 国家自然科学基金资助项目（61271007）；云南省人才培养项目（14118844）
作者简介:
张文(1993-), 男, 硕士生.研究领域为信息融合.
赵宣植(1981-), 男, 博士, 讲师.研究领域为信息融合, 阵列信号处理.
刘增力(1966-), 男, 硕士, 教授.研究领域为现代信号处理及其应用, 水声信号处理.
Sparse Gauss-Hermite PHD Maneuvering Multi-target Tracking Algorithm
School of information engineering and automation, Kunming University of Science and Technology, Kunming 650500, China
0 引言
随着各种先进传感技术的飞速发展与日益完善，多目标跟踪技术已经在军事和民用领域进行了广泛的应用.传统的多目标跟踪方法主要有联合概率数据关联(JPDA)[1]和多假设跟踪(MHT)[2]，但在杂波环境且目标数目未知、时变的情况下这些方法计算代价很高.近年来，Mahler等学者提出基于随机集的多目标跟踪方法，通过随机有限集描述多目标状态集和传感器观测集，利用多目标贝叶斯滤波算法递推得到后验概率密度，从而估计目标数目和运动状态，成功地避免了复杂的数据关联，又具有较严格的数学理论基础，越来越受到国内外学者的重视[3-5].
该方法主要通过递推形式求解多目标后验概率密度的一阶矩，即概率假设密度(probability hypothesis density，PHD)，最后从PHD中提取出目标的个数和状态.文[6-7]提出采用蒙特卡洛(SMC)方法来实现PHD滤波器，即粒子PHD滤波.但在滤波过程中，该算法需要大量粒子来近似PHD且粒子存在“退化”的现象，在状态与目标的估计中，需要采用聚类技术，增加了算法的复杂度，同时也很难保证目标状态提取的精度.针对以上问题，Vo和Ma[8]提出在线性情况下通过混合高斯的形式来实现PHD滤波器，从而通过每个高斯分量的权值和分布情况估计目标的数目和状态.尽管GM-PHD滤波器的基本假设比粒子PHD滤波器更加严格，但它却能为PHD提供真正的解析形式的代数解，具有计算量小、易于实现等优点.但非线性情况下GM-PHD没有解析表达式，而工程应用中遇到的大部分问题都是非线性的，因此，在PHD框架下，国内外学者相继提出将扩展卡尔曼、容积卡尔曼和不敏卡尔曼与概率假设密度滤波相结合的算法[9-11].但此类算法的性能都受限于与其结合的非线性滤波算法，扩展卡尔曼采用非线性函数的一阶泰勒展开式近似非线性函数，对于强非线性系统效果较差；容积卡尔曼和不敏卡尔曼通过Sigma点近似高斯概率密度函数，克服了EKF在强非线性系统中的局限性，但滤波精度最多只能达到三阶. Ito等[12]在2001年首次提出将高斯厄米特数值积分方法与高斯滤波器结合的高斯厄米特求积分滤波器(Gauss-Hermite quadrature Filter，GHF)，文[13]详细分析了GHF的滤波精度和算法复杂度，通过增加积分点数可以使GHF达到任意阶泰勒精度，并将GHF应用到混合高斯滤波器中.但随着状态矢量的维度增加，所需Sigma点的数目将呈指数级增加，导致“维度灾难”，当维度大于4时很少使用GHF.为解决“维度灾难”问题，文[14]提出通过结合稀疏网格取点规则来减少GHF的Sigma点数，生成一种稀疏网格求积滤波器(Sparse Grid quadrature Filter，SGF)，相比GHF，该方法适用于高维情况，且对于高度非线性情况，可使用高阶高斯厄米特多项式来得到SGF，在保证滤波精度的前提下减少了计算成本.
现实场景中运动目标的运行规律是复杂而多变的，采用单一运动模式对运动目标进行模型化描述往往并不适用，在对机动目标进行跟踪时，单一运动模型容易造成滤波精度低、跟踪目标丢失等现象.因此，对机动目标跟踪一般采取多模型的方法，多模型算法假设任意时刻目标运动模型与有限个已知模型中的一个匹配，目标的真实运动模型会从一个到另一个随机跳变.目前主要的多模型算法主要有自主多模型算法、协作多模型算法和变结构多模型算法[15].在这些多模型算法里，IMM是最有效的跟踪算法之一，它通过合理的假设估计，表现出了良好的跟踪性能，并得到了广泛的应用[16-17].
针对以上问题，本文提出一种基于稀疏高斯厄米特滤波的交互式多模型混合高斯概率假设滤波器.首先将稀疏高斯厄米特与PHD滤波器相结合，实现非线性系统下的PHD滤波，再通过IMM算法，完成运动模型的判定和目标状态的估计融合.仿真结果表明，本文提出的算法具有较高的滤波精度且数目估计更准确.
1 概率假设密度滤波器
1.1 多目标跟踪随机有限集理论
由于多目标贝叶斯递推过程中包含多重积分，一般情况下常采用多目标后验概率密度的一阶矩即概率假设密度(probabilityhypothesisdensity，PHD)来代替多目标后验概率密度来进行计算. PHD滤波器的形式如下：
PHD预测器方程为
|
(1) |
其中，PHD的“伪马尔科夫转移密度”为
|
(2) |
PHD的校正方程为
|
(3) |
对于任意观测集Z，PHD的伪似然LZ(x)=LZ(x|Z(k))为
|
(4) |
其中，
|
(5) |
式中，Dk+1|k+1(·)表示多目标随机集的PHD，bk+1|k(x)表示新生目标强度，bk+1|k(x|x′)表示衍生目标强度，fk+1|k(x|x′)为单目标马尔可夫转移密度，Ps(x′)为目标存活概率，PD(x)为目标检测概率，λc(z)为杂波强度，Lz(x)为单目标似然函数.
1.2 混合高斯PHD
PHD滤波器极大地缓解了完全多目标贝叶斯滤波器的计算负载问题.但是式(1)和式(4)中涉及到多维积分，无法得到解析解，混合高斯近似PHD能有效地解决此问题[18].假定在k时刻，多目标先验PHD函数可以表示为混合高斯形式：
|
(6) |
式(6)中，Jk表示k时刻高斯分量个数，xk|ki、Pk|ki、wk|ki分别表示第i个高斯分量的均值、方差和权值.
在k+1时刻，GM-PHD滤波器预测和更新结果表示为
|
(7) |
|
(8) |
式(7)中Jk|k+1表示k+1时刻预测得到的高斯分量个数，xk+1|ki、Pk+1|ki、wk+1|ki分别表示预测得到的高斯分量的均值、方差和权值. (8)式中PD为目标检测概率，xk+1|k+1i、Pk+1|k+1i、wk+1|k+1i分别表示更新后的高斯分量的均值、方差和权值.
2 基于交互式多模型的高斯厄米特PHD算法
2.1 高斯厄米特滤波
高斯厄米特滤波是一种基于高斯厄米特数值积分方法的递推贝叶斯滤波，它通过确定性点采样方法来提高非线性系统的后验估计精度[13].对于一维函数f(x)，高斯厄米特数值求积形式如下：
|
(9) |
式中，xi和ωi为第i个积分点和它对应的权值，其选取规则为：构造一个主对角线上元素全为0的对称矩阵J，矩阵中元素
，i=1，…，s-1，则一维高斯厄米特Sigma点
，λi为矩阵J的第i个特征值，权值wi=(vj)12，其中(vj)1为第j个归一化特征向量中的首个元素.
对于n维情况x~N(x；0，In)，通过张量积可得n维高斯厄米特求积形式[19]：
|
(10) |
其中，χi =[χi1，χi2，…，χin]T为第i个Sigma点，wi= 
为第i个Sigma点的权值，s为用于逼近的高斯厄米特多项式的阶数.
当目标的状态转移方程和量测方程已知，系统噪声和量测噪声都是均值都为零的高斯噪声，方差分别为R和Q，k时刻目标的状态估计均值和方差为xk|k和Pk|k，通过高斯厄米特数值积分方法生成Sigma点，得
|
(11) |
高斯厄米特状态预测方程为
|
(12) |
|
(13) |
观测方程的预测为
|
(14) |
|
(15) |
|
(16) |
|
(17) |
2.2 稀疏网格取点
为使得高斯厄米特积分更易于实现，结合Smolyak[20]提出的稀疏网格的取点方法，高斯厄米特积分形式变为
|
(18) |
In，L(f)表示对n维向量的L级精度积分，Iij表示精度为ij的单变量积分逼近，ij∈Φ，Φ为n个自然数的精度序列，Nqn为精度序列集合，定义为
Xij为单维积分Iij的积分点集，式(18)中Iij具体形式为
多维积分点集为
权值：
即对每一个q值，将所有可以组合的单维积分点任意组合，得到n维积分点.
对比GHF，SGF仅对阶数和不超过2L-1的项进行匹配.随着状态矢量维度的增加，GHF所需要的积分点数呈指数增加，而SGF使得积分点数呈多项式增加，大大减少了积分点数.
2.3 基于交互式多模型的高斯厄米特PHD滤波
在多目标跟踪过程中，针对目标机动采用交互式多模型算法(IMM)实现对目标运动模型的有效识别，首先通过IMM输入交互过程完成模型交互，再通过Sparse Gauss Hermite PHD(以下简称SGH-PHD)滤波器实现不同运动模型下多目标的跟踪滤波，同时计算每个目标的模型似然和模型概率，实现目标的估计融合，最后对高斯分量集实行剪枝和合并，得出目标的状态和数目估计. IMM-SGH-PHD算法流程如下
初始PHD选择下述混合高斯形式：
|
(19) |
预测 假设已知先验PHD为
|
(20) |
由式(1)、式(2)、式(17)可知k+1时刻预测PHD为
|
(21) |
其中，bk+1|k(x)表示新生目标PHD，且有
|
(22) |
DS，k+1|k(x)和Dβ，k+1|k(x)表示目标存活集和衍生集的预测PHD，由于目标运动模式不确定，需要采用IMM算法进行模型交互，由于目标衍生集预测机理与目标存活集相同，此处不再说明
|
(23) |
其中，μk|ki，pq为模型混合概率，πpq为模型概率转移矩阵，
，对于q=1，2，…，r存活目标状态从新初始化的均值、协方差和权值为
|
(24) |
|
(25) |
|
(26) |
利用精度L=3的稀疏高斯厄米特滤波对混合估计得到的均值和协方差进行预测
|
(27) |
|
(28) |
|
(29) |
|
(30) |
|
(31) |
|
(32) |
|
(33) |
|
(34) |
校正 预测的高斯分量个数为nk+1|k，k+1观测集的势为|Zk+1|=mk+1，更新的PHD可表示为
|
(35) |
其中，根据量测信息更新高斯分量为
|
(36) |
|
(37) |
|
(38) |
计算模型似然和模型概率，并进行估计融合
|
(39) |
|
(40) |
|
(41) |
|
(42) |
|
(43) |
进行权值归一化处理
|
(44) |
枝剪与合并 假设当前的高斯分量集为(wk|ki，xk|ki，Pk|ki)，i=1，…，nk+1.首先剔除权值小于截断阈值S的分量，对分布距离小于合并阈值U的分量进行合并，Clark、Panta和Vo[21]提出的合并方法如下，他们采用下式定义的合并距离：di，j=(xk|ki-xk|kj)T(Pk|ki)-1(xk|ki-xk|kj)，如果该距离足够小，则将两个分量(wk|ki，xk|ki，Pk|ki)与(wk|kj，xk|kj，Pk|kj)合并为一个单个分量(wi，j，xi，j，Pi，j)，且wi，j=wk|ki+wk|kj，
，
+(xk|ki-xk|kj)(xk|ki-xk|kj)T.
经过对高斯分量的剪枝和合并后，目标数目可估计为所有分量的权值之和：
，目标状态一般取权值大于0.5的高斯分量的均值.
3 仿真结果及分析
考虑一个二维雷达多机动目标跟踪，目标运动方程和量测方程为
其中，
包含目标的位置和速度，目标匀速运动和匀速转弯运动的状态模型分别为F1= 
，
，vθk~N(0，σθ2)，vrk~N(0，σr2)，取σw=0.2 m，σθ=0.1 rad，σr=0.3 m，角速度ω=0.3 rad/s.目标存活概率PS=0.98，目标检测概率PD=0.99，剪枝阈值S=1e-5，合并阈值U=5，IMM算法模型初始概率为
，概率转移矩阵为
，采样周期T=1 s，共采取100个时刻，杂波在观测空间上服从均匀分布，数目服从均值λ=5的泊松分布.目标真实运动参数如表 1所示，其中目标3为目标2在50 s时的衍生目标.
表 1 目标真实运动参数
Tab.1 Real motion parameters of targets
| 目标 |
起始时刻/s |
消失时刻/s |
目标初始状态 |
| 1 |
1 |
60 |
(-40 m, 1.5 m/s, 40 m, -1.5 m/s) |
| 2 |
22 |
80 |
(-40 m, 1.5 m/s, 60 m, -1.5 m/s) |
| 3 |
50 |
80 |
目标2在50 s时的衍生 |
| 4 |
38 |
100 |
(60 m, 1.5 m/s, -40 m, -1.5 m/s) |
| 5 |
50 |
100 |
(-40 m, 1.5 m/s, -40 m, -1.5 m/s) |
仿真选取最优子模式分配(optimal subpattern assignment，OSPA)脱靶距离[22]为算法估计精度的评价指标，其中参数p=2，c=8.为了验证本算法的稳定性和有效性，仿真给出了IMM-SGH-GMPHD跟IMM-UK-GMPHD和IMM-EK-GMPHD的对比结果，其中IMM-UK-GMPHD和IMM-EK-GMPHD算法框架与IMM-SGH-GMPHD相同. 图 1为目标真实轨迹和观测数据，图 2~图 4为各种算法滤波轨迹和真实估计的比较图.
分析图 2~图 4，从滤波精度来讲3种算法对目标跟踪都出现了较小偏差，但都能一直跟上目标，没有出现发散的情况.从算法稳定性来讲，IMM-EK-GMPHD和IMM-UK-GMPHD出现了目标丢失的情况，IMM-SGH-GMMPHD滤波过程比较稳定，未出现目标丢失. 图 5为2种算法目标数目估计对比.
从图 5可以看出：除个别时间点上杂波观测太接近目标导致3种算法数目估计都出现错误外，IMM-SGHGMPHD对目标数目估计更加稳定，更接近真实目标数.
图 6为60次蒙特卡洛仿真的平均OSPA距离，表 2为3种算法的平均运行时间和OSPA距离比较. OSPA距离代表状态矢量和目标个数估计的综合误差，OSPA距离越小，算法效果越好.可以看出，IMM-SGH-GMPHD虽然在滤波时间略高于其它两种算法，但却提供了更高的跟踪精度和数目估计准确度.
表 2 算法运行时间和OSPA距离和比较
Tab.2 Comparison of algorithm time and OSPA total distance
| 滤波算法 |
平均运行时间/s |
平均OSPA距离和/m |
| IMM-EK-GMPHD |
2.342 6 |
55.520 1 |
| IMM-UK-GMPHD |
1.964 5 |
46.143 7 |
| IMM-SGH-GMPHD |
2.923 1 |
27.017 3 |
4 结论
滤波算法的设计和优化是目标定位与跟踪的基础.本文针对多机动目标跟踪精度低，目标跟踪丢失等问题提出一种基于交互式多模型的稀疏高斯厄米特概率假设滤波器.经典的概率假设滤波器在抗杂波和跟踪目标数目未知或时变方面表现出了优异的性能，但缺乏处理非线性系统和目标机动的能力.本文通过稀疏高斯厄米特滤波将GM-PHD滤波器扩展到非线性系统，并有效地结合IMM算法，实现了运动模型的判断，实现目标状态的估计融合. 60次蒙特卡洛仿真结果表明相比于现有算法，本文提出的算法以合理的计算成本提供了更高的状态估计精度和数目估计准确度.
