近几年来，随着雷达分辨率的不断提高，针对扩展目标跟踪问题的研究已经引起国内外学者的广泛关注^[1-8]。2009年，Mahler推导出了随机有限集框架下扩展目标概率假设密度(probability hypothesis density, PHD)滤波器的递推方程^[3]。随后，Granstrom等学者在线性高斯条件的假设下，给出了扩展目标PHD滤波器的高斯混合(Gaussian mixture, GM)实现方式^[4]。2012年，Granstrom等进一步对扩展目标GMPHD滤波算法作了详细论述，并用实测数据对算法进行了验证^[5]。然而，该算法存在2个问题:①随着扩展目标个数的增多，量测划分数也随之增多，运算量大幅上升；②在扩展目标发生交叉运动的时刻，目标数目易产生漏估。针对算法所存在的问题，Zhang等提出一种基于快速模糊自适应谐振理论(adaptive resonance theory, ART)的扩展目标量测集划分方法^[6]，实现了对量测集的快速划分，提高了算法的实时性，然而，在目标较为密集和杂波条件下，模糊ART算法易出现“饱和”问题导致分类错误^[7]。孔云波等提出一种基于网格密度分布和谱聚类的扩展目标量测集划分算法^[7]，与传统的距离划分相比，运算量降低了38%。刘风梅等采用mean shift对量测集进行划分，使算法的运算量进一步降低^[8]。鉴于mean shift算法在量测划分中表现出的优良性能，本文采用该算法对扩展目标量测集进行划分。但以上算法均未能解决扩展目标交叉时刻的漏估问题。

针对上述问题，本文提出一种基于mean shift和图结构(graph structure, GS)的扩展目标GMPHD滤波算法。采用mean shift算法对扩展目标的量测集进行划分，并依据图结构更新后反馈回的信息判断是否需要对量测集进行子划分；使用图理论对扩展目标间的关系进行建模，搭建图结构，并使用滤波值的一步预测更新图结构，反馈更新后的图结构信息至量测划分步用以指导下一时刻的量测划分。仿真实验验证了算法的有效性和准确性。

1 mean shift算法

mean shift算法本质上是一种基于密度梯度估计的非参数方法，被广泛应用于视频跟踪、图像分割和聚类分析等领域。本文使用mean shift算法对扩展目标量测集进行划分，降低算法的运算量。

mean shift算法是从量测集出发对密度函数进行估计，然后从密度函数梯度的非参数估计中获得。其中，核密度估计是最常用的密度估计方法，它依据预先选定的核函数κ(·)对扩展目标量测集中每个量测的密度函数进行计算。

假定在k时刻从雷达屏幕上所获得的扩展目标的量测集为Z_k, 选取的核函数用κ(·)表示, 则量测z_k的核密度估计值可计算如下:

(1)

本文采用较为常用的高斯核函数^[9]:

(2)

式中, d表示量测维度, 本文取2; , 表示核函数的带宽^[10], α≥1。核函数κ(·)主要用于度量每个样本点z_kⁱ与核中心点z_k之间的相似性, 带宽h则决定了核函数所能够影响的范围。经推导可得最终的mean shift迭代公式为：

(3)

mean shift算法具体步骤如下:

① 计算k时刻扩展目标量测集中每个量测的核密度估计值。

② 获取收敛点。首先, 从给定的量测集中选择任意一个量测, 沿核密度估计梯度∇的方向不断进行搜索, 直到其收敛于局部密度极大值点; 然后, 将沿途的所有访问过的量测划分到对应的类别中; 最后, 在没有访问过的量测中重新选择一个量测作为起始点再次进行递归搜索, 直到量测集中所有的量测都被访问过为止, 便可获得所有的收敛点。

③ 收敛点的合并。计算各收敛点间的马氏距离, 然后将距离小于预设门限的收敛点合并, 同时将属于被合并收敛点的量测集也合并为一类。

④ 输出最终的划分结果。

2 扩展目标的图结构

当扩展目标发生交叉时, 目标产生的量测相距较近, 此时由于量测划分不准确导致目标数目产生漏估。针对该问题, 本文使用图理论对扩展目标间的关系进行建模, 通过搭建图结构动态地获取目标间的距离关系, 并将图结构信息反馈回量测划分步以判断是否需要进行子划分。

定义图结构如下:图结构由节点集合N_o和边的集合E_d构成, 记为G=(N_o, E_d), 其中, N_o表示节点的有限非空集合, E_d表示边的集合。将每个目标视为图结构中的一个节点, 用目标的状态表示节点的状态。同时本文使用连通图的概念来描述相互邻近的扩展目标, 当节点相距较近时, 用边将之连接, 构成连通图。这样, 当不同节点可以构成连通图时, 表示节点对应的目标相距较近, 否则, 对应目标相距较远。

记k时刻的图结构为G_k, 表示k时刻节点(目标)的状态估计值, 其中i=1, …, N_k, N_k为k时刻节点个数, 则k时刻由各节点状态估计值构成的向量集可表示为。计算该时刻各节点之间的马氏距离d_k(i, j), 将之与预先选定的门限ε进行比较, 若d_k(i, j)〈ε, 则认为在k时刻时节点与相距较近, 用边将之连接, 构成连通图。

本文使用邻接矩阵A_k表示k时刻的图结构, 其计算方法如下:

(4)

式中, a_k(i, j)的计算方法如下:

(5)

邻接矩阵A_k中为1的位置表示对应的2个节点构成连通图, 对应目标相距较近, 存在交叉或靠近运动。

4个扩展目标T1、T2、T3和T4对应的图结构如图 1所示, 点迹“·”分别表示4个扩展目标节点, 箭头表示不同节点的运动方向, 椭圆表示预先选定的门限ε, 不同节点之间的连线表示构成连通图的边。该图结构对应的邻接矩阵为:

3 基于mean shift和图结构的GMPHD扩展目标跟踪 3.1 扩展目标运动模型和量测模型

假设扩展目标都是线性运动的, 且不同扩展目标之间的运动状态相互独立, 其状态方程和量测方程分别为:

(6)

(7)

式中, x_k表示扩展目标在k时刻的运动状态, F_k为其对应的状态转移矩阵; z_k表示扩展目标在该时刻产生的量测, H_k为其对应的量测矩阵。ω_k表示过程噪声, 是一个均值为0, 协方差矩阵为Q_k的高斯白噪声; v_k表示量测噪声, 也是一个均值为0的高斯白噪声, 其协方差矩阵用R_k表示。

3.2 算法流程

1) 初始化

根据所给目标初始状态集X₀, 按公式(4)与公式(5)计算得到初始化邻接矩阵A₀, 进而得到初始化的图结构G₀。

2) 量测划分

使用高斯核函数计算量测集中每个量测点的核密度估计值, 将其中小于所设门限的量测点视为杂波，直接剔除。然后，对扩展目标量测集进行mean shift划分，并根据图结构更新后反馈回的目标之间的关系信息判断是否需要进行子划分。如果需要, 则对mean shift聚类所得划分结果的每个子集W_j中的扩展目标个数进行最大似然估计, 如果最大似然估计值, 则使用k-means算法对该集合进行子划分; 否则, 直接进行下一步滤波。

假设各扩展目标产生量测的泊松数为γ, 且相互独立, 同时，子集W_j中杂波可忽略不计。那么，W_j中所包含的扩展目标个数的最大似然估计值_j可计算如下:

(8)

式中，pois(·)表示泊松分布函数，|W_j|表示子集W_j中的量测个数。

3) 扩展目标GMPHD滤波

本文重点在于mean shift聚类和图结构两部分内容, 因此这里仅对扩展目标GMPHD滤波器进行简单介绍, 滤波器的其他内容如预测、高斯项的修剪合并以及扩展目标状态提取等可参见文献[4-5]。假设k时刻扩展目标预测强度函数的高斯混合形式表示为:

(9)

式中, J_k|k-1表示用于拟合预测强度的高斯项个数, w_k|k-1^(j), m_k|k-1^(j), P_k|k-1^(j)分别表示第j个高斯项对应的权值、均值和协方差矩阵。那么, k时刻更新后的扩展目标强度函数可计算如下:

(10)

式中, (1-p_d(x))D_k|k-1(x)表示对漏检部分的更新，表示对检测到部分的更新。J_k|k表示用于拟合更新后扩展目标强度的高斯项个数, w_k|k^(j), m_k|k^(j)与P_k|k^(j)分别代表更新后的第j个高斯项的权值、均值和协方差矩阵。

4) 图结构更新

依据状态方程对滤波后的扩展目标状态进行一步预测, 得到下一时刻各扩展目标的大致位置, 并根据公式(4)与公式(5)计算下一时刻的邻接矩阵, 更新图结构。最后, 将代表下一时刻图结构信息的邻接矩阵反馈回量测划分步以指导下一时刻扩展目标的量测划分。

4 仿真实验与分析

为验证本文所提基于mean shift和图结构的扩展目标GMPHD滤波算法的性能, 本文通过仿真对比所提算法(MSGS-GMPHD)与基于mean shift的扩展目标GMPHD算法(MS-GMPHD)以及基于距离划分的扩展目标GMPHD算法(DP-GMPHD)在同一仿真场景中的跟踪效果来验证所提算法的性能。matlab仿真所使用的PC机硬件平台是3.00GHz Inter(R) Pentium(R), 内存4.00 GB。

观测区域大小为[-500, 500]m×[-500, 500]m, 考虑存在交叉情况的4个扩展目标, 整个过程持续40 s。目标存活时间为:目标1, 1~40时刻; 目标2, 1~40时刻; 目标3, 25~40时刻; 目标4, 33~40时刻。将扩展目标运动模型建模为CV模型, 采样时间间隔为T_s=1 s, 各目标产生量测的泊松率为γ(x)=15, 目标存活概率为p_s=0.99, 检测概率为p_d=0.99。杂波平均数为5, 均匀分布在观测区域内。本文采用OSPA距离^[11]作为算法性能的评价指标, 其参数设置为p=2, c=200, 高斯项的修剪门限T=1×10^-5, 合并门限U=4, 最大高斯项数设置为J_max=100。过程噪声标准差设置为σ_w=2, 量测噪声标准差为σ_v=10, 状态转移矩阵F_k和量测矩阵H_k分别为:

新生目标个数为J_b=4, 初始权重为w_{b, k}^(j)=0.1, 方差为P_{b, k}^(j)=diag(10, 10, 5, 5), 各目标初始状态分别为:

单次蒙特卡罗仿真的结果如图 2所示, 从中可以看出, 3种算法均可实现对扩展目标的稳定跟踪, 但在目标交叉时刻, MS-GMPHD与DP-GMPHD产生了漏估, 而本文所提MSGS-GMPHD算法可以给出准确的估计。

经100次蒙特卡罗仿真后的平均目标数估计如图 3所示, 图 4为其对应的OSPA距离, 可以看出, 在目标交叉时刻, 本文所提算法可以给出更为准确的估计, 而MS-GMPHD与DP-GMPHD均存在漏估问题。

运行时间上, 100次蒙特卡罗仿真的平均划分数如图 5所示, 各时刻的平均运行时间如图 6所示。可以看出, 传统的DP-GMPHD算法随目标数增多, 划分数在不断增多, 运算量也随之增大, 运算时间也在不断增加。由于mean shift能够直接给出扩展目标量测集的正确划分, 因而总是只有一种划分结果, 运算量大幅下降。为了更清楚地看出同样基于mean shift的2种算法中量测划分数和运行时间方面所表现出的异同, 将图 5和图 6中DP-GMPHD算法结果去除, 结果如图 7所示, 2种算法划分数相同, 而本文所提算法MSGS-GMPHD由于在交叉时刻要进行子划分, 所以在交叉时刻运算时间多于MS-GMPHD, 但本文算法解决了交叉时刻扩展目标数目的漏估问题。

下面给出3种算法运行情况的定量比较, 结果如表 1中。表中给出了3种算法单次蒙特卡罗仿真的平均量测划分数和平均运行时间。

从表 1可以看出，相比于传统的DP-GMPHD扩展目标跟踪算法，MS-GMPHD和MSGS-GMPHD算法的量测划分数大大减少，运算效率大幅提升。其中，MS-GMPHD算法的运算时间降低为传统的DP-GMPHD滤波算法的5.68%，而本文所提MSGS-GMPHD算法则降低为DP-GMPHD算法的6.00%。本文所提MSGS-GMPHD算法的运行时间略高于MS-GMPHD算法，但是本文算法解决了扩展目标交叉时刻的漏估问题，以极小的时间代价获得了更准确的估计结果。

5 结论

本文针对传统的基于距离划分的扩展目标GMPHD滤波算法量测划分数多、运算量大的问题以及交叉时刻产生的目标数漏估问题，提出一种基于mean shift和图结构的扩展目标GMPHD滤波算法。采用mean shift聚类代替传统的距离划分，有效地降低了算法的运算量。使用图理论对不同扩展目标间的关系进行建模，实时动态地获取目标间的距离关系，并将该信息反馈回量测划分步以判断是否进行子划分，图结构的引入很好地解决了目标交叉时刻的漏估问题。仿真实验证明，本文所提MSGS-GMPHD算法拥有比传统算法更好的估计性能。