CN110376581A - 基于高斯混合概率假设密度滤波器的显式多目标跟踪方法 - Google Patents

Abstract

本发明公开了一种基于高斯混合概率假设密度滤波器的显式多目标跟踪方法，其核心技术在于利用查寻矩阵，基于更新高斯分量的标签和权值屏蔽目标高先验密度区域杂波的干扰，实现各目标的后验信息、量测及其状态估计的同标签化，从而得到无需额外关联程序的显式多目标跟踪。并且，保持与基本的GM‑PHD滤波器相似的实时性。工程上易于实现，且很好地改善多目标跟踪的精度。

Description

基于高斯混合概率假设密度滤波器的显式多目标跟踪方法

技术领域

本发明涉及一种基于高斯混合概率假设密度滤波器的显式多目标跟踪方法，属于信号处理技术领域。

背景技术

多目标跟踪是从一系列量测中估计出目标的数量和状态，其广泛应用于军事和民用领域，例如，雷达多目标跟踪。随着人们对雷达功能的需求不断提高，其应用场景日趋复杂。如信噪比、低信杂波以及目标密集等现象，严重影响雷达的性能，增加了虚警概率并导致目标的检测概率降低。这些将直接影响目标状态－量测关联的准确性，从而降低状态估计的精度。关于屏蔽多目标跟踪中各类干扰、实现精确的多目标状态估计及航迹的相关课题，已成为国内外本领域研究的热点。目前，应用最多的是联合概率数据关联、多假设关联和随机有限集(RFS)。

由于RFS避免了传统的数据关联，因此在多目标跟踪应用中得到了广泛的关注。在贝叶斯滤波框架下基于RFS，已提出了多种贝叶斯多目标滤波器的近似实现，主要包括概率假设密度(PHD)滤波器、势化概率假设密度滤波器和多伯努力滤波器。虽然这些滤波器不能提供多目标的航迹，但仍得到了广泛的应用。近些年，提出了一些贴标签的、基于RFS的标签化多伯努力(GLMB)滤波器，不仅可以显著改善多目标状态提取的精度，而且提供显式的多目标航迹。然而，此类滤波器融合了数据关联技术，使得计算负担增加。它们难以应用于实时性要求高的场景。因此，有必须设计一种计算量相对较小、且可提供显式的多目标航迹的滤波器。

不同于全贝叶斯滤波器，PHD滤波器只迭代多目标后验概率的一阶矩。从实际应用角度来讲，PHD滤波器特别适用于实时性要求高的场景。但其不提供显式的多目标航迹，而且对多目标数量的估计不够稳定。近年来，国际上对此问题的研究取得了很大的进展。基于PHD滤波器的序贯蒙特卡罗(SMC)近似实现，已有一种显式的多目标跟踪算法提出。基于PHD滤波器的高斯混合(GM)近似实现，杂波密集环境下的多目标状态提取的精度已得到很好的提高。然而，当杂波密集、目标紧邻或漏检时，仍难以得到精确的显式多目标航迹。与SMC-PHD滤波器相比，GM-PHD滤波器具有状态提取简单、计算量低的特点。因此，基于GM-PHD去设计可提供显式多目标航迹的滤波器是首选。

发明内容

基于GM-PHD滤波器的多目标跟踪技术尚未解决复杂环境下的多目标航迹管理，针对这个技术问题，本发明的目的在于提出一种基于GM-PHD滤波器的显式多目标跟踪算法，其核心技术在于利用查寻矩阵，基于更新高斯分量的标签和权值屏蔽目标高先验密度区域杂波的干扰，实现各目标的后验信息、量测及其状态估计的同标签化，从而得到无需额外关联程序的显式多目标跟踪。并且，保持与基本的GM-PHD滤波器相似的实时性。工程上易于实现，且很好地改善多目标跟踪的精度。

本发明为解决上述技术问题采用以下技术方案：

本发明提供一种基于高斯混合概率假设密度滤波器的显式多目标跟踪方法，基于多目标先验信息已知的假设，提出适用于高斯混合概率假设密度GM-PHD滤波器的两个“一对一”原则：某簇相同标签的存活高斯分量只对应于一个量测，此量测也只对应此簇存活高斯分量。

该方法包括以下步骤：

步骤1、高斯分量初始化

初始化k＝0时刻多目标的后验强度D_0|0(x)，得到近似多目标后验强度D_0|0(x)的J_γ,k个高斯分量的参数集其中，x为单目标的状态，J_γ,k为先验的新生目标位置数，表示高斯密度，分别为第i个新生目标分量的权值、均值和协方差矩阵，为第i个新生目标分量的标签，V₁是一个远大于场景中可能出现的所有目标总数的值；

步骤2、高斯分量预测

已知k-1时刻得到的近似多目标后验强度的高斯分量的参数集对k时刻的多目标信息进行预测，得到k时刻的多目标先验强度得到近似k时刻的多目标先验强度D_k|k-1(x)的J_k|k-1个高斯分量的参数集其中，J_k-1是k-1时刻近似多目标后验强度的高斯分量数目，和分别是k-1时刻第j5个分量的权值、均值、协方差矩阵和标签，J_k|k-1＝J_k-1+J_γ,k，和分别是第j个预测分量的权值、均值、协方差矩阵和标签；

步骤3、剔除杂波

根据雷达获得的k时刻的量测集以及步骤2得到的高斯分量参数集基于适用于GM-PHD滤波器的波门法剔除杂波，得到k时刻的有效量测集其中，是雷达获得的第ze1个量测，是采用波门法后得到的第ze个有效量测，|·|为集合的势；

步骤4、更新高斯分量

根据获得的Z_k,ef和更新k时刻的多目标后验强度p_D,k为k时刻单目标的检测概率，(1-p_D,k)D_k|k-1(x)由J_k|k-1个更新高斯分量的参数集近似；其中，和分别为根据有效量测z_k,ef和第j个预测分量得到的更新高斯分量的权值、均值、协方差矩阵和标签，λ为k时刻雷达获得的平均杂波数，c(z_k)为雷达获得的任意量测z_k的空间分布的概率密度，H_k为量测矩阵，R_k为量测噪声协方差矩阵； I为单位矩阵，得到近似k时刻的多目标后验强度D_k|k(x)的J_k|k-1(1+|Z_k,ef|)个更新高斯分量的参数集和分别是第j1个更新高斯分量的权值、均值、协方差矩阵和标签；

步骤5、多目标状态提取

从更新的高斯分量中提取多个带身份标识的单目标状态估计，具体包括：

(5a)设四个查询矩阵U_w、U_m、U_P、U_l为J_k|k-1(1+|Z_k,ef|)维的矩阵，步骤4中得到的J_k|k-1(1+|Z_k,ef|)个更新高斯分量的权值、均值、协方差矩阵和标签分别存储在U_w、U_m、U_P和U_l中，其中，中的元素依次存储至U_w的第一列，中的元素依次存储至U_m的第一列，中的元素依次存储至U_P的第一列，中的元素依次存储至U_l的第一列；按ze＝1,...,|Z_k,ef|，依次存储由更新的高斯分量的参数至U_w、U_m、U_P和U_l中，具体为：中的元素依次存储至U_w的第ze+1列，中的元素依次存储至U_m的第ze+1列，中的元素依次存储至U_P的第ze+1列，中的元素依次存储至U_l的第ze+1列；

(5b)令

(5c)设k时刻的状态估计集航迹标签集发生时刻集有效量测在U_w中的列号集

(5d)取U_w中第二列至最后一列的所有元素中的最大值，如果此最大值小于状态提取的基本权值阈值w₁，则直接跳至步骤(5h)。否则，根据此最大值所在的不同位置，从U_l中取这些不同位置处的标签，接着对这些标签进行处理，去掉其中重复的元素，得到不重复的标签集

(5e)按依次对所有标签为的高斯分量进行处理，具体包括：

(5e.1)根据标签在U_l中的行号，从U_w的相应行中取最大值，记这个最大值的列号为记这个最大值在U_m中相应位置处的元素为m⁽ⁱ¹⁾，U_l中所有标签为的行序号集记作

(5e.2)判断标签的属性，根据其属性对其进行相应的操作，存在以下三种情况：

A.为新生分量的标签，具体操作步骤包括：

①从矩阵U_w中所有列在行中的元素中取大于或等于w₂的元素所在列的不重复的列序号构成集合同时将这些列序号存至CI中，其中，w₂为大权值分量的阈值；

②若依次按执行以下操作：

首先，记U_w中第列所在行的元素中最大值所在行的行号为提取U_m的第列在第行中的元素为单目标状态估计，将其存储至状态估计集中，并建立新的确定航迹r_max＝r_max+1，且将r_max和k分别存储至航迹标签集和发生时刻集中；

接着，分别对和进行修改，为第列元素所依据的大权值量测标出其后验信息，具体操作为：

其中，记中行在所有列中的元素为记中第行在第列中的元素为记中行在所有列中的元素为记的最后一行下面添加的那一行元素为记的第列在最后一行的元素为记U_w中第行在第列中的元素为记的最后一行下面添加的那一行元素为记U_m中第行在所有列中的元素为记的最后一行下面添加的那一行元素为记U_P中第行在所有列中的元素为记的最后一行下面添加的那一行元素为

③取中权值大于或等于w₁、且小于w₂的分量所在列的不重复的列序号集，记为从中剔除包含的元素，得到同时将这些列序号存至CI中；

④若依次按执行以下操作：

首先，记U_w中第列所在行的元素中最大值所在行的行号为建立新的暂态航迹r_maxt＝r_maxt+1；

接着，分别对和进行修改，为第列元素所依据的小权值量测标出其后验信息，具体操作为：

其中，记的第行在第列的元素为记U_w的第行在第列的元素为记的最后一行元素在第列的元素为记U_m和U_P的第行在所有列的元素分别为和

B.为确定分量的标签，从此类分量中提取单目标状态估计，具体步骤包括：

①如果对应的最大值小于状态预处理的权值阈值w₃，记U_w的第1列在行中的元素中最大值所在的行号为依次执行以下状态估计预处理步骤：

其中，||b1-b2||表示b1与b2之间的欧氏距离；U_w和U_m的第行在第1列的元素分别记为和w₃满足w₁≤w₃＜w₂；已知零均值的量测噪声协方差矩阵 和分别为x轴和y轴的量测方差，令σ＝[σ_x σ_y]^T，d(1)＝||σ-0||；

②将m⁽ⁱ¹⁾作为单目标状态估计，将其存储至状态估计集中，并将和k分别存储至航迹标签集和发生时刻集中，同时，将存储至CI；

③高斯分量权值校正

先逐一计算m⁽ⁱ¹⁾至同标签分量之间的欧氏距离：

再根据条件校正权值：

其中，U_w和U_m的第行在第ze1列的元素分别记为和

C.为暂态分量的标签，具体步骤包括：

①当时，或当且时，取m⁽ⁱ¹⁾为单目标状态估计，将其和k分别存储至状态估计集和发生时刻集中；建立新的确定航迹r_max＝r_max+1，将r_max存储至航迹标签集修改和中元素，置以及将存储至其中，为k-1时刻的小权值暂态分量集，为k-2时刻的小权值暂态分量集；

②否则，将存储在k时刻的小权值暂态分量集中；

(5e.3)根据两个“一对一”原则，将U_w中标签为的所有元素置零；

(5f)执行完循环之后，以CI的元素为列号集，根据两个“一对一”原则，将这些列在U_w的任意行的元素全部置零；

(5g)返回至第(5d)步；

(5h)记为和矩阵的行数，分别将和中的个元素逐列转换为维的数组，得到修正后的更新高斯分量的参数集

步骤6、高斯分量剪枝，具体包括：

(6a)剪去连续n_nost次没有状态估计的确定分量或暂态分量，具体步骤包括：

(6a.1)取中不重复的标签集，从中删除和中包含的标签，得到k时刻未获得状态估计的确定分量和暂态分量的标签集

(6a.2)按依次判别在前n_nost-1个时刻未获得状态估计的确定分量和暂态分量的标签集中的出现情况：如果在k1＝k-n_nost+1,...,k-1时刻的每个中都出现过，则置中所有标签为的分量权值为零；

(6b)将所有权值大于小权值阈值w_et的高斯分量的序号记录在I1中；

(6c)令依此类推，得到其中，和分别是第j3个高斯分量的权值、均值、协方差矩阵和标签；

步骤7、基于步骤6剪枝后得到高斯分量的参数集进行高斯分量合并，令为高斯分量的序号集，且J_k＝0，具体步骤包括：

(7a)当I2为空时，跳至步骤(7e)，否则，J_k＝J_k+1，得到标签

(7b)取与同标签、且可合并的分量的序号集为I3，其中，d_Th为合并的距离阈值；

(7c)合并序号集I3中的所有高斯分量参数， 得到一个新的高斯分量；

(7d)I2＝I2-I3，并返回步骤(7a)；

(7e)得到

步骤8、将不同时刻具有相同身份标识的状态估计相连，即得到显式的多目标航迹；

步骤9、k＝k+1，返回步骤2。

进一步地，步骤3中的波门法的具体步骤为：

(3a)已知零均值的量测噪声协方差矩阵 和分别为x轴和y轴的量测方差，令σ＝[σ_x σ_y]^T，以此为依据计算波门阈值d(a)＝a||σ-0||。其中，a为置信系数。

(3b)设有效量测集

(3c)为k时刻获得的最新量测集，按ze1＝1,...,|Z_k|，依次判断每个的属性。判断时，按j＝1,...,J_k|k-1，逐一计算至每个高斯分量均值的距离如果d^(ze1,j)≤d(a)，则判定附近存在高斯分量，视为有效量测，将其加入到有效量测集Z_k,ef中，接着跳出此循环，执行对的判断。

进一步地，若标签l≤V₂，则为确定分量的标签；若V₂＜l≤V₁，则为暂态分量的标签；若l＞V₁，则为新生分量的标签。

进一步地，k＝0时刻，确定分量和暂态分量的标签分别以r_max＝0和r_maxt＝V₂为基数，每新建一个确定航迹，r_max加1；每新建一个暂态航迹，r_maxt加1。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：

(1)将高斯分量的标签分为三类，并且设定确定高斯分量的标签在整个滤波过程中不变。为实现显式的多目标航迹提供基础；

(2)建立查寻矩阵来存储更新高斯分量的四类参数，基于高斯分量的权值和标签，将各目标的后验信息、量测及其状态估计同标签化，不仅屏蔽了近距离杂波的干扰、提高了多目标状态提取的精度，而且得到了无需额外关联程序的显式多目标跟踪；

(3)只有相同标签的高斯分量才可以合并，可避免因漏检或量测噪声导致的权值相对小的目标的高斯分量被较大权值的目标的高斯分量合并，解决近距离目标相互间的干扰问题。

(4)利用波门技术，得到有效量测；基于有效量测，更新高斯分量；在同标签化的过程中，没有额外的计算处理。这些操作保证了本发明具有与基本的GM-PHD滤波器相似的实时性。

附图说明

图1是本发明的方法结构图；

图2是本发明用于剔除杂波的波门算法流程图；

图3是本发明查寻矩阵结构图；

图4是本发明目标后验信息的有效区域示意图；

图5(a)是杂波数为100、检测概率为0.9时单次实验得到的各目标在x轴和y轴的航迹；

图5(b)是杂波数为100、检测概率为0.9时单次实验得到的各目标随时间变化的x轴的航迹；

图5(c)是杂波数为100、检测概率为0.9时单次实验得到的各目标随时间变化的y轴的航迹；

图6(a)是杂波数为100、检测概率为0.9时100次实验得到的各滤波器的平均OSPA距离；

图6(b)是杂波数为100、检测概率为0.9时100次实验得到的各滤波器的平均计算时间；

图7(a)是检测概率为0.95、杂波数从1至100时100次实验得到的各滤波器的平均OSPA距离；

图7(b)是检测概率为0.95、杂波数从1至100时100次实验得到的各滤波器的平均计算时间；

图8(a)是杂波数为100、检测概率从0.8至1时100次实验得到的各滤波器的平均OSPA距离；

图8(b)是杂波数为100、检测概率从0.8至1时100次实验得到的各滤波器的平均计算时间。

具体实施方式

下面结合附图对本发明的技术方案做进一步的详细说明：

基于GM-PHD滤波器的多目标跟踪技术尚未解决复杂环境下的多目标航迹管理，针对这个技术问题，本发明的目的在于提出一种基于GM-PHD滤波器的显式多目标跟踪算法，其核心技术在于利用查寻矩阵，基于更新高斯分量的标签和权值屏蔽目标高先验密度区域杂波的干扰，实现各目标的后验信息、量测及其状态估计的同标签化，从而得到无需额外关联程序的显式多目标跟踪。并且，保持与基本的GM-PHD滤波器相似的实时性。本发明工程上易于实现，且很好地改善多目标跟踪的精度。

为达到上述技术目的，本发明采用以下技术方案予以实现。

本发明假设多目标先验信息已知，提出适用于GM-PHD滤波器的两个“一对一”原则：某簇相同标签的存活高斯分量只对应于一个量测，此量测也只对应此簇存活高斯分量。

基于提出的这两个“一对一”原则，本发明设计了一种基于GM-PHD滤波器的显式多目标跟踪算法，如图1所示，包括以下步骤：

步骤1、高斯分量初始化

初始化k＝0时刻多目标的后验强度D_0|0(x)，得到近似多目标后验强度的高斯分量，具体包括：

(1a)已知k时刻新生目标的先验强度表示为其中，x为单目标的状态，J_γ,k为先验的新生目标位置数，表示高斯密度，分别为第i个新生目标分量的权值、均值和协方差矩阵，为第i个新生目标分量的标签，V₁是一个远大于场景中可能出现的所有目标总数的值。本发明中特取V₁＝400。k＝0时刻的多目标后验强度为由高斯分量的参数集近似。

(1b)对高斯分量l的标签进行分类，分为三类。若l≤V₂，为确定分量的标签；若V₂＜l≤V₁，为暂态分量的标签；若l＞V₁，新生分量的标签。滤波初始步，确定分量和暂态分量的标签分别以r_max＝0和r_maxt＝V₂为基数，本发明中特取V₂＝200。每新建一个确定航迹，r_max加1；每新建一个暂态航迹，r_maxt加1。

步骤2、高斯分量预测

已知k-1时刻得到的近似多目标后验强度的高斯分量参数集对k时刻的多目标信息进行预测，得到k时刻多目标的先验强度由J_k|k-1个高斯分量的参数集近似。其中，J_k-1是k-1时刻近似多目标后验强度的高斯分量数目，和分别是k-1时刻第j5个分量的权值、均值、协方差矩阵和标签。J_k|k-1＝J_k-1+J_γ,k，和分别是第j个预测分量的权值、均值、协方差矩阵和标签。具体步骤为：令j＝1对应j5＝1， 依此类推，直至j＝J_k-1对应j5＝J_k-1；令j＝J_k-1+1对应i＝1，依此类推，直至j＝J_k|k-1对应i＝J_γ,k。其中，F_k-1和Q_k-1分别为单目标状态转移矩阵和过程噪声协方差矩阵。

步骤3、剔除杂波

根据雷达获得的k时刻的最新量测集以及步骤2得到的高斯分量参数集基于适用于GM-PHD滤波器的波门法剔除杂波，得到k时刻的有效量测集其中，是雷达获得的第ze1个量测，是采用波门法后得到的第ze个有效量测，|·|为集合的势。

如图2所示，具体包括：

(3a)已知零均值的量测噪声协方差矩阵 和分别为x轴和y轴的量测方差，令σ＝[σ_x σ_y]^T，以此为依据计算波门阈值d(a)＝a||σ-0||。其中，a为置信系数。

(3b)设有效量测集

(3c)为k时刻获得的最新量测集，按ze1＝1,...,|Z_k|，依次判断每个的属性。判断时，按j＝1,...,J_k|k-1，逐一计算至每个高斯分量均值的距离如果d^(ze1,j)≤d(a)，则判定附近存在高斯分量，视为有效量测，将其加入到有效量测集Z_k,ef中，接着跳出此循环，执行对的判断。

步骤4、更新高斯分量

根据获得的有效量测集Z_k,ef和k时刻预测的高斯分量集更新k时刻的多目标后验强度D_k|k(x)。已知p_D,k为k时刻单目标的检测概率，具体为：

其中，

其中，(1-p_D,k)D_k|k-1(x)由高斯分量参数集近似。D_D,k(x；z_k,ef)由J_k|k-1个高斯分量的参数集近似， 和分别为根据有效量测z_k,ef和第j个预测高斯分量得到的更新高斯分量的权值、均值、协方差矩阵和标签，称第j个预测高斯分量为此更新高斯分量的父分量。λ为k时刻雷达获得的平均杂波数，c(z_k)为雷达获得的任意量测z_k的空间分布的概率密度，H_k为量测矩阵，R_k为量测噪声协方差矩阵，I为单位矩阵。

据此，k时刻的多目标后验强度由J_k|k-1(1+|Z_k,ef|)个高斯分量的参数集近似。其中，和分别是第j1个更新高斯分量的权值、均值、协方差矩阵和标签。

步骤5、多目标状态提取

从更新的高斯分量中提取多个带身份标识的单目标状态估计，具体包括：

(5a)设四个查询矩阵U_w、U_m、U_P、U_l为J_k|k-1(1+|Z_k,ef|)维的矩阵，步骤4中得到的J_k|k-1(1+|Z_k,ef|)个更新高斯分量的权值、均值、协方差矩阵和标签分别存储在U_w、U_m、U_P和U_l中。图3给出了U_w、U_m、U_P和U_l矩阵的结构图。其中，中的元素依次存储至U_w的第一列，中的元素依次存储至U_m的第一列，中的元素依次存储至U_P的第一列，中的元素依次存储至U_l的第一列。按ze＝1,...,|Z_k,ef|，依次存储由更新的高斯分量的参数至U_w、U_m、U_P和U_l中，具体为：中的元素依次存储至U_w的第ze+1列，中的元素依次存储至U_m的第ze+1列，中的元素依次存储至U_P的第ze+1列，中的元素依次存储至U_l的第ze+1列。显然，U_P中每一行的元素相同，U_l中每一行的元素相同。

将更新的高斯分量存储至U_w、U_m、U_P和U_l中，便于查寻各更新分量的父分量及其对应的量测。

(5b)创建四个临时矩阵和

(5c)设k时刻的状态估计集航迹标签集发生时刻集有效量测在U_w中的列号集

(5d)记矩阵U_w的所有行在第二列至最后一列的元素为从中取最大值，记作U_w中可能有多个当跳至步骤(5h)。否则，当根据在U_w中不同位置处的行号与列号，从U_l中取这些不同位置处的标签，接着对这些标签进行处理，去掉其中重复的元素，得到不重复的标签集其中，w₁为状态提取的基本权值阈值，取w₁＝0.02。

(5e)按依次对所有标签为的高斯分量进行处理，具体包括：

(5e.1)根据标签在U_l中的行号，从U_w的相应行中取最大值，记这个最大值的列号为记这个最大值在U_m中相应位置处的元素为m⁽ⁱ¹⁾，U_l中所有标签为的行序号集记作

(5e.2)判断标签的属性，根据其属性对其进行相应的操作。存在以下三种情况：

A.新生分量的标签。具体步骤包括：

①记矩阵U_w中所有列在行中的元素为从取大于或等于w₂的元素所在列的不重复的列序号集，记为同时将这些列序号存至CI中。具体步骤为：

式中，uni(·)是取集合中不重复元素的函数。w₂为大权值分量的阈值，令w₂＝0.4。

②若依次按执行以下操作：

首先，记U_w中第列所在行的元素为取中最大值所在行的行号为记U_m的第列在第行中的元素为提取为单目标状态估计，将其存储至状态估计集并建立新的确定航迹r_max＝r_max+1，且执行和为相应的状态估计贴上标签r_max和时间序号k。

接着，分别对和进行修改，为第列元素所依据的大权值量测标出其后验信息，具体操作为：

式中，记中行在所有列中的元素为记中第行在第列中的元素为记中行在所有列中的元素为记的最后一行下面添加的那一行元素为记的第列在最后一行的元素为记U_w中第行在第列中的元素为记的最后一行下面添加的那一行元素为记U_m中第行在所有列中的元素为记的最后一行下面添加的那一行元素为记U_P中第行在所有列中的元素为记的最后一行下面添加的那一行元素为

式中，当c1＞1时，即有两个或以上的大权值新生分量存在、且之前已提取过新生目标的状态估计，则在已有的和的最后一行下面分别增加一行，为第列元素所依据的大权值量测标出其后验信息。

③取中权值大于或等于w₁、且小于w₂的分量所在列的不重复的列序号集，记为从中剔除包含的元素，得到同时将这些列序号存至CI中，具体操作为：

④若依次按执行以下操作：

首先，记U_w中第列所在行的元素为取中最大值所在行的行号为建立新的暂态航迹r_maxt＝r_maxt+1。

接着，分别对和进行修改，为第列元素所依据的小权值量测标出其后验信息，具体操作为：

式中，记的第行在第列的元素为记U_w的第行在第列的元素为记的最后一行元素在第列的元素为记U_m和U_P的第行在所有列的元素分别为和

式中，当没有大权值新生分量存在且第一次为小权值量测标后验信息时，即c2＝1时，先置接着修改当没有大权值新生分量存在且c2＞1时，或者时，则分别在和的最后一行下面分别增加一行，为第列元素所依据的小权值量测标出其后验信息。

B.确定分量的标签。从此类分量中提取单目标状态估计。

当目标紧邻处存在杂波，又或量测噪声大导致的量测偏移真实位置，导致此目标后验信息(即高斯分量)的分布发散，使得其量测的权值较小。或者，目标在k-1时刻漏检，导致后验分布的信息量小，使得k时刻其量测的权值较小。因此，需要先对权值较小的高斯分量进行状态估计预处理，即修正状态估计及其权值。具体步骤包括：

①如果记U_w的第1列在行中的元素为取中最大值所在的行号为依次执行以下状态估计预处理步骤：

其中，||b1-b2||表示b1与b2之间的欧氏距离；U_w和U_m的第行在第1列的元素分别记为和w₃满足w₁≤w₃＜w₂；已知零均值的量测噪声协方差矩阵 和分别为x轴和y轴的量测方差，令σ＝[σ_x σ_y]^T，d(1)＝||σ-0||。此举可尽量消除各类干扰产生的偏差。

②取单目标状态估计贴标签记录时间序号同时，将存储至

③高斯分量权值校正。值越大，单目标的后验分布越集中，有效区域半径越小；值越小，单目标的后验分布越分散，有效区域半径应越大。用图4来示例此现象。

为此，先逐一计算m⁽ⁱ¹⁾至同标签分量之间的距离d^(i2,ze1)，

再根据条件校正权值，

其中，U_w和U_m的第行在第ze1列的元素分别记为和在有效区域内，高斯分量的权值不变；区域外的高斯分量权值都削弱至原来的a1倍。本发明中，a1＝1/3。根据权值来划定后验分布的有效区域，在状态估计预测处理的基础上，可进一步屏蔽各类干扰。

C.暂态分量的标签。具体步骤包括：

①当时，或当且(为k-1时刻的小权值暂态分量集)时，取m⁽ⁱ¹⁾为单目标状态估计，并执行和建立新的确定航迹r_max＝r_max+1，状态估计贴标签接着，修改和中元素，置以及并且，将存储至

②否则，将存储在k时刻的小权值暂态分量集中，即

(5e.3)根据两个“一对一”原则，将U_w中标签为的所有元素置零：

(5f)执行完循环之后，以CI的元素为列号集，记这些列在U_w的任意行的元素为根据两个“一对一”原则，将这些列的元素全部置零：

(5g)返回至第(5d)步；

(5h)记为和矩阵的行数，分别将和中的个元素逐列转换为维的数组，得到修正后的更新高斯分量的参数集

步骤6、高斯分量剪枝。具体步骤包括：

(6a)剪去连续n_nost次没有状态估计的确定分量或暂态分量，n_nost＝6。具体步骤包括：

(6a.1)取中不重复的标签集，从中删除和中包含的标签，得到k时刻未获得状态估计的确定分量和暂态分量的标签集

(6a.2)按依次判别在前n_nost-1个时刻未获得状态估计的确定分量和暂态分量的标签集中的出现情况：如果在k1＝k-n_nost+1,...,k-1时刻的每个中都出现过，则执行置高斯分量参数集中所有标签为的分量权值为零。

(6b)根据小权值阈值w_et＝10^-5，将所有权值大于w_et的高斯分量的序号记录在I1中，即剪去小权值分量后的序号集I1为

(6c)令J_k|k＝|I1|，依此类推，得到其中，和分别是第j3个高斯分量的权值、均值、协方差矩阵和标签。

步骤7、基于剪枝后得到高斯分量的参数集进行高斯分量合并。。多目标紧邻或交叉时，相互间的后验信息在合并时容易产生干扰，因漏检或量测噪声导致的权值相对小的目标的高斯分量容易被较大权值的目标的高斯分量合并，导致小权值目标的部分真实信息丢失。因此，只有标签相同的分量才可合并。令为高斯分量的序号集，且J_k＝0，具体步骤包括：

(7a)当I2为空时，跳至步骤(7e)，否则，J_k＝J_k+1，得到标签

(7b)取与同标签、且可合并的分量的序号集为I3，

其中，d_Th为合并的距离阈值，一般d_Th＝4。

(7c)合并序号集I3中的所有高斯分量参数，

从而得到一个新的高斯分量。

(7d)I2＝I2-I3，并返回步骤(7a)。

(7e)得到

步骤8、将不同时刻具有相同身份标识的状态估计相连，即得到显式的多目标航迹。

步骤9、k＝k+1，返回步骤2。

至此，本发明的一种基于高斯混合概率假设密度滤波器的显式多目标跟踪算法结束。

实施例

通过以下仿真实现对本发明效果作进一步验证说明。

1.实验条件：

在一个二维场景[-1000 1000]m×[-1000 1000]m中，每个目标的运动方程为：

其中，x_k＝[x_1,k x_2,k x_3,k x_4,k]^T，[x_1,k x_3,k]^T为k时刻目标的位置，[x_2,k x_4,k]^T为k时刻目标的速度。Δ＝1s为采样间隔，σ_ω＝5m/s²。目标在场景中可以在任意时刻出现或消失，存活概率p_S,k＝0.95。新生目标出现的强度函数为其中，

目标的量测方程为

其中，υ_x,k和υ_y,k为相互独立的零均值高斯白噪声，均方差分别为σ_x＝10m、σ_y＝10m。杂波均匀地分布于监测区域[-1000 1000]m×[-1000 1000]m，每帧的杂波数为λ。基本的GM-PHD滤波器中，杂波强度为κ＝λ/2000²。本发明提出的滤波器与基本的GM-PHD滤波器比较。由于本发明提出的滤波器在高斯分量更新前使用了波门法去杂波，为了公平，基本的GM-PHD滤波器也采用了波门法去杂波。其基于所有预测高斯分量，量测门概率p₀＝0.999。高斯分量的最大数目为J_max＝100，截取阈值w_et＝10^-5，合并阈值d_Th＝4。

OSPA距离用来评估多目标跟踪性能，c＝100，p＝1。基于相同的目标轨迹，运行100次MC实验得到平均性能，每次实验中的量测相互独立。

2.仿真内容

仿真1，杂波数目λ＝100、检测概率p_D,k＝0.9时，图5给出了本发明所提方法在单次实验确认出的状态估计与航迹。图5中，相同灰度的皆属于同一条航迹的。图6(a)是本发明所提方法在100次实验的平均OSPA距离，图6(b)是相应的计算时间。

仿真2，检测概率p_D,k＝0.95保持不变、杂波数目从1增加到100时，图7(a)是本发明所提方法在100次实验的平均OSPA距离，图7(b)是相应的计算时间。

仿真3，杂波数目λ＝100保持不变、检测概率从0.8增加到1时，图8(a)是本发明所提方法在100次实验的平均OSPA距离，图8(b)是相应的计算时间。

3.仿真结果分析：

从图5可以看出，本发明所设计的滤波器可以实现密集杂波下的复杂多目标跟踪，给出正确的多目标航迹。从图6可以看出，本发明所设计的滤波器的跟踪精度较基本的GM-PHD滤波器有明显的提高，实时性却相似；与GLMB滤波器相比，跟踪精度略低，实时性却完胜。

从图7可以看出，在p_D,k＝0.95、杂波数目从1增加到100时，本发明所设计的滤波器的跟踪精度较基本的GM-PHD滤波器有明显的提高，实时性相似。

从图8可以看出，杂波数目λ＝100保持不变、检测概率从0.8增加到1时，本发明所设计的滤波器的跟踪精度较基本的GM-PHD滤波器有明显的提高，实时性相似。

综上所述，仿真实验验证了本发明的正确性、有效性和可靠性。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围；这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (4)

1.基于高斯混合概率假设密度滤波器的显式多目标跟踪方法，其特征在于，基于多目标先验信息已知的假设，提出适用于高斯混合概率假设密度GM-PHD滤波器的两个“一对一”原则：某簇相同标签的存活高斯分量只对应于一个量测，此量测也只对应此簇存活高斯分量；

该方法包括以下步骤：

步骤1、高斯分量初始化

初始化k＝0时刻多目标的后验强度D_0|0(x)，得到近似多目标后验强度D_0|0(x)的J_γ,k个高斯分量的参数集其中，x为单目标的状态，J_γ,k为先验的新生目标位置数，表示高斯密度，分别为第i个新生目标分量的权值、均值和协方差矩阵，为第i个新生目标分量的标签，V₁是一个远大于场景中可能出现的所有目标总数的值；

步骤2、高斯分量预测

已知k-1时刻得到的近似多目标后验强度的高斯分量的参数集对k时刻的多目标信息进行预测，得到k时刻的多目标先验强度得到近似k时刻的多目标先验强度D_k|k-1(x)的J_k|k-1个高斯分量的参数集其中，J_k-1是k-1时刻近似多目标后验强度的高斯分量数目，和分别是k-1时刻第j5个分量的权值、均值、协方差矩阵和标签，J_k|k-1＝J_k-1+J_γ,k，和分别是第j个预测分量的权值、均值、协方差矩阵和标签；

步骤3、剔除杂波

根据雷达获得的k时刻的量测集以及步骤2得到的高斯分量参数集基于适用于GM-PHD滤波器的波门法剔除杂波，得到k时刻的有效量测集其中，是雷达获得的第ze1个量测，是采用波门法后得到的第ze个有效量测，|·|为集合的势；

步骤4、更新高斯分量

根据获得的Z_k,ef和更新k时刻的多目标后验强度p_D,k为k时刻单目标的检测概率，(1-p_D,k)D_k|k-1(x)由J_k|k-1个更新高斯分量的参数集近似；其中，和分别为根据有效量测z_k,ef和第j个预测分量得到的更新高斯分量的权值、均值、协方差矩阵和标签，λ为k时刻雷达获得的平均杂波数，c(z_k)为雷达获得的任意量测z_k的空间分布的概率密度，H_k为量测矩阵，R_k为量测噪声协方差矩阵； I为单位矩阵，得到近似k时刻的多目标后验强度D_k|k(x)的J_k|k-1(1+|Z_k,ef|)个更新高斯分量的参数集和分别是第j1个更新高斯分量的权值、均值、协方差矩阵和标签；

步骤5、多目标状态提取

从更新的高斯分量中提取多个带身份标识的单目标状态估计，具体包括：

(5a)设四个查询矩阵U_w、U_m、U_P、U_l为J_k|k-1(1+|Z_k,ef|)维的矩阵，步骤4中得到的J_k|k-1(1+|Z_k,ef|)个更新高斯分量的权值、均值、协方差矩阵和标签分别存储在U_w、U_m、U_P和U_l中，其中，中的元素依次存储至U_w的第一列，中的元素依次存储至U_m的第一列，中的元素依次存储至U_P的第一列，中的元素依次存储至U_l的第一列；按ze＝1,...,|Z_k,ef|，依次存储由更新的高斯分量的参数至U_w、U_m、U_P和U_l中，具体为：中的元素依次存储至U_w的第ze+1列，中的元素依次存储至U_m的第ze+1列，中的元素依次存储至U_P的第ze+1列，中的元素依次存储至U_l的第ze+1列；

(5b)令

(5c)设k时刻的状态估计集航迹标签集发生时刻集有效量测在U_w中的列号集

(5d)取U_w中第二列至最后一列的所有元素中的最大值，如果此最大值小于状态提取的基本权值阈值w₁，则直接跳至步骤(5h)；否则，根据此最大值所在的不同位置，从U_l中取这些不同位置处的标签，接着对这些标签进行处理，去掉其中重复的元素，得到不重复的标签集

(5e)按依次对所有标签为的高斯分量进行处理，具体包括：

(5e.1)根据标签在U_l中的行号，从U_w的相应行中取最大值，记这个最大值的列号为记这个最大值在U_m中相应位置处的元素为m⁽ⁱ¹⁾，U_l中所有标签为的行序号集记作

(5e.2)判断标签的属性，根据其属性对其进行相应的操作，存在以下三种情况：

A.为新生分量的标签，具体操作步骤包括：

①从矩阵U_w中所有列在行中的元素中取大于或等于w₂的元素所在列的不重复的列序号构成集合同时将这些列序号存至CI中，其中，w₂为大权值分量的阈值；

②若依次按执行以下操作：

首先，记U_w中第列所在行的元素中最大值所在行的行号为提取U_m的第列在第行中的元素为单目标状态估计，将其存储至状态估计集中，并建立新的确定航迹r_max＝r_max+1，且将r_max和k分别存储至航迹标签集和发生时刻集中；

接着，分别对和进行修改，为第列元素所依据的大权值量测标出其后验信息，具体操作为：

其中，记中行在所有列中的元素为记中第行在第列中的元素为记中行在所有列中的元素为记的最后一行下面添加的那一行元素为记的第列在最后一行的元素为记U_w中第行在第列中的元素为记的最后一行下面添加的那一行元素为记U_m中第行在所有列中的元素为记的最后一行下面添加的那一行元素为记U_P中第行在所有列中的元素为记的最后一行下面添加的那一行元素为

③取中权值大于或等于w₁、且小于w₂的分量所在列的不重复的列序号集，记为从中剔除包含的元素，得到同时将这些列序号存至CI中；

④若依次按执行以下操作：

首先，记U_w中第列所在行的元素中最大值所在行的行号为建立新的暂态航迹r_maxt＝r_maxt+1；

接着，分别对和进行修改，为第列元素所依据的小权值量测标出其后验信息，具体操作为：

其中，记的第行在第列的元素为记U_w的第行在第列的元素为记的最后一行元素在第列的元素为记U_m和U_P的第行在所有列的元素分别为和

B.为确定分量的标签，从此类分量中提取单目标状态估计，具体步骤包括：

①如果对应的最大值小于状态预处理的权值阈值w₃，记U_w的第1列在行中的元素中最大值所在的行号为依次执行以下状态估计预处理步骤：

其中，||b1-b2||表示b1与b2之间的欧氏距离；U_w和U_m的第行在第1列的元素分别记为和w₃满足w₁≤w₃＜w₂；已知零均值的量测噪声协方差矩阵 和分别为x轴和y轴的量测方差，令σ＝[σ_x σ_y]^T，d(1)＝||σ-0||；

②将m⁽ⁱ¹⁾作为单目标状态估计，将其存储至状态估计集中，并将和k分别存储至航迹标签集和发生时刻集中，同时，将存储至CI；

③高斯分量权值校正

先逐一计算m⁽ⁱ¹⁾至同标签分量之间的欧氏距离：

再根据条件校正权值：

其中，U_w和U_m的第行在第ze1列的元素分别记为和

C.为暂态分量的标签，具体步骤包括：

①当时，或当且时，取m⁽ⁱ¹⁾为单目标状态估计，将其和k分别存储至状态估计集和发生时刻集中；建立新的确定航迹r_max＝r_max+1，将r_max存储至航迹标签集修改和中元素，置以及将存储至其中，为k-1时刻的小权值暂态分量集，为k-2时刻的小权值暂态分量集；

②否则，将存储在k时刻的小权值暂态分量集中；

(5e.3)根据两个“一对一”原则，将U_w中标签为的所有元素置零；

(5f)执行完循环之后，以CI的元素为列号集，根据两个“一对一”原则，将这些列在U_w的任意行的元素全部置零；

(5g)返回至第(5d)步；

(5h)记为和矩阵的行数，分别将和中的个元素逐列转换为维的数组，得到修正后的更新高斯分量的参数集

步骤6、高斯分量剪枝，具体包括：

(6a)剪去连续n_nost次没有状态估计的确定分量或暂态分量，具体步骤包括：

(6a.1)取中不重复的标签集，从中删除和中包含的标签，得到k时刻未获得状态估计的确定分量和暂态分量的标签集

(6a.2)按依次判别在前n_nost-1个时刻未获得状态估计的确定分量和暂态分量的标签集中的出现情况：如果在k1＝k-n_nost+1,...,k-1时刻的每个中都出现过，则置中所有标签为的分量权值为零；

(6b)将所有权值大于小权值阈值w_et的高斯分量的序号记录在I1中；

(6c)令J_k|k＝|I1|，依此类推，得到其中，和分别是第j3个高斯分量的权值、均值、协方差矩阵和标签；

步骤7、基于步骤6剪枝后得到高斯分量的参数集进行高斯分量合并，令为高斯分量的序号集，且J_k＝0，具体步骤包括：

(7a)当I2为空时，跳至步骤(7e)，否则，J_k＝J_k+1，得到标签(7b)取与同标签、且可合并的分量的序号集为I3，其中，d_Th为合并的距离阈值；

(7c)合并序号集I3中的所有高斯分量参数， 得到一个新的高斯分量；

(7d)I2＝I2-I3，并返回步骤(7a)；

(7e)得到

步骤8、将不同时刻具有相同身份标识的状态估计相连，即得到显式的多目标航迹；

步骤9、k＝k+1，返回步骤2。

2.如权利要求1所述的基于高斯混合概率假设密度滤波器的显式多目标跟踪方法，其特征在于，步骤3中的波门法的具体步骤为：

(3a)已知零均值的量测噪声协方差矩阵 和分别为x轴和y轴的量测方差，令σ＝[σ_x σ_y]^T，以此为依据计算波门阈值d(a)＝a||σ-0||。其中，a为置信系数；

(3b)设有效量测集

(3c)为k时刻获得的最新量测集，按ze1＝1,...,|Z_k|，依次判断每个的属性。判断时，按j＝1,...,J_k|k-1，逐一计算至每个高斯分量均值的距离如果d^(ze1,j)≤d(a)，则判定附近存在高斯分量，视为有效量测，将其加入到有效量测集Z_k,ef中，接着跳出此循环，执行对的判断。

3.如权利要求1所述的基于高斯混合概率假设密度滤波器的显式多目标跟踪方法，其特征在于，若标签l≤V₂，则为确定分量的标签；若V₂＜l≤V₁，则为暂态分量的标签；若l＞V₁，则为新生分量的标签。

4.如权利要求3所述的基于高斯混合概率假设密度滤波器的显式多目标跟踪方法，其特征在于，k＝0时刻，确定分量和暂态分量的标签分别以r_max＝0和r_maxt＝V₂为基数，每新建一个确定航迹，r_max加1；每新建一个暂态航迹，r_maxt加1。