-
目标跟踪在军事和民用领域均具有广泛的应用[1-3], 而在现实工程中, 例如, 在海上目标跟踪[4]以及空中目标跟踪[5]等领域, 观测设备在获取与目标相关的真实量测同时会获取与目标不相关的虚警数据. 虚警数据也被称为杂波[6], 由于无法在混有杂波信息的数据集中直接识别真实量测, 基于真实量测设计的滤波器[7]无法对目标进行跟踪. 此外密集的杂波信号会影响真实目标的回波信号, 导致雷达探测能力下降进而影响目标的跟踪精度. 因此如何提高杂波环境下目标跟踪精度是当前目标跟踪领域亟需解决的问题之一.
许多研究者均对杂波环境下的目标跟踪问题展开研究, 并取得了不错的研究成果. 对于单目标跟踪系统而言, 最近邻算法[8-10]以及概率数据关联算法(PDA)[11-14]是较为经典的解决方法. 最近邻算法通过构建距离指标函数选取与预测中心最近的数据作为真实量测, 然而当杂波距离预测中心较近时, 其往往会错误的将杂波选取为真实量测, 导致跟踪精度下降. PDA算法本质上是一种全邻算法, 其认为关联波门内的每一个数据均可作为真实量测, 并通过新息函数为每一个数据赋予相应的权重. 由于该方法相比于最近邻算法可降低杂波对状态更新的影响, 因此具有更高的目标跟踪精度. PDA算法的设计思想同样被应用于解决杂波环境下机动目标跟踪问题, 相关的算法包括交互多模型概率数据关联算法[15] (IMM-PDA)、联合交互多模型概率数据关联算法[16]以及联合交互多模型距离加权概率数据关联算法[17]等.
从贝叶斯统计原理的角度来看, 杂波环境下目标真实状态后验概率密度函数应为真实量测作为观测信息情况下的状态后验概率密度函数. 由于无法直接识别真实量测, 因此目标真实状态后验概率密度函数无法获取并且每一个数据作为真实量测而获取的状态后验概率密度函数均可能为真实状态后验概率密度函数. 此外, 当检测概率不为1时, 所有数据与目标不相关时获取的状态后验概率密度函数也可能成为真实状态后验概率密度函数. 代理概率密度函数表示某一变量可能的概率分布, 所有数据关联形式下获取的状态后验概率密度函数可构成一个状态代理概率密度函数集合. 基于该集合的统计特征以实现状态近似后验概率密度函数的获取是一个合理的选择. PDA算法通过提取每一个数据-状态后验概率密度函数的均值与协方差并基于加权融合的方式以实现状态近似后验概率密度函数的获取. 这类近似方法虽然简单且易于实现, 但是在数据数目较多的情形下会存在计算效率低的缺陷. 本文在变分贝叶斯框架[18]下设计了一种新的概率数据关联算法(VB-PDA), 该算法首先将关联事件视为随机变量并利用多项分布对两者进行建模, 随后基于量测集、目标状态、关联事件的联合概率密度函数求取关联事件的后验概率密度函数, 最后将关联事件的后验概率密度函数引入变分贝叶斯框架中以获取状态近似后验概率密度函数. VB-PDA算法本质上是一种在PDA建模框架下通过权重KL平均[19]准则来实现状态近似后验概率密度函数获取的方法. 相比于PDA算法, VB-PDA算法在提高计算效率的同时, 获取的近似后验概率密度函数也更能体现状态代理概率密度函数集的统计特征. 相关仿真实验结果验证了算法的有效性.
本文内容安排如下: 第1节对相关问题进行了描述; 第2节对PDA算法进行了介绍以及对VB-PDA算法进行了推导; 第3节从权重KL平均角度对VB-PDA算法进行了分析; 第4节对算法计算复杂度进行了分析; 第5节通过相关仿真实验对提出的算法进行验证; 第6节对本文进行了总结.
-
假设线性离散系统状态空间模型的状态方程和量测方程分别为:
$${{{x}}_{k + 1}} = {{{F}}_k}{{{x}}_k} + {{{w}}_k}$$ (1) $${{{z}}_{k + 1}} = {{{H}}_{k + 1}}{{{x}}_{k + 1}} + {{{v}}_{k + 1}}$$ (2) 其中,
${{{x}}_{k + 1}}$ 为目标在k+1时刻的状态信息,${{{F}}_k}$ 为状态转移矩阵,${{{z}}_{k + 1}}$ 为真实量测,${{{H}}_{k + 1}}$ 为量测矩阵,${{{x}}_{k + 1}}$ 和${{{z}}_{k + 1}}$ 的维数分别为$r$ 和$m$ ,$r > m$ .${{{w}}_k}$ 和${{{v}}_{k + 1}}$ 为离散系统的过程噪声与量测噪声,${{{w}}_k}$ 和${{{v}}_{k + 1}}$ 均服从高斯分布.$$p({{{w}}_k}) = {\cal{N}}({{{w}}_k};{0^{r \times 1}},{{{Q}}_k})$$ (3) $$p({{{v}}_{k + 1}}) = {\cal{N}}({{{v}}_{k + 1}};{0^{m \times 1}},{{{R}}_{k + 1}})$$ (4) 其中,
$p( \bullet )$ 表示概率密度函数,${\cal{N}}( \bullet ;{{d}},\sum )$ 表示高斯分布,${{d}}$ 和$\sum $ 表示变量的均值以及协方差.${{{Q}}_k}$ 为过程噪声协方差,${{{R}}_{k + 1}}$ 为量测噪声协方差. 考虑系统$k + 1$ 时刻获取的数据集为${{{Z}}_{k + 1}} = {\big\{} {{{z}}_{k + 1,1}}, {{{z}}_{k + 1,2}},{{{z}}_{k + 1,3}}, \cdots {{{z}}_{k + 1,{n_{k + 1}}}} {\big\}}$ ,${n_{k + 1}}$ 为数据集元素数量. 数据集由目标量测和与目标不相关的杂波数据组成. 假设与目标相关联的数据个数至多为1, 状态的后验概率密度函数$p({{{x}}_{k + 1}}\left| {{{{Z}}_{1:k + 1}}} \right.)$ 可表示为:$$\begin{array}{l} p({{{x}}_{k + 1}}\left| {{{{Z}}_{1:k + 1}}} \right.)=\\ \left\{ \begin{array}{l} p({{{x}}_{k + 1}}\left| {{{{Z}}_{1:k}}} \right.){\rm{ }}{\text{所有数据均为杂波}}\\ p({x_{k + 1}}\left| {{{{Z}}_{1:k}},{{{z}}_{k + 1,i}}} \right.){\rm{ }}{{{z}}_{k + 1,i}}{\text{与目标关联}}{\rm{ }} \end{array} \right. \end{array}$$ (5) 其中,
$p({{{x}}_{k + 1}}\left| {{{{Z}}_{1:k}},{{{z}}_{k + 1,i}}} \right.)$ 表示${{{z}}_{k + 1,i}}$ 为真实量测情况下目标状态后验概率密度函数. 根据贝叶斯公式,$p({{{x}}_{k + 1}}\left| {{{{Z}}_{1:k}},{{{z}}_{k + 1,i}}} \right.)$ 可表示为$$\begin{split} &p({{{x}}_{k + 1}}\left| {{{{Z}}_{1:k}},{{{z}}_{k + 1,i}}} \right.) = \\ &\dfrac{{p({{{z}}_{k + 1,i}}\left| {{{{x}}_{k + 1}},{{{Z}}_{1:k}}} \right.)p({{{x}}_{k + 1}}\left| {{{{Z}}_{1:k}}} \right.)}}{{p({{{z}}_{k + 1,i}}\left| {{{{Z}}_{1:k}}} \right.)}}\end{split}$$ (6) 假设
$p({{{x}}_k}\left| {{{{Z}}_{1:k}}} \right.) = {\cal{N}}({{{x}}_k};{{{x}}_{k\left| k \right.}},{{{P}}_{k\left| k \right.}})$ , 由于系统的状态空间模型为隐马尔可夫模型, 因此,$p({{{z}}_{k + 1,i}}\left| {{{{x}}_{k + 1}},{{{Z}}_{1:k}}} \right.)$ 和$p({{{x}}_{k + 1}}\left| {{{{Z}}_{1:k}}} \right.)$ 分别为$$\begin{split} & p({{{z}}_{k + 1,i}}\left| {{{{x}}_{k + 1}},{{{Z}}_{1:k}}} \right.) = \\ & p({{{z}}_{k + 1,i}}\left| {{{{x}}_{k + 1}}} \right.) \propto {\cal{N}}({{{z}}_{k + 1,i}};{{{H}}_{k + 1}}{{{x}}_{k + 1}},{{{R}}_{k + 1}}) \end{split} $$ (7) $$p({{{x}}_{k + 1}}\left| {{{{Z}}_{1:k}}} \right.) = {\cal{N}}({{{x}}_{k{\rm{ + 1}}}};{{{x}}_{k + 1\left| k \right.}},{{{P}}_{k + 1\left| k \right.}})$$ (8) $${{{x}}_{k + 1\left| k \right.}} = {{{F}}_k}{{{x}}_{k\left| k \right.}}$$ (9) $${{{P}}_{k + 1\left| k \right.}} = {{{F}}_k}{{{P}}_{k\left| k \right.}}{{F}}_k^{\rm{T}} + {{{Q}}_k}$$ (10) 根据公式(7)和(8),
$p({{{x}}_{k + 1}}\left| {{{{Z}}_{1:k}},{{{z}}_{k + 1,i}}} \right.)$ 可被计算为:$$ \begin{split} &p({{{x}}_{k + 1}}\left| {{{{Z}}_{1:k}},{{{z}}_{k + 1,i}}} \right.) =\\ &{\cal{N}}({{{x}}_{k{\rm{ + 1}}}};{{{x}}_{k + 1\left| {k + 1,i} \right.}},{{{P}}_{k + 1\left| {k + 1,i} \right.}}) \end{split}$$ (11) $$\begin{split} {{{x}}_{k + 1\left| {k + 1,i} \right.}} =& {{{x}}_{k + 1\left| k \right.}} +\\ &{{{K}}_{k + 1}}({{{z}}_{k + 1,i}} - {{{H}}_{k + 1}}{{{x}}_{k + 1\left| k \right.}})\end{split}$$ (12) $${{{P}}_{k + 1\left| {k + 1,i} \right.}} = {{{P}}_{k + 1\left| k \right.}} - {{{K}}_{k + 1}}{{{H}}_{k + 1}}{{{P}}_{k + 1\left| k \right.}}$$ (13) $${{{K}}_{k + 1}} = {{{P}}_{k + 1\left| k \right.}}{{H}}_{k + 1}^{\rm{T}}{({{{H}}_{k + 1}}{{{P}}_{k + 1\left| k \right.}}{{H}}_{k + 1}^{\rm{T}} + {{{R}}_{k + 1}})^{ - 1}}$$ (14) 由于无法直接区分真实量测与杂波, 状态真实后验概率密度函数无法获取. 如何获取合适的状态近似后验概率密度函数是本文的研究目的.
-
不同的杂波分布会导致算法具备不同的数学形式. 假设目标的检测概率为
${p_d}$ , 门概率为${p_g}$ . 系统视场${{\varTheta }}$ 被视为关联门, 杂波在${{\varTheta }}$ 内服从均匀分布, 该区域体积为${V_{k + 1}}$ . 杂波数目服从泊松分布, 杂波密度为${\lambda _{k + 1}}$ .${V_{k + 1}}$ ,${p_d}$ ,${\lambda _{k + 1}}$ 被认为是已知的.${p_g}$ 表示真实量测在视场内的概率. 其可由如下公式获取.$$\begin{split} {p_g} =& \int_{{{{z}}_{k + 1}} \in {{\varTheta }},{{{x}}_{k{\rm{ + 1}}}} \in {{R}}} {{\cal{N}}({{{z}}_{k + 1}};{{{H}}_{k + 1}}{{{x}}_{k + 1}},{{{R}}_{k + 1}})} \times \\ &{\cal{N}}({{{x}}_{k{\rm{ + 1}}}};{{{x}}_{k + 1\left| k \right.}},{{{P}}_{k + 1\left| k \right.}}){\rm{d}}{{{z}}_{k + 1}}{\rm{d}}{{{x}}_{k{\rm{ + 1}}}} = \\ & \int_{{z_{k + 1}} \in {{\varTheta }}} {{\cal{N}}({{{z}}_{k + 1}};{{{H}}_{k + 1}}{{{x}}_{k + 1\left| k \right.}},}\\ & {{{H}}_{k + 1}}{{{P}}_{k + 1\left| k \right.}}{{H}}_{k + 1}^T+ {{{R}}_{k + 1}}){\rm{d}}{{{z}}_{k + 1}} \\[-10pt] \end{split} $$ (15) 其中
${{R}}$ 表示实数域. 本文介绍的PDA算法以及设计的VB-PDA算法均是基于上述杂波环境给出的. -
PDA算法基于高斯分布
${\cal{N}}({{{x}}_{k{\rm{ + 1}}}};{{{x}}_{k + 1\left| {k + 1} \right.}}, {{{P}}_{k + 1\left| {k + 1} \right.}})$ 获取目标状态近似后验概率密度函数. 均值${{{x}}_{k + 1\left| {k + 1} \right.}}$ 与协方差${{{P}}_{k + 1\left| {k + 1} \right.}}$ 从状态代理概率密度函数集中获取, 相应的公式分别为:$$\begin{split} {{{x}}_{k + 1\left| k \right.{\rm{ + 1}}}} =& {{{x}}_{k + 1\left| k \right.}} + \left(\sum\limits_{i = 0}^{{n_{k + 1}}} {\Pr ({{{z}}_{k + 1}} = {{{z}}_{k + 1,i}})} \right)\\ &{{{K}}_{k + 1}}\left(\frac{{\sum\limits_{i = 1}^{{n_{k + 1}}} {\Pr ({{{z}}_{k + 1}} = {{{z}}_{k + 1,i}}){{{z}}_{k + 1,i}}} }}{{\sum\limits_{i = 1}^{{n_{k + 1}}} {\Pr ({{{z}}_{k + 1}} = {{{z}}_{k + 1,i}})} }} -\right.\\ & {{{H}}_{k + 1}}{{{x}}_{k + 1\left| k \right.}})\\[-10pt]\end{split}$$ (16) $$\begin{split} & {{{P}}_{k + 1\left| k \right.{\rm{ + 1}}}} = {{{P}}_{k + 1\left| k \right.}} - \left(\sum\limits_{i = 1}^{{n_{k + 1}}} {\Pr ({{{z}}_{k + 1}} = {{{z}}_{k + 1,i}})} \right)\\ & {{{K}}_{k + 1}}{{{H}}_{k + 1}}{{{P}}_{k + 1\left| k \right.}} + \sum\limits_{i = {\rm{0}}}^{{n_{k + 1}}} {\big(\Pr ({{{z}}_{k + 1}} = {{{z}}_{k + 1,i}})}\\ &({{{x}}_{k + 1\left| k \right.{\rm{ + 1}},i}} - {{{x}}_{k + 1\left| k \right.{\rm{ + 1}}}}){({{{x}}_{k + 1\left| k \right.{\rm{ + 1}},i}} - {{{x}}_{k + 1\left| k \right.{\rm{ + 1}}}})^{\rm{T}}}\big) \end{split} $$ (17) 其中
${{{x}}_{k + 1\left| k \right.{\rm{ + 1}},0}} = {{{x}}_{k + 1\left| k \right.}}$ .$\Pr ({{{z}}_{k + 1}} = {{{z}}_{k + 1,i}}){\rm{ }}$ ${\rm{(}}i = {\rm{0}}, 1 \cdots {n_{k + 1}}{\rm{)}}$ 为相应状态代理概率密度函数权重, 当$i = {\rm{0}}$ 时表示所有数据均为杂波.$\Pr ({{{z}}_{k + 1}} = {{{z}}_{k + 1,0}})$ 和$\Pr ({{{z}}_{k + 1}} = {{{z}}_{k + 1,i}})$ ${\rm{ (}}i = 1,2 \cdots {n_{k + 1}}{\rm{)}}$ 的计算公式如下.$$\Pr ({{{z}}_{k + 1}} = {{{z}}_{k + 1,0}}) = \frac{{1 - {p_d}{p_g}}}{{1 - {p_d}{p_g} + \displaystyle\sum\limits_{j = 1}^{{n_{k + 1}}} {{l_{k + 1,j}}} }}$$ (18) $$\Pr ({{{z}}_{k + 1}} = {{{z}}_{k + 1,i}}){\rm{ }} = \frac{{{l_{k + 1,i}}}}{{1 - {p_d}{p_g} + \displaystyle\sum\limits_{j = 1}^{{n_{k + 1}}} {{l_{k + 1,j}}} }}$$ (19) -
在VB-PDA算法中, 量测集中任意一个数据为真实量测的概率相同被视为先验信息. 此外,
${p_d}$ 、${p_g}$ 以及杂波相关信息同样被视为已知的先验信息并被使用构建相应的概率密度函数. 将关联事件视为随机变量, 由于关联事件的数目是有限的, 因此可用多项分布进行建模. 此外, 关联事件数目的确定意味着量测数目同样确定. 因此基于已有的先验信息, 关联事件的先验概率密度函数$p({\eta _{k + 1,0}}, {\eta _{k + 1,}} 1, \cdots {\eta _{k + 1,{n_{k + 1}}}}\left| {{{{Z}}_{1:k}}} \right.)$ 可被建模为$$\begin{split} & p\left({\eta _{k + 1,0}},{\eta _{k + 1,1}}, \cdots {\eta _{k + 1,{n_{k + 1}}}}\left| {{{{Z}}_{1:k}}} \right.\right) \propto \\ & {\left((1 - {p_d}{p_g}){e^{ - {\lambda _{k + 1}}{V_{k + 1}}}}\dfrac{{{{({\lambda _{k + 1}}{V_{k + 1}})}^{{n_{k + 1}}}}}}{{{n_{k + 1}}!}}\right)^{{\eta _{k + 1,0}}}}\times \\ & \prod\limits_{i = 1}^{{n_{k + 1}}} {{{\left(\dfrac{{{p_d}{p_g}}}{{{n_{k + 1}}}}{e^{ - {\lambda _{k + 1}}{V_{k + 1}}}}\dfrac{{{{({\lambda _{k + 1}}{V_{k + 1}})}^{{n_{k + 1}} - 1}}}}{{({n_{k + 1}} - 1)!}}\right)}^{{\eta _{k + 1,i}}}}} \propto \\ & {((1 - {p_d}{p_g}))^{{\eta _{k + 1,0}}}}\prod\limits_{i = 1}^{{n_{k + 1}}} {{{\left(\dfrac{{{p_d}{p_g}}}{{{n_{k + 1}}{\lambda _{k + 1}}{V_{k + 1}}}}\right)}^{{\eta _{k + 1,i}}}}} \\[-10pt] \end{split} \tag{21}$$ (21) $${l_{k + 1,i}} = \frac{{{p_d}{\cal{N}}({{{z}}_{k{\rm{ + 1}},i}};{{{H}}_{k + 1}}{{{x}}_{k + 1\left| k \right.}},{{{H}}_{k + 1}}{{{P}}_{k + 1\left| k \right.}}{{H}}_{k + 1}^T + {{{R}}_{k + 1}})}}{{{\lambda _{k + 1}}}}\tag{20}$$ (20) 其中,
${\eta _{k + 1,i}} \in \left\{ {{\rm{0,1}}} \right\}$ ,$\sum\nolimits_{i = 0}^{{n_{k + 1}}} {{\eta _{k + 1,i}}} = 1$ .${\eta _{k + 1,0}}$ 表示所有数据均为杂波情况下的关联事件.${\eta _{k + 1,i}}$ 表示${{{z}}_{k + 1,i}}$ 为真实量测情况下的关联事件. 根据状态空间模型的量测方程以及门概率信息, 似然函数$p({{{z}}_{k + 1,i}}\left| {{{{x}}_{k + 1}},{\eta _{k + 1,i}} = 1} \right.)$ ${\rm{(}}i > 0{\rm{)}}$ 的表达式为$$\begin{split} &p({{{z}}_{k + 1,i}}\left| {{{{x}}_{k + 1}},{\eta _{k + 1,i}} = 1} \right.){\rm{ }} = \\ &\dfrac{{{\cal{N}}({{{z}}_{k + 1,i}};{{{H}}_{k + 1}}{{{x}}_{k + 1}},{{{R}}_{k + 1}})}}{{{p_g}}}\end{split} $$ (22) 根据(22)以及杂波分布, 量测似然函数
$p ({{{Z}}_{k + 1}}\left| {{{{x}}_{k + 1}},{\eta _{k + 1,0}},{\eta _{k + 1,1}}, \cdots {\eta _{k + 1,n}}} \right.)$ 表达式为$$\begin{split} & p\left({{{Z}}_{k + 1}}\left| {{{{x}}_{k + 1}},{\eta _{k + 1,0}},{\eta _{k + 1,1}}, \cdots {\eta _{k + 1,n}}} \right.\right) = \\ & {\left(\dfrac{1}{{V_{k + 1}^{{n_{k + 1}}}}}\right)^{{\eta _{k + 1,0}}}}\prod\limits_{i = 1}^{{n_{k + 1}}} \Bigg(\dfrac{1}{{{p_g}V_{k + 1}^{{n_{k + 1}} - 1}}}\\ &{\cal{N}}\left({{{z}}_{k + 1,i}};{{{H}}_{k + 1}}{{{x}}_{k + 1}},{{{R}}_{k + 1}}\right)\bigg)^{{\eta _{k + 1,i}}} \end{split} $$ (23) 根据(8), (21), (23)以及条件独立原则,
$p({{{x}}_{k + 1}}, {\eta _{k + 1,0}},{\eta _{k + 1,1}}, \cdots {\eta _{k + 1,{n_{k + 1}}}},{{{Z}}_{1:k + 1}})$ 的表达式为$$\begin{split} & p({{{x}}_{k + 1}},{\eta _{k + 1,0}},{\eta _{k + 1,1}}, \cdots {\eta _{k + 1,n}},{{{Z}}_{1:k + 1}}) = \\ & p({{{Z}}_{k + 1}}\!\left| {{{{x}}_{k + 1}},{\eta _{k + 1,0}},{\eta _{k + 1,1}}, \cdots {\eta _{k + 1,n}}} \right.)p({{{x}}_{k + 1}}\!\left| {{{{Z}}_{1:k}}}\! \right.) \!\times\\ & p({\eta _{k + 1,0}},{\eta _{k + 1,1}}, \cdots {\eta _{k + 1,n}}\left| {{{{Z}}_{1:k}}} \right.)p({{{Z}}_{1:k}})\propto \\ & {((1 - {p_d}{p_g}){\cal{N}}({{{x}}_{k + 1}};{{{x}}_{k + 1\left| k \right.}},{{{P}}_{k + 1\left| k \right.}}))^{{\eta _{k + 1,0}}}} \times \\ & \prod\limits_{i = 1}^{{n_{k + 1}}} {\left(\dfrac{{{p_d}}}{{{\lambda _{k + 1}}}}{\cal{N}}({{{x}}_{k + 1}};{{{x}}_{k + 1\left| k \right.}},{{{P}}_{k + 1\left| k \right.}}\right)} \times \\ & {\cal{N}}({{{z}}_{k{\rm{ + 1}},i}};{{{H}}_{k + 1}}{{{x}}_{k + 1}},{{{R}}_{k + 1}}){)^{{\eta _{k + 1,i}}}} \\[-10pt] \end{split} $$ (24) 根据公式(24),
$\left\{ {{\eta _{k + 1,0}},{\eta _{k + 1,1}}, \cdots {\eta _{k + 1,{n_{k + 1}}}}} \right\}$ 的后验概率密度函数$p({\eta _{k + 1,0}},{\eta _{k + 1,1}}, \cdots \eta _{k + 1, {n_{k + 1}}}\left| {{{{Z}}_{1:k + 1}}} \right.)$ 表达式为$$\begin{split} & p({\eta _{k + 1,0}},{\eta _{k + 1,1}}, \cdots {\eta _{k + 1,{n_{k + 1}}}}\left| {{{{Z}}_{1:k + 1}}} \right.) \propto \\ & p({\eta _{k + 1,0}},{\eta _{k + 1,1}}, \cdots {\eta _{k + 1,{n_{k + 1}}}},{{{Z}}_{1:k + 1}}) =\\ & \int p({\eta _{k + 1,0}},{\eta _{k + 1,1}}, \cdots {\eta _{k + 1,{n_{k + 1}}}},{{{Z}}_{1:k + 1}},{{{x}}_{k + 1}})\\ & d{{{x}}_{k + 1}} \propto \prod\limits_{i = 1}^{{n_{k + 1}}} \bigg(\frac{{{p_d}}}{{{\lambda _{k + 1}}}}{\cal{N}}\big({{{z}}_{k{\rm{ + 1}},i}};{{{H}}_{k + 1}}{{{x}}_{k + 1\left| k \right.}},\\ & {{{H}}_{k + 1}}{{{P}}_{k + 1\left| k \right.}}{{H}}_{k + 1}^T + {{{R}}_{k + 1}}\big)\!{\bigg)^{{\eta _{k + 1,i}}}} \!\! \times\! {(1 \!\!-\!\! {p_d}{p_g})^{{\eta _{k + 1,0}}}} \Rightarrow \\ & p({\eta _{k + 1,0}},{\eta _{k + 1,1}}, \cdots {\eta _{k + 1,{n_{k + 1}}}}\left| {{{{Z}}_{1:k + 1}}} \right.) = \\ & {(\Pr ({{{z}}_{k + 1}}\! =\! {{{z}}_{k + 1,0}}))^{{\eta _{k + 1,0}}}}\prod\limits_{i = 1}^{{n_{k + 1}}} {{{(\Pr ({{{z}}_{k + 1}} \!=\! {{{z}}_{k + 1,i}}))}^{{\eta _{k + 1,i}}}}} \end{split} $$ (25) 可以发现基于本文建模框架下推导的关联事件权重模型与PDA算法相同. 利用VB算法[18]对状态近似后验密度函数
${q_{{\rm{VB}}}}({{{x}}_{k + 1}})$ 进行求解, 其求解公式如下$$\begin{split} & {q_{{\rm{VB}}}}({{{x}}_{k + 1}}) \propto \\ & \exp \!\! {\int \!\!{\ln p({{{x}}_{k + 1}},{\eta _{k + 1,0}},{\eta _{k + 1,1}}, \cdots {\eta _{k + 1,{n_{k + 1}}}},{{{Z}}_{1:k + 1}})} } \times\\ & {q_{{\rm{VB}}}}({\eta _{k + 1,0}},{\eta _{k + 1,1}}, \cdots {\eta _{k + 1,{n_{k + 1}}}}) \!\times \\ & \left. {{\rm{d}}\left\{ {{\eta _{k + 1,0}},{\eta _{k + 1,1}}, \cdots {\eta _{k + 1,{n_{k + 1}}}}} \right\}} \right\} \\[-10pt] \end{split} $$ (26) 其中,
${q_{{\rm{VB}}}}( \cdot )$ 表示变量在VB框架下获取的近似后验概率密度函数.${q_{{\rm{VB}}}}({\eta _{k + 1,0}},{\eta _{k + 1,1}}, \cdots {\eta _{k + 1,{n_{k + 1}}}})$ 满足:$$\begin{split} &{q_{{\rm{VB}}}}({\eta _{k + 1,0}},{\eta _{k + 1,1}}, \cdots {\eta _{k + 1,{n_{k + 1}}}})= \\ & p({\eta _{k + 1,0}},{\eta _{k + 1,1}}, \cdots {\eta _{k + 1,{n_{k + 1}}}}\left| {{{{Z}}_{1:k + 1}}} \right.) \end{split} $$ (27) 根据公式(26)和(27),
${q_{{\rm{VB}}}}({{{x}}_{k + 1}})$ 的表达式为$$\begin{split} & {q_{{\rm{VB}}}}({{{x}}_{k + 1}}) \propto \\ & \exp { - ({{{x}}_{k + 1}} - {{{x}}_{k + 1\left| k \right.}}){{({{{P}}_{k + 1\left| k \right.}})}^{ - 1}}{{({{{x}}_{k + 1}} \!-\! {{{x}}_{k + 1\left| k \right.}})}^{\rm{T}}}} \!- \\ & \sum\limits_{i = 1}^{{n_{k + 1}}} {\Pr ({{{z}}_{k + 1}} \!=\! {{{z}}_{k + 1,i}})({{{z}}_{k{\rm{ + 1}},i}} \!-\! {{{H}}_{k + 1}}{{{x}}_{k + 1}}){{({{{R}}_{k + 1}})}^{ \!-\! 1}}} \times\\ & { {{({{{z}}_{k{\rm{ + 1}},i}} - {{{H}}_{k + 1}}{{{x}}_{k + 1}})}^{\rm{T}}}} \propto \\ & {\cal{N}}\Bigg(\frac{{\sum\limits_{i = 1}^{{n_{k + 1}}} {\Pr ({{{z}}_{k + 1}} = {{{z}}_{k + 1,i}}){{{z}}_{k + 1,i}}} }}{{\sum\limits_{i = 1}^{{n_{k + 1}}} {\Pr ({{{z}}_{k + 1}} = {{{z}}_{k + 1,i}})} }};\\ & {{{H}}_{k + 1}}{{{x}}_{k + 1}},\frac{{{{{R}}_{k + 1}}}}{{\sum\limits_{i = 1}^{{n_{k + 1}}} {\Pr ({{{z}}_{k + 1}} = {{{z}}_{k + 1,i}})} }}\Bigg) \times \\ & {\cal{N}}({{{x}}_{k + 1}};{{{x}}_{k + 1\left| k \right.}},{{{P}}_{k + 1\left| k \right.}}) \Rightarrow\\ &{q_{{\rm{VB}}}}({{{x}}_{k + 1}}) = {\cal{N}}({{{x}}_{k + 1}};{{{x}}_{{\rm{VB}},k + 1\left| {k{\rm{ + 1}}} \right.}},{{{P}}_{{\rm{VB}},k + 1\left| {k{\rm{ + 1}}} \right.}}) \\[-10pt] \end{split} $$ (28) $$\begin{split} & {{{x}}_{{\rm{VB}},k + 1\left| {k{\rm{ + 1}}} \right.}} = {{{x}}_{k + 1\left| k \right.}} + {{{K}}_{{\rm{VB}},k + 1}}\\ &\left(\frac{{\sum\limits_{i = 1}^{{n_{k + 1}}} {\Pr ({{{z}}_{k + 1}} = {{{z}}_{k + 1,i}}){{{z}}_{k + 1,i}}} }}{{\sum\limits_{i = 1}^{{n_{k + 1}}} {\Pr ({{{z}}_{k + 1}} = {{{z}}_{k + 1,i}})} }} - {{{H}}_{k + 1}}{{{x}}_{k + 1\left| k \right.}}\right) \\ \end{split} $$ (29) $$\begin{split} {{{K}}_{{\rm{VB}},k + 1}} =& {{{P}}_{k + 1\left| k \right.}}{{H}}_{k + 1}^{\rm{T}}({{{H}}_{k + 1}}{{{P}}_{k + 1\left| k \right.}}{{H}}_{k + 1}^{\rm{T}} +\\ &\frac{{{{{R}}_{k + 1}}}}{{\sum\limits_{i = 1}^{{n_{k + 1}}} {\Pr ({{{z}}_{k + 1}} = {{{z}}_{k + 1,i}})} }}\Bigg)^{ - 1}\end{split}$$ (30) $${{{P}}_{{\rm{VB}},k + 1\left| {k{\rm{ + 1}}} \right.}} = {{{P}}_{k + 1\left| k \right.}} - {{{K}}_{{\rm{VB}},k + 1}}{{{H}}_{k + 1}}{{{P}}_{k + 1\left| k \right.}}$$ (31) 实际上, 基于推导过程可以发现, 状态后验概率密度函数在贝叶斯框架下可直接被获取, 并且该概率密度函数服从高斯混合分布. 然而随着状态更新的过程不断进行, 该后验概率密度函数的高斯分量会呈几何式增长, 导致其计算复杂度过大. 考虑到实时性的影响, 采用VB算法获取状态近似后验概率密度函数. VB-PDA算法获取的状态近似后验概率密度函数服从高斯分布, 这样的结果会为算法在后续的状态更新过程中提供便利. 此外, VB-PDA算法获取的状态协方差小于状态预测协方差以及PDA算法获取的状态协方差, 因此该算法可有效避免奇异现象的产生.
-
本部分从代理概率密度函数集的角度说明VB-PDA算法的优势. 权重KL平均
$\bar p$ 被认为是最能体现代理概率密度函数集统计特征的概率密度函数[19], 权重KL平均$\bar p(\Psi )$ 定义为:$$\bar p(\Psi ) = \arg \mathop {\min }\limits_{\bar p(\Psi )} \sum\limits_{i = 1}^l {{u_i}(\Psi ){\rm{KL}}(\bar p(\Psi )\left| {{p^i}(\Psi )} \right.)} $$ (32) $${\rm{KL}}(\bar p(\Psi )\left| {{p^i}} \right.(\Psi )) = \int {\bar p(\Psi )\ln \frac{{\bar p(\Psi )}}{{{p^i}(\Psi )}}} {\rm{d}}\Psi $$ (33) 其中
${\rm{KL(}} \cdot {\rm{)}}$ 为KL散度函数,${p^i}(\Psi )$ 为第i个代理概率密度函数,${u_i}(\Psi )$ 为第i个代理概率密度函数的权重,$\Psi $ 为待估计参数集,$l$ 为代理概率密度函数数目. 对于PDA算法而言,$l = {n_{k + 1}} + 1$ , 其待估计参数集为${{{x}}_{k + 1}}$ .${p^i}({{{x}}_{k + 1}})$ 和${u_i}({{{x}}_{k + 1}})$ 的表达式分别为$${p^i}({{{x}}_{k + 1}}) \!=\! \left\{ \begin{array}{l} \!\!p({{{x}}_{k + 1}}\left| {{{{Z}}_{1:k}}} \right.)\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;i = 0 \\ \!\! p({{{x}}_{k + 1}}\left| {{{{Z}}_{1:k}},{{{z}}_{k + 1,i}}} \right.)\;\;\;\;\;i = 1,2 \cdots {n_{k + 1}} \\ \end{array} \right.$$ (34) $${u_i}({{{x}}_{k + 1}}) = \Pr ({{{z}}_{k + 1}} = {{{z}}_{k + 1,i}})$$ (35) 根据公式(32)-(35), 基于 PDA算法建模框架下的目标状态权重KL平均
${\bar p_{{\rm{PDA}}}}({{{x}}_{k + 1}})$ 为$$\begin{split} & {{\bar p}_{{\rm{PDA}}}}({{{x}}_{k + 1}}) = \arg \mathop {\min }\limits_{\bar p({{{x}}_{k + 1}})} \left\{ {\sum\limits_{i = 0}^{{n_{k + 1}}} {\Pr ({{{z}}_{k + 1}} = {{{z}}_{k + 1,i}})} } \right. \times \\ &\left. { {\rm{KL}}(\bar p({{{x}}_{k + 1}})\left| {p({{{x}}_{k + 1}}\left| {{{{Z}}_{1:k}},{{{z}}_{k + 1,i}}} \right.)} \right.)} \right\} \\[-15pt] \end{split} $$ (36) 对于VB-PDA算法而言, 其待估计参数集
${\Psi _{{\rm{VB}}}}$ 中的元素为${{{x}}_{k + 1}}$ 和$\left\{ {{\eta _{k + 1,0}},{\eta _{k + 1,1}}, \cdots {\eta _{k + 1,{n_{k + 1}}}}} \right\}$ . 从VB框架[18]来看, VB-PDA算法获取的参数集近似后验概率密度函数${q_{{\rm{VB}}}}({\Psi _{{\rm{VB}}}})$ 的等价形式为$$\begin{split} & {q_{{\rm{VB}}}}({\Psi _{{\rm{VB}}}}) = \arg \mathop {\max }\limits_{\bar p({\Psi _{{\rm{VB}}}})} {\rm{L}}(\bar p({\Psi _{{\rm{VB}}}})\left| {p({\Psi _{{\rm{VB}}}},{{{Z}}_{1:k + 1}})} \right.) = \\ & \arg \mathop {\min }\limits_{\bar p({\Psi _{{\rm{VB}}}})} {\rm{KL}}(\bar p({\Psi _{{\rm{VB}}}})\left| {p({\Psi _{{\rm{VB}}}},{{{Z}}_{1:k + 1}})} \right.) \\[-15pt] \end{split} $$ (37) $$\begin{split} &{\rm{L}}({q_{{\rm{VB}}}}({\Psi _{{\rm{VB}}}})\left| {p({\Psi _{{\rm{VB}}}},{{{Z}}_{1:k + 1}})} \right.) = \\ & \int {{q_{{\rm{VB}}}}({\Psi _{{\rm{VB}}}})\ln \frac{{p({\Psi _{{\rm{VB}}}},{{{Z}}_{1:k + 1}})}}{{{q_{{\rm{VB}}}}({\Psi _{{\rm{VB}}}})}}} {\rm{d}}{\Psi _{{\rm{VB}}}} \end{split} $$ (38) 其中
${\rm{L}}( \cdot )$ 表示下界函数, VB算法是一种基于平均场理论实现参数近似后验概率密度函数获取的方法,${q_{{\rm{VB}}}}({\Psi _{{\rm{VB}}}})$ 可被分解为$$ \begin{split} {q_{{\rm{VB}}}}({\Psi _{{\rm{VB}}}}) =& {q_{{\rm{VB}}}}({{{x}}_{k + 1}})\\ {q_{{\rm{VB}}}}({\eta _{k + 1,0}},{\eta _{k + 1,1}}, \cdots {\eta _{k + 1,{n_{k + 1}}}}) \end{split}$$ (39) 根据公式(25),(27)和(39),公式(37)的等价形式为
$$\begin{split} {q_{{\rm{VB}}}}({{{x}}_{k + 1}}) = &\arg \mathop {\min }\limits_{\bar p({{{x}}_{k + 1}})} {\rm{KL}}(\bar p({{{x}}_{k + 1}})\\ &\prod\limits_{i = 0}^{{n_{k + 1}}} {\Pr {{({{{z}}_{k + 1}} = {{{z}}_{k + 1,i}})}^{{\eta _{k + 1,i}}}}} \\ &{\rm{ }}\left| {p({\Psi _{{\rm{VB}}}}\left| {{{{Z}}_{1:k + 1}}} \right.)} \right.) \end{split} $$ (40) 根据公式(24)和(25),
$p({{{x}}_{k + 1}}\left| {\eta _{k + 1,0}},{\eta _{k + 1,1}}, \cdots {\eta _{k + 1,{n_{k + 1}}}},{{{Z}}_{1:k + 1}} \right.)$ 的表达式为$$\begin{split} & p({{{x}}_{k + 1}}\left| {{\eta _{k + 1,0}},{\eta _{k + 1,1}}, \cdots {\eta _{k + 1,{n_{k + 1}}}},{{{Z}}_{1:k + 1}}} \right.) = \\ & \dfrac{{p({{{x}}_{k + 1}},{\eta _{k + 1,0}},{\eta _{k + 1,1}}, \cdots {\eta _{k + 1,{n_{k + 1}}}},{{{Z}}_{1:k + 1}})}}{{p({\eta _{k + 1,0}},{\eta _{k + 1,1}}, \cdots {\eta _{k + 1,{n_{k + 1}}}}\left| {{{{Z}}_{1:k + 1}}} \right.)p({{{Z}}_{1:k + 1}})}} \propto \\ & {\cal{N}}{({{{x}}_{k + 1}};{{{x}}_{k + 1\left| k \right.}},{{{P}}_{k + 1\left| k \right.}})^{{\eta _{k + 1,0}}}}\\ & \prod\limits_{i = 1}^{{n_{k + 1}}} {{{({\cal{N}}({{{x}}_{k + 1}};{{{x}}_{k + 1\left| {k + 1,i} \right.}},{{{P}}_{k + 1\left| {k + 1,i} \right.}}))}^{{\eta _{k + 1,i}}}}} \Rightarrow\\ &p({{{x}}_{k + 1}}\left| {{\eta _{k + 1,0}},{\eta _{k + 1,1}}, \cdots {\eta _{k + 1,{n_{k + 1}}}},{{{Z}}_{1:k + 1}}} \right.)= \\ & {\cal{N}}{({{{x}}_{k + 1}};{{{x}}_{k + 1\left| k \right.}},{{{P}}_{k + 1\left| k \right.}})^{{\eta _{k + 1,0}}}}\\ &\prod\limits_{i = 1}^{{n_{k + 1}}} {{{({\cal{N}}({{{x}}_{k + 1}};{{{x}}_{k + 1\left| {k + 1,i} \right.}},{{{P}}_{k + 1\left| {k + 1,i} \right.}}))}^{{\eta _{k + 1,i}}}}} = \\ & p{({{{x}}_{k + 1}}\left| {{{{Z}}_{1:k}}} \right.)^{{\eta _{k + 1,0}}}}\\ &\prod\limits_{i = 1}^{{n_{k + 1}}} {{{(p({{{x}}_{k + 1}}\left| {{{{Z}}_{1:k}},{{{z}}_{k + 1,i}}} \right.))}^{{\eta _{k + 1,i}}}}} \\[-16pt] \end{split} $$ (41) 根据公式(39)和(41),
${\rm{KL}}(\bar p({{{x}}_{k + 1}})\prod\limits_{i = 0}^{{n_{k + 1}}} \Pr ({{{z}}_{k + 1}} = {{{z}}_{k + 1,i}})^{{\eta _{k + 1,i}}} \left| {p({\Psi _{{\rm{VB}}}}\left| {{{{Z}}_{1:k + 1}}} \right.)} \right.)$ 的表达式为$$\begin{split} & {\rm{KL}}\!(\bar p({{{x}}_{k + 1}})\!\!\!\prod\limits_{i = 0}^{{n_{k + 1}}} \!\!{\Pr \!{{({{{z}}_{k + 1}}\! \!=\!\! {{{z}}_{k + 1,i}}\!)}^{{\eta _{k + 1,i}}}}} \!\left| {p({\Psi _{{\rm{VB}}}}\left| {{{{Z}}_{1:k \!+\! 1}}}\! \right.)} \!\right.) \!\!= \\ & \int {\bar p({{{x}}_{k + 1}})(\prod\limits_{i = 1}^{{n_{k + 1}}} {\Pr {{({{{z}}_{k + 1}} = {{{z}}_{k + 1,i}})}^{{\eta _{k + 1,i}}}}} )}\times \\ & \ln \!\dfrac{{\bar p({{{x}}_{k + 1}})}}{{p({{{x}}_{k + 1}}\!\left| {{\eta _{k + 1,0}},{\eta _{k + 1,1}}, \cdots {\eta _{k + 1,{n_{k + 1}}}},{{{Z}}_{1:k \!+\! 1}}} \right.)}}{\rm{d}}{{{x}}_{k \!+\! 1}}\!\times \\ &{\rm{d}}\left\{ {{\eta _{k + 1,0}},{\eta _{k + 1,1}}, \cdots {\eta _{k + 1,{n_{k + 1}}}}} \right\}= \\ & \int \bar p({{{x}}_{k + 1}})\ln \bar p({{{x}}_{k + 1}}) - \bar p({{{x}}_{k + 1}})\\ &\sum\limits_{i = 1}^{{n_{k + 1}}} {\Pr ({{{z}}_{k + 1}} = {{{z}}_{k + 1,i}})} \times \\ &p({{{x}}_{k + 1}}\left| {{{{Z}}_{1:k}},{{{z}}_{k + 1,i}}} \right.){\rm{d}}{{{x}}_{k + 1}} = \\ & \sum\limits_{i = 0}^{{n_{k + 1}}} \!\!\Pr ({{{z}}_{k + 1}} \!\!=\!\! {{{z}}_{k + 1,i}}){\rm{KL}}(\bar p({{{x}}_{k + 1}})\!\left| {p({{{x}}_{k + 1}}\left| {{{{z}}_{1:k}},{{{z}}_{k \!+\! 1,i}}}\! \right.)} \!\right.) \end{split} $$ (42) 根据公式(36)和(42)可知, VB-PDA算法获取的参数集近似后验概率密度函数即为PDA算法建模框架下的状态代理概率密度函数集的权重KL平均, 因此该算法获取的状态近似后验概率密度函数比PDA算法获取的状态近似后验概率密度函数更能体现代理概率密度函数集统计特征.
-
算法的计算复杂度同样是影响算法跟踪性能的重要指标. 在这里我们将加减运算与乘除运算次数作为指标对四种算法的计算复杂度进行分析. 表1展示了PDA算法、距离加权概率数据关联算法[17](DW-PDA)以及VB-PDA算法在一步状态更新过程中所需的加减运算与乘除运算次数. 可以看出, VB-PDA算法比PDA算法以及DW-PDA算法拥有更低的算法复杂度. 当杂波数目较大多时, VB-PDA算法计算效率优势更加明显.
表 1 一步状态更新过程中所需的加减运算与乘除运算次数
Table 1. The number of addition and subtraction operations and multiplication and division operations required in the process of one-step state update
算法 加减法运算次数 乘除法运算次数 PDA $\begin{aligned} &{r^3}{n_{k + 1} } + {r^2}(m{n_{k + 1} } - {n_{k + 1} } + 2) + \\ & r({n_{k + 1} } + 1 + 2m) + 2{n_{k + 1} } - 1 +\\ & m{n_{k + 1} } - m\end{aligned}$ $\begin{aligned} & {r^3} + {r^2}(2{n_{k + 1} } + m + 3) + \\ & r(2m + 1) + m{n_{k + 1} } + 1 \end{aligned}$ DW-PDA $\begin{aligned} & {r^3}{n_{k + 1} } + {r^2}(m{n_{k + 1} } - {n_{k + 1} } + 2) + \\ & r({n_{k + 1} } + 1 + 2m) + 4{n_{k + 1} } - 3 + \\ & m{n_{k + 1} } - m\end{aligned}$ $\begin{aligned} &{r^3} + {r^2}(2{n_{k + 1} } + m + 3)+ \\ & r(2m + 1) + m{n_{k + 1} } + 4{n_{k + 1} } + 1\end{aligned}$ VB-PDA $\begin{aligned} & {r^3}{n_{k + 1} } + {r^2}(m{n_{k + 1} } - 2{n_{k + 1} } + 1) + \\ & 2mr + 2{n_{k + 1} } - 1 + m{n_{k + 1} } - m \end{aligned}$ $\begin{aligned} & {r^3} + {r^2}m + r(2m + 1)+ \\ & m{n_{k + 1} } + 1 + {m^2} \end{aligned}$ -
利用目标跟踪案例对算法的跟踪性能进行测试. 假设目标在二维平面进行匀速直线运动, 利用匀速模型(CV)对该目标进行跟踪并且量测设备可以获取目标的位置信息. 因此, 目标的状态信息
${{{x}}_k}$ , 状态转移矩阵${{{F}}_k}$ , 量测矩阵${{{H}}_{k + 1}}$ 分别为:$${{{x}}_k} = {\left[ {\begin{array}{*{20}{c}} {{a_k}}&{{b_k}}&{{{\dot a}_k}}&{{{\dot b}_k}} \end{array}} \right]^T}$$ (43) $${{{F}}_k} = \left[ {\begin{array}{*{20}{c}} 1&0&t&0 \\ 0&1&0&t \\ 0&0&1&0 \\ 0&0&0&1 \end{array}} \right]$$ (44) $${{{H}}_{k + 1}} = \left[ {\begin{array}{*{20}{l}} 1&0&0&0 \\ 0&1&0&0 \end{array}} \right]$$ (45) 其中,
${a_k}$ 和${b_k}$ 表示目标位置信息,${\dot a_k}$ 和${\dot b_k}$ 表示目标速度信息,$t = 1{\rm{s}}$ 表示采样时间. 系统过程噪声协方差和系统噪声量测协方差为:$${{{Q}}_{k + 1}} = 0.25\left[ {\begin{array}{*{20}{c}} {{1 / 3}}&{{1 / 2}}&0&0 \\ {{1 / 2}}&1&0&0 \\ 0&0&{{1 / 3}}&{{1 / 2}} \\ 0&0&{{1 / 2}}&1 \end{array}} \right]$$ (46) $${{{R}}_{k + 1}} = \left[ {\begin{array}{*{20}{c}} {1{\rm{00}}{{\rm{m}}^2}}&0 \\ 0&{1{\rm{00}}{{\rm{m}}^2}} \end{array}} \right]$$ (47) 目标真实初始状态信息
${{{x}}_{\rm{0}}}$ 为$[{\rm{100\;m}},{\rm{100\;m}},{\rm{1\;m/s}}, {\rm{1\;m/s}}]$ , 初始状态协方差${{{P}}_{{\rm{0}}\left| 0 \right.}}$ 为${\rm{diag(}}[{\rm{100}}\;{{\rm{m}}^2},{\rm{100}}\;{{\rm{m}}^2}, {\rm{1}}\;{{\rm{m}}^2}{\rm{/}}{{\rm{s}}^2},{\rm{1}}\;{{\rm{m}}^2}{\rm{/}}{{\rm{s}}^2}])$ , 目标初始状态信息${{{x}}_{0\left| {\rm{0}} \right.}}$ 为${x_{\rm{0}}}$ 与协方差为${{{P}}_{{\rm{0}}\left| 0 \right.}}$ 高斯白噪声的和[20], 该高斯白噪声由matlab程序随机产生. 因此在每一次蒙特卡洛仿真实验中, 目标初始状态信息是不同的.在仿真实验中, 我们将均方根误差(RSME)和平均时间均方根误差(TRSME)作为指标来比较各算法的估计性能, 目标位置和速度的RMSE和TRSME计算公式如下:
$$\begin{split} &RMS{E_{pos}} = \\ &\sqrt {\frac{1}{M}\sum\limits_{s = 1}^M {\left( {{{(\hat a_{k + 1}^s - a_{k + 1}^s)}^2} + {{(\hat b_{k + 1}^s - b_{k + 1}^s)}^2}} \right)} } \end{split}$$ (48) $$\begin{split} &RMS{E_{{\rm{vel}}}} =\\ &\sqrt {\frac{1}{M}\sum\limits_{s = 1}^M {\left( {{{(\hat {\dot{ a}}_{k + 1}^s - \dot a_{k + 1}^s)}^2} + {{{(\hat {\dot b}}_{k + 1}^s - \dot b_{k + 1}^s)}^2}} \right)} } \end{split}$$ (49) $$\begin{split} &TRMS{E_{pos}} = \frac{1}{L}\sum\limits_{k = 0}^{L - 1}\\ &{\sqrt {\frac{1}{M}\sum\limits_{s = 1}^M {\left( {{{(\hat a_{k + 1}^s - a_{k + 1}^s)}^2} + {{(\hat b_{k + 1}^s - b_{k + 1}^s)}^2}} \right)} } } \end{split}$$ (50) $$\begin{split} &TRMS{E_{{\rm{vel}}}} = \frac{1}{L}\sum\limits_{k = 0}^{L - 1} \\ &{\sqrt {\frac{1}{M}\sum\limits_{s = 1}^M {\left( {{{{(\hat {\dot a}}_{k + 1}^s - \dot a_{k + 1}^s)}^2} + {{{(\hat {\dot b}}_{k + 1}^s - \dot b_{k + 1}^s)}^2}} \right)} } } \end{split}$$ (51) 其中,
$M = 100$ 表示蒙特卡洛实验的次数,$L = 100{\rm{s}}$ 表示仿真实验总时长,$\hat a_{k + 1}^s$ 和$\hat b_{k + 1}^s$ 表示第$s$ 次蒙特卡洛实验目标位置估计信息,$a_{k + 1}^s$ 和$b_{k + 1}^s$ 表示第$s$ 次蒙特卡洛实验目标位置真实信息.${\hat {\dot a}}_{k + 1}^s$ 和${\hat {\dot b}}_{k + 1}^s$ 表示第$s$ 次蒙特卡洛实验目标速度估计信息.$\dot a_{k + 1}^s$ 和$\dot b_{k + 1}^s$ 表示第$s$ 次蒙特卡洛实验目标速度真实信息. 为了验证算法的有效性, 我们对PDA算法、DW-PDA算法和VB-PDA算法进行测试.假设杂波分布在
$[0 \sim 400\;{\rm{m}},{\rm{ }}0 \sim 400\;{\rm{m}}]$ 区域里服从均匀分布, 因此${V_{k + 1}} = 160000\;{{\rm{m}}^{\rm{2}}}$ . 杂波数目服从泊松分布. 检测概率为0.95, 杂波数目统计期望为20. 图1-2展出了三种算法的位置与速度RMSE, 表2展示了三种算法的位置与速度TRMSE. 可以发现,VB-PDA算法的位置估计精度在三种算法中最好, 而DW-PDA算法的速度估计精度在三种算法中最好. 从TRMSE数据可以发现, 相比于PDA算法, VB-PDA算法在位置和速度估计精度上分别提高了2.18%和0.11%. DW-PDA算法是在权重计算的过程中对PDA算法的改进, 而VB-PDA算法是从状态后验概率密度函数获取的过程中对PDA算法的改进. 由于两种算法的改进方式不同, 出现DW-PDA算法在速度估计精度优于VB-PDA算法这一现象属于正常情形. 表3展示了四种算法一次蒙特卡洛实验所需的计算时间, 可以看出VB-PDA算法的所需的计算时间最少, 这也验证了第五部分分析的结果. 从数据中我们可以发现, 相比于PDA算法, VB-PDA算法在计算效率上提高了17.37%.
图 1 场景1下三种算法的位置RSME
Figure 1. The RMSE of position from three algorithms in scenario 1
图 2 仿真场景1下三种算法的速度RSME
Figure 2. The RMSE of velocity from three algorithms in scenario 1
表 2 场景1下三种算法的TRSME
Table 2. The TRSME of three algorithms in scenario 1
算法 位置TRMSE/m 速度TRMSE/(m/s) PDA 6.892 0.873 DW-PDA 6.792 0.839 VB-PDA 6.742 0.872 表 3 场景1下一次蒙特卡洛仿真实验所需的计算时间
Table 3. The computational time of three algorithms at one Monte Carlo simulation experiment in scenario 1
算法 计算时间/ms PDA 56.52 DW-PDA 58.72 VB-PDA 46.70 -
考虑基于机动目标跟踪案例对提出算法的跟踪性能进行测试, 利用IMM算法[22]与上述三种算法结合对目标进行跟踪. IMM算法中采用两个带不同转弯速率匀速转弯模型(CT)对该目标进行跟踪并且量测设备可以获取目标的位置信息. 因此, 目标的状态信息
${x_k}$ , 子模型状态转移矩阵${{{F}}_{k,j}}{\rm{ }}(j = 1,2)$ , 量测矩阵${{{H}}_{k + 1,j}}{\rm{ }}(j = 1,2)$ 分别为:$${{{x}}_k} = {\left[ {\begin{array}{*{20}{c}} {{a_k}}&{{b_k}}&{{{\dot a}_k}}&{{{\dot b}_k}} \end{array}} \right]^T}$$ (52) $${{{F}}_{k,j}} = \left[ {\begin{array}{*{20}{c}} 1&0&{\dfrac{{\sin ({w_j}t)}}{w}}&{\dfrac{{ - (1 - \cos ({w_j}t))}}{{{w_j}}}} \\ 0&1&{\dfrac{{(1 - \cos ({w_j}t))}}{{{w_j}}}}&{\dfrac{{\sin ({w_j}t)}}{w}} \\ 0&0&{\cos ({w_j}t)}&{ - \sin ({w_j}t)} \\ 0&0&{\sin ({w_j}t)}&{\cos ({w_j}t)} \end{array}} \right]$$ (53) $${{{H}}_{k + 1,j}} = \left[ {\begin{array}{*{20}{l}} 1&0&0&0 \\ 0&1&0&0 \end{array}} \right]$$ (54) 其中,
${a_k}$ 和${b_k}$ 表示目标位置信息,${\dot a_k}$ 和${\dot b_k}$ 表示目标速度信息,$t = 1{\rm{s}}$ 表示采样时间. 子模型的转弯速率分别为${w_1} = - \dfrac{\pi }{{20}}{\rm{rad/s}}$ ,${w_2} = \dfrac{\pi }{{20}}{\rm{rad/s}}$ . 系统过程噪声协方差${{{Q}}_{k,j}}$ 和量测噪声协方差${{{R}}_{k + 1,j}}$ 设置为:$${{{Q}}_{k,j}} = 0.25\left[ {\begin{array}{*{20}{c}} {{1 / 3}}&{{1 / 2}}&0&0 \\ {{1 / 2}}&1&0&0 \\ 0&0&{{1 / 3}}&{{1 / 2}} \\ 0&0&{{1 / 2}}&1 \end{array}} \right]$$ (55) $$.{{{R}}_{k + 1,j}} = \left[ {\begin{array}{*{20}{c}} {1{\rm{00}}{{\rm{m}}^2}}&0 \\ 0&{1{\rm{00}}{{\rm{m}}^2}} \end{array}} \right]$$ (56) 目标真实初始状态信息
${{{x}}_{\rm{0}}}$ 为$[{\rm{100\;m}},{\rm{100\;m}},{\rm{3\;m/s}}, {\rm{3\;m/s}}]$ , 目标在1-50 s做转弯速率为$\dfrac{\pi }{{100}}{\rm{rad/s}}$ 的匀速转弯运动, 在51-100 s做转弯速率为$- \dfrac{\pi }{{100}}{\rm{rad/s}}$ 的匀速转弯运动. 假设两个子模型初始权重均为0.5, 马尔可夫转移矩阵为$\left[ {\begin{aligned} {0.9}&{0.1} \\ {0.1}&{0.9} \end{aligned}} \right]$ . 子模型初始状态协方差为${\rm{diag(}}[{\rm{100}}\;{{\rm{m}}^2},{\rm{100}}\;{{\rm{m}}^2},{\rm{1}}\;{{\rm{m}}^2}{\rm{/}}{{\rm{s}}^2},{\rm{1}}\;{{\rm{m}}^2}{\rm{/}}{{\rm{s}}^2}])$ . 为了验证算法的有效性, 我们对IMM-PDA算法、IMM-DW-PDA算法和IMM-VB-PDA算法进行测试.假设杂波分布在
$[0 \sim 400\;{\rm{m}},{\rm{ }}0 \sim 400\;{\rm{m}}]$ 区域里服从均匀分布, 因此${V_{k + 1}} = 160000\;{{\rm{m}}^{\rm{2}}}$ . 杂波数目服从泊松分布. 检测概率为0.9, 杂波数目统计期望为10. 图3-4给出了三种算法获取的位置与速度RMSE. 表4给出了三种算法获取的位置与速度TRSME. 可以看到在此类场景下, IMM-VB-PDA算法的估计精度最佳. 从TRMSE数据可以发现, 相比于IMM-PDA算法, IMM-VB-PDA算法在位置和速度估计精度上分别提高了5.26%和4.8%.
图 3 仿真场景2下三种算法的位置RSME
Figure 3. The RMSE of position from three algorithms in scenario 2
图 4 仿真场景2下三种算法的速度RSME
Figure 4. The RMSE of velocity from three algorithms in scenario 2
表 4 仿真场景2下三种算法的TRSME
Table 4. The TRSME of three algorithms in scenario 2
算法 位置TRMSE/m 速度TRMSE/(m/s) IMM-PDA 8.409 1.874 IMM-DW-PDA 8.101 1.818 IMM-VB-PDA 7.967 1.784 -
PDA算法采用从状态代理概率密度函数集里提取均值信息与协方差信息的方式获取状态近似后验概率密度函数. 该方法虽然便捷, 但在数据数目较多的情形下存在计算效率低的缺陷. 针对这一问题, 本文提出了一种VB-PDA算法, 该算法本质上是在PDA算法建模框架下提取相应的权重KL平均以获取状态近似后验概率密度函数. 相比PDA算法, VB-PDA算法在提高计算效率的同时, 获取了近似程度更高的状态后验概率密度函数. 相关仿真实验结果验证了VB-PDA算法的有效性. 由于本文提出的算法是针对单目标跟踪问题设计的, 如何将类似的思想拓展到多目标跟踪领域是未来的研究目标之一.
Variational Bayesian Probabilistic Data Association Algorithm
More Information-
摘要: 针对杂波环境下的目标跟踪问题, 提出了一种基于变分贝叶斯的概率数据关联算法(VB-PDA). 该算法首先将关联事件视为一个随机变量并利用多项分布对其进行建模, 随后基于数据集、目标状态、关联事件的联合概率密度函数求取关联事件的后验概率密度函数, 最后将关联事件的后验概率密度函数引入变分贝叶斯框架中以获取状态近似后验概率密度函数. 相比于概率数据关联算法, VB-PDA算法在提高算法实时性的同时在权重KL平均准则下获取了近似程度更高的状态后验概率密度函数. 相关仿真实验对提出算法的有效性进行了验证.Abstract: Aiming at the problem of target tracking in clutter, this paper proposes a variational Bayesian based probabilistic data association algorithm (VB-PDA). Firstly, associated events are regarded as a random variable and modelled by the multinomial distribution. Then, the joint probability density function of data set, target state and associated events is constructed and the posterior probability density function of associated events is obtained by using this joint probability density function. Finally, the posterior probability density function of associated events is introduced into the framework of variational Bayesian to obtain the approximate posterior probability density function of state. Compared with the probabilistic data association algorithm, the VB-PDA algorithm obtains a state posterior probability density function with higher approximation degree based on the weight KL average criterion while improving real-time performance. The simulation experiments verify the effectiveness of proposed algorithm.
-
表 1 一步状态更新过程中所需的加减运算与乘除运算次数
Table 1 The number of addition and subtraction operations and multiplication and division operations required in the process of one-step state update
算法 加减法运算次数 乘除法运算次数 PDA $\begin{aligned} &{r^3}{n_{k + 1} } + {r^2}(m{n_{k + 1} } - {n_{k + 1} } + 2) + \\ & r({n_{k + 1} } + 1 + 2m) + 2{n_{k + 1} } - 1 +\\ & m{n_{k + 1} } - m\end{aligned}$ $\begin{aligned} & {r^3} + {r^2}(2{n_{k + 1} } + m + 3) + \\ & r(2m + 1) + m{n_{k + 1} } + 1 \end{aligned}$ DW-PDA $\begin{aligned} & {r^3}{n_{k + 1} } + {r^2}(m{n_{k + 1} } - {n_{k + 1} } + 2) + \\ & r({n_{k + 1} } + 1 + 2m) + 4{n_{k + 1} } - 3 + \\ & m{n_{k + 1} } - m\end{aligned}$ $\begin{aligned} &{r^3} + {r^2}(2{n_{k + 1} } + m + 3)+ \\ & r(2m + 1) + m{n_{k + 1} } + 4{n_{k + 1} } + 1\end{aligned}$ VB-PDA $\begin{aligned} & {r^3}{n_{k + 1} } + {r^2}(m{n_{k + 1} } - 2{n_{k + 1} } + 1) + \\ & 2mr + 2{n_{k + 1} } - 1 + m{n_{k + 1} } - m \end{aligned}$ $\begin{aligned} & {r^3} + {r^2}m + r(2m + 1)+ \\ & m{n_{k + 1} } + 1 + {m^2} \end{aligned}$ 
下载: 导出CSV
表 2 场景1下三种算法的TRSME
Table 2 The TRSME of three algorithms in scenario 1
算法 位置TRMSE/m 速度TRMSE/(m/s) PDA 6.892 0.873 DW-PDA 6.792 0.839 VB-PDA 6.742 0.872
下载: 导出CSV
表 3 场景1下一次蒙特卡洛仿真实验所需的计算时间
Table 3 The computational time of three algorithms at one Monte Carlo simulation experiment in scenario 1
算法 计算时间/ms PDA 56.52 DW-PDA 58.72 VB-PDA 46.70
下载: 导出CSV
七星彩规则表 4 仿真场景2下三种算法的TRSME
Table 4 The TRSME of three algorithms in scenario 2
算法 位置TRMSE/m 速度TRMSE/(m/s) IMM-PDA 8.409 1.874 IMM-DW-PDA 8.101 1.818 IMM-VB-PDA 7.967 1.784
下载: 导出CSV
-
[1] 孟琭, 杨旭. 目标跟踪算法综述. 自动化学报, 2019, 45(7): 1244?1260 MENG Lu, YANG Xu. A Survey of Object Tracking Algorithms. Acta Automatica Sinica, 2019, 45(7): 1244?1260 [2] Sobhani B, Paolini E, Giorgetti A. Target Tracking for UWB Multistatic Radar Sensor Networks. IEEE Journal of Selected Topics in Signal Processing, 2017, 8(1): 125?136 [3] 甘林海, 王刚, 刘进忙, 李松. 群目标跟踪技术综述. 自动化学报, 2020, 46(3): 411?426 GAN Lin-Hai, WANG Gang, LIU Jin-Mang, LI Song. An Overview of Group Target Tracking. Acta Automatica Sinica, 2020, 46(3): 411?426 [4] 何友, 黄勇, 关键, 陈小龙. 海杂波中的雷达目标检测技术综述. 现代雷达, 2014, 36(12): 1?9 HE You, HUANG Yong, GUAN Jian, CHEN Xiao -Long. An overview on radar target detection in sea clutter. Modern radar, 2014, 36(12): 1?9 [5] 陈一梅, 刘伟峰, 孔明鑫, 张桂林. 基于GLMB滤波和Gibbs采样的多扩展目标有限混合建模与跟踪算法. 自动化学报, 2020, 46(7): 1445?1456 CHEN Yi-Mei, LIU Wei-Feng, KONG Ming-Xin, ZHANG Gui-Lin. A Modeling and Tracking Algorithm of Finite Mixture Models for Multiple Extended Target Based on the GLMB Filter and Gibbs Sampler. Acta Automatica Sinica, 2020, 46(7): 1445?1456 [6] Kirubarajan T, Bar-Shalom Y. Probabilistic data association techniques for target tracking in clutter. Proceedings of the IEEE, 2004, 92(3): 536?557 doi: 10.1109/JPROC.2003.823149 [7] Wu P L, Zhou Y, Li X X. Space-based passive tracking of non-cooperative space target using robust filtering algorithm. Proceedings of the Institution of Mechanical Engineers, Part I: Journal of Systems and Control Engineering, 2016, 230(6): 551?561 doi: 10.1177/0959651816637770 [8] Song T L, L D G, R J. A probabilistic nearest neighbor filter algorithm for tracking in a clutter environment. Signal Processing, 2005, 85(10): 2044?2053 doi: 10.1016/j.sigpro.2005.01.016 [9] Sinha A, Ding Z, Kirubarajan T, Farooq, M. Track Quality Based Multitarget Tracking Approach for Global Nearest-Neighbor Association. IEEE Transactions on Aerospace & Electronic Systems, 2012, 48(2): 1179?1191 [10] Az iz, Ashraf M. A new nearest-neighbor association approach based on fuzzy clustering. Aerospace Science & Technology, 2013, 26(1): 87?97 [11] Bar-Shalom Y, Tse E. Tracking in a cluttered environment with probabilistic data association. Automatica, 1975, 11(5): 451?460 doi: 10.1016/0005-1098(75)90021-7 [12] Wu P L, Li X X, Kong J S, Liu J L. Heterogeneous multiple sensors joint tracking of maneuvering target in clutter. Sensors, 2015, 15(7): 17350?17365 doi: 10.3390/s150717350 [13] 程婷, 何子述, 李亚星. 一种具有自适应关联门的杂波中机动目标跟踪算法. 电子与信息学报, 2012, 34(4): 865?870 CHENG Ting, HE Zi-Shu, LI Ya-Xing. A maneuvering target tracking algorithm in clutter with adaptive correlation gate. Journal of Electronics and Information, 2012, 34(4): 865?870 [14] Barshalom Y, Daum F, Huang J. The Probabilistic Data Association Filter. IEEE Control Systems Magazine, 2012, 29(6): 82?100 [15] Kirubarajan T, Bar-Shalom Y, Blair W D, Watson G A. IMMPDAF for radar management and tracking benchmark with ECM. IEEE Transactions on Aerospace & Electronic Systems, 1998, 34(4): 1115?1134 [16] 潘泉, 刘刚, 戴冠中, 张洪才. 联合交互式多模型概率数据关联算法. 航空学报, 1999, 20(3): 234?238 PAN Quan, LIU Gang, DAI Guan-Zhong, ZHANG Hong-Cai. Combined interacting multiple models probabilistic data association algorithm. Acta Aeronautica ET Astronautica Sinica, 1999, 20(3): 234?238 [17] 陈晓, 李亚安, 李余兴, 蔚婧. 基于距离加权的概率数据关联机动目标跟踪算法. 上海交通大学学报, 2018, 52(4): 100?105 CHEN Xiao, LI Ya-An, Li Yu-Xing, WEI Jing. Maneuvering Target Tracking Algorithm Based on Weighted Distance of Probability Data Association. Journal of Shanghai Jiao Tong University, 2018, 52(4): 100?105 [18] Yun P, Wu P L, He S. Pearson Type VⅡ Distribution-Based Robust Kalman Filter under Outliers Interference. IET Radar, Sonar & Navigation, 2019, 13(8): 1389?1399 [19] Battistelli G, Chisci L. Kullback–Leibler average, consensus on probability densities, and distributed state estimation with guaranteed stability. Automatica, 2014, 50(3): 707?718 doi: 10.1016/j.automatica.2013.11.042 [20] Huang Y L, Zhang Y G, Li N, Wu Z M, Chambers J. A novel robust student's t-based Kalman filter. IEEE Transactions on Aerospace and Electronic Systems, 2017, 53(3): 1545?1554 doi: 10.1109/TAES.2017.2651684 [21] Panta K, Ba-Ngu V, Singh S. Novel data association schemes for the probability hypothesis density filter. IEEE Transactions on Aerospace and Electronic Systems, 2007, 43(2): 556?570 doi: 10.1109/TAES.2007.4285353 [22] Blom H A P, Bar-Shalom Y. The interacting multiple model algorithm for systems with Markovian switching coefficients. IEEE Transactions on Automatic Control, 1988, 33(8): 780?783 doi: 10.1109/9.1299 -
点击查看大图
计量
- 文章访问数: 14
- HTML全文浏览量: 7
- 被引次数: 0
