CN106556818A - 一种用于单目标跟踪的低计算复杂度贝努利滤波器 - Google Patents

Abstract

本发明公开了一种用于单目标跟踪的低计算复杂度贝努利滤波器，简化了整体的贝叶斯推导过程，优化了目标状态转移函数与目标量测似然函数。通过对目标量测与状态独立杂波的统计特性进行分析，推导出改进的序贯蒙特卡罗实现过程，可有效提取出杂波环境下的真实目标。本发明的滤波器可有效跟踪杂波环境下的机动单目标，并具有较低的计算复杂度与较高的跟踪精度。

Description

一种用于单目标跟踪的低计算复杂度贝努利滤波器

技术领域

本发明涉及一种滤波器，具体涉及一种用于单目标跟踪的低计算复杂度贝努利滤波器。

背景技术

单目标跟踪是对杂波环境下移动目标运动状态的实时估计过程。随着随机有限集理论的不断发展，单目标跟踪技术有效地摆脱了复杂数据关联操作的束缚，广泛应用于民用与军用等各类监控领域。

近年来，国内外学者深入研究了单目标跟踪技术并取得了大量的优秀成果，许多基于随机有限集理论的目标跟踪学术论文陆续出版在重要的国际期刊上。贝努利滤波器主要是假定目标运动状态为贝努利随机有限集。为了得到近似解，常采用序贯蒙特卡罗方法予以实现。在现有文献中，B.Ristic等学者首先提出了贝努利滤波的数学框架，由采样粒子的运动规律推导出复杂动态系统的预测与更新公式，并应用于目标定位、检测与跟踪等领域。随后，现有文献讨论了一种用于杂波环境下纯方位角跟踪的贝努利滤波器，并结合信息理论判据给出了对应的序贯蒙特卡罗实现。现有文献提出了一种用于不确定量测条件下海事辐射源跟踪的贝努利滤波器，适用于对信源工作切换时间内的目标跟踪。现有文献描述了一种用于漏检条件下的贝努利滤波器，可用于杂波密集环境下的扩展目标检测与跟踪。目前，现有文献提出了一种新型的贝努利滤波器及其序贯蒙特卡罗实现。首先通过改进采样粒子建议密度分布函数降低预测方差，然后采用辅助粒子滤波策略，智能抽取鲁棒性强的粒子及标签来实现目标跟踪。然而，由于粒子滤波的固有缺陷，现有贝努利滤波器的序贯蒙特卡罗实现在预测步与更新步仍需执行很多繁杂的数学操作，存在较高的计算复杂度，在很大程度上限制了对实际目标的跟踪性能。因此，如何提高传统贝努利滤波器的跟踪性能已成为实际应用中的重要课题。

目前，传统的贝努利滤波器特征为：

假定k时刻的目标状态集存在于状态空间中，量测集存在于量测空间中，那么随机动态系统可表示为：

x_k＝F_k|k-1(x_k-1)+v_k-1 (1)

z_k＝h_k(x_k)+u_k (2)

其中，F_k|k-1(·)为基于目标当前运动状态x_k所演化的非线性转移函数，h_k(·)表示当前量测z_k与运动状态x_k之间的非线性关系，v_k-1和u_k分别为系统的过程噪声和量测噪声。同时，定义π_k|k-1(x_k|x_k-1)为由前一时刻运动状态x_k-1过渡到当前运动状态x_k的转移概率密度，g_k(z_k|x_k)为单目标量测的似然函数。

由随机有限集理论，贝努利滤波器可并行计算后验的目标运动状态空间概率密度函数s_k(x_k)＝Pr(x_k|Z_1：k)和后验的目标存在概率p_k＝Pr{|X_k|＝1|Z_1：k}。假定状态集X_k可由贝努利集对(p_k，s_k(x_k))描述，那么对应的后验的概率生成函数记为：

其中，第一个条件表明当前检测场景中仅有一个状态为x_k的目标，第二个条件表明当前检测场景中无目标存在。

于是，利用幸存目标概率p_s，k|k-1(x_k-1)、新生目标概率p_b，k|k-1和新生目标密度b_k|k-1(x_k)可定义如下的状态转移概率生成函数：

其中，X_k-1为k-1时刻的目标运动状态集，它的四个限制条件与式(3)所描述的相类似。

因此，预测的目标存在概率和空间密度函数可分别表示为：

p_k|k-1＝p_b，k|k-1(1-p_k-1)+p_s，k|k-1(x_k-1)p_k-1 (5)

其中，表示变量与ζ的内积。由式(6)可以得出，预测的空间密度函数s_k|k-1(x_k)主要由预测的目标存在概率p_k|k-1决定。

考虑到杂波随机有限集Z_k的势分布满足独立同分布，于是定义对应的概率生成函数为：

在式(7)中，当集合Z_k服从于泊松分布时，对应的势分布为：

其中，λ为杂波率。

另一方面，用杂波分布函数c_k(z_k|x_k)表示式(7)中的标准概率密度函数p(z_k)：

p(z_k)＝c_k(z_k|x_k) (9)

那么，可以得出：

其中，λc_k(z_k|x_k)为k时刻杂波过程。

假定p_D，k(x_k)为传感器的实际检测概率，由式(10)可得出量测似然函数的概率生成函数为：

其中，第一个公式中的左项表示杂波过程，右项则表示X_k＝{x_k}条件下的目标运动状态分量；而第二个公式仅含有条件下的杂波过程。

因此，更新的目标存在概率和空间密度函数可表示为：

其中，增量Δ_k定义为：

由式(12)和式(13)可知，在更新过程中，函数s_k(x_k)与p_k随时间推移相互耦合，其递推都依赖于变量Δ_k。因此，在传统的贝努利滤波器中，函数s_k(x_k)具有较高的计算复杂度。

发明内容

为解决上述问题，本发明提供了一种用于单目标跟踪的低计算复杂度贝努利滤波器。

为实现上述目的，本发明采取的技术方案为：

一种用于单目标跟踪的低计算复杂度贝努利滤波器，分别采用三个随机有限集来表示k时刻的量测方程，即主要目标产生的量测T_k(x_k)、可疑目标产生的量测S_k(x_k)和状态独立杂波产生的量测C_k；

Z_k＝T_k(x_k)∪S_k(x_k)∪C_k (1)

其中，T_k(x_k)项定义如下：

可以看出，由状态x_k产生的主要量测所对应的概率为1-p_D，k(x_k)。

为了表述方便，再将非主要量测z_k所对应的集合S_k(x_k)与集合C_k统一为：

K_k(x_k)＝S_k(x_k)∪C_k (3)

这里，K_k(x_k)为两个统计独立的随机有限集的并集，对应的强度函数为：

v_K，k(z_k|x_k)＝v_S，k(z_k)+v_C，k(z_k|x_k) (4)

其中，v_S，k(·|x_k)和V_C，k(·)分别为S_k(x_k)和C_k的强度函数。

对于状态的独立x_k，每个量测z_k的概率密度都服从于以下独立同分布：

优选地，贝努利滤波器的序贯蒙特卡罗实现过程包括如下步骤：

预测步：假定k-1时刻后验的目标存在概率密度可由权重粒子集近似：

其中，δ(·)表示狄拉克函数，L_k-1为所需的粒子数，为第i个粒子的运动状态，为该粒子的归一化权值，且满足条件

更新步：假定q_k(x_k|x_k-1，Z_k)为k时刻采样粒子的建议概率密度：

那么，目标存在概率密度p_k(x_k|Z_1：k)可由一组新的粒子集近似：

其中，L_k|k-1为更新的粒子数，对应粒子的更新权值为：

为提升运算效率，这里采用建议概率函数q_k(x_k|x_k-1，Z_k)用于表示状态转移函数于是，式(9)可简化为：

重采样与状态更新步：为解决序贯蒙特卡罗方法中的粒子退化问题，在粒子集中重采样L_k个粒子，其中第i个粒子的权值为那么新得到的粒子可表示为最后，可得出目标的状态估计：

本发明具有以下有益效果：

本发明的滤波器可有效跟踪杂波环境下的机动单目标，并具有较低的计算复杂度与较高的跟踪精度。

附图说明

图1为本发明实施例中目标运动轨迹与量测示意图。

图2为本发明实施例中x坐标下的目标运动轨迹估计图。

图3为本发明实施例中y坐标下的目标运动轨迹估计图。

图4为本发明实施例中目标数目估计图。

图5为本发明实施例1阶优化子模式分配距离图。

图6为本发明实施例中计算复杂度。

具体实施方式

为了使本发明的目的及优点更加清楚明白，以下结合实施例对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本发明实施例提供了一种用于单目标跟踪的低计算复杂度贝努利滤波器，其特征在于，分别采用三个随机有限集来表示k时刻的量测方程，即主要目标产生的量测T_k(x_k)、可疑目标产生的量测S_k(x_k)和状态独立杂波产生的量测C_k；

Z_k＝T_k(x_k)∪S_k(x_k)∪C_k (1)

其中，T_k(x_k)项定义如下：

可以看出，由状态x_k产生的主要量测所对应的概率为1-p_D，k(x_k)。

为了表述方便，再将非主要量测z_k所对应的集合S_k(x_k)与集合C_k统一为：

K_k(x_k)＝S_k(x_k)∪C_k (3)

这里，K_k(x_k)为两个统计独立的随机有限集的并集，对应的强度函数为：

v_K，k(z_k|x_k)＝v_S，k(z_k)+V_C，k(z_k|x_k) (4)

其中，v_S，k(·|x_k)和V_C，k(·)分别为S_k(x_k)和C_k的强度函数。

对于状态的独立x_k，每个量测z_k的概率密度都服从于以下独立同分布：

贝努利滤波器的序贯蒙特卡罗实现过程包括如下步骤：

预测步：假定k-1时刻后验的目标存在概率密度可由权重粒子集近似：

其中，δ(·)表示狄拉克函数，L_k-1为所需的粒子数，为第i个粒子的运动状态，为该粒子的归一化权值，且满足条件

更新步：假定q_k(x_k|x_k-1，Z_k)为k时刻采样粒子的建议概率密度：

那么，目标存在概率密度p_k(x_k|Z_1：k)可由一组新的粒子集近似：

其中，L_k|k-1为更新的粒子数，对应粒子的更新权值为：

为提升运算效率，这里采用建议概率函数q_k(x_k|x_k-1，Z_k)用于表示状态转移函数于是，式(9)可简化为：

重采样与状态更新步：为解决序贯蒙特卡罗方法中的粒子退化问题，在粒子集中重采样L_k个粒子，其中第i个粒子的权值为那么新得到的粒子可表示为最后，可得出目标的状态估计：

注释：假定无可疑目标产生的量测，即当S_k(x_k)≡0或v_S，k(z_k)＝0条件满足时，式(3)和(4)可简化成单目标贝叶斯公式，表示当前场景中无新生目标出现或无原始目标消失，只有随机动态模型所描述的当前目标。这时，式(4)可记为只含有杂波分量的简化形式v_K，k(z_k|x_k)＝v_C，k(z_k|x_k)，式(5)可以简化成标准的归一化公式。

根据贝叶斯滤波原理，在预测步求得：

s_k|k-1(x_k|Z_1：k-1)＝<π_k|k-1(x_k|x_k-1)，s_k-1(x_k-1|Z_1：k-1)> (12)

其中，s_k-1(x_k-1|Z_1：k-1)表示k-1时刻的后验概率密度函数。

在更新步，利用当前量测集Z_k可以得出更新的后验概率密度函数：

其中，若后验的概率分布函数s_k-1(x_k-1|Z_1：k-1)已知，那么如何寻求优化的量测似然函数φ_k(Z_k|x_k)是首要解决的问题。

假定贝努利推导过程至多适用于一个目标产生的量测，那么有p_b，k|k-1＝0，于是式(5)可简化为：

p_k|k-1＝p_s，k|k-1(x_k)p_k-1 (14)

可以看出，式(14)中只含有幸存目标的存在项，而无新生目标的出生项。

通常，在实际的目标检测过程中，较低的检测概率是不期望的。那么，当检测概率为较为理想并假定p_D(x_k)→1时，漏检概率则为1-p_D(x_k)→0。这时，利用可化简s_k(x_k)为：

同时，Δ_k可改写为：

接下来，用项分别乘以式(15)中的分子与分母，并由式(13)得出如下公式：

其中，φ_k(Z_k|x_k)为给定状态x_k与量测集Z_k下概率分布的牛顿-莱布尼兹推导：

当存在漏检分量时，其完整表达式为：

序贯蒙特卡罗实现

根据上述滤波原理，可推导出改进贝努利滤波器的序贯蒙特卡罗实现过程：

预测步：假定k-1时刻后验的目标存在概率密度可由权重粒子集近似：

其中，δ(·)表示狄拉克函数，L_k-1为所需的粒子数，为第i个粒子的运动状态，为该粒子的归一化权值，且满足条件

更新步：假定q_k(x_k|x_k-1，Z_k)为k时刻采样粒子的建议概率密度：

那么，目标存在概率密度p_k(x_k|Z_1：k)可由一组新的粒子集近似：

其中，L_k|k-1为更新的粒子数，对应粒子的更新权值为：

为提升运算效率，这里采用建议概率函数q_k(x_k|x_k-1，Z_k)用于表示状态转移函数于是，式(23)可简化为：

重采样与状态更新步：为解决序贯蒙特卡罗方法中的粒子退化问题，在粒子集中重采样L_k个粒子，其中第i个粒子的权值为那么新得到的粒子可表示为最后，可得出目标的状态估计：

注释2：由贝叶斯滤波原理，可以看出，式(24)表示了更新粒子权值的核估计。当检测概率为趋近于上界为1的固定常数时，可继续化简式(18)并用于计算于是有：

计算复杂度分析

假定主要目标产生量测当前时刻运动状态x_k和前一时刻运动状态x_k-1已知，那么函数和的计算复杂度则分别为O(α)和O(β)。由于单个传感器能处理检测场景中的n个目标，若每个目标需要L_k个粒子，那么传统贝努利滤波器的计算复杂度可表示为：

O(n(α+n！β)L_k) (27)

对改进的贝努利滤波器来说，只存在幸存目标项，而不存在新生目标项。因此，函数的计算复杂度降低了一半。接下来，对函数来说，由于主要目标量测的数目为1，其计算复杂度降至nβ/|Z_k|。因此，总计算复杂度可记为：

可以看出，随着λ取值的逐渐增大，改进贝努利滤波器的总计算复杂度趋于O(αL_k/2)，而同样条件下传统滤波器的计算复杂度则为O((α+β)L_k)。因此，改进贝努利滤波器的总计算效率得到了明显的提升。

实验仿真与结果分析

为了验证低计算复杂度贝努利滤波器的实际跟踪性能，本节讨论并分析了一个典型的机动单目标跟踪实验。

实验环境

在[-2000，2000]×[0，2000]m²的半圆形监测区域内跟踪机动单目标的实际运动状态。其中，被动传感器定位在坐标原点(0，0)m，检测概率为98.5％，检测周期为60s，采样时间为1s。假定机动目标由起始点(10，10)m出发在整个检测周期内执行逆时针协同转弯运动，转弯速率为0.01rad s^-1，速度为(-10，30)m s^-1，那么随机动态方程可表示为：

其中，噪声向量v_k-1和u_k的标准差分别为diag(100，10，100，10，1)和diag(2π/180，10)，diag(·)表示对角矩阵。此外，设定现有的目标幸存概率为99％，杂波率为10。

为了评估跟踪性能，采用1阶优化子模式分配距离由1000次蒙特卡罗实验比较改进滤波器与传统滤波器。若和分别为原始状态集和估计状态集，为内的排列集，c(c＞0)为截断系数，并由更新粒子权值决定势分布误差与本地误差：

为了有效评价检测过程中的目标过估与漏检现象，本实验中设置c＝100m，达到有效区分集合X与是否匹配的目的。

结果分析

具体的实验结果如图1-图6所示。首先，图1描述了检测区域内该目标的实际运动轨迹与实际量测。可以看出，在含有随机杂波的实际环境下，曲线运动轨迹表明该目标正在执行协同转弯运动。

图2和图3分别给出了x和y坐标下的目标真实运动轨迹、实际量测和两类滤波器的估计结果。可以看出，尽管两类滤波器都能估计目标的位置，但是本文提出的滤波器却给出了更为稳定的位置估计。相比而言，传统滤波器给出的位置估计在整个检测周期内具有较大的偏差，并且在第40s时发生了虚警。

图4绘出了两类滤波器对目标数的估计曲线。可以看出，本文提出的滤波器产生了与真实目标数相一致的势估计。然而，由于将目标运动轨迹附近的杂波误认为真实目标，传统滤波器在第40s时过估了一个目标。

图5比较了两类滤波器的1阶优化子模式分配距离。可以看出，传统的滤波器具有较大的距离误差，而本文的滤波器在整个检测周期内，估计位置与目标真实位置相接近，其误差低于传统滤波器。尤其是在强峰值出现的第40s时，发生了错误的势估计。

图6描绘了不同杂波率条件下两类滤波器的计算复杂度。可以看出，随着杂波率λ的增大，本文提出的滤波器计算复杂度曲线具有单调递减特性。相反地，传统滤波器的计算复杂度较大。例如，在杂波率为λ＝10的条件下，传统滤波器采用80个采样粒子所需的复杂度就达到了O(3818.18)。对于同样的复杂度，本文提出的滤波器可应用到200个粒子，即计算复杂度减小了近150％。

由以上实验结果与数值分析可以得出，相比于传统的贝努利滤波器，本文提出的滤波器具有良好的综合跟踪性能，无论是目标跟踪精度还是实际计算效率都得到了明显的改善以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以作出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (2)

1.一种用于单目标跟踪的低计算复杂度贝努利滤波器，其特征在于，分别采用三个随机有限集来表示k时刻的量测方程，即主要目标产生的量测T_k(x_k)、可疑目标产生的量测S_k(x_k)和状态独立杂波产生的量测C_k；

Z_k＝T_k(x_k)∪S_k(x_k)∪C_k (1)

其中，T_k(x_k)项定义如下：

可以看出，由状态x_k产生的主要量测所对应的概率为1-p_D，k(x_k)。

为了表述方便，再将非主要量测z_k所对应的集合S_k(x_k)与集合C_k统一为：

K_k(x_k)＝S_k(x_k)∪C_k (3)

这里，K_k(x_k)为两个统计独立的随机有限集的并集，对应的强度函数为：

v_K，k(z_k|x_k)＝v_S，k(z_k)+v_C，k(z_k|x_k) (4)

其中，v_S，k(·|x_k)和v_C，k(·)分别为S_k(x_k)和C_k的强度函数。

而对于状态的独立x_k，每个量测z_k的概率密度都服从于以下独立同分布：

$c_k\left(z_k\right)=c_k\left(z_k\right|x_k)=\frac{v_{K,k}\left(z_k\right|x_k)}{<v_{K,k}\left(z_k\right|x_k),1>}---\left(5\right)$

2.如权利要求1所述的一种用于单目标跟踪的低计算复杂度贝努利滤波器，其特征在于，贝努利滤波器的序贯蒙特卡罗实现过程包括如下步骤：

预测步：假定k-1时刻后验的目标存在概率密度可由权重粒子集近似：

$p_{k-1}\left(x_{k-1}\right|Z_{1:k-1})\stackrel{\&CenterDot;}{=}\underset{i=1}{\overset{L_{k-1}}{\&Sigma;}}{w}_{k-1}^{\left(i\right)}\mathrm{\&delta;}(x_{k-1}-x_{k-1}^{\left(i\right)})---\left(6\right)$

其中，δ(·)表示狄拉克函数，L_k-1为所需的粒子数，为第i个粒子的运动状态，为该粒子的归一化权值，且满足条件

更新步：假定q_k(x_k|x_k-1，Z_k)为k时刻采样粒子的建议概率密度：

$x_k^{\left(i\right)}~q_k\left(x_k\right|x_{k-1},Z_k)---\left(7\right)$

那么，目标存在概率密度p_k(x_k|Z_1:k)可由一组新的粒子集近似：

$p_k\left(x_k\right|Z_{1:k})\stackrel{\&CenterDot;}{=}\underset{i=1}{\overset{L_{k|k-1}}{\&Sigma;}}{\tilde{w}}_k^{\left(i\right)}\mathrm{\&delta;}(x_k-x_k^{\left(i\right)})---\left(8\right)$

其中，L_k|k-1为更新的粒子数，对应粒子的更新权值为：

${\tilde{w}}_k^{\left(i\right)}=\frac{\frac{{\mathrm{\&phi;}}_k\left(Z_k\right|x_k^{\left(i\right)}){\mathrm{\&pi;}}_{k|k-1}\left(x_k^{\left(i\right)}\right|x_{k-1}^{\left(i\right)})}{q_k\left(x_k^{\left(i\right)}\right|x_k^{\left(i\right)},Z_k)}{w}_{k-1}^{\left(i\right)}}{\underset{i=1}{\overset{L_{k|k-1}}{\&Sigma;}}\frac{{\mathrm{\&phi;}}_k\left(Z_k\right|x_k^{\left(i\right)}){\mathrm{\&pi;}}_{k|k-1}\left(x_k^{\left(i\right)}\right|x_{k-1}^{\left(i\right)})}{q_k\left(x_k^{\left(i\right)}\right|x_k^{\left(i\right)},Z_k)}{w}_{k-1}^{\left(i\right)}}---\left(9\right)$

为提升运算效率，这里采用建议概率函数q_k(x_k|x_k-1，Z_k)用于表示状态转移函数于是，式(9)可简化为：

${\tilde{w}}_k^{\left(i\right)}=\frac{{\mathrm{\&phi;}}_k\left(Z_k\right|x_k^{\left(i\right)})w_{k-1}^{\left(i\right)}}{\underset{i=1}{\overset{L_{k|k-1}}{\&Sigma;}}{\mathrm{\&phi;}}_k\left(Z_k\right|x_k^{\left(i\right)})w_{k-1}^{\left(i\right)}}---\left(10\right)$

重采样与状态更新步：为解决序贯蒙特卡罗方法中的粒子退化问题，在粒子集中重采样L_k个粒子，其中第i个粒子的权值为那么新得到的粒子可表示为最后，可得出目标的状态估计：

${\hat{x}}_k=\underset{i=1}{\overset{L_k}{\&Sigma;}}{w}_k^{\left(i\right)}{x}_k^{\left(i\right)}---\left(11\right)$