CN104156984B - 一种不均匀杂波环境下多目标跟踪的概率假设密度方法 - Google Patents

Abstract

本发明实施例提供了一种不均匀杂波环境下多目标跟踪的概率假设密度方法，涉及智能信息处理领域，可以在计算量减小的情况下，提高计算结果的准确性。所述方法包括：监视区域为稀疏杂波区域的，直接进行标准PHD预测及更新，获得估计的目标数以及目标状态；监视区域为密集杂波区域的则寻找凸包，确定杂波区；并在逐时刻进行PHD预测时选择未包含在杂波区的回波进行计算，在进行PHD更新时，加入真实目标的量测刚好落入所述杂波区的回波；将所述后验强度中的高斯成分进行剪枝与合并，获得目标强度；然后计算所述目标强度中所有高斯成分的权值和，得到监视区域内的估计目标数；提取目标强度中权值大于τ的高斯成分，作为估计的目标状态。

Description

一种不均匀杂波环境下多目标跟踪的概率假设密度方法

技术领域

本发明涉及智能信息处理领域，尤其涉及一种不均匀杂波环境下多目标跟踪的概率假设密度方法。

背景技术

多目标跟踪是研究从目标量测和杂波中估计目标数目与各个目标状态的方法。传统的处理多目标跟踪问题的方法大都是基于数据关联的，需要建立量测和目标的对应关系，如最近邻法(NN)，联合概率数据关联算法(JPDA)，多假设跟踪算法(MHT)。近年来，基于随机有限集(Random Finite Set，RFS)理论的多目标跟踪算法得到很大的关注。Mahler提出的概率假设密度(Probability Hypothesis Density，PHD)滤波算法就是一种在随机有限集框架下，相比完全多目标贝叶斯滤波器，运算复杂度得到有效降低的算法。

在传统的PHD滤波器中，通常假设在可观测的范围内杂波是均匀分布的，但在实际环境中杂波的空间分布在观测区域是不均匀的，杂波区可分为稀疏杂波区和密集杂波区，在密集杂波区，虚警的数目远远大于真实目标的量测，这种情况下进行多目标跟踪是一个困难的课题。因此，当杂波分布不均匀时，传统的PHD滤波器的跟踪性能会急剧下降。针对此问题，现有技术中有一种方法是在PHD粒子滤波的过程中，通过设定门限，将位于门限之外的观测值滤除，仅利用位于门限内的观测值进行粒子权重的更新。还有一种方法是提出一种应用目标起始和维持规则辅助的SMC-PHD滤波方法，每个粒子除了状态矢量和权值外，还增加了两个辅助变量，一个用来表示粒子的存活年龄，一个用来表示最近N个周期粒子受到观测值更新的情况。

上述这些方法，都是针对密集杂波环境下目标跟踪的PHD方法，而在既有稀疏杂波又有密集杂波的复杂的杂波环境下应用这些方法，就会导致在稀疏杂波区和密集杂波区计算量都很大，不能实现快速运算，然而在稀疏杂波区计算量本不用那么大，应用这些方法反而导致计算量大，而且得到的结果也不一定理想。

发明内容

本发明的实施例提供一种不均匀杂波环境下多目标跟踪的概率假设密度方法，可以在计算量减小的情况下，提高计算结果的准确性。

为达到上述目的，本发明的实施例采用如下技术方案：

一种不均匀杂波环境下多目标跟踪的概率假设密度方法，包括：

判断监视区域内n帧累积的回波数是否大于门限值ε，其中，所述n，ε为预设数值；

若否，则直接进行标准PHD预测及PHD更新，获得估计的目标数以及目标状态；

若是，则寻找凸包，确定杂波区；

初始化，获得初始时刻的后验强度，所述初始时刻为0时刻；

针对k时刻测得的回波值，选择未包含在所述杂波区的回波进行PHD预测计算，根据k-1时刻的后验强度，获得k时刻的预测强度，其中，所述k为大于等于1的整数；

根据所述k时刻的预测强度，用更新的回波进行PHD更新，获得k时刻的后验强度，其中，所述更新的回波包括k时刻测得的回波值中所述未包含在所述杂波区的回波，以及k时刻真实目标的量测刚好落入所述杂波区的回波；所述后验强度为高斯混合形式，包括各个高斯成分的权值、均值和协方差；

将所述后验强度中权值小于T的高斯成分滤除，将高斯成分间的距离小于U的高斯成分进行合并，获得目标强度；其中，U≤10，T≤10^-5；

计算所述目标强度中所有高斯成分的权值和，得到所述监视区域内的估计目标数；提取所述目标强度中权值大于τ的高斯成分，作为估计的目标状态，其中，所述τ≥0.5。

可选的，所述寻找凸包，确定杂波区，包括：

应用近似传播(AP)聚类算法对所述监视区域内的n帧累积的回波进行聚类，获得聚类结果；

对于所述聚类结果中聚类点数大于η的区域，用Graham scan方法寻找凸包，确定杂波区，所述η根据仿真场景而定。

上述技术方案提供的不均匀杂波环境下多目标跟踪的概率假设密度方法，利用AP聚类算法，首先对监视区域内多帧累积的回波进行聚类，用寻找凸包的方法确定杂波区，然后再逐观测时刻进行PHD预测和更新，在PHD预测时，不用杂波区内的回波，但当进行PHD更新时，需要将发现的真实目标的量测刚好落入杂波区中的回波取出进行PHD更新，这样可以保证算法的高精度，使估计结果更准确。另外，当监视区域内的回波数大于某个门限时才进行聚类和确定杂波区的运算，否则不作此运算直接进行标准的PHD预测和更新，这样计算过程简单，对于不均匀杂波区可以实现快速运算，从而实现不均匀杂波环境下PHD方法跟踪性能的提升。

附图说明

图1为本发明实施例提供的一种不均匀杂波环境下多目标跟踪的概率假设密度方法的流程示意图；

图2为本发明实施例提供的一种AP聚类算法数据点之间消失传递示意图；

图3为本发明实施例提供的一种应用Graham scan算法寻找凸包的示意图；

图4为本发明实施例提供的一种不均匀杂波环境下进行多目标跟踪的仿真场景图；

图5为本发明实施例提供的一种应用本发明提供的方法仿真出来的状态估计图；

图6为本发明实施例提供的一种应用本发明提供的方法仿真出来的目标数估计图；

图7为本发明实施例提供的一种应用高斯混合PHD方法仿真出来的状态估计图；

图8为本发明实施例提供的一种应用高斯混合PHD方法仿真出来的目标数估计图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明实施例提供了一种不均匀杂波环境下多目标跟踪的概率假设密度方法，如图1所示，所述方法包括以下步骤：

101、判断监视区域内n帧累积的回波数是否大于ε。

其中，所述n，ε为预设数值；所述ε是依据监视区域内的回波数而由使用者预设的，所述n也是由使用者自己预设的。

若监视区域内n帧累积的回波数大于ε，则表明所述监视区域为密集杂波区域，进行步骤103-107；若监视区域内n帧累积的回波数小于ε，则表明所述监视区域为稀疏杂波区域，进行步骤102。

102、直接进行标准PHD预测及PHD更新，获得估计的目标数以及目标状态。

在这里所述的标准PHD预测及PHD更新，即为现有技术中的传统的PHD预测及PHD更新，本领域技术人员都清楚了解其计算过程，故在此不再进行赘述。

103、寻找凸包，确定杂波区。

可选的，步骤103包括以下步骤：

A、应用AP聚类算法对所述监视区域内的n帧累积的回波进行聚类，获得聚类结果。

AP(Affinity Propagation，近似传播)聚类算法是根据N个数据点之间的相似度进行聚类，这些相似度可以是对称的，即两个数据点互相之间的相似度一样(如欧式距离)；也可以是不对称的，即两个数据点互相之间的相似度不等。这些相似度组成N×N的相似度矩阵S(其中N为有N个数据点)。AP算法不需要事先指定聚类数目，相反它将所有的数据点都作为潜在的聚类中心，称之为exemplar。

AP聚类算法的主要思想是：以S矩阵的对角线上的数值S(k′，k′)作为k′点能否成为聚类中心的评判标准，这意味着该值越大，这个点成为聚类中心的可能性也就越大，这个值便是参考度P。聚类的数量受到参考度P的影响，如果认为每个数据点都有可能作为聚类中心，那么P就应取相同的值。如果取输入的相似度的均值作为P的值，得到聚类数量是中等的。如果取最小值，得到类数较少的聚类。AP聚类算法中传递两种类型的消息，即Responsibility和Availability。R(i′，k′)表示从点i′发送到候选聚类中心k′的数值消息，反映k′点是否适合作为i′点的聚类中心。A(i′，k′)则表示从候选聚类中心k′发送到i′的数值消息，反映i′点是否选择k′作为其聚类中心。R(i′，k′)与A(i′，k′)越强，则k′点作为聚类中心的可能性就越大，并且i′点隶属于以k′点为聚类中心的聚类的可能性也越大。AP聚类算法通过迭代过程不断更新每一个点的吸引度和归属度值，直到产生M个高质量的exemplar，同时将其余的数据点分配到相应的聚类中。

在这里需要说明的是：

1、exemplar：指的是聚类中心。

2、similarity：数据点i′和点j′的相似度记为S(i′，j′)。是指点j′作为点i′的聚类中心的相似度。一般使用欧氏距离来计算。在该算法中，所有点与点的相似度全部取为负值。这样，相似度值越大说明点与点的距离越近，便于后续计算。

3、preference：数据点i′的参考度称为P(i′)或S(i′，i′)。是指点i′作为聚类中心的参考度。一般取S相似度值的中值。

4、Responsibility：R(i′，k′)用来描述点k′适合作为数据点i′的聚类中心的程度。

Availability：A(i′，k′)用来描述点i′选择点k′作为其聚类中心的适合程度。

如图2所示，为Responsibility(图2中记为R)和Availability(图2中记为A)两者之间的关系。

R(i′，k′)、A(i′，k′)与R(k′，k′)的计算公式为：

R(i′，k′)＝S(i′，k′)-max{A(i′，j′)+S(i′，j′)}(j′∈{1，2，…，N}&j′≠k′) (1)

R(k′，k′)＝P(k′)-max{A(k′，j′)+S(k′，j′)}(j′∈{1，2，…，N}&j′≠k′) (3)

上面可以看出，当P(k′)较大使得R(k′，k′)较大时，A(i′，k′)也较大，从而类代表k′作为最终聚类中心的可能性较大；同样，当越多的P(i′)较大时，越多的类代表倾向于成为中心的聚类中心。因此，增大或减小P可以增加或减少AP输出的聚类数目。

5、Damping factor(阻尼系数)λ：主要是起收敛作用的。每次迭代，吸引度R_i′和归属度A_i′要与上一次的R_i′-1和A_i′-1进行加权更新。公式如下：

R_i′＝(1-λ)*R_i′+λ*R_i′-1 (4)

A_i′＝(1-λ)*A_i′+λ*A_i′-1 (5)

其中，λ∈[0.5，1)。

AP聚类算法的具体步骤如下：

(1)计算N个点之间的相似度值，将值放在S矩阵中，再选取P值(一般取S的中值)。

(2)设置一个最大迭代次数，迭代过程开始后，计算每一次的R值和A值，根据R(k′，k′)+A(k′，k′)的值来判断是否为聚类中心，当迭代次数超过最大值或者当聚类中心连续多少次迭代不发生改变时终止计算，获得聚类结果。

B、对于所述聚类结果中聚类点数大于η的区域，用Graham scan方法寻找凸包，确定杂波区。

所述η为预设值，可以根据仿真场景而定。

其中，Graham scan算法是一种计算一组平面有限点的凸包的算法。其求解凸包的具体步骤如下：

(1)在所有点中选取y坐标最小的一点H，当作基点。如果存在多个点的y坐标都为最小值，则选取x坐标最小的一点。坐标相同的点应排除。然后按照其它各点P和基点构成的向量＜H，P＞与x轴的夹角进行排序，夹角由大至小进行顺时针扫描，反之则进行逆时针扫描。实现中无需求得夹角，只需根据向量的内积公式求出向量的模即可。以图3为例，基点为H，根据夹角由小至大排序后依次为H，K，C，D，L，F，G，E，I，B，A，J。下面进行逆时针扫描。

(2)线段＜H，K＞一定在凸包上，接着加入C。假设线段＜H，C＞也在凸包上，因为就H，K，C三点而言，它们的凸包就是由此三点所组成。但是接下来加入D时会发现，线段＜K，D＞才会在凸包上，所以将线段＜K，C＞排除，C点不可能是凸包。即当加入一点时，必须考虑到前面的线段是否会出现在凸包上。

(3)从基点开始，凸包上每条相临的线段的旋转方向应该一致，并与扫描的方向相反。如果发现新加的点使得新线段与上线段的旋转方向发生变化，则可判定上一点必然不在凸包上。实现时可用向量叉积进行判断，设新加入的点为P_n+1，上一点为P_n，再上一点为P_n-1。顺时针扫描时，如果向量＜P_n-1，P_n＞与＜P_n，P_n+1＞的叉积为正(逆时针扫描判断是否为负)，则将上一点删除。删除过程需要回溯，将之前所有叉积符号相反的点都删除，然后将新点加入凸包。在图3中，加入K点时，由于线段＜H，K＞相对于＜H，C＞为顺时针旋转，所以C点不在凸包上，应该删除，保留K点。接着加入D点，由于线段＜K，D＞相对＜H，K＞为逆时针旋转，故D点保留。按照上述步骤进行扫描，直到点集中所有的点都遍历完成，即得到如图3所示的由点H、J、B、G、D以及K连线形成的凸包，确定出杂波区。

当然，寻找凸包，确定杂波区还可以有现有技术中的其他方法，在这里不再详述。

1041、初始化，获得初始时刻的后验强度，所述初始时刻为0时刻。

在进入仿真场景进行PHD预测时，应当假设目标运动模型和观测模型都满足线性高斯条件；假设目标的存活概率和探测概率相互独立，且与目标状态无关；假设新生目标随机集和衍生目标随机集的强度均为高斯混合形式；用J₀个高斯成分的加权和进行初始化。

假设目标运动模型和观测模型都满足线性高斯条件：

目标运动模型：f_k|k-1(x|ζ)＝N(x；F_k-1ζ，Q_k-1) (6)

观测模型：g_k(z|x)＝N(z；H_kx，R_k) (7)

其中，N(·；m，P)表示其密度均值为m，方差为P的高斯分布；F_k-1为状态转移矩阵，Q_k-1为***噪声协方差矩阵；H_k为观测矩阵，R_k为测量噪声协方差矩阵。

假设目标的存活概率和探测概率相互独立，且与目标状态无关：

存活概率：p_S，k(x)＝p_S，k (8)

探测概率：p_D，k(x)＝p_D，k (9)

假设新生目标随机集和衍生目标随机集的强度均为高斯混合形式：

新生目标随机集的强度：

衍生目标随机集的强度：

其中，J_γ，k，i＝1，2，...，J_γ，k等参数决定了新生目标随机集的强度，为新生目标随机集强度的第i个高斯成分的权值、均值和协方差，J_γ，k是k时刻新生目标高斯成分的个数。同理：J_β，k，j＝1，2，...，J_β，k等参数决定了由目标ζ衍生的目标随机集的强度，为新生目标随机集强度的第i个高斯成分的权值、均值和协方差，J_β，k是k时刻新生目标高斯成分的个数。在这里所述的各个参数的在初始时刻即0时刻时的值都是已知的。在后续时刻的值可以根据初始时刻即0时刻时的值推导计算出来。

用J₀个高斯成分的加权和进行初始化，获得初始时刻的后验强度：

其中N(·；m，P)是均值为m，方差为P的高斯分布。其权值的和是期望的初始目标数：公式(12)中的其他相关参数J₀，在用户应用本发明进行仿真时，由用户根据具体仿真环境设定自己设定。

1042、针对k时刻得到的回波值，选择未包含在所述杂波区的回波进行PHD预测计算，根据k-1时刻的后验强度，获得k时刻的预测强度。

本发明实施例中的k时刻、k-1时刻为观测时刻，其中k-1时刻为k时刻的前一观测时刻。k-1＝0时刻为初始时刻，此时，k＝1。

在逐时刻进行PHD预测时，对于每个观测时刻得到的回波值，若其中有回波在步骤103中的杂波区，则去掉那些回波，再进行PHD预测。

假定k-1时刻的后验强度具有如下混合高斯形式：

在k-1＝0时，公式(13)可以用公式(12)表示：

则针对k时刻得到的回波值，选择未包含在所述杂波区的回波进行PHD预测计算，根据k-1时刻的后验强度，获得k时刻的预测强度，k时刻的预测强度也为高斯混合形式，如下给出：

v_k|k-1(x)＝v_S，k|k-1(x)+v_β，k|k-1(x)+γ_k(x) (14)

其中，γ_k(x)在公式(10)中已给出。

由公式(13)中含有的参数以及J_k-1可以计算出公式(10)以及公式(15)-(20)，然后通过公式(10)以及公式(15)-(20)就可以计算出公式(14)中k时刻的预测强度v_k|k-1(x)。

在这里，可以将由公式(10)，以及(15)-(20)，计算得到k时刻的预测强度记为如下高斯混合形式：

105、根据所述k时刻的预测强度，用更新的回波进行PHD更新，获得k时刻的后验强度，其中，所述更新的回波包括k时刻所得的回波值中所述未包含在所述杂波区的回波，以及k时刻真实目标的量测刚好落入所述杂波区的回波。

所述后验强度为高斯混合形式，包括各个高斯成分的权值、均值和协方差。

在进行PHD更新时，需要对杂波区中的回波进行自适应判定，若确定杂波区里的都是杂波，而无真实目标的量测，则杂波区内的回波均不用；但若发现真实目标的量测刚好落入杂波区中，则将其取出放入PHD更新要用的回波中，改变杂波区的大小。判断是否落入杂波区的方法为：用目标状态计算量测，将得出的量测与杂波区的回波比对，若发现量测落入杂波区，则将量测取出，放入PHD更新需用的回波中；反之，则进行PHD更新时用不包括杂波区的回波。

由步骤1042可以获得k时刻的预测强度为如下高斯混合形式：

则k时刻的后验强度也为高斯混合形式：

其中，

由公式(21)中含有的参数以及J_k|k-1可以计算出公式(23)-(27)，然后通过公式(23)-(27)就可以计算出公式(22)中k时刻的后验强度v_k(x)。

106、将所述后验强度中权值小于T的高斯成分滤除，将高斯成分间的距离小于U的高斯成分进行合并，获得目标强度；其中，U≤10，T≤10^-5。

设定一个阈值T(T≤10^-5)，当一些高斯成分的权值小于T时，将这些高斯成分滤除掉，具体公式可以表示如下：

其中，i＝N_p+1，...，J_k是大于阈值的权值。

设定一个阈值U(U＜10)，当一些高斯成分间的距离小于U时，将这些高斯成分进行合并。给出一个最大允许的高斯成分数目J_max。

令l＝0，

重复

l：＝l+1 (29)

I：＝I\L (35)

直到

如果l＞J_max，用J_max个权值较大的高斯成分代替得到的即为合并后的高斯成分。

其具体合并过程本领域技术人员都清楚了解，在此不再详述。

进行高斯合并后的所述目标强度也是高斯混合形式，可以如下所示：

107、计算所述目标强度中所有高斯成分的权值和，得到所述监视区域内的估计目标数；提取所述目标强度中权值大于τ的高斯成分，作为估计的目标状态，其中，所述τ≥0.5。

获取估计目标数：通过混合高斯项的权值和获取监视区域内的目标数，即

获取目标状态：提取权值大于τ(τ≥0.5)的高斯成分，即

为了验证本发明方案，下面给出面向不均匀杂波环境下多目标跟踪的自适应概率假设密度方法仿真结果。

仿真场景：考虑不均匀杂波环境下监视区域内有三个目标的场景。

目标1在k＝1时刻出现，初始状态为[-25m，4m/s，10m，2m/s]^T，在k＝38时刻消失；目标2在k＝10时刻出现，初始状态为[550m，-15m/s，100m，4m/s]^T，在k＝60时刻消失；目标3在k＝22时刻出现，初始状态为[-270m，18m/s，-190m，6m/s]^T，在k＝51时刻消失。

目标状态方程为：

X_k+1＝F_kX_k+v_k (39)

其中，状态向量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时刻的目标速度。采样间隔T＝1s。状态转移矩阵F_k、过程噪声的协方差分别是式中，过程噪声是零均值的高斯白噪声，且和量测噪声序列相互独立。目标存活概率为p_S＝0.975。

量测方程为：

Z_k＝H_kX_k+ω_k (40)

式中，量测噪声是零均值的白噪声。量测矩阵和量测噪声协方差分别是目标探测概率为p_D＝0.96。

杂波的空间分布是非均匀的，在某几个区域杂波很密集。整个监视空域是[-400m，600m]×[-400m，800m]，在整个区域每个观测时刻的平均杂波数是5。密集杂波区域为：在[-200m，-150m]×[-200m，-150m]区域，每个观测时刻的平均杂波数是20；在[50m，100m]×[0m，100m]区域，每个观测时刻的平均杂波数是20；在[300m，350m]×[100m，150m]区域，每个观测时刻的平均杂波数是20。

仿真场景见图4。

仿真时，n＝5，ε＝200，η＝80。

在仿真实验中，本发明方法与传统的高斯混合PHD(简称为：GM-PHD)方法(标准PHD的一种实现方式)进行对比。本发明仿真结果见图5及图6，高斯混合PHD的仿真结果见图7及图8。

对比两种方法，可以看出，在不均匀杂波环境下，即在既有稀疏杂波区又有密集杂波区的情况下，本发明(图5和图6显示)可以更好的估计每个时刻的目标状态及目标数，从图中可看出，相比高斯混合PHD，本发明估计每个时刻目标状态及目标数的精度有显著的提高，达到了提升目标跟踪性能的目的。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应所述以权利要求的保护范围为准。

Claims (2)

1.一种不均匀杂波环境下多目标跟踪的概率假设密度方法，其特征在于，包括：

判断监视区域内n帧累积的回波数是否大于门限值ε，其中，所述n，ε为预设数值；

若否，则直接进行标准PHD预测及PHD更新，获得估计的目标数以及目标状态；

若是，则寻找凸包，确定杂波区；

初始化，获得初始时刻的后验强度，所述初始时刻为0时刻；

针对k时刻测得的回波值，选择未包含在所述杂波区的回波进行PHD预测计算，根据k-1时刻的后验强度，获得k时刻的预测强度，其中，所述k为大于等于1的整数；

根据所述k时刻的预测强度，用更新的回波进行PHD更新，获得k时刻的后验强度，其中，所述更新的回波包括k时刻测得的回波值中所述未包含在所述杂波区的回波，以及k时刻真实目标的量测刚好落入所述杂波区的回波；所述后验强度为高斯混合形式，包括各个高斯成分的权值、均值和协方差；

将所述后验强度中权值小于T的高斯成分滤除，将高斯成分间的距离小于U的高斯成分进行合并，获得目标强度；其中，U≤10，T≤10^-5；

计算所述目标强度中所有高斯成分的权值和，得到所述监视区域内的估计目标数；提取所述目标强度中权值大于τ的高斯成分，作为估计的目标状态，其中，所述τ≥0.5。

2.根据权利要求1所述的方法，其特征在于，所述寻找凸包，确定杂波区，包括：

应用近似传播(AP)聚类算法对所述监视区域内的n帧累积的回波进行聚类，获得聚类结果；

对于所述聚类结果中聚类点数大于η的区域，用Graham scan方法寻找凸包，确定杂波区，所述η根据仿真场景而定。