自动化学报  2017, Vol. 43 Issue (2): 238-247   PDF    
均方根嵌入式容积粒子PHD多目标跟踪方法
熊志刚, 黄树彩, 赵炜, 苑智玮, 徐晨洋     
空军工程大学防空反导学院 西安 710051
摘要: 针对基于概率假设密度算法(Probability hypothesis density,PHD)的非线性多目标跟踪精度低、滤波发散等问题,提出了一种新的PHD算法——改进的均方根嵌入式容积粒子PHD算法(Advanced square-root imbedded cubature particle PHD,ASRICP-PHD).新的算法在初始化采样时将整个采样区域等概率划分为若干个区域,然后利用既定的准则从每个区域抽取粒子,并利用均方根嵌入式容积滤波方法对每个粒子进行滤波,来拟合重要密度函数,预测和更新多目标状态的PHD.仿真结果表明该算法能对多目标进行有效跟踪,相比拟随机采样法和伪随机采样,等概率采样的方法在多目标位置估计和数目估计上有更高的精度.
关键词: 多目标跟踪     概率假设密度     均方根嵌入式容积滤波     等概率采样    
Square-root Imbedded Cubature Particle PHD Multi-target Tracking Algorithm
XIONG Zhi-Gang, HUANG Shu-Cai, ZHAO Wei, YUAN Zhi-Wei, XU Chen-Yang     
Air and Missile Defense College, Air Force Engineering University, Xi'an 710051
Received: 2015-12-30, Accepted: 2016-06-12.
Foundation Item: Supported by National Natural Science Foundation of China (61503408, 61573374), Natural Science Basic Research Plan in Shaanxi Province of China (2012JM8020), and Aeronautical Science Foundation of China (20130196004)
Author brief: HUANG Shu-Cai Professor at the Air and Missile Defense College, Air Force Engineering University. He received his Ph. D. degree from Air Force Engineering University in 2005. His research interest covers aerospace cooperative target tracking and interception cueing;
ZHAO Wei Ph. D. candidate at the Air and Missile Defense College, Air Force Engineering University. He received his master degree from Air Force Engineering University in 2014. His research interest covers system modeling and simulation, target tracking and detection;
YUAN Zhi-Wei Master student at the Air and Missile Defense College, Air Force Engineering University. He received his bachelor degree from Changchun University of Science and Technology in 2014. His research interest covers infrared target tracking and detection, and system identification;
XU Chen-Yang Ph. D. candidate at the Air and Missile Defense College, Air Force Engineering University. He received his master degree from Air Force Engineering University in 2016. His research interest covers interception system modeling and simulation
Corresponding author. XIONG Zhi-Gang Ph. D. candidate at the Air and Missile Defense College, Air Force Engineering University. He received his master degree from Air Force Engineering University in 2016. His main research interest is aerospace cooperative target tracking. Corresponding author of this paper
Recommended by Associate Editor LAI Jian-Huang
Abstract: Considering the low accuracy, filter divergence and other problems of nonlinear multi-target tracking based on probability hypothesis density (PHD), a new filter named advanced square-root imbedded cubature particle PHD (ASRICP-PHD) is proposed. ASRICP-PHD divides the whole particle sampling area into several parts of equal probability, then uses a special rule to obtain particles from each part, and matches the important density function with square-root imbedded cubature particle filter, and therefore predicts and updates PHD. Simulation shows that ASRICP-PHD is able to track multiple targets effectively. Moreover, compared with quasi random sampling, the method of particle sampling based on probability has higher accuracy in terms of multi-target positions and number's estimations.
Key words: Multi-target tracking     probability hypothesis density (PHD)     square-root imbedded cubature filter     sampling with equal probability    

概率假设密度算法(Probability hypothesis density, PHD) 是公认的多目标跟踪的有效手段, 与其他的多目标跟踪算法相比, PHD在目标数目未知情况下性能更加优越, 有很好的发展应用前景, 因此是现在研究的热点[1].然而现有的PHD算法在进行非线性多目标跟踪时仍存在精度低、实时性差的问题.基于蒙特卡罗的粒子PHD算法[2] (Monte Carlo particle PHD, MCP-PHD) 和基于拟蒙特卡罗的粒子PHD算法[3-4] (Quasi-Monte Carlo particle PHD, QMCP-PHD) 需要大量的粒子, 在提高估计精度的同时也增加了时间成本.高斯厄米特粒子PHD算法[5] (Gauss Hermite particle PHD, GHP-PHD) 和容积粒子PHD算法[6] (Cubature particle PHD, CP-PHD) 在PHD预测阶段减少了粒子的开支, 但是却引入了各自不同的问题: GHP-PHD高斯厄米特积分点数目随目标状态的维度增加呈指数增长, 导致“维数灾难”; CP-PHD有效解决了“维数灾难”的问题, 但是容积积分点偏差随着目标状态维度增加而增大, 导致高维下粒子溢出有效积分区域, 即“粒子溢出”[7].基于此, 本文引入嵌入式容积准则求解积分点, 得到均方根嵌入式容积粒子PHD算法(Square-root imbedded cubature particle PHD, SRICP-PHD), 嵌入式容积准则可通过自由变量控制积分点的偏差, 能有效防止“粒子溢出”, 并且, 在某些条件下, 嵌入式容积积分近似精度高于容积准则[8].

另外, 引入高斯混合[9]的手段实现非高斯噪声估计, 改善滤波效果.具体做法是初始化采样得到粒子集, 并将粒子视为子目标, 估计其噪声特性, 继而混合得到目标非高斯噪声.该过程中, 采样方法至关重要, 它将决定粒子使用效率, 影响算法精度和稳定性.传统的伪随机采样(Pseudo random sampling, PRS) 会出现“粒子贫化”、“滤波发散”的问题, 可增加粒子数目, 但同时会增加时间成本, 而且增加粒子数并不能解决本质问题, 粒子数不可能无休止地增大.另一途径是寻找新的采样手段, 如拟随机采样(Quasi random sampling, QRS), 然而拟随机采样效果有限[10].基于此, 本文提出了等概率采样(Sampling with equal probability, EPS) 准则, 将整体采样区域等概率划分, 并按既定的规则从各个子区域选取粒子.该采样手段能保证区域之间无交集, 区域内部粒子不会聚集在一起, 因此能有效避免粒子“聚集”的现象.由于等概率区域以多维球面为边界, 因此可以通过容积准则中的球规则获取对称分布的积分点, 如此便得到改进的均方根嵌入式容积粒子PHD算法(Advanced square-root imbedded cubature particle PHD, ASRICP-PHD), 解决了粒子“聚集”和积分点获取的问题, 改善了算法的精度.

1 均方根嵌入式容积粒子PHD

PHD是多目标状态随机集的一阶矩[11].与单目标状态估计一样, 多目标跟踪旨在估计多目标状态的后验概率密度:

$ {f_{k\left| {k - 1} \right.}}\;({{\pmb{X}}_k}\left| {{{\pmb{Z}}_{1:k - 1}}} \right.) = \int {{f_{k\left| {k - 1} \right.}}({{\pmb{X}}_k}\left| {{{\pmb{X}}_{k - 1}}} \right.)} \times \nonumber \\ \qquad \quad {f_{k - 1\left| {k - 1} \right.}}({{\pmb{X}}_{k - 1}}\left| {{{\pmb{Z}}_{1:k - 1}}} \right.){\rm d}{{\pmb{X}}_{k - 1}} $ (1)
$ {f_{k\left| k \right.}}({{\pmb{X}}_k}\left| {{{\pmb{Z}}_{1:k}}} \right.) = \nonumber \\ \quad \qquad \frac{{{f_k}({{\pmb{Z}}_k}\left| {{{\pmb{X}}_k}} \right.){f_{k\left| {k - 1} \right.}}({{\pmb{X}}_k}\left| {{{\pmb{Z}}_{1:k - 1}}} \right.)}}{{\int {{f_k}({{\pmb{Z}}_k}\left| {{{\pmb{X}}_k}} \right.) {f_{k\left| {k - 1} \right.}}({{\pmb{X}}_k}\left| {{{\pmb{Z}}_{1:k - 1}}} \right.){\rm d}{{\pmb{X}}_k}} }} $ (2)

其中, ${f_{k\left| {k - 1} \right.}}({{\pmb{X}}_k}\left| {{{\pmb{Z}}_{1:k - 1}}} \right.)$为多目标状态集的预测概率密度; ${{f_{k\left| {k - 1} \right.}}({{\pmb{X}}_k}\left| {{{\pmb{X}}_{k - 1}}} \right.)}$, ${{f_k}({{\pmb{Z}}_k}\left| {{{\pmb{X}}_k}} \right.)}$分别为状态转移函数和量测似然函数; ${f_{k\left| k \right.}}({{\pmb{X}}_k}\left| {{{\pmb{Z}}_{1:k}}} \right.)$为多目标状态集的后验概率密度.

难点在于如何表示数目变化的多目标状态以非线性积分计算. Mahler在总结前人研究成果的基础上, 利用随机有限集(Random finite sets, RFS) 来表示多目标状态.假设k时刻目标数和量测量数目分别为$N_k$$M_k$, 则多目标状态RFS $\pmb{X}_k$和量测RFS $\pmb{Z}_k$可表示为

$ {{\pmb{X}}_k} = \left( {{{\pmb{x}}_{k,1}}, \cdots, {{\pmb{x}}_{k,{N_k}}}} \right) $ (3)
$ {{\pmb{Z}}_k} = \left( {{{\pmb{z}}_{k,1}}, \cdots, {{\pmb{z}}_{k,{M_k}}}} \right) $ (4)

由于影响目标数目变化的因素有“目标的衍生”及“新生目标”出现, 因此式(3) 和(4) 可表示为

$ {{\pmb{X}}_k} =\left[ {\mathop \cup \limits_{{\pmb{\varsigma }} \in {{\pmb{X}}_{k - 1}}} {{\pmb{S}}_{k\left| {k - 1} \right.}}({\pmb{\varsigma }})} \right] \cup\nonumber\\ \qquad \left[ {\mathop \cup \limits_{{\pmb{\varsigma }} \in {{\pmb{X}}_{k - 1}}} {{\pmb{B}}_{k\left| {k - 1} \right.}}({\pmb{\varsigma }})} \right] \cup {{\pmb{\Gamma }}_k} $ (5)
$ {{\pmb{Z}}_k} = {{\pmb{K}}_k} \cup \left [ {\mathop \cup \limits_{{\pmb{x}} \in {{\pmb{X}}_{k - 1}}} {{\pmb{\Theta }}_k}({\pmb{x}})} \right] $ (6)

其中, ${{\pmb{S}}_{k\left| {k - 1} \right.}}( \cdot )$表示存活目标状态转移; ${{\pmb{B}}_{k\left| {k - 1} \right.}}( \cdot )$表示目标的衍生; ${{\pmb{\Gamma }}_k}$, ${{\pmb{\Theta }}_k}( \cdot )$为目标对应的量测; ${{\pmb{K}}_k}$为虚警杂波.

针对式(1) 和(2) 中的非线性积分, Mahler提出利用多目标状态一阶矩近似求解.假设k时刻先验PHD和后验PHD分别为${\upsilon _{k\left| {k - 1} \right.}}({\pmb{x}})$${\upsilon _{k\left| k \right.}}({\pmb{x}})$, 则:

$ {\upsilon _{k\left| {k - 1} \right.}}({\pmb{x}}) =\nonumber\\ \quad \int {{P_{S,k\left| {k - 1} \right.}}({\pmb{x}}'){f_{k\left| {k - 1} \right.}}({\pmb{x}}\left| {{\pmb{x}}'} \right.)} {\upsilon _{k - 1}}({\pmb{x}}'){\rm d}{\pmb{x}}' +\nonumber\\ \quad \int {{\beta _{k\left| {k - 1} \right.}}({\pmb{x}}\left| {{\pmb{x}}'} \right.)} {\upsilon _{k - 1}}({\pmb{x}'}){\rm d}{\pmb{x}}' + {\gamma _k}({\pmb{x}}) $ (7)
$ {\upsilon _{k\left| k \right.}}({\pmb{x}}) = (1 - {P_{D,k}}({\pmb{x}})){\upsilon _{k\left| {k - 1} \right.}}({\pmb{x}}) +\nonumber\\ \quad \sum\limits_{{\pmb{z}} \in {{\pmb{Z}}_k}} {\frac{{{P_{D,k}}({\pmb{x}}){f_{k\left| k \right.}} ({\pmb{z}}\left| {\pmb{x}} \right.){\upsilon _{k\left| {k - 1} \right.}}({\pmb{x}})}}{{{\kappa _k}({\pmb{z}}) + \int {{P_{D,k}}({\pmb{x}}){f_{k\left| k \right.}}({\pmb{z}}\left| {\pmb{x}} \right.){\upsilon _{k\left| {k - 1} \right.}}({\pmb{x}})} {\rm d}{\pmb{x}}}}} $ (8)

其中, ${\gamma _k}( \cdot )$为新生目标状态集PHD; ${\beta _{k\left| {k - 1} \right.}}( \cdot )$为衍生目标状态集PHD; ${P_{S, k\left| {k - 1} \right.}}( \cdot )$为单目标的生存概率; ${P_{D, k}}( \cdot )$为当前时刻对状态为$\pmb{x}$的目标的探测概率; ${{\kappa _k}({\pmb{z}})}={{\lambda _k}{c_k}({\pmb{z}})}$为杂波集的PHD, 杂波服从均值为${{\lambda _k}}$的泊松分布, ${{c_k}({\pmb{z}})}$为杂波概率密度.

1.1 高斯混合PHD

式(7) 和(8) 中的积分无法得到解析解, 可通过高斯混合的方式来解决[12], 先验的PHD为

$ {\upsilon _{k - 1}} = \sum\limits_{i = 1}^{{J_{k - 1}}} {w_{k - 1}^{(i)}} {\rm N}({\pmb{x}};{\pmb{m}}_{k - 1}^{(i)},{\pmb{P}}_{k - 1}^{(i)}) $ (9)

式(9) 中, ${{J_{k - 1}}}$为先验的高斯分量个数; ${\pmb{m}}_{k - 1}^{(i)}$${\pmb{P}}_{k - 1}^{(i)}$分别为前一时刻估计得到的目标状态值和状态协方差值; ${w_{k - 1}^{(i)}}$为高斯分量对应的权重, 其和为前一时刻估计的目标数目.估计得到$k-1$时刻的PHD, 然后进行多目标贝叶斯估计的递推:

$ {\upsilon _{\left. k \right|k - 1}} ({\pmb{x}}) =\nonumber\\ \quad \sum\limits_{i = 1}^{{J_{\left. k \right|k - 1}}} {w_{k\left| {k - 1} \right.}^{(i)}} {\rm N}\left( {{\pmb{x}}; {\pmb{m}}_{k\left| {k - 1} \right.}^{(i)},{\pmb{P}}_{k\left| {k - 1} \right.}^{(i)}} \right) $ (10)
$ {\upsilon _k}({\pmb{x}}) = (1 - {P_{D,k}}){\upsilon _{k\left| {k - 1} \right.}}({\pmb{x}}) +\nonumber\\ \qquad \sum\limits_{{\pmb{z}} \in {{\pmb{Z}}_k}} {\sum\limits_{i = 1}^{{J_{\left. k \right|k - 1}}} {w_k^{(i)}({\pmb{z}})} {\rm N}\left( {{\pmb{x}};{\pmb{m}}_k^{(i)}({\pmb{z}}),{\pmb{P}}_k^{(i)}({\pmb{z}})} \right)} $ (11)

式(10) 中, ${{J_{\left. k \right|k - 1}}}$预测得到的高斯分量个数; ${w_{k\left| {k - 1} \right.}^{(i)}}$为预测高斯分量的权值, 等于${w_{k - 1}^{(i)}}$; ${{\pmb{m}}_{k\left| {k - 1} \right.}^{(i)}}$, ${{\pmb{P}}_{k\left| {k - 1} \right.}^{(i)}}$分别为预测得到的高斯分量均值和协方差.式(11) 中, ${w_k^{(i)}({\pmb{z}})}$为更新的权重; ${{\pmb{m}}_k^{(i)}({\pmb{z}})}$, ${{\pmb{P}}_k^{(i)}({\pmb{z}})}$为状态和方差的更新值.

1.2 SRICP-PHD及其主要步骤

三阶嵌入式容积准则通过自由变量$\delta$来控制积分点偏差, 积分点数目取决于目标状态维数[13], 优点是高维下积分点不会溢出“有效的积分区域”.利用嵌入式容积准则产生积分点, 估计多目标状态及误差平方根, 从而得到一种新的PHD滤波算法SRICP-PHD, 可改善PHD算法性能, 并解决因舍入造成的误差矩阵奇异的问题.简便起见, 对权值为$w_j$, 均值为$\pmb{m}$的变量$\pmb{x}_j$, 定义$\Theta$运算:

$ \Theta \left( {\left\{ {{w_j},{{\pmb{x}}_j}} \right\},{\pmb{m}}} \right) = {w_j}({{\pmb{x}}_j} - {\pmb{m}}){({{\pmb{x}}_j} - {\pmb{m}})^{\rm T}} $ (12)

假设得到式(9) 所示的先验PHD, 采样得到粒子集$\left\{ {{\pmb{x}}_{k - 1}^{(i)(j)}} \right\}_{j = 1}^N \sim {\rm N}\left( {{\pmb{x}};{\pmb{m}}_{k - 1}^{(i)}, {\pmb{P}}_{k - 1}^{(i)}} \right)$, 同时根据Tria变换求解对应状态协方差的平方根${\pmb{S}}_{k - 1}^{(i)}$, Tria定义见文献[7], 则PHD更新计算如下:

1) 积分点采样

$ {\pmb{x}}_{k - 1}^{(i)(j)(l)} = {\pmb{S}}_{k - 1}^{(i)}{{\pmb{\xi }}_l} + {\pmb{x}}_{k - 1}^{(i)(j)},\quad l = 1, \cdots ,m $ (13)

式中, $m=2^n+1$, n为状态变量维数; ${\pmb{\delta }} = \left[ {\delta, \cdots, \delta } \right]$为自由变量, 其中${\delta ^2} \ge 1$.

$ {{\pmb{\xi }}_l} = \left\{ \begin{array}{ll} \bf{0} \\ {{{\left[ {\pmb{\delta }} \right]}_l}}\\ \end{array},\right.\ {\omega _l} = \left\{ \begin{array}{ll} {1 - \dfrac{1}{{{\delta ^2}}}}, & {l = 1} \\ {\dfrac{1}{{{2^n}{\delta ^2}}}}, & {l = 2, \cdots ,m} \\ \end{array} \right.\notag $

2) 时间更新, 状态转移函数为$f\left( \cdot \right)$

$ {\pmb{x}}_{\left. k \right|k - 1}^{(i)(j)(l)} = f({\pmb{x}}_{k - 1}^{(i)(j)(l)}) $ (14)
$ {{\hat x}}_{\left. k \right|k - 1}^{(i)(j)} = \sum\limits_{l = 1}^m {{\omega _l}{\pmb{x}}_{\left. k \right|k - 1}^{(i)(j)(l)}} $ (15)
$ {\pmb{S}}_{\left. k \right|k - 1}^{(i)(j)} ={\rm Tria}([\begin{array}{*{20}{c}} {{\pmb{\zeta }}_{\left. k \right|k - 1}^{(i)(j)}}{{{\pmb{S}}_{\pmb{Q}}}} \end{array}]) $ (16)
$ {\pmb{\zeta }}_{\left. k \right|k - 1}^{(i)(j)} =[{\pmb{x}}_{\left. k \right|k - 1}^{(i)(j)(1)} - {\pmb{x}}_{k - 1}^{(i)(j)}, \cdots, \nonumber\\ \quad {\pmb{x}}_{\left. k \right|k - 1}^{(i)(j)(m)} - {\pmb{x}}_{k - 1}^{(i)(j)}]{\pmb{W}} $ (17)
$ {\pmb{W}} = {\rm diag}\left\{\sqrt {{\omega _1}} , \cdots ,\sqrt {{\omega _m}} \right\} $ (18)

式(15) 中, ${{{\pmb{S}}_{\pmb{Q}}}}$为过程噪声方差平方根.

3) 状态更新

$ {{\pmb x}'}_{\left. k \right|k - 1}^{(i)(j)(l)} = {\pmb{S}}_{\left. k \right|k - 1}^{(i)(j)}{{{\xi }}_l} + {{\hat x}}_{\left. k \right|k - 1}^{(i)(j)},l = 1, \cdots ,m $ (19)
$ {\pmb{z}}_{\left. k \right|k - 1}^{(i)(j)(l)} = h({{\pmb x}'}_{\left. k \right|k - 1}^{(i)(j)(l)}) $ (20)
$ {\pmb{z}}_{\left. k \right|k - 1}^{(i)(j)} = \sum\limits_{l = 1}^m {{\omega _l}{\pmb{z}}_{\left. k \right|k - 1}^{(i)(j)(l)}} $ (21)
$ {\pmb{x}}_{\left. k \right|k}^{(i)(j)} = {{\hat x}}_{\left. k \right|k - 1}^{(i)(j)} + {{\pmb{K}}_k}({{\pmb{z}}^j} - {\pmb{z}}_{\left. k \right|k - 1}^{(i)(j)}) $ (22)
$ {{\pmb{K}}_k} = \dfrac {{\pmb{P}}_{\pmb{xz}}} {\frac{{\left( {{\pmb{S}}_{\left. k \right|k}^{(i)(j)} }\right)^{\rm T}}} {{\left( {{\pmb{S}}_{\left. k \right|k}^{(i)(j)} }\right)}}} $ (23)
$ {\pmb{\zeta }}_{\left. {{\pmb{z}},k} \right|k - 1}^{(i)(j)} = [{\pmb{z}}_{\left. k \right|k - 1}^{(i)(j)(1)} - {\pmb{z}}_{\left. k \right|k - 1}^{(i)(j)}, \cdots,\nonumber\\ \qquad {\pmb{z}}_{\left. k \right|k - 1}^{(i)(j)(m)} - {\pmb{z}}_{\left. k \right|k - 1}^{(i)(j)}]{\pmb{W}} $ (24)
$ {{\pmb{P}}_{{\pmb{xz}}}} = {\pmb{\zeta }}_{\left. k \right|k - 1}^{(i)(j)}{\left( {{\pmb{\zeta }}_{\left. {{\pmb{z}},k} \right|k - 1}^{(i)(j)}} \right)^{\rm T}} $ (25)
$ {\pmb{S}}_{\left. k \right|k}^{(i)(j)} = {\rm Tria}([\begin{array}{*{20}{c}} {{\pmb{\zeta }}_{\left. {{\pmb{z}},k} \right|k - 1}^{(i)(j)}}{{{\pmb{S}}_{\pmb{R}}}}\end{array}]) $ (26)

式(19) 中, $h\left( \cdot \right)$为单个目标的量测函数.式(23) 中, ${{\pmb{S}}_{\pmb{R}}}$为量测噪声方差平方根.式(21) 中, ${{\pmb{z}}_j}$是由最近邻域法得到的与积分点对应的观测值.

4) PHD预测和更新

$ {\pmb{m}}_{k\left| {k - 1} \right.}^{(i)} = \sum\limits_{j = 1}^N {{\omega _j}{\pmb{x}}_{\left. k \right|k}^{(i)(j)}} $ (27)
$ {\pmb{P}}_{k\left| {k - 1} \right.}^{(i)} = \sum\limits_{j = 1}^N {\Theta \left( {\left\{ {{\omega _j},{\pmb{x}}_{\left. k \right|k}^{(i)(j)}} \right\},{{\hat x}}_{k\left| {k - 1} \right.}^i} \right)} $ (28)

基于重要密度函数采样, 即从${\rm N}\left( {{\pmb{x}};{\pmb{m}}_{k|k - 1}^{(i)}, {\pmb{P}}_{k|k - 1}^{(i)}} \right)$采样得到$\left\{ {{\pmb{x}}_k^{(i)(j)}} \right\}_{j = 1}^N, i = 1, \cdots, {J_{k|k - 1}}$.然后对于任意的${\pmb{z}} \in {{\pmb{Z}}_k}$, 对目标的状态进行更新计算.

$ {w'}_k^{(i)(j)}(z)=\dfrac{{{N_k}\Bigg({{\pmb{z}};h({\pmb{x}}_k^{(i)(j)}), { R}} \Bigg)}} {{\pi_k^{(i)}({\pmb{x}}_k^{(i)(j)}|{{\pmb{x}}_{k-1}^{(i)},{\pmb Z}_{1:k}})}} \times\nonumber\\ \qquad \qquad\quad {\rm N} \left( {{\pmb{x}}_k^{(i)(j)};{\pmb{m}}_{k\left| {k - 1} \right.}^{(i)},{\pmb{P}}_{k\left| {k - 1} \right.}^{(i)}} \right) $ (29)
$ w_k^{(i)}({\pmb{z}}) = \frac{{{P_{D,k}}w_{\left. k \right|k - 1}^{(i)}\sum\limits_{j = 1}^N {{w'}_k^{(i)(j)}({\pmb{z}})} }}{{N{\kappa _k}({\pmb{z}}) + {P_{D,k}}\sum\limits_{i = 1}^{{J_{k\left| {k - 1} \right.}}} {w_{\left. k \right|k - 1}^{(i)}\sum\limits_{j = 1}^N {{w'}_k^{(i)(j)}({\pmb{z}})} }}} $ (30)
$ \qquad \quad m_k^{(i)}({\pmb{z}}) = \sum\limits_{j = 1}^N \bar {w'}_k^{(i)(j)}({\pmb{z}}){\pmb{x}}_k^{(i)(j)} $ (31)
$ {\pmb{P}}_k^{(i)}({\pmb{z}}) = \sum\limits_{j = 1}^N {\Theta \left( {\left\{ {\bar {w'}_k^{(i)(j)}({\pmb{z}}),{\pmb{x}}_k^{(i)(j)}} \right\},{\pmb{m}}_k^{(i)}({\pmb{z}})} \right)} $ (32)

其中, $\bar {w'}_k^{(i)(j)}({\pmb{z}}) = {{w'}_k^{(i)(j)}({\pmb{z}})} /{\sum\nolimits_{j = 1}^N {{w'}_k^{(i)(j)}({\pmb{z}})} }$.

2 ASRICP-PHD

低维度下, SRICP-PHD算法改善的效果不明显.鉴于初始化采样决定粒子使用效率, 影响滤波精度和稳定性, 因此本文针对采样手段提出了EPS, 由此得到ASRICP-PHD算法.

2.1 等概率采样

初始化采样目的在于通过高斯组合实现非高斯噪声估计. PRS无法避免粒子“聚集”的现象, QRS可以解决“粒子聚集”的问题[14], 然而, 得到的粒子在统计特性上更接近均匀分布, 会导致多数粒子无效.基于此, 本文提出等概率粒子采样法, 其本质是整体采样区域的划分和若干子区域粒子选取, 子区域粒子采样借鉴三阶容积准则[7].

n维的标准高斯分布${\pmb{x}} \sim {\rm N}\left( {{\pmb{x}};{\bf{0}}, {{{I}}_n}} \right)$中采样得到粒子集$\left\{ {{{\pmb{x}}_i}} \right\}_{i = 1}^N$, 其中$\bf{0}$为0向量, ${ {I}_n}$为单位矩阵, N$2n$的倍数.粒子集应能反映变量的统计特性, 关键在于采样粒子的概率分布情况, 将整体采样区域等概率划分为M个区域(对应的等概率区间为$[0, {P_1}), [{P_1}, {P_2}), \cdots, [{P_{M-1}}, 1]$, 则理论上各区域粒子采样数目应相同.首要工作是建立一维概率区间与多维采样区域之间的联系, 对于标准的高斯分布, 任取粒子${{{\pmb{x}}_i}}$, 其概率密度值为

$ {\rho _i} = \frac{1}{{\left( {{{\left( {2\pi } \right)}^n}\det ({{ {I}}_n})} \right)}^{\frac{1}{2}}} \exp \left( {\frac{{{{\pmb{x}}_i}{\pmb{x}}_i^{\rm T}}}{2}} \right) $ (33)

由式(33) 可知, n维的标准的高斯分布的等概率面是n维的球面, 球面上任意一点的2范数值是相等的, 如此便建立了一维概率区间与n维粒子分布区域之间的联系.为方便之后的计算, 现给出如下定义: 1) 等概率粒子采样区域称之为${\Omega _1}, {\Omega _2}, \cdots, {\Omega _M}$; 2) 任意点${{\pmb{x}}_i}$的2范数值为${r_i}$; 3) 半径为${r_i}$球面对应的概率值为$\phi \left( {{r_i}} \right) = {P_i}$; 4) 由$P_i$得到$r_i$的运算定义为${r_i} = {\phi ^{ - 1}}\left( {{P_i}} \right)$.

对任意的概率区间$[{P_{i-1}}, {P_i})$ ($i=2, \cdots, M-1$), 可用$[{r_{i-1}}, {r_i})$来表示其对应的粒子分布区域${\Omega _i}$, ${\Omega _i}$是由两个球面所形成的封闭区域; 概率区间$[0, {P_1})$对应的粒子分布区域${\Omega _1}$是一个球体, 球体的半径为${r_1}$; 概率区间$[{P_{M-1}}, {P_M})$对应的粒子分布区域${\Omega _M}$是没有外边界的, 内边界为$r_{M-1}$对应的球面.求解得到等概率粒子分布区域之后就是如何选取粒子的问题, 即粒子采样规则.

取任意的粒子分布区域${\Omega _i}$ ($i=2, \cdots, M-1$), 从中任意的选取一个粒子${\pmb{x}}_i^j$, 满足${r_{i- 1}} \le \left\| {{\pmb{x}}_i^j} \right\| \le {r_i}$, $\left\| \cdot \right\|$为2范数值.可以将${\pmb{x}}_i^j$表示为模与单位向量的乘积, 即${\pmb{x}}_i^j = r_i^j{\pmb{y}}$, 则$r_i^j = \left\| {{\pmb{x}}_i^j} \right\|$, $\pmb{y}$为任意方向的单位向量, 这两个参数是粒子采样的关键.可以发现, 上述的粒子分布区域与容积准则的积分区域很相似, 因此, 对于n维的粒子, $\pmb{y}$可以选取为${\left[{\bf{1}} \right]_j}\left( {1 \le j \le 2n} \right)$, ${\left[{\bf{1}} \right]_j} \in {{\bf {R}}^n}$, 且

$ {\left[ \bf{1} \right]_j} \in \left\{ {{{\pmb{e}}_i}, - {{\pmb{e}}_i}} \right\}_{i = 1}^n $ (34)

其中, $\pmb{e}_i$i个元素为1, 余下为零.

关于$r_i^j$, 考虑到状态估计误差逐渐减小的特点, 可令$r_i^j = \alpha {r_i} + \left( {1 - \alpha } \right){r_{i - 1}}$, 同理, 对于粒子分布区域${\Omega _1}$${\Omega _M}$, 由于区域的半封闭性质, $r_1^j$$r_M^j$选取其边界值即可.据上分析可以确定概率区间的个数, 即$M = N/2n$.由此可知积分点选取为

$ {\pmb{x}} = {\pmb{S}}_{k - 1}^{(i)}{\pmb{\xi }} + {\pmb{m}}_{k - 1}^{(i)} $ (35)

式中, ${\pmb{\xi }}$服从标准高斯分布.如此便建立了基于概率的粒子采样规则.

2.2 EPS近似误差分析

对于分布函数为$p\left( {\pmb{x}} \right)$的变量, 定义非线性变换函数为$g\left( {\pmb{x}} \right)$, 则积分近似的原理可表示为

$ { {I}}\left[ {g({\pmb{x}})} \right] = \int {g({\pmb{x}})p({\pmb{x}}){\rm d}{\pmb{x}}} \approx \sum\limits_{i = 1}^n {g({{\pmb{x}}_i})\omega ({{\pmb{x}}_i})} $ (36)

假设变量服从高斯分布N$\left( {{\pmb{x}};{\bar{\pmb x}}, {\pmb{S}}{{\pmb{S}}^{\rm T}}} \right)$, 则基于PRS和EPS的积分近似可分别表示为

$ { {I}}\left[ {g({\pmb{x}})} \right] \approx \sum\limits_{i = 1}^n {g({\pmb{S}}{ {I}_i} + {\bar{\pmb x}}){\omega _i}} $ (37)
$ { {I}}\left[ {g({\pmb{x}})} \right] \approx \sum\limits_{r = {r_1}}^{{r_m}} {\sum\limits_{i = 1}^n {\frac{1}{m}g(r{\pmb{S}}{{\bf 1}_i} + {\bar{\pmb x}})\omega ({{\bf 1}_i})} } $ (38)

式(37) 中$ {I}_i$n维标准高斯分布随机数.

考虑简单的一维高斯分布N$\left( {x;\bar x, \sigma } \right)$及非线性变换函数$g(x) = {x^2} + 2$, 真实积分值为

$ I\;\left[ {g(x)} \right] = \int {g(x)p(x){\rm d}x} =\nonumber\\ \quad \int {({x^2} + 2)} \frac{1}{{\sqrt {2\pi } \sigma }}\exp \left( {\frac{{ - {{(x - \bar x)}^2}}}{{2{\sigma ^2}}}} \right){\rm d}x $ (39)

${x = \sqrt 2 \sigma t + \bar x}$则上式可转换为

$ I\left[ {g(x)} \right] = \frac{1}{{\sqrt \pi }}\int 2 \exp \left( { - {{(t)}^2}} \right){\rm d}t + \nonumber\\ \qquad \frac{1}{{\sqrt \pi }}\int {(2{\sigma ^2}{t^2} + {{\bar x}^2})} \exp \left( { - {{(t)}^2}} \right){\rm d}t $ (40)

上式的计算需要借助于gamma函数的性质$\Gamma \left( {1/2} \right) = \sqrt \pi $, 因此上式可得:

$ \int {g(x)p(x){\rm d}x} = 2 + {\bar x^2} + {\sigma ^2} $ (41)

图 3图 4是不同统计特性下PRS与EPS积分误差比较, 误差用自定义变量$er$表示.对比可知EPS稳定性强, 收敛性好, 优势明显.

图 1 误差比较($N = 1 000$) Figure 1 Comparison of error ($N = 1 000$)
图 2 误差比较($N = 20 000$) Figure 2 Comparison of error ($N = 20 000$)
图 3 目标数目估计比较($q = 0.1,{\lambda _k} = 5; {\rm CV}$) Figure 3 Comparison of number estimation ($q = 0.1,{\lambda _k} = 5; {\rm CV}$
图 4 OSPA比较($q = 0.1,{\lambda _k} = 5; {\rm CV}$) Figure 4 Comparison of OSPA ($q = 0.1,{\lambda _k} = 5; {\rm CV}$
2.3 ASRICP-PHD算法流程

步骤 1. 初始化采样. $k-1$时刻, 依据EPS对每一个高斯分量${\rm N}({\pmb{x}};{\pmb{m}}_{k - 1}^{(i)}, {\pmb{P}}_{k - 1}^{(i)})$进行粒子采样得到粒子集$\left\{ {{\pmb{x}}_{k - 1}^{(i)(j)}} \right\}_{j = 1}^N$, 采样方式如第2.1节所述.

步骤 2. PHD预测.根据第1.2节中所述的SRICP计算存活目标的重要密度函数$\pi _k^{(i)}({{\pmb{x}}_k}\left| {{\pmb{x}}_{k - 1}^{(i)}, {{\pmb{Z}}_{1:k}}} \right.)$, 预测存活目标的状态, 根据新生目标随机集PHD预测新生目标状态.

步骤 3. PHD更新.根据观测集$\pmb{Z}_k$更新目标的PHD, 具体计算过程如第1.2节所述.

步骤 4. 筛选和融合处理.设置阈值滤除权值小的高斯分量, 融合接近的高斯分量[12].

2.4 算法复杂度分析

根据算法流程计算关键步运算复杂度(统计浮点操作数flops[15]).算法复杂度由初始化采样手段和非线性滤波算法积分点数目决定. PRS复杂度为O$( {2N{n^2}} )$, EPS存在粒子区域划分以及偏差的选择, 复杂度为O$( {3N{n^2}+{Nn}} )$.算法中, 嵌入式容积准则积分数目为${2^{n+1}}$, 容积准则积分点数目为$2n$, 无迹变换积分点个数为${2n}+1$, 高斯厄米特积分点个数为${m_1^n}$ (${m_1}>2$为整数).而滤波步中, 单个粒子积分点采样的复杂度为O$({7/6}{n^3}m + 3{n^2}m)$, 积分点预测复杂度为O$( {4{n^2}m - 2mn} )$, 粒子状态预测复杂度为O$( {mn} )$, 粒子平方根误差预测复杂度为O$( {{n^3} + 3{n^2} + 2mn} )$.积分点量测预测复杂度为O$( {4m{l^2} - 2ml} )$, 粒子量测预测复杂度为O$( {ml})$, 粒子量测平方根误差估计复杂度为O$( {{l^2} + 3{l^2} + 2ml} )$, 状态与量测互协方差估计复杂度为O$( {{l^2} + {l^2} + 2ml} )$.由此可知, 由粒子数决定的复杂度为$N{\rm O}( {7/6}{n^3}m + 7{n^2}m + mn + {n^3} + 3{n^2} + 4m{l^2} +3ml + 2{l^3} + 4{l^2}) $.分析可知高斯厄米特算法最为复杂, 嵌入式容积准则缺点在于积分点数目较多, EPS复杂度在可控范围内.粒子数相同时, ASRICP-PHD复杂度较高, 减少EPS所需的粒子数能降低算法复杂度, 前提是保障滤波估计精度.

3 仿真校验

考虑二维平面内的多目标跟踪.目标的状态变量取为${{\pmb{x}}_k} = {[{x_k}, {\dot x_k}, {y_k}, {\dot y_k}]^{\rm T}}$.令${P_{S, k}} = 0.99$, ${P_{D, k}} = 0.98$; 定义采样时间间隔${T = 1 {\rm s}}$; 杂波在监控区域内均匀分布, 杂波数目服从泊松分布, 期望值为${\lambda _k}$.以最优子模式分配(Optimal subpattern assignment, OSPA) 脱靶距离为指标比较算法性能[16], OSPA中, 参数p取为2, c取为100.在多次试验下, 定义“跟踪失败率”评价指标${P_f}$, 单次实验跟踪失败定义为:任意时刻OSPA大于设定阈值$offset$即为跟踪失败.各目标的运动情况如表 1所示, 其中${t_0}$${t_f}$分别为监测时间范围内目标运动起止时间.

表 1 初始化目标运动参数 Table 1 Motion parameters initialization

目标的运动模型和量测模型:

$ \left\{ {\begin{array}{*{20}{l}} {{{\pmb{x}}_k} = f\left( {{{\pmb{x}}_{k - 1}}} \right) + {{\pmb{w}}_k}}\\ {{{\pmb{z}}_k} = h\left( {{{\pmb{x}}_k}} \right) + {{\pmb{v}}_k}} \end{array}} \right. $ (42)

式(42) 中过程噪声${\pmb{w}_k}$与量测噪声${\pmb{v}_k}$均为零均值高斯白噪声, 误差协方差矩阵分别为为${\pmb{Q}_k}$${\pmb{R}_k}$.

实验中采用不同的运动模型及量测方程, 不考虑目标的衍生, 新生目标(泊松分布) 状态随机集为

$ {\gamma _k}({\pmb{x}}) = 0.1{\rm N}({\pmb{x}};{\pmb{m}}_{\pmb{\gamma }}^{(1)},{{\pmb{P}}_{\pmb{\gamma }}}) + 0.1{\rm N}({\pmb{x}};{\pmb{m}}_{\pmb{\gamma }}^{(2)},{{\pmb{P}}_{\pmb{\gamma }}}) $ (43)

式(43) 中, ${\pmb{m}}_{\pmb{\gamma }}^{(1)}$, ${\pmb{m}}_{\pmb{\gamma }}^{(2)}$取为表 1中目标3、4的状态值, 误差矩阵为${{\pmb{P}}_{\pmb{\gamma }}} = {\rm diag}\left\{ {{{\left[{100, 25, 100, 25} \right]}^{\rm T}}} \right\}$.

实验中产生真实轨迹时加上非高斯噪声误差, 定义零均值高斯噪声$\left\{ {{{\pmb{\nu }}_i}} \right\}_{i = 1}^3$, 非高斯噪声为

$ {{\pmb{\nu }}_{ nG}} =\sum\limits_{i = 1}^3 {{{\pmb{\nu }}_i}} + \underbrace {\sum {{{\pmb{\nu }}_i}{{\pmb{\nu }}_j}} }_{1 \le i \ne j \le 3} + \,{\pmb\nu }_1{{\pmb{\nu }}_2}{{\pmb{\nu }}_3} $ (44)
3.1 实验 1

匀速(Constant velocity, CV) 运动多目标跟踪.监控区域为$\left[{-300 \rm {m}, 300 \rm {m}} \right] \times \left[{-300 \rm {m}, 300 \rm {m}} \right]$, 直接量测量为目标xy方向位置信息.线性状态转移矩阵及过程噪声协方差阵参考文献[7], 过程噪声协方差阵相关参数q设置多变.

量测噪声${{\pmb{v}}_k}$同样为零均值高斯白噪声, 误差协方差矩阵为${ {R}} = {\rm diag}\left\{ {\left[{\sigma _x^2, \sigma _y^2} \right]} \right\}, {\sigma _x} = {\sigma _y} = 5$ m.实验1在均方根嵌入式容积粒子PHD中比较了EPS、QRS以及PRS的效果, 结果取60次的平均值, 表 2是算法“跟踪失败率”比较, 值越小表示算法性能越好, 囿于布局, 没有给出$q=1$时的$P_f$. 图 3图 6是算法多目标数目估计结果, 用自定义变量$Tn$表示, 图中未标示的直线代表真实的目标数目. 图 4图 7是算法OSPA比较, 其中针对EPS进行了减少采样粒子数的实验.为比较各采样手段的复杂度, 图 5图 8给出了不同采样手段下的单步运行时间, 用自定义变量$St$表示.对比可知, EPS能提高PHD算法的稳定性, 相比于另外两种采样手段能更好地改善算法性能.通过理论分析, EPS复杂度高于PRS, 然而由时间对比结果可知, 在相同粒子数下, 由EPS运算复杂度增加的时间成本相对于整体算法的时间成本较小, 图 5图 8中部分时刻的时间对比甚至呈现相反的结果, 原因是QRS和PRS的缺点会导致目标粒子权值较小, 杂波被识别为目标, 因此时间成本稍高.但无论何种采样手段, 相互之间时间成本差异不大.此外, 实验结果表明减少EPS粒子数的实验是成功的, 表明在保障精度的前提下, 合理减少EPS粒子数PHD算法性能仍能得到较好的改善.

表 2 $P_f$分析($offset = 40\,{\rm {m}}; {\rm CV}$) Table 2 Analysis of $P_f$ ($offset = 40\,{\rm {m}}; {\rm CV}$)
图 5 单步运行时间比较($q = 0.1,{\lambda _k} = 5; {\rm CV}$) Figure 5 Comparison of step time ($q = 0.1,{\lambda _k} = 5; {\rm CV}$)
图 6 目标数目估计比较($q = 1,{\lambda _k} = 5; {\rm CV}$) Figure 6 Comparison of number estimation ($q = 1,{\lambda _k} = 5; {\rm CV}$)
图 7 OSPA比较($q = 1,{\lambda _k} = 5; {\rm CV}$) Figure 7 Comparison of OSPA ($q = 1,{\lambda _k} = 5; {\rm CV}$)
图 8 单步运行时间比较($q = 1,{\lambda _k} = 5; {\rm CV}$) Figure 8 Comparison of step time ($q = 1,{\lambda _k} = 5; {\rm CV}$)
3.2 实验 2

比较了ASRICP-PHD与EPS支持下的ASRCP-PHD、PRS支持下的SRGHP-PHD以及QRS支持下平方根无迹粒子PHD (Square-root unscented particle PHD with Halton points, HSRUP-PHD) 对匀转弯(Constant turning, CT) 多目标跟踪效果.监视区域为$\left[{-400 \rm {m}, 200 \rm {m}} \right] \times \left[{-400 \rm {m}, 400 \rm {m}} \right]$, 非线性观测方程为

$ {\pmb{z}}_k^i = \left[ {\begin{array}{*{20}{c}} {\sqrt {x_k^2 + y_k^2} }\\ {\arctan \left(\dfrac {{y_k}} {{x_k}} \right)} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{\nu _{r,k}}}\\ {{\nu _{a,k}}} \end{array}} \right] $ (45)

式(45) 是单个探测器的量测方程, 其中量测噪声协方差为$R = {\rm diag}\{ [ \sigma _a^2, \sigma _r^2] \}$, $\sigma _a =0.0075 {\rm rad/s}$, ${\sigma _r} = 5$ m.部署两个探测器, 位置坐标分别为$( - 200 \rm {m}, -250 \rm {m} )$, $( - 250 {\rm m}, 200 {\rm m} )$.杂波强度多变.状态转移线性矩阵和过程噪声误差协方差阵与$\omega $直接相关[7] ($\omega $为转弯速率, 取为$\pi /36 {\rm rad/s}$).

图 9图 14以及表 3是不同噪声以及不同杂波强度下实验结果(粒子数取为200, 实验结果取60次平均值).分析结果可知EPS适应性强, 可改善多种算法性能. SRGHP-PHD在非高斯噪声环境下效果差, 原因是PRS会出现“粒子聚集”以至粒子使用效率低、杂波干扰大的问题, ASRGHP-PHD性能有所改善, 但效果依然差, 这表明$m_1=3$时, 基础算法GHP-PHD稳定性差, 增加积分点数目可能会提高精度, 然而对比时间结果可知相应时间成本高、计算复杂, 即使在“漏跟”情况下, SRGHP-PHD及ASRGHP-PHD时间成本依然高, 不具实时性. HSRUP-PHD在本次实验中效果较差, 调节尺度因子后改进的ASRUP-PHD效果相对较好, 表明EPS优于QRS, 后者过于注重粒子集的均匀分布而导致多数粒子无效.对比之下, ASRICP-PHD可相对较好地实现多目标跟踪.对比ASRICP-PHD与ASRCP-PHD结果, 两者OSPA相差较小, 原因是低维度下, 容积准则积分点偏差较小, 因而跟踪效果较好, 表明在低维度下, 仅依靠嵌入式容积准则并不能很好地改善PHD算法的性能.然而, 对比$P_f$值, 可知ASRICP-PHD性能上优于ASRCP-PHD.

图 9 目标数目估计比较($q = 3,{\lambda _k} = 10; {\rm CT}$) Figure 9 Comparison of number estimation ($q = 3,{\lambda _k} = 10; {\rm CT}$)
图 10 OSPA比较($q = 3,{\lambda _k} = 10; {\rm CT}$) Figure 10 Comparison of OSPA ($q = 3,{\lambda _k} = 10; {\rm CT}$)
图 11 单步运行时间比较($q = 3,{\lambda _k} = 10; {\rm CT}$) Figure 11 Comparison of step time ($q = 3,{\lambda _k} = 10; {\rm CT}$)
图 12 目标数目估计比较($q = 3,{\lambda _k} = 10; {\rm CT}$) Figure 12 Comparison of number estimation ($q = 3,{\lambda _k} = 10; {\rm CT}$)
图 13 OSPA比较($q = 3,{\lambda _k} = 10; {\rm CT}$) Figure 13 Comparison of OSPA ($q = 3,{\lambda _k} = 10; {\rm CT}$)
图 14 单步运行时间比较($q = 3,{\lambda _k} = 10; {\rm CT}$) Figure 14 Comparison of step time ($q = 3,{\lambda _k} = 10; {\rm CT}$)
表 3 $P_f$分析($offset = 100\,{\rm {m}}; {\rm CT}$) Table 3 Analysis of $P_f$ ($offset = 100\,{\rm {m}}; {\rm CT}$)

综上分析可知, 针对邻近多目标跟踪, ASRICP-PHD算法能够提高目标状态估计的精度, 等概率采样手段可改善PHD算法性能, 同时在满足精度的需求下, 可合理减少粒子数以降低运算复杂度.

4 结论

针对非线性多目标跟踪提出了ASRICP-PHD算法.新算法引入了嵌入式容积粒子滤波算法对多目标状态集的PHD进行估计, 解决了容积粒子PHD “粒子溢出"的问题, 在一定程度上提高了PHD算法的精度.同时, 针对传统伪随机采样和拟随机采样的不足, 提出了等概率采样方法, 使得采样得到的粒子充满了有效的积分空间, 且具有较强的对称分布特性, 明显提高了ASRICP-PHD算法的精度.最后证明了ASRICP-PHD算法在保障精度的前提下, 可以通过减少粒子数目降低时间成本.后续继续对PHD更新算法进行研究, 以减少不必要运算, 提高时效性.同时, 需深入研究未知信噪比、未知噪声以及未知杂波强度条件下多目标跟踪技术, 以解决实际的目标跟踪难题.

参考文献
1 Li B, Pang F W. Improved probability hypothesis density filter for multitarget tracking. Nonlinear Dynamics, 2014, 76 (1): 367–376. DOI:10.1007/s11071-013-1131-1
2 Zhan Rong-Hui, Liu Sheng-Qi, Ou Jian-Ping, Zhang Jun. Improved multitarget track before detect algorithm using the sequential Monte Carlo probability hypothesis density filter. Journal of Electronics & Information Technology, 2014, 36 (11): 2593–2599.
( 占荣辉, 刘盛启, 欧建平, 张军. 基于序贯蒙特卡罗概率假设密度滤波的多目标检测前跟踪改进算法. 电子与信息学报, 2014, 36 (11): 2593–2599. )
3 Xu L L, Ökten G. High-performance financial simulation using randomized quasi-Monte Carlo methods. Quantitative Finance, 2015, 15 (8): 1425–1436. DOI:10.1080/14697688.2015.1032549
4 Guo D, Wang X D. Quasi-Monte Carlo filtering in nonlinear dynamic systems. IEEE Transactions on Signal Processing, 2006, 54 (6): 2087–2098. DOI:10.1109/TSP.2006.873585
5 Yang Jin-Long, Ji Hong-Bin, Liu Jin-Mang. Gauss-Hermite particle PHD filter for bearings-only multi-target tracking. Systems Engineering and Electronics, 2013, 35 (3): 457–462.
( 杨金龙, 姬红兵, 刘进忙. 高斯厄米特粒子PHD被动测角多目标跟踪算法. 系统工程与电子技术, 2013, 35 (3): 457–462. )
6 Wang Hua-Jian, Jing Zhan-Rong. Probability hypothesis density filter based on cubature rule and its application to multi-target tracking. Transactions of Beijing Institute of Technology, 2014, 34 (12): 1304–1309.
( 王华剑, 景占荣. 基于容积原则的概率假设密度滤波算法. 北京理工大学学报, 2014, 34 (12): 1304–1309. )
7 Zhan Rong-Hui, Zhang Jun. Nonlinear Filtering Theory with Target Tracking Application. Beijing: National Defend Industry Press, 2013.
( 占荣辉, 张军. 非线性滤波理论与目标跟踪应用. 北京: 国防工业出版社, 2013. )
8 Zhang X C. Cubature information filters using high-degree and embedded cubature rules. Circuits, Systems, and Signal Processing, 2014, 33 (6): 1799–1818. DOI:10.1007/s00034-013-9730-0
9 Lian Feng, Han Chong-Zhao, Li Chen. Multiple-model GM-CBMeMBer filter and track continuity. Acta Automatica Sinica, 2014, 40 (2): 336–347.
( 连峰, 韩崇昭, 李晨. 多模型GM-CBMeMBer滤波器及航迹形成. 自动化学报, 2014, 40 (2): 336–347. )
10 Li Cui-Yun, Jiang Zhou, Ji Hong-Bin, Cao Xiao-Nan. GMP-PHD filter based on quasi-Monte Carlo in unknown clutter. Control and Decision, 2014, 29 (11): 1997–2001.
( 李翠芸, 江舟, 姬红兵, 曹潇男. 基于拟蒙特卡罗的未知杂波GMP-PHD滤波器. 控制与决策, 2014, 29 (11): 1997–2001. )
11 Mahler R P S. Multitarget Bayes filtering via first-order multitarget moments. IEEE Transactions on Aerospace and Electronic Systems, 2003, 39 (4): 1152–1178. DOI:10.1109/TAES.2003.1261119
12 Wang Y, Meng H D, Liu Y M, Wang X Q. Collaborative penalized Gaussian mixture PHD tracker for close target tracking. Signal Processing, 2014, 102 : 1–15. DOI:10.1016/j.sigpro.2014.01.034
13 Zhang Xin-Chun, Guo Cheng-Jun. Square-root imbedded cubature Kalman filtering. Control Theory & Applications, 2013, 30 (9): 1116–1121.
( 张鑫春, 郭承军. 均方根嵌入式容积卡尔曼滤波. 控制理论与应用, 2013, 30 (9): 1116–1121. )
14 Wang X, Hickernell F J. Randomized Halton sequences. Mathematical and Computer Modelling, 2000, 32 (7-8): 887–899. DOI:10.1016/S0895-7177(00)00178-3
15 Grewal M S, Andrew A P. Kalman Filtering:Theory and Practice Using MATLAB. Hoboken, N.J.: John Wiley & Sons, 2008.
16 Yang Feng, Wang Yong-Qi, Liang Yan, Pan Quan. A survey of PHD filter based multi-target tracking. Acta Automatica Sinica, 2013, 39 (11): 1944–1956.
( 杨峰, 王永齐, 梁彦, 潘泉. 基于概率假设密度滤波方法的多目标跟踪技术综述. 自动化学报, 2013, 39 (11): 1944–1956. DOI:10.3724/SP.J.1004.2013.01944 )