-
摘要: 针对基于概率假设密度算法(Probability hypothesis density,PHD)的非线性多目标跟踪精度低、滤波发散等问题,提出了一种新的PHD算法——改进的均方根嵌入式容积粒子PHD算法(Advanced square-root imbedded cubature particle PHD,ASRICP-PHD).新的算法在初始化采样时将整个采样区域等概率划分为若干个区域,然后利用既定的准则从每个区域抽取粒子,并利用均方根嵌入式容积滤波方法对每个粒子进行滤波,来拟合重要密度函数,预测和更新多目标状态的PHD.仿真结果表明该算法能对多目标进行有效跟踪,相比拟随机采样法和伪随机采样,等概率采样的方法在多目标位置估计和数目估计上有更高的精度.
-
关键词:
- 多目标跟踪 /
- 概率假设密度 /
- 均方根嵌入式容积滤波 /
- 等概率采样
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. -
概率假设密度算法(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) 其中, 为多目标状态集的预测概率密度; , 分别为状态转移函数和量测似然函数; 为多目标状态集的后验概率密度.
难点在于如何表示数目变化的多目标状态以非线性积分计算. 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{\Gamma }}_k}$ , 为目标对应的量测; ${{\pmb{K}}_k}$ 为虚警杂波.
针对式(1) 和(2) 中的非线性积分, Mahler提出利用多目标状态一阶矩近似求解.假设k时刻先验PHD和后验PHD分别为、, 则:
$ {\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; 为衍生目标状态集PHD; 为单目标的生存概率; 为当前时刻对状态为 $\pmb{x}$ 的目标的探测概率; 为杂波集的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{P}}_{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)}}$ ; , 分别为预测得到的高斯分量均值和协方差.式(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, 采样得到粒子集, 同时根据Tria变换求解对应状态协方差的平方根, 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{\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) 基于重要密度函数采样, 即从采样得到.然后对于任意的, 对目标的状态进行更新计算.
$ {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) 其中, .
2. ASRICP-PHD
低维度下, SRICP-PHD算法改善的效果不明显.鉴于初始化采样决定粒子使用效率, 影响滤波精度和稳定性, 因此本文针对采样手段提出了EPS, 由此得到ASRICP-PHD算法.
2.1 等概率采样
初始化采样目的在于通过高斯组合实现非高斯噪声估计. PRS无法避免粒子“聚集”的现象, QRS可以解决“粒子聚集”的问题[14], 然而, 得到的粒子在统计特性上更接近均匀分布, 会导致多数粒子无效.基于此, 本文提出等概率粒子采样法, 其本质是整体采样区域的划分和若干子区域粒子选取, 子区域粒子采样借鉴三阶容积准则[7].
从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) 等概率粒子采样区域称之为; 2) 任意点 ${{\pmb{x}}_i}$ 的2范数值为 ${r_i}$ ; 3) 半径为 ${r_i}$ 球面对应的概率值为; 4) 由 $P_i$ 得到 $r_i$ 的运算定义为.
对任意的概率区间 $[{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$ , 满足, 为2范数值.可以将 ${\pmb{x}}_i^j$ 表示为模与单位向量的乘积, 即, 则 $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$ 选取其边界值即可.据上分析可以确定概率区间的个数, 即.由此可知积分点选取为
$ {\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, 则基于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及非线性变换函数 $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稳定性强, 收敛性好, 优势明显.
2.3 ASRICP-PHD算法流程
步骤 1. 初始化采样. $k-1$ 时刻, 依据EPS对每一个高斯分量进行粒子采样得到粒子集, 采样方式如第2.1节所述.
步骤 2. PHD预测.根据第1.2节中所述的SRICP计算存活目标的重要密度函数, 预测存活目标的状态, 根据新生目标随机集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} )$ .由此可知, 由粒子数决定的复杂度为.分析可知高斯厄米特算法最为复杂, 嵌入式容积准则缺点在于积分点数目较多, EPS复杂度在可控范围内.粒子数相同时, ASRICP-PHD复杂度较高, 减少EPS所需的粒子数能降低算法复杂度, 前提是保障滤波估计精度.
3. 仿真校验
考虑二维平面内的多目标跟踪.目标的状态变量取为.令, ${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目标 位置(m) 速度(m·s-1) t0 (s) tf (s) 1 (-250, 250) (5, -5) 1 45 2 (-250, -250) (5, 5) 1 50 3 (-160, 160) (7, -9) 17 50 4 (-160, -160) (9, 7) 27 50 目标的运动模型和量测模型:
$ \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{\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]$ , 直接量测量为目标x、y方向位置信息.线性状态转移矩阵及过程噪声协方差阵参考文献[7], 过程噪声协方差阵相关参数q设置多变.
量测噪声 ${{\pmb{v}}_k}$ 同样为零均值高斯白噪声, 误差协方差矩阵为 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$ 分析()Table 2 Analysis of $P_f$ ()采样手段 失败次数 仿真次数 Pf QRS 20 60 0.334 EPS (N=96) 8 60 0.134 EPS (N=200) 3 60 0.05 PRS 35 60 0.584 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) 多目标跟踪效果.监视区域为, 非线性观测方程为
$ {\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) 是单个探测器的量测方程, 其中量测噪声协方差为, $\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 $ 为转弯速率, 取为).
图 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.
表 3 $P_f$ 分析()Table 3 Analysis of $P_f$ ()算法 失败次数 仿真次数 Pf ASRICP-PHD 17 60 0.28 ASRCP-PHD 32 60 0.502 HSRUP-PHD 50 60 0.83 SRGHP-PHD 55 60 0.916 综上分析可知, 针对邻近多目标跟踪, ASRICP-PHD算法能够提高目标状态估计的精度, 等概率采样手段可改善PHD算法性能, 同时在满足精度的需求下, 可合理减少粒子数以降低运算复杂度.
4. 结论
针对非线性多目标跟踪提出了ASRICP-PHD算法.新算法引入了嵌入式容积粒子滤波算法对多目标状态集的PHD进行估计, 解决了容积粒子PHD “粒子溢出"的问题, 在一定程度上提高了PHD算法的精度.同时, 针对传统伪随机采样和拟随机采样的不足, 提出了等概率采样方法, 使得采样得到的粒子充满了有效的积分空间, 且具有较强的对称分布特性, 明显提高了ASRICP-PHD算法的精度.最后证明了ASRICP-PHD算法在保障精度的前提下, 可以通过减少粒子数目降低时间成本.后续继续对PHD更新算法进行研究, 以减少不必要运算, 提高时效性.同时, 需深入研究未知信噪比、未知噪声以及未知杂波强度条件下多目标跟踪技术, 以解决实际的目标跟踪难题.
-
表 1 初始化目标运动参数
Table 1 Motion parameters initialization
目标 位置(m) 速度(m·s-1) t0 (s) tf (s) 1 (-250, 250) (5, -5) 1 45 2 (-250, -250) (5, 5) 1 50 3 (-160, 160) (7, -9) 17 50 4 (-160, -160) (9, 7) 27 50 表 2 $P_f$ 分析( $offset = 40\,{\rm {m}}; {\rm CV}$ )
Table 2 Analysis of $P_f$ ( $offset = 40\,{\rm {m}}; {\rm CV}$ )
采样手段 失败次数 仿真次数 Pf QRS 20 60 0.334 EPS (N=96) 8 60 0.134 EPS (N=200) 3 60 0.05 PRS 35 60 0.584 表 3 $P_f$ 分析( $offset = 100\,{\rm {m}}; {\rm CT}$ )
Table 3 Analysis of $P_f$ ( $offset = 100\,{\rm {m}}; {\rm CT}$ )
算法 失败次数 仿真次数 Pf ASRICP-PHD 17 60 0.28 ASRCP-PHD 32 60 0.502 HSRUP-PHD 50 60 0.83 SRGHP-PHD 55 60 0.916 -
[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] 占荣辉, 刘盛启, 欧建平, 张军.基于序贯蒙特卡罗概率假设密度滤波的多目标检测前跟踪改进算法.电子与信息学报, 2014, 36(11):2593-2599 http://www.cnki.com.cn/Article/CJFDTOTAL-DZYX201411009.htmZhan 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 http://www.cnki.com.cn/Article/CJFDTOTAL-DZYX201411009.htm [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] 杨金龙, 姬红兵, 刘进忙.高斯厄米特粒子PHD被动测角多目标跟踪算法.系统工程与电子技术, 2013, 35(3):457-462 http://www.cnki.com.cn/Article/CJFDTOTAL-XTYD201303003.htmYang 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 http://www.cnki.com.cn/Article/CJFDTOTAL-XTYD201303003.htm [6] 王华剑, 景占荣.基于容积原则的概率假设密度滤波算法.北京理工大学学报, 2014, 34(12):1304-1309 http://www.cnki.com.cn/Article/CJFDTOTAL-BJLG201412018.htmWang 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 http://www.cnki.com.cn/Article/CJFDTOTAL-BJLG201412018.htm [7] 占荣辉, 张军.非线性滤波理论与目标跟踪应用.北京:国防工业出版社, 2013.Zhan Rong-Hui, Zhang Jun. Nonlinear Filtering Theory with Target Tracking Application. Beijing:National Defend Industry Press, 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] 连峰, 韩崇昭, 李晨.多模型GM-CBMeMBer滤波器及航迹形成.自动化学报, 2014, 40(2):336-347 http://www.aas.net.cn/CN/abstract/abstract18295.shtmlLian Feng, Han Chong-Zhao, Li Chen. Multiple-model GM-CBMeMBer filter and track continuity. Acta Automatica Sinica, 2014, 40(2):336-347 http://www.aas.net.cn/CN/abstract/abstract18295.shtml [10] 李翠芸, 江舟, 姬红兵, 曹潇男.基于拟蒙特卡罗的未知杂波GMP-PHD滤波器.控制与决策, 2014, 29(11):1997-2001 http://www.cnki.com.cn/Article/CJFDTOTAL-KZYC201411015.htmLi 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 http://www.cnki.com.cn/Article/CJFDTOTAL-KZYC201411015.htm [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] 张鑫春, 郭承军.均方根嵌入式容积卡尔曼滤波.控制理论与应用, 2013, 30(9):1116-1121 http://www.cnki.com.cn/Article/CJFDTOTAL-KZLY201309007.htmZhang Xin-Chun, Guo Cheng-Jun. Square-root imbedded cubature Kalman filtering. Control Theory & Applications, 2013, 30(9):1116-1121 http://www.cnki.com.cn/Article/CJFDTOTAL-KZLY201309007.htm [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] 杨峰, 王永齐, 梁彦, 潘泉.基于概率假设密度滤波方法的多目标跟踪技术综述.自动化学报, 2013, 39(11):1944-1956 doi: 10.3724/SP.J.1004.2013.01944Yang 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 doi: 10.3724/SP.J.1004.2013.01944 期刊类型引用(3)
1. 王建华,冉煜琨. 基于最优几何匹配的时间连贯3D动画重建. 电子设计工程. 2021(20): 36-42 . 百度学术
2. 张煌军,徐雪松,张文清,刘瑞. 基于SCKF的4旋翼无人机的姿态估计. 江西师范大学学报(自然科学版). 2019(02): 154-159 . 百度学术
3. 郑婷婷,杨旭升,张文安,俞立. 一种高斯渐进滤波框架下的目标跟踪方法. 自动化学报. 2018(12): 2250-2258 . 本站查看
其他类型引用(2)
-