目标跟踪中的传感器控制 ^[1] (亦称为传感器管理)的核心思想是依据一定的最优准则, 选择传感器的工作方式及运行参数以确保能够从全局使目标跟踪的性能达到最优. 它本质上属于最优非线性控制问题 ^[2]. 对于一个自主可控的单传感器, 它的控制运行轨迹会影响到目标跟踪系统的性能, 而控制它的轨迹的主要目的是依据某种评价准则使得多目标跟踪系统整体跟踪性能达到最优. 这类问题的解决通常在部分可观测马尔科夫决策过程(Partially observed Markov decision processes, POMDPs) 的理论框架下进行 ^[3]. 但由于目标运动参数未知(状态空间上的不确定性), 目标往往又是机动运行的(模型的不确定性), 跟踪任务中不可避免的会同时出现多个目标, 加上传感器量测误差、传感器漏检和杂波的共存带来的量测起源的不确定性, 传感器对各个目标的最优运行轨迹相互制约, 综合以上这些因素使得求解机动多目标的传感器轨迹控制的优化策略其实是非常困难的.

基于随机有限集(Random finite set, RFS)的多目标跟踪算法由于提供了一种无需进行数据关联即可解决多目标跟踪问题的理论框架而广受关注. 该方法将多目标状态和量测利用RFS进行建模, 从集值估计的整体角度描述和解决多目标跟踪问题. 这种集值估计的理论体系为传感器控制的研究提供了很大的便利, 可以从集值随机变量所描述的多目标跟踪的整体状态信息对传感器控制的决策做指导. Mahler已经利用有限集统计(Finite set statistics, FISST)理论 ^[1]将基于信息论的传感器控制方法推广到更一般的情况, 使得解决更为复杂的传感器系统控制成为可能. 近几年, 国外已有一些学者基于FISST的多目标递推贝叶斯滤波器求解并最终实现传感器的控制 ^[4-7]. 但是, 目前对基于FISST方法的机动多目标的传感器控制并没有得到系统性研究.

本文的主要内容是采用信息论方法研究机动多目标的传感器控制问题, 对机动多目标概率密度提出了一种目标导向的多伯努利(Target-oriented multi-Bernoulli, TOMB)粒子采样方法, 借助该采样方法去完成相应控制方案评价函数的求取. 文中利用多模型(Multi-model, MM)方法去近似描述机动目标运动模式的不确定性, 结合多模型多伯努利滤波器去递推描述机动多目标状态的不确定性, 研究利用Bhattacharyya距离 ^[8] (Bhattacharyya distance)作为评价函数, 利用TOMB采样方法求解评价函数并对传感器的运动做决策规划. 在实验部分将对提出的传感器运行控制和机动多目标跟踪策略进行详细的分析研究.

1.   问题描述

假定 k 时刻有$n_x(k)$个目标, 它们的状态分别为${ x}_{1, k}, \cdots, { x}_{n_x(k), k}$, 且每一个状态取值在目标状态空间$\mathcal{X}\subseteq R^{n_x}$. 同时, 传感器在k 时刻接收到$n_z(k)$个量测为$ z_{1, k}, \cdots, z_{n_z(k), k}$, 每一个量测取值在观测空间$\mathcal Z \subseteq R^{{n_z}}$. 那么, k 时刻多目标状态${{X}_{k}}=\{{ x}_{1, k}, \cdots, { x}_{n_x(k), k}\} \in \mathcal{F}(X) $和多目标量测${{Z}_{k}}=\{{ z}_{1, k}, \cdots, { z}_{n_{ z}(k), k}\} \in \mathcal{F}(Z) $都分别构成一个RFS, 其中, $\mathcal{F}(X)$代表所有$\mathcal{X}$ 的有限子集组成的集合, $\mathcal{F}(Z)$代表所有$\mathcal{Z}$ 的有限子集组成的集合. 考虑到目标的存活、消亡和新生, k 时刻目标数目$n_x(k)$是时变的. 直角坐标系下的单目标状态向量${ {x}_{k}} \in {{X}_{k}}$, 它可表示为

假设 $k-1$ 时刻, 状态$ x_{k-1} \in X_{k-1}$的目标以存活概率$p_{S, k}( x_{k-1})$继续生存在 k 时刻, 若不考虑衍生目标, 则多目标动态系统可建模为 ^[9-10]

其中, ${{S_{k|k -1}}}( x_{k -1})$是从 $k-1$ 时刻状态$ x_{k -1}$演化而来的 k 时刻存活目标状态的RFS, ${\Gamma _k}$是 k 时刻新生目标状态的RFS. 假设所有存活和新生目标的状态演变方程为

其中, ${{f}_{k}}$ 是系统的状态转移矩阵, ${ w}_{k}$是过程噪声向量. 对于一个运动模式未知的机动目标, 显然${{f}_{k}}$是不易获取的. 本文研究利用MM方法近似机动目标的运动模式, 并在第3节详细推导并给出序贯蒙特卡罗多模型多目标多伯努利滤波器的一般实现形式.

另外, 对于 k 时刻目标状态$ x_{k} \in X_{k}$, 传感器以检测概率${{p}_{D, k}}( x_k)$ 检测到它对应的量测${ {z}_{k}}$, ${ {z}_{k}}\in\mathcal{Z}_k$, 若量测噪声为${{ v}_{k}}$, 假设量测方程具有如下形式:

其中, 传感器在 k 时刻的位置${ x}_{s, k}=[{{x}_{s, k}}, {{y}_{s, k}}]^{\rm T}$是精确可测量的. 传感器对应目标量测是一个RFS, 可表示为$\Theta _k( x_k)$. 假设传感器在 k 时刻接收到的杂波表示为 ${{K}_{k}}$, 该集合为已知强度为${{\lambda }_{C}}$的泊松RFS. 则传感器在 k 时刻所接收到的多目标量测就可建模为 ^[9-10]

本文传感器控制的目的是, 基于连续的决策过程能够在线控制传感器的实时位置${ x}_{s, k}$, 对应该位置传感器会收集到对应时刻的量测集${Z_k}$. 假定过去时刻所有的传感器控制决策和量测集已知, 以下将详细研究如何利用信息论方法对传感器下一时刻的运动(或位置${ x}_{s, k+1}$)做决策.

2.   基于信息论的传感器控制

本文的传感器控制的研究在POMDP理论框架下进行. POMDP包含几个要素. 首先是当前状态信息的不确定性描述, RFS框架下即可用k时刻的多目标后验概率密度函数${{p}_{k}}(X_k |{{Z}_{1:k}})$表示. 其次, POMDP包括一套可实现的传感器控制集合${U}_k$, 每一个传感器控制${{\nu}}\in {U}_k$决定传感器将来的位置. 而最后一个关键的因素就是针对每一个传感器运动 ${\nu}$ 给定一个相应评价函数$\mathcal{R}(\nu)$. 基于信息论方法, 本文选择将整体信息增益的测度作为传感器控制的评价指标, 即将传感器控制后的多目标后验概率密度和多目标先验概率密度之间信息增量作为评价.

将${U}_{{k+1}:{k+H}}$表示为向后H步总的控制方案的集合, $H\geq 1$. 最优控制序列可简写为$ u_{{k+1}:{k+H}} \mathop = \limits^{\text{abbr.}} u_H$, 则该序列按以下准则确定:

其中, ${Z_{k + 1:k + H}}({\nu _H})$是由控制方案${\nu}_{H}$决定的量测集.

本文选用Bhattacharyya距离 ^[8]作为多目标概率密度的信息增益的评价函数, 距离越大说明传感器控制所获取信息越大. 评价函数表示为

其中, $BC({{p}_{k+H|k+H}}, {{p}_{k+H|k}})$是对应分布间的Bhattacharyya系数, $p_{k+H|k}(X)$是H步预测多目标概率密度, $p_{k+H|k+H}(X;{\nu}_H)$是采用控制${\nu}_H$ 后的后验多目标概率密度. 这两个概率密度的计算可连续递推H步的最优多目标贝叶斯递推滤波公式进行计算. 最优多目标贝叶斯递推滤波公式为

其中, ${{f_{k|k -1}}} ({X_k}|X_{k-1})$和${g_k}({Z_k}|{X_k})$分别是多目标状态转移密度和多目标似然函数.

为做出最优传感器控制决策, 对于式(7)应考虑$H \to \infty $时的多目标概率密度, 但是, 精确的目标运动模型难以获取, 过程噪声的存在使得较高的H就会显著影响到多步预测的多目标概率密度的准确性, 而机动目标的模型本身存在着更大的不确定性. 另外, 控制集合${U}_{{k+1}:{k+H}}$的势会伴随着H 的增大呈指数增长. 为了简化计算, 我们选择一步控制方法, 即$H=1$, 该方案又被称为“Myopic” 方案 ^[11].

式(7)评价函数的计算必须要求解式(8)和(9), 但它们并不存在闭式解. 文献^[4]利用多目标状态粒子$\{X_k^i\}_{i=1}^N$ 的采样去近似相应的多目标概率密度. 但对于随机集合的采样, 不仅包含目标个数的采样, 而且包含目标状态的采样, 运算处理非常复杂. 本文将在第4节提出一种基于目标导向的多伯努利(TOMB) 的粒子采样方法, 直接在单目标状态空间上进行采样, 逼近相应的多目标概率密度.

另外, 需要考虑如何在"Myopic"方案下根据 k 时刻每一个可能的控制${\nu}$生成量测集合$Z_k(\nu)$. 依据检测概率p_D和杂波强度${{\lambda }_{C}}$ 按照式(4)和(5)生成量测是一种直观的做法, 但是这样会带来极大的计算负担. 我们借鉴文献^{[4, 11]}的做法, 在不考虑杂波和量测噪声, 以及检测概率$p_D=1$ 的理想情况下, 对每一个${\nu}$仅生成一个量测集, 这个理想量测集合可表示为

其中, $ x_{s, k}(\nu)$是控制${\nu}$对应的传感器位置, $\hat{X}_{k|k-1}$是预测的多目标状态集合的估计.

3.   机动多目标跟踪的多伯努利滤波器的实现

3.1   多伯努利随机有限集

若X是状态空间$\mathcal{X}$上的(单)伯努利RFS, 它可用单目标存在概率r和单目标状态分布p来联合表示, 而X的势分布是一个参数为r的伯努利分布. 若$\emptyset$表示空集, 伯努利RFS的概率密度为 ^[12]

若X是状态空间$\mathcal{X}$上的多伯努利RFS, 它是一个确定数目且相互独立的伯努利RFS的集合. 组成X的第i个伯努利RFS表示为 ${{X}}^{(i)}$, 它的存在概率为${{r}}^{(i)}$, 概率密度为${{p}}^{(i)}$, $i=1$, M是伯努利RFS的索引, 那么$X={{\cup }_{i=1}}^{M}{{X}^{(i)}}$. 显然, 它的势平均为$\sum\limits_{i=1}^{M}{{{r}^{(i)}}}$. 则X的概率密度$\pi$可表示为$\pi (\varnothing )=\prod\limits_{j=1}^{M}{(1-{{r}^{(j)}})}$, 且

为描述方便, 可以将上述密度简写为参数集表达形式 ^[13], 即$\pi {\rm{ = }}\{ ({r^{(i)}},{\rm{ }}{p^{(i)}})\} _{i = 1}^M$.

3.2   势均衡多目标多伯努利滤波器

区别于PHD^[14]和CPHD^[15]多目标矩递推滤波器, 多目标多伯努利(Multi-target multi-Bernoulli, MeMBer)滤波器直接近似递推了多目标状态的后验概率分布, 使得多目标跟踪问题的求解及其状态的递推估计显得更为直观. Vo等已证明MeMBer滤波器存在势的过估计 ^[13], 且基于后验多目标概率密度的概率生成泛函(Probability generating functional, PGFl)给出一种势修正策略, 提出势均衡多目标多伯努利(Cardinality balanced MeMBer, CBMeMBer)滤波器. 相对于传统的SMC-PHD和传统的SMC-CPHD滤波器, SMC-CBMeMBer滤波器可以更方便可靠地进行状态的提取和估计 ^[13]. 以下先给出CBMeMBer滤波器的递推公式.

1) 预测步

假设$k-1$ 时刻后验多目标多伯努利密度表示为

则预测的多目标密度也是一个多伯努利密度

其中, $\{ ({r_{\Gamma ,k}}^{(i)},{p_{\Gamma ,k}}^{(i)})\} _{i = 1}^{{M_{\Gamma ,k}}}$是k时刻新生多伯努利密度. 而存活目标预测伯努利密度参数

其中, ${f_{k|k -1}}(\cdot |{\zeta})$是k时刻先前状态为${\zeta}$条件下的单目标状态转移密度, ${p_{S, k}}(\zeta)$是k时刻先前状态为${\zeta}$条件下的目标存活概率.

2) 更新步

令CBMeMBer滤波器预测多目标多伯努利密度为

那么后验多目标密度可用多伯努利密度近似如下

其中, ${Z_k}$是k时刻量测集, $\left\{ \left( {{r}_{L, k}}^{(i)}, p_{L, k}^{(i)} \right) \right\}_{i=1}^{{{M}_{k|k-1}}}$ 是继承航迹(漏检)的多伯努利密度参数集, 若检测概率表示为$p_{D, k}( x)$, 则

$\{({{r}_{U, k}}(z), {{p}_{U, k}}(x;z))\}z\in {{Z}_{k}}$ 是量测更新的多伯努利密度的参数集, 有:

机动多目标滤波器的实现} CBMeMBer滤波器无法直接求取闭式解, 并且考虑到机动多目标跟踪的要求, 以下给出序贯蒙特卡罗多模型CBMeMBer (Sequential Monte Carlo MMCBMeMBer, SMC-MMCBMeMBer)滤波器的实现形式.

对于一个机动目标, 通常可用多个目标模型的混合去近似它的运动模型, 即多模型(MM)方法. 模型参数定义为$\{{{d}_{k}}\}\in \mathcal M$, $\mathcal M$ 是所有可能的目标模型构成的集合, 考虑一个跳变的马尔科夫状态空间系统(Jump Markov state space systems, JMSS), 目标模型$\{{{d}_{k}}\}$根据一个确定的马尔科夫链随时间进行演变. 首先, 扩维目标的新状态

那么对于式(13), 在状态扩维的情况下, 假设$k-1$ 时刻概率密度为

其中, ${{\delta}}({x})$是关于扩维状态${x}$的狄拉克$\delta$函数, 扩维状态粒子

1) 预测步

将式(25)带入式(15), 得到伯努利RFS预测的存在概率

JMSS的状态演化模型可表示为

其中, $p(d_{k}|{d_{k-1}})$是模型转移概率的马尔科夫链.

先来采样模型粒子

根据状态采样的建议分布$q{{_{k}^{(i)(x}}_{k}}|{{x}_{k-1}}, {{d}_{k}}^{(i, j)})$, 采样存活目标状态粒子

新粒子可表示为$\widetilde{x}_{P, k|k-1}^{(i, j)}={{[x_{P, k|k-1}^{(i, j)\text{T}}, d_{k}^{(i, j)}]}^{\text{T}}}$, 那么预测概率密度可近似表达为

其中

对于k时刻新生多伯努利密度, 第i个伯努利RFS的概率密度${{p}_{\Gamma , k}}^{(i)}(\tilde{x})$表示为

其中, $p_{\Gamma, k}^{(i)}$是对应目标模型的新生概率.

新生伯努利密度的采样过程为: 采样模型粒子${{d}_{\Gamma , k}}^{(i, j)}\tilde{\ }{{p}_{\Gamma , k}}^{(i)}(d)$, 然后选择新生状态密度的建议分布$q_{\Gamma , k}^{(i)(x}|{{d}_{\Gamma , k}}^{(i, j)})$采样新生状态粒子

新生粒子权重

那么, 对应的新生概率密度

其中, $\widetilde{x}_{\Gamma , k}^{(i, j)}={{[x_{\Gamma , k}^{(i, j)\text{T}}, {{d}_{\Gamma , k}}^{(i, j)}]}^{\text{T}}}$.

2) 更新步

假设k时刻联合预测和新生的多伯努利密度为${{\pi }_{k|k-1}}=\left( {{r}_{k|k-1}}^{(i)},p_{k|k-1}^{(i)} \right)$, 根据前述的粒子采样, 则概率密度$p_{k|k-1}^{(i)}$ 可表示为

那么, 更新后的多目标密度可表示为

其中, Z_k是执行传感器控制方案${ u}_k \in {U}_k$ 后, 传感器在相应的控制位置所接收到的量测.

将式(39)代入式(19)和(20), 继承航迹(漏检部分)所对应的多伯努利密度的参数

其中

量测更新航迹多伯努利密度的参数

其中

3) 重采样 为减少粒子退化对滤波器估计性能的影响, 对每一个假设航迹的粒子集进行重采样. 新粒子按照更新后的粒子权重的大小进行重采样, 对应航迹的重采样规模$L_k^{(i)}=L_k\times r_k^{(i)}$, L_k是单目标平均采样个数, $r_k^{(i)}$是更新后的第i个航迹的存在概率. 为保证每个假设航迹的粒子采样个数, 可设置每个航迹的粒子采样规模$L_k^{(i)}$不小于$L_{\min}$.

4.   传感器控制的求解

4.1   目标导向的多目标概率密度的采样方法

如第2节所述, 我们采用$H=1$的“Myopic”控制方案. 那么, 对所有${\nu}\in {U}_k$的控制方案计算相应的评价函数, 然后依据式(6)确定k时刻的最终控制方案. 评价函数的求解必须近似求解相应的多目标概率密度, 可采用多目标SMC方法, 即复杂的集合采样方式 ^[4]. 本文以下将基于多伯努利密度, 提出一种目标导向的多伯努利(TOMB)密度的粒子采样方法, 直接在单目标状态空间上进行采样去逼近多目标概率密度.

式(12)利用多伯努利概率密度去近似多目标概率密度. 当从多目标状态空间上采样时, 由于多伯努利RFS 是相互独立的(单)伯努利RFS 的集合 ^[13], ${r}^{(i)}$表示第i个伯努利RFS (航迹)的存在概率, 因此, 我们将多目标密度粒子采样过程分两步进行. 首先, 根据${r}^{(i)}$的取值采样航迹; 其次, 再对已采样航迹的概率密度${p}^{(i)}$进行采样. 对于第 l 次采样, $1\le l\le {{L}_{k}}\times \sum\limits_{i=1}^{M}{{{r}^{(i)}}}$, 这种采样方式可以表示为

其中, $q(r)$是伯努利RFS (航迹)采样的分布函数. ${q}(\tilde{ x}|r)$是采样航迹条件下状态(含模型参数)采样的分布函数.

这种采样方法将使采样过程集中于${r}^{(i)}$较大的RFS的概率密度, 即以事实上的目标航迹的采样为主导. 那么经过多次采样以后, 所有粒子的加权和将逼近于多目标多伯努利RFS的概率密度$\pi$.

1) 伯努利RFS (航迹)的采样 ${r}^{(i)}$代表第i个伯努利RFS的存在概率, 则设计航迹的采样分布

2) 状态的采样 假设第i个伯努利RFS的概率密度$p^{(i)}$可用加权粒子集近似表达, 那么状态采样的分布

其中, $w^{(i, j)}$是对应状态粒子的权重.

4.2   传感器控制方案的决策

我们利用TOMB的粒子采样方法近似预测的多目标概率密度, 则预测的多目标概率密度

其中, ${{{\hat{L}}}_{k|k-1}}={{L}_{k}}\times \sum\limits_{j=1}^{L_{k-1}^{(i)}}{{{w}_{k-1}}^{(i, j)}}r_{k|k-1}^{(i)}, \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{w}_{k|k-1}^{(j)}=1/{{{\hat{L}}}_{k|k-1}}$.

按照预测多伯努利密度, 提取相应的目标状态集合${{{\hat{X}}}_{k|k-1}}=\{{{{\hat{x}}}_{k|k-1}}^{(i)}\}_{i=1}^{{{{\hat{N}}}_{k|k-1}}}$. 并对每个${\nu}\in {U}_k$确定传感器位置$\tilde{ x}_{s, k}(\nu)$, 然后按式(10) 生成对应控制方案${\nu}$的量测集合.

如第2节所分析, 控制方案的理想量测集合不考虑杂波, 且检测概率$p_D=1$. 在这种理想情况下, 不存在继承航迹(漏检)的伯努利RFS, 而对于量测更新的伯努利RFS, 即$ z \in Z_k({\nu})$, 参考式(21), 此时

这就意味着, 每个理想量测的伯努利RFS (航迹)的存在概率都恒为1. 在这种情况下, 如果根据提出的TOMB 粒子采样方法, 对于各个量测更新RFS (航迹)的采样是等概率的. 参考第3节的采样过程, SMC-MMCBMeMBer 滤波器的更新多伯努利概率密度和预测多伯努利概率密度对应同样的粒子集, 只是粒子权重不同. 所以, 根据后验密度${{p}_{U, k}}(\tilde{x};z)$ 的表达式, 后验多目标密度可近似为

其中, $j^*$代表采样粒子对应的预测伯努利RFS (航迹)的索引.

一步控制方案(Myopic)的评价函数为

利用提出的TOMB粒子采样方法可以近似多目标概率密度, 将式(55)和(57)代入式(64), 多项式相乘后狄拉克$\delta$函数的正交项全为0, 则

将上式的被积函数展开, 再次利用$\delta$函数的正交项为0消除多余项, 评价函数最终可近似为

4.3   总体算法程序的伪码

为了说明提出算法的总体流程, 列出算法伪码如下:

输入. $\{{{r}_{k-1}}^{(i)}, {{p}_{k-}}{{1}^{(i)}}\}_{i=1}^{{{M}_{k-1}}}$, 传感器控制位置$ x_{s, k-1}$, ${{p}_{k-1}}^{(i)}(x)=\sum\limits_{j=1}^{L_{k-1}^{(i)}}{{{w}_{k-1}}^{(i, j)}}{{\delta }_{\widetilde{x}_{k-1}^{(i, j)}}}(x).$

步骤1. 预测和新生

for $i = 1:M_{k-1}$

for $j = 1:L_{k-1}^{(i)}$

依式(27)计算$r_{P, k|k-1}^{(i)}$, 按式(29)和(30)分别采样模型粒子和状态粒子, 按式(32)和(33)计算$\tilde{w}_{P, k|k-1}^{(i, j)}$.

end for

end for

给定新生多伯努利密度为$({{r}_{\Gamma , k}}^{(i)}, p_{\Gamma , k}^{(i)})_{i=1}^{{{M}_{\Gamma , k}}}$, 依据 式(35) $\sim$ (38)采样新生状态粒子, 包括模型粒子, 并计算新生粒子权重.

步骤 2. 传感器控制

按照TOMB采样方法近似估计$p_{k|k-1}(X)$

按预测多伯努利密度提取目标状态集合$\hat{X}_{k|k-1}$.

确定所有可能的控制集合${U}_k$, 并对每个${\nu} \in {U}_k$ 确定传感器

位置$x_{s, k}(\nu)$,

按$\hat{X}_{k|k-1}$生成相应量测集合$Z_k(\nu)$.

按$Z_k(\nu)$, 根据式(57) $\sim$ (63)近似估计$p_{k|k}(X)$.\按式(66)计算${\mathcal{R(\nu) }}$, 并确定最终控制方案${ u}_k$.

传感器在新控制位置$x_{s, k}( {u}_k)$接收实际量测$Z_k({u}_k)$.

步骤3. 更新

$M_{k|k-1}=M_{k-1}+M_{\Gamma , k}$.

for $i = 1:M_{k|k-1}$

计算$\varrho _{L, k}^{(i)}, r_{L, k}^{(i)}, p_{L, k}^{(i)}$.

end for

for $m = 1:{\rm length}(Z_k({u}_k))$

for $i = 1:M_{k|k-1}$

初始化${\varrho _{U,k}}^{(i)}(z)$;

for $j = 1:L_{k|k-1}^{(i)}$

计算${{\psi }_{k, z}}({{\widetilde{x}}_{k}}^{(i, j)})$.${{\varrho }_{U, k}}^{(i)}(z)={{\varrho }_{U, k}}^{(i)}(z)+{{w}_{k|k-1}}^{(i, j)}{{\psi }_{k, z}}({{\widetilde{x}}_{k}}^{(i, j)})$.\计算${{w}^{*}}{{_{U, k}}^{(i, j)}}(z)$.

end for

end for

计算$r_{U, k}({z}), p_{U, k}({ x}; {z})$, 权重归一化.

end for

步骤 4. 航迹删减和重采样

$M_k=M_{k|k-1}+{\rm length}(Z_k( {u}_k))$.

for $ i = 1:M_k$

if ${{r}_{k}}^{(i)}<$航迹删除阈值

删除第i个航迹;

else

按$r_k^{(i)}$求重采样粒子个数, 并进行重采样.

end if

end for

步骤 5. 状态提取

${\rm target\_index} = {\rm find}(r_k^{(i)}>0.5)$.

求目标个数, 即$\hat{N} = {\rm length(target\_inde}x)$.

粒子加权求目标状态集合${{\hat{X}}_{k}}={{\hat{x}}_{k}}^{(i)_{i=1}^{{\hat{N}}}}$.

输出. 目标状态集合$\hat{X}_k$, 多伯努利参数集$\{{{r}_{k}}^{(i)}, {{p}_{k}}^{(i)}\}_{i=1}^{{{M}_{k}}}$, 传感器位置$ x_{s, k}$, 其中

${{p}_{k}}^{(i)}(x)=\sum\limits_{j=1}^{L_{k-1}^{(i)}}{{{w}_{k-1}}^{(i, j)}}{{\delta }_{\widetilde{x}_{k}^{(i, j)}}}(x)$

5.   仿真分析

5.1   多目标跟踪性能评价

为了能联合评价多目标势估计和多目标状态估计的综合性能, 本文采用OSPA (Optimal subpattern assignment)距离 ^[16]来评估多目标跟踪的质量. 先给出OSPA的定义: 设多目标状态的真值集合为$X=\{{x}_1, \cdots, {x}_m\}$, 相应的状态估计集合$\hat X=\{\hat{ x}_1, \cdots, \hat{ x}_n\}$, 若$m\leq n$, 则OSPA距离为

其中, ${\bar d^{(c)}}(x\hat ,x) = \min (c,\left\| {x - \hat x} \right\|)$, ${\prod _k}$表示所有$\{ 1, \cdot \cdot \cdot , k\}$ 的排列构成的集合, 距离阶次$ p \ge 1$, 截断系数 $c > 0$. 如果$m>n$, 则$\bar{d}_{p}^{(c)}(X, \hat{X})=\bar{d}_{p}^{(c)}(\hat{X}, X)$.

5.2   机动多目标跟踪场景的设计

为了跟踪机动目标, 以下采用常速 (Constant velocity, CV)模型和协同转弯 (Coordinated turn, CT)模型 ^[17]的组合多模型对机动目标进行跟踪. CT模型对应可变的转角速度$\omega$, 模型演变的马尔科夫链

CV模型的状态向量参考式(1), 状态转移密度

其中, ${{F}_{k}}=\ \ \left[ \begin{matrix} {{I}_{2}}&T{{I}_{2}} \\ {{\text{0}}_{\text{2}}}&{{I}_{2}} \\ \end{matrix} \right]\ , {{Q}_{k}}={{\sigma }_{v}}^{2}\left[ \begin{array}{*{35}{l}} \frac{{{T}^{4}}}{4}{{I}_{2}}\ \ \frac{{{T}^{\text{3}}}}{\text{2}}{{I}_{2}}\ \\ frac{{T}^{\text{3}}}\text{2}{{I}_{2}}\ \ {{T}^{\text{2}}}{{I}_{2}} \\ \end{array} \right]$; ${I_n}$代表$n\times n$的单位阵, ${\sigma _v}=5 {\rm m/s}^2$.

CT模型的状态${ {{y}}}=[{ x_k}^{\text{T}}, \omega _k]^{\text{T}}$, 它的状态演化模型为

其中, $\tilde { w}_{k}\sim {\rm N}(\cdot;0_{2\times 1}, \sigma _w^2 I_2)$, $u_{k}\sim {\rm N}(\cdot;0, \sigma _u^2 )$, $0_{n\times m}$是${n\times m}$的零矩阵, $\sigma _w=5 {\rm m/s}^2$, $\sigma _u=(5 \pi /180){\rm rad/s}$, 状态转移矩阵

噪声矩阵

仿真区域内均伴随着目标的新生和消亡, 对应着可变的目标数目(势), 所有目标轨迹均采用随机机动模式, 在CV 和CT之间随机切换, 目标运动模式随机切换的概率满足设定的马尔科夫链. 若目标为CT运动, 其对应的随机转角速度$\omega \sim {\rm N} (\omega ; 0, (5\pi/180)^2)$.

5.3   新生多伯努利密度的设计

参考式(34), 初始模型由$p_{\Gamma, k}^{(i)}(d)$决定, 在目标初始的机动类型没有明确先验知识的情况下, 我们认为在模型间的选取是等概率的. 利用量测RFS建立新生多伯努利密度, 假设目标到传感器的新生先验距离为$R_\gamma$, 先验距离的误差是均值为零的高斯噪声, 其标准差为${\sigma _{R_\gamma}}$. 假设角量测$\theta$的噪声是独立于先验距离误差且均值为零的高斯噪声, $\sigma _\theta$是其标准差. 基于修正的量测转换方法(Modified unbiased converted measurement, MUCM)^[18]将量测转换到直角坐标系, 表示补偿因子$\lambda _\theta={{\rm e}^{ -{\sigma _\theta }^2/2}}$, 有转换量测

对应的转换后的协方差

其中

假设k时刻的多目标量测RFS为${{Z}_{k}}( {u}_k)$, 势为$n_z(k)$, 对于一个新接收的量测$z_{k}^{(i)}\in {{Z}_{k}}({{u}_{k}}), i=1, \cdots , {{n}_{z}}(k)$, 由它确定一个独立的伯努利RFS, 其密度为$\pi _{b, k}^{(i)}=\{({{r}_{b, k}}^{(i)}, {{p}_{b, k}}^{(i)})\}, \pi _{b, k}^{(i)}$ 的预测密度作为$k+1$ 时刻的新生伯努利密度. 其中$r_{b, k}^{(i)}=0.01$, 密度

其中, ${w_{b, k}}^{(i, j)}=1/L_{b, k}^{(i)}$, j 是采样粒子索引, ${L_{b, k}}^{(i)}=L_{\rm min}$. 初始状态粒子$\tilde{ {x}}_{b, k}^{(i, j)}$由位置、速度和模型粒子构成, 模型粒子在两个模型间均匀采样 , 而位置由分布${\rm N}( x; { m}_{\rm MUCM}( z_k^{(i)}), P_{\rm MUCM}( z_k^{(i)}))$采样. 新生速度分量根据分布${\rm N}(\cdot; 0_{2\times1}, {v^2_{\rm max}}\times I_2)$采样. $v_{\rm max}$是目标的最大速度, 设为50 m/s.

5.4   传感器控制集合

传感器控制位置理论上是连续且存在无限多个, 但这显然是没有必要的. 我们采用启发式方法, 在传感器运动容许的范围内选择典型的控制方案. 如果k时刻传感器的实际控制位置为${ x}_{s, k}=[{{x}_{s, k}}, {{y}_{s, k}}]$, 那么$k+1$时刻传感器控制集合为${U}_{k+1}$, 则选择所有可能的传感器位置为

式中, 选择${N_{\theta}}=8$, ${N_R}=2$, 则总共对应17种控制方案(包含传感器不动). $v_{s, c}$是传感器自身的容许控制速度, 设为50 m/s.

5.5   机动多目标距离方位跟踪

考虑机动多目标距离方位跟踪(Range-bearing tracking, RBT), 模型轨迹如图 1所示, 传感器初始位置位于坐标原点, 多目标的新生和消亡的时刻如表 1所示.

图 1  实际目标的轨迹

Fig. 1  Actual target trajectories

下载: 全尺寸图片 幻灯片

表 1  目标存活周期

Table 1  Survival periods of targets

	目标1	目标2	目标3	目标4	目标5
新生时刻(s)	1	1	10	20	20
消亡时刻(s)	50	30	50	50	40

下载: 导出CSV 

| 显示表格

对于距离方位跟踪(RBT), 式(4)有如下形式:

量测噪声${ v_k}\sim {\rm N} (\cdot, 0, R_k)$, $R_k={\rm diag}\{[\sigma _r^2, \sigma _\theta^2]^{\text{T}}\}$, $\sigma _\theta=(\pi/180) {\rm rad}$, $\sigma _r=5$ m. 利用MUCM 建立新生伯努利密度, $R_\gamma$ 和${\sigma _{R_\gamma}}$ 可分别用距离量测和$\sigma _r$代替.

量测采样周期$T=1$ s, 总共采样50次. 杂波是一个泊松RFS, 且在观测区域$[0 , 2\pi]\times [0, 5 000 {\rm m}]$内均匀分布, 每周期的杂波平均数为10个. 目标存活概率为${{p}_{S, k}}=0.98$, 检测概率${{p}_{D, k}}=0.95$. 选取截断系数\emph{c} = 50 m, 距离阶次\emph{p} = 1. 单航迹平均采样数$L_k=500$, 单航迹最小采样数$L_{\rm min}=100$. 航迹删除的存在概率阈值为$10^{-4}$. 仿真硬软件环境为: Matlab R2010a, Win 7 SP1 64-bit, Intel Core i5-4570 CPU 3.20 GHz, RAM 4.00 GB.

利用第3节的SMC-MMCBMeMBer滤波器跟踪场景中的机动目标. 给定4种传感器控制方案来对比验证本文控制算法对机动多目标的跟踪效果. 其中, 方案一是"Stationary", 代表静默方案, 即传感器静止在原点. 方案二是"Prior zigzag", 是预先设定好的传感器控制方案, 如图 2所示, 传感器在实验内以恒速率运动并经历多个航向的变化产生机动, 这是一个纯方位跟踪(Bearings-only tracking, BOT)中常用的"Zigzag"控制轨迹, 为可靠保证多目标BOT的可观测性. 方案三是"Random control", 代表每个时刻的传感器控制方案在控制集合中随机选取. 方案四是"Proposed control", 即文中提出的传感器控制方案.

图 2  先验控制方案

Fig. 2  Prior sensor control

下载: 全尺寸图片 幻灯片

图 3是单次实验对机动多目标RBT跟踪中的传感器轨迹控制, 看出轨迹呈复杂的"Zigzag"形状, 在整体多目标跟踪的过程中, 传感器会根据当前的滤波结果, 根据提出的控制方案不断调整自身的位置以获取最大的观测信息量, 保证传感器形成和所有目标全局最佳的观测位置, 每个周期控制的目的都使得多目标后验概率密度相对于多目标先验概率密度的信息增量最大化, 并伴随着目标机动、目标的新生和消亡都会呈现出显著的控制作用.

图 3  距离方位跟踪中的传感器控制轨迹

Fig. 3  Sensor control trajectory in range-bearing tracking

下载: 全尺寸图片 幻灯片

为了验证提出的机动多目标传感器控制策略在RBT中整体的有效性, 做200次Monte Carlo (MC)仿真实验. 所有控制方案的多目标位置估计的OSPA距离均值如图 4所示, 目标个数(势)估计效果如图 5所示.

图 4  距离方位跟踪中多目标位置估计OSPA 的比较

Fig. 4  Position OSPA distances for di®erent control strategies in range-bearing tracking

下载: 全尺寸图片 幻灯片

图 5  距离方位跟踪中目标个数估计的比较

Fig. 5  Cardinality statistics for di®erent control strategies in range-bearing tracking

下载: 全尺寸图片 幻灯片

通过图 4多目标跟踪的OSPA距离评价, 可以分析不同传感器控制策略下机动多目标RBT的整体性能, 虽然RBT具有完备的量测信息使得4种控制方案都有不错的跟踪性能, 但本文所提出的控制策略相对于其他控制方案在多目标整体跟踪性能上具有明显的优势, 显著提升了多目标跟踪的质量. 图 5是多目标势估计效果, 各控制方案的势估计均值都接近于实际的目标势(图 5 (a)), 但从表 2 仍然可以清晰看出本文的算法具有相对较小的势估计误差, 而且从图 5 (b)也可以看出, 本文算法具有相对稳定的势估计. 注意到, 方案一(Stationary) 在这个实验中的跟踪效果整体要好于方案二(Prior zigzag), 这也说明RBT的控制要求不同于传统的BOT, BOT更多的是需要考虑可观测因素, 所以受传感器机动能力的影响很大. 而RBT具有完备的测量信息, 为了获取更好的跟踪效果, 在传感器的实际控制上需要参考不同传感器位置的信息增量. 本文正是结合给出的SMC-MMCBMeMBer滤波器, 对机动多目标的多伯努利后验概率密度进行估计, 以此为基础, 依据提出的TOMB采样方法和对评价函数的具体求解, 基于信息增益的最大化准则, 总能为不同的机动多目标跟踪问题给出一个适合于当前条件下, 使多目标后验概率密度信息增量最大化的控制策略. 而因为这种基于信息测度的评价完全基于机动多目标状态(含模型参数)的后验概率密度, 由它构造多目标跟踪的全局代价函数, 使得控制策略整体考虑了影响到多目标跟踪质量的各种因素.

表 2  单步势估计误差均值的绝对值

Table 2  Absolute values of step-averaged cardinality error

	方案一	方案二	方案三	方案四
势误差Ne	0.1319	0.1348	0.1342	0.1258

下载: 导出CSV 

| 显示表格

图 6显示了机动多目标跟踪的MC仿真中所有实验的传感器控制位置(轨迹云). 由于杂波个数、分布以及噪声的随机性, 显然每次独立实验的传感器轨迹都不大可能一致, 但是该轨迹云仍能够充分呈现出该场景下实现机动多目标RBT跟踪的最优轨迹控制的大致趋势.

图 6  距离方位跟踪蒙特卡罗仿真中所有传感器控制位置

Fig. 6  Total position for sensor control in range-bearing tracking

下载: 全尺寸图片 幻灯片

为了评价各个控制算法的计算效率, 我们统计了MC仿真中程序单步运行的平均时间, 如表 3所示. 从结构上看, 整体控制算法分为多目标跟踪部分和传感器控制部分. 由于4种控制方案在多目标跟踪部分都利用本文给出的SMC-MMCBMeMBer滤波器跟踪机动多目标, 所以运行时间的实际差距在于传感器控制的计算. 方案三仅需要在控制集合中随机选取, 所以前三种方案的程序运行时间基本都由机动多目标跟踪本身所花费的时间所主导. 而方案四相对其他控制方案所多出的计算时间则是计算后验和先验多目标概率密度的Bhattacharyya距离和最终的传感器控制求解所花费的时间, 以换取多目标跟踪性能的显著提升.

表 3  各方案单步平均运行时间

Table 3  Step-averaged run time for di®erent control strategies

	方案一	方案二	方案三	方案四
时间(s)	1.2894	1.2988	1.3025	4.5737

下载: 导出CSV 

| 显示表格

6.   结论与展望

本文主要的工作和创新是提出了一种机动多目标传感器控制策略, 该方法在POMDP理论框架下进行. 文中基于信息论方法利用Bhattacharyya距离作为传感器控制的评价指标, 结合所研究的机动多目标跟踪问题给出了SMC-MMCBMeMBer滤波器的一般实现形式. 此外, 为求解最终的控制方案, 提出了一种TOMB粒子采样方式去近似多目标概率密度, 并基于该方法对最终控制方案进行求解. 通过最终典型非线性多目标跟踪问题的仿真实验验证了提出控制策略的有效性. 对于今后的工作, 可以进一步开展目标跟踪中多传感器控制问题的研究, 包括多传感器资源的合理分配和优化协同调度, 并研究多传感器控制的最优准则. 另外, 本文虽然以RBT 为例, 但本文研究方法仍然适用于其他跟踪问题, 尤其是一些对传感器轨迹控制有迫切需求的典型跟踪问题, 例如纯角度跟踪(Angle-only tracking, AOT)问题、纯距离跟踪问题(Range-only tracking, ROT)等.