基于核密度估计的GM-AI-PHD多目标跟踪方法

基于核密度估计的gm-ai-phd多目标跟踪方法
技术领域
1.本发明属于单传感器多目标跟踪领域，涉及一种基于核密度估计的gm-ai-phd多目标跟踪方法，能够在复杂环境下实时且准确地实现对多目标的定位与跟踪；


背景技术：

2.多目标跟踪指的是：从传感器获取的复杂数据中估计出目标的数量、状态等信息；单传感器多目标跟踪方法主要分为两类：
①
基于数据关联(data association,da)的方法；
②
基于随机有限集(random finite set,rfs)的方法；其中，基于da的跟踪方法在面临传感器质量下降，杂波强度变大等诸多更为复杂的目标跟踪问题时，往往因为其计算量复杂而存在着组合爆炸的风险，在一些复杂环境中，难以满足应用的需要；近年来，基于rfs的跟踪方法引起了广泛的关注，一些基于rfs的滤波器被陆续提出，例如概率假设密度(probability hypothesis density,phd)滤波器、势phd(cphd)滤波器和幅度信息(amplitude information,ai)辅助的phd(ai-phd)滤波器等；
3.以下是基于rfs跟踪方法的部分代表性工作：
4.2003年，r.mahler提出基于rfs的phd滤波器，将状态变量的概率空间建模简化为状态变量的概率假设空间建模，但是没有得到封闭式的解；随后，澳大利亚学者b.n.vo提出了可实现的高斯混合phd(gm-phd)滤波器和序贯蒙特卡罗phd(smc-phd)滤波器；2007年，r.mahler 提出cphd滤波器；2008年，clark等人首次将幅度信息引入phd滤波器，其中目标和杂波幅度采用为rayleigh分布；2009年，edmund f.brekke等人将幅度信息假设为k分布，并引入多目标跟踪领域；随后，li等人针对高分辨雷达和低擦地角引起杂波重拖尾的情况，假设幅度服从威布尔分布，提出一种ai辅助的phd滤波器；
5.上述利用ai辅助目标跟踪的方法中，均假设目标和杂波的幅值服从一种已知统计分布，具体地，该方法将统计分布对应的概率密度函数引入跟踪滤波器中，实现目标量测和杂波量测的快速有效区分，从而提高滤波器的跟踪性能；但是当真实的统计分布未知时，利用上述方法是不可行的，因此本发明提出利用核密度估计方法实时估计目标和杂波的幅度概率密度函数，使得在统计分布未知场景下，也能够实现对多目标的准确跟踪；


技术实现要素：

6.本发明的目的是针对现有技术的不足，提出基于核密度估计(kernel density estimation,kde) 的gm-ai-phd多目标跟踪方法，简称gm-kde-ai-phd，具体内容如下：
7.步骤(1)、构建单传感器多目标跟踪场景，设置相关参数，包括目标运动的过程噪声和传感器的观测噪声；
8.步骤(2)、对多目标的状态、观测进行建模；
9.步骤(3)、构建gm-kde-ai-phd滤波器；
10.步骤(4)、利用建好的gm-kde-ai-phd滤波器实现多目标跟。
11.本发明的第二个目的是提供基于核密度估计的gm-ai-phd多目标跟踪装置，包括
gm
‑ꢀ
kde-ai-phd滤波器。
12.本发明的第三个目的是提供一种计算机可读存储介质，其上存储有计算机程序，当所述计算机程序在计算机中执行时，令计算机执行所述的方法。
13.本发明的第四个目的是提供一种计算设备，包括存储器和处理器，所述存储器中存储有可执行代码，所述处理器执行所述可执行代码时，实现所述的方法。
14.本发明的有益效果是：
15.本发明提出了一套完整的多目标跟踪方案，利用核密度估计方法对目标、杂波幅度分布的概率密度函数进行估计，并将该估计结果嵌入ai-phd滤波器过程中，最后采用线性高斯模型实现滤波器的求解；本发明可以在杂波数多、幅度信息统计分布未知的复杂场景下实时在线地实现对多目标的跟踪，较原始的phd滤波器跟踪性能有巨大提升，较已知统计分布的ai-phd滤波器跟踪性能也有巨大提升，可以有效提高对多目标跟踪的效率、精度和鲁棒性；
附图说明
16.图1是本发明方法gm-kde-ai-phd核心部分具体实现的流程图；
17.图2是多目标运动轨迹图；
18.图3是在200次蒙特卡洛下，本发明方法与真实分布gm-ai-phd滤波器、瑞利分布gm
‑ꢀ
ai-phd滤波器、无幅度gm-phd滤波器得出的目标估计个数平均值的比较图；
19.图4是在200次蒙特卡洛下，本发明方法与真实分布gm-ai-phd滤波器、瑞利分布gm
‑ꢀ
ai-phd滤波器、无幅度gm-phd滤波器得出的ospa平均值的比较图；
具体实施方式
20.下面结合具体实施方式对本发明做进一步的分析；
21.基于核密度估计的gm-ai-phd的多目标跟踪方法，包括以下步骤：
22.步骤(1)、构建单传感器多目标跟踪场景，设置相关参数，包括目标运动的过程噪声和传感器的观测噪声，具体是：
23.目标在二维平面内运动，其状态表示为x＝[p
x
,v
x
,py,vy]
t
，其中p
x
,py分别为目标在x、y方向的位置，v
x
、vy分别为目标在x、y方向的速度，t表示转置；
[0024]
设置目标运动的过程噪声协方差为σ
x
(k)2表示k时刻环境对目标在x 速度的噪声方差，σy(k)2表示k时刻环境对目标在y速度的噪声方差；
[0025]
设置传感器的观测噪声协方差为其中δ
x2
、δ
y2
分别表示传感器对目标x位置与y位置的噪声方差；
[0026]
步骤(2)、对多目标的状态、观测进行建模，具体是：
[0027]
2-1目标状态建模
[0028]
设k时刻的目标个数为nk，目标状态依次为k时刻目标状态可由随机有限集(rfs)xk表示，即
[0029][0030]
其中χ为状态空间，是χ的所有有限子集集合；
[0031]
目标从k-1时刻的某一状态x
k-1
转移到k时刻的某一状态xk称为状态转移，假设每个目标状态转移概率密度f
k|k-1
(xk|x
k-1
)符合线性高斯模型，表示为：
[0032]fk|k-1
(xk|x
k-1
)＝n(xk；f
k-1
x
k-1
,q
k-1
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
[0033]
其中n(xk；f
k-1
x
k-1
,q
k-1
)表示目标状态xk的空间分布，其服从状态均值为f
k-1
x
k-1
，协方差为 q
k-1
的高斯分布，f
k-1
为状态转移矩阵，q
k-1
为过程噪声协方差矩阵；
[0034]
2-2目标运动学观测建模
[0035]
设k时刻接收的目标量测个数为mk，目标量测依次为k时刻传感器量测矢量表示为zk＝[z
x
，zy]
t
，k时刻目标量测可由有限集zk表示，即
[0036][0037]
其中，z
x
,zy分别为目标在x,y方向传感器测得的位置，t表示转置，为观测空间，是的所有有限子集集合；
[0038]
假设k时刻传感器测量模型g
k|k-1
(zk|xk)符合线性高斯模型，表示为：
[0039]gk|k-1
(zk|xk)＝n(zk；hkxk,rk)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0040]
其中n(zk；hkxk,rk)表示k时刻某一量测zk的空间分布，其服从均值为hkxk，协方差为rk的高斯分布，hk为观测矩阵，rk为测量噪声协方差矩阵；
[0041]
2-3杂波运动学观测建模
[0042]
设k时刻接收的杂波个数服从参数为λ的泊松分布，杂波在量测空间中的位置分布服从均匀分布，即
[0043][0044]
其中，nc代表k时刻杂波期望数，ρ(nc)表示参数为λ的泊松分布的概率函数，表示观测到杂波j的观测取值为的概率密度，v表示观测空间的体积；
[0045]
那么，k时刻杂波的强度表示为：
[0046][0047]
2-4、目标幅度观测建模：
[0048]
为构建目标观测幅度数据，设本发明的目标幅度服从威布尔分布，其概率密度表示为：
[0049][0050]
其中，a表示幅度，u是比例参数，v是形状参数，本发明中，取u＝25，v＝1.5；
[0051]
用ga(a)表示a》0的目标幅度概率密度函数，即：
[0052][0053]
设幅值阈值为τ，当目标量测幅度大于τ时，认为有概率pd被传感器检测为目标，该概率称为检测概率，此时目标幅度的概率密度函数用表示：
[0054][0055][0056]
2-5杂波幅度观测建模
[0057]
为构建杂波观测幅度数据，设本发明的杂波幅度服从威布尔分布，其概率密度如式(7)所示，其中，取u＝3,v＝1；用ca(a)表示a》0的杂波幅度概率密度函数，即：
[0058][0059]
当杂波幅度大于τ时，认为有概率p
fa
被传感器检测为目标，该概率称为虚警概率，此时杂波幅度的概率密度函数用表示：
[0060][0061][0062]
一般地，p
fa
为恒定值，因此，由式(12)可确定唯一的幅度阈值τ，与便不难确定；
[0063]
步骤(3)、构建gm-kde-ai-phd滤波器，具体是：
[0064]
3-1phd滤波器
[0065]
贝叶斯滤波器中，利用随机有限集(rfs)模型，多目标预测概率密度和更新概率密度如式 (14)-(15)所示：
[0066]
p
k|k-1
(xk|z
1:k-1
)＝∫f
k|k-1
(xk|x
k-1
)p
k-1
(x
k-1
|z
1:k-1
)μs(dx
k-1
)
ꢀꢀꢀꢀꢀꢀ
(14)
[0067][0068]
其中，z
1:k
表示从1到k时刻的所有观测的集合，p
k|k-1
(xk|z
1:k-1
)表示k时刻多目标预测概率密度函数，f
k|k-1
(xk|x
k-1
)表示从k-1时刻状态x
k-1
到k时刻状态xk的转移概率密度函数， p
k-1
(x
k-1
|z
1:k-1
)表示k-1时刻多目标更新概率密度函数，pk(xk|z
1:k
)表示k时刻多目标更新概率密度函数,gk(zk|xk)表示k时刻似然函数，μs是在上的适当参考量，dx
k-1
表示x
k-1
的导数；
[0069]
递归方程(14)-(15)涉及到空间上的多个积分，这在计算上很难处理；phd滤波器是为解决多目标贝叶斯滤波器计算难解性而开发的一种近似方法，phd滤波器没有及时传播多目标更新密度，而是传播更新强度；预测强度和更新强度分别如式(16)-(17)所示：
[0070]vk|k-1
(xk)＝∫p
s,kfk|k-1
(xk|x
k-1
)v
k-1
(x
k-1
)dx
k-1
γk(xk)
ꢀꢀꢀꢀꢀꢀ
(16)
[0071][0072]
其中v
k|k-1
(xk)和vk(xk)分别表示递归式(14)(15)中多目标预测密度p
k|k-1
和多目标更新密度pk的强度，γk(xk)为k时刻新出现目标rfs的强度，p
s,k
为目标在k时刻仍保存先前状态的概率即生存概率，f
k|k-1
(xk|x
k-1
)为从状态x
k-1
到状态xk的转移概率密度函数，p
d,k
(xk)为k时刻的检测状态为xk的概率即检测概率，ck(zk)为k时刻杂波rfs的强度，gk(zk|xk)为状态xk到zk的似然函数；
[0073]
3-2ai-phd滤波器
[0074]
将幅度信息引入phd滤波器，建立ai-phd的递归方程；
[0075]
①
预测强度：ai-phd滤波器状态矢量与phd滤波器完全相同，因此预测步公式与phd滤波器一致，预测强度如下所示：
[0076]dk|k-1
(xk)＝∫p
s,kfk|k-1
(xk|x
k-1
)d
k-1
(x
k-1
)dx
k-1
γk(xk)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(18)
[0077]
②
更新强度：引入目标和杂波的幅度信息后，更新强度发生变化，更新强度如下所示：
[0078][0079][0080][0081]
其中，a表示该量测对应的幅度大小，为k时刻的带幅度量测有限集，表示带幅度杂波rfs强度，表示带幅度似然函数；表示加入幅度信息后的检测概率，且gk(zk|xk)和ck(zk)为式(17)的似然函数和杂波rfs强度， ga(a)和ca(a)分别为目标和杂波幅度的概率密度函数；
[0082]
将式(20)(21)带入式(19)，更新强度的表达式由式(19)变为式(22)：
[0083][0084]
3-3 kde-ai-phd(基于核密度估计的幅度信息辅助的phd)滤波器
[0085]
利用kde方法结合观测幅度数据对目标幅度概率密度函数ga(a)杂波幅度概率密度函数 ca(a)进行估计，核密度估计的公式如式(23)所示：
[0086][0087]
其中，α表示随机变量，f(α)为核密度估计的概率密度函数，n表示已知样本的数
量，h表示带宽，αi表示α的测量值，表示核函数；
[0088]
在开始滤波之前，可采集待滤波场景的离线观测幅度数据对目标离线幅度概率密度函数和杂波离线幅度概率密度函数进行估计，但是如果将其直接代替滤波过程式(22) 中的ga(a)和ca(a)，则不具有准确性和实时性，故本发明将kde估计嵌入滤波过程中结合观测幅度数据进行实时估计，保证滤波器的准确性和实时性；具体实现方法如下所示：
[0089]
设滤波前采集的离线幅度数据为第0时刻数据a0，且滤波过程中第k 时刻接收的幅度数据为其中k≥1，a表示辐度，n0表示第0时刻采集的幅度数据个数，nk表示第k时刻采集的幅度数据个数，一般n0》nk，kde-ai-phd滤波器更新强度计算如下：
[0090]
step1：数据处理
[0091]
为保证实时性，用于滤波的第k时刻的幅度数据表示为前k时刻幅度数据的累加，
[0092]
step2：核密度估计
[0093]
步骤
①
：采用合适的阈值τ(τ一般不超过5)，对离线数据进行区分，如果幅度a大于τ则认为是目标幅度，构成目标幅度数据集反之则认为是杂波幅度，构成杂波幅度数据集上标p和q分别表示集合a
t,k
和a
c,k
的数量；
[0094]
步骤
②
：选择核函数和带宽h，可选择高斯函数，h可取值为1.8，即：
[0095][0096]
步骤
③
：对离线样本a
t
和离线样本ac进行采用式(23)进行kde估计，得出相应的和为：
[0097][0098][0099]
step3：将step2估计的和代替滤波过程式(22)中的ga(a)和ca(a)，得出实时的更新强度dk(xk)，但预测强度不变；
[0100][0101]
3-4 gm-kde-ai-phd滤波器
[0102]
当系统的状态方程和观测方程满足式(2)、(4)的线性高斯模型时，kde-ai-phd的预测和更新可以由高斯和的闭合形式描述，k时刻预测和更新步骤的高斯实现如下：
[0103]
假设在k-1时刻，多目标更新强度表示为如下的高斯混合形式：
[0104][0105]
其中，j
k-1
为k-1时刻更新强度高斯项的个数，和分别为k-1时刻第i个更新高斯项的权重、均值和协方差矩阵；
[0106]
新出现目标rfs的phd可以表示为式(29)：
[0107][0108]
其中，j
γ,k
表示k时刻新生目标高斯项的个数，分别表示k时刻第j个新生目标高斯项的权重、均值和协方差矩阵；
[0109]
①
预测强度
[0110][0111]
其中，γk(xk)由式(29)给出，j
k-1
为k-1时刻更新强度高斯项的个数，式(30)中高斯项的权重、均值和协方差表示为：
[0112][0113]
其中，p
s,k
表示k时刻目标的生存概率，f
k-1
表示状态转移矩阵，q
k-1
表示过程噪声协方差矩阵；
[0114]
预测强度改写为式(32)：
[0115][0116]jk|k-1
＝j
k-1
j
γ,k
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(33)
[0117][0118]
其中，和分别表示k时刻第l个预测目标高斯项的权重、均值和协方差矩阵
[0119]
②
更新强度
[0120][0121]
其中，和分别表示k时刻第l个更新高斯项的权重、均值和协方差矩阵；
[0122][0123][0124][0125][0126][0127][0128][0129]
其中，表示kde实时估计的k时刻目标幅度概率密度函数，表示kde实时估计的k 时刻杂波幅度概率密度函数，hk为观测矩阵，rk为测量噪声协方差矩阵，t表示转置；
[0130]
步骤(4)：利用建好的gm-kde-ai-phd滤波器实现多目标跟踪；
[0131]
本发明目标跟踪效果如图2、图3、图4所示；图2为跟踪场景，该场景包含了3个目标的真实轨迹和滤波轨迹；图3、图4是在200次蒙特卡洛仿真后得出的目标估计个数平均值比较图和ospa平均值比较图，由图3可以看出，本发明方法的目标估计个数较瑞利分布gm-ai-phd 滤波器(基于参数估计的gm-ai-phd)、gm-phd滤波器更加准确，且与真实分布gm-ai-phd 滤波器相差较小，由图4可以看出，本发明方法的ospa平均值比瑞利分布gm-ai-phd滤波器、 gm-phd滤波器的ospa平均值低，且与真实分布gm-ai-phd滤波器相差较小；因此，本发明算法与基于参数估计的gm-ai-phd、无幅度gm-phd相比，能够对多目标进行实时准确跟踪。