0 引　言

概率假设密度（PHD）滤波是一种基于随机有限集理论的多目标贝叶斯滤波，对多目标跟踪过程中的目标随机出现、分裂、消失、漏检和虚警杂波等现象有着严格的数学描述，它代表着多目标跟踪技术的一个新方向^[1]。通常，实现递推传递概率假设密度函数的多目标跟踪算法比较困难，Vo等利用粒子滤波技术给出了PHD滤波器的近似实现形式^[2]，简称为PF-PHD滤波算法，这就加快了PHD滤波的工程应用。

文献[3]中将势分布概率假设密度滤波（CPHD）应用于水下多目标跟踪，并针对传感器量测的不确定性，在CPHD滤波的框架下引入量测的幅值信息，增加数据关联的可靠性。文献[4]探究了基于多波束前视声呐数据的2种不同多目标跟踪方法，一种方法是给每个目标分配卡尔曼滤波器，通过设定量测门限和数据关联技术进行目标跟踪，另一种是应用PF-PHD滤波算法结合数据关联技术的目标状态估计。仿真表明PHD滤波器比传统卡尔曼滤波方法更优越。文献[5]利用声信号得到目标的距离测量信息，采用蒙特卡罗方法实现的PHD滤波器成功地跟踪3个水下自主航行器（AUV）。文献[6]通过单一的水下传感器获得目标位置量测信息，针对杂波存在时，通常的跟踪算法仅依赖当前的位置观测，故使用高斯混合概率假设密度（PHD）前向-后向平滑算法，它是利用多个量测数据对滤波值进行平滑，进而减小目标估计的位置误差。结果表明在高密度杂波的水下环境中该算法表现良好。

目前大部分文献侧重研究目标的主动跟踪，实际上，被动跟踪在某些方面更具应用价值。本文采用PF-PHD滤波算法对监测区域内目标进行跟踪，并利用双被动声呐^[7]来获得目标的方位角信息。

1 双被动声呐系统观测模型

假设二维监测区域内有多个目标,采用匀速运动（CV）模型模拟每个目标的位置变化情况。采样时间为 $T$ ，其中一个目标 $k$ 时刻的状态矢量表示为 ${x_k} = {[{\rm{x}},{v_x},{\rm{y}},{v_y}]^{\rm T}}$ ， ${x_k}$ 包括目标分别在 $x,y$ 方向上位置和速度。那么，任意目标的状态方程可表示为：

${x_k} = { F}{x_{k - 1}} + { G}{w_k}\text{，}$

(1)

其中， ${ F} = \left[ {\left. {\begin{array}{*{20}{c}} 1&T&0&0 \\ 0&1&0&0 \\ 0&0&1&T \\ 0&0&0&1 \end{array}} \right]} \right.$ ， ${ G} = \left[ {\left. {\begin{array}{*{20}{c}} {{T^2}/2}&0 \\ T&0 \\ 0&{{T^2}/2} \\ 0&T \end{array}} \right]} \right.\text{，}$

${ F}$ 表示目标状态转移矩阵， ${ G}$ 为状态噪声强度的输入矩阵。通常假设过程噪声 ${w_k}$ 服从均值为0、协方差矩阵为 ${ Q}$ 的高斯分布，并且假设过程噪声和量测噪声相互独立。

假设观测平台上装备有仅得到目标角度观测信息的2部被动声呐基阵，等效于雷达跟踪中的2个测向传感器，那么这里传感器个数为2，接下来文中描述的单个传感器均指单被动声呐基阵，并且假设被跟踪目标与声呐设备在一个平面上，以便研究算法特性。令 $({a_s},{b_s})$ 为第 $s$ 个传感器的位置。则对于第 $s$ 个被动传感器，可以得到目标 ${x_k}$ 的方位角：

${\beta _{k,s}} = {\rm arc}\tan \frac{{{\rm{x}} - {a_s}}}{{{\rm{y}} - {b_s}}}\text{，}$

(2)

则对于第 $s$ 个传感器量测方程：

${z_{k,s}} = {h_s}({x_k}) + v_k^s = {\beta _{k,s}} + v_k^s\text{，}$

(3)

式中： $v_k^s$ 为第 $s$ 个传感器的量测噪声矢量，本文采用同类型传感器，假设各传感器量测噪声相互独立，且服从同一高斯分布： $v_k^s \sim N(0,\sigma _v^2)$ ， $s = 1,{\rm{2}}$ ，其中， $\sigma _v^2$ 为量测角度方差。那么，所有传感器对目标 ${x_k}$ 的集中式量测可表示为 ${z_k} = {[{z_{k,1}},{z_{k,{\rm{2}}}}]^{\rm T}}$ 。

2 概率假设密度（PHD）滤波 2.1 随机集（RFS）的多目标模型

多目标跟踪中当每一时刻目标数、杂波或虚警个数不同时，状态空间和量测空间的维数也会随之变化，此时目标跟踪的模型应表示为各个目标状态和量测的集合。基于随机集理论的跟踪模型有着坚实的数学基础，能准确地描述多目标跟踪中目标产生、消失、衍生等现象。

假设在采样时刻 $k$ ，多目标状态集为 ${X_k}\!\!=\!\!\{ {x_{k,1}},\cdots\!,{x_{k,{N_k}}}\} $ ， ${N_k}$ 表示在 $k$ 时刻的目标数， ${x_{i,k}}$ $(i = 1,\cdots,{N_k})$ 表示目标 $i$ 的状态向量。多目标量测集为 ${Z_k} = \{ z_k^1,\cdots,z_k^{{M_k}}\} $ ，其中 ${M_k}$ 表示 $k$ 时刻得到的量测数目， $z_k^j$ $(j = 1,\cdots,{M_k})$ 表示第 $j$ 个量测向量，部分量测可能来源于杂波或虚警。那么， $k$ 时刻的目标状态集 ${X_k}$ 可表示为：

$ {X_k} \!\!= [\mathop \cup \limits_{{x_{k - 1}} \in {X_{k - 1}}} \!\!\!\!\!\!\!{S_{k\left| {k - 1} \right.}}({x_{k - 1}})] \cup [\mathop \cup \limits_{{x_{k - 1}} \in {X_{k - 1}}} \!\!\!\!\!\!\!\!{B_{k\left| {k - 1} \right.}}({x_{k - 1}})] \cup {\varGamma _k}\text{。}\!\!\!\!\! $

(4)

其中： ${X_{k - 1}}$ 为 $k - 1$ 时刻的目标状态集； ${S_{k\left| {k - 1} \right.}}$ 为 $k - 1$ 时刻到 $k$ 时刻的仍存活的目标状态集； ${B_{k\left| {k - 1} \right.}}$ 为从 $k - 1$ 时刻所衍生的目标状态集； ${\Gamma _k}$ 为在 $k$ 时刻新生目标随机集。类似的，在 $k$ 时刻多目标的量测有限集可描述为：

${Z_k} = {K_k} \cup [\mathop \cup \limits_{{x_{k - 1}} \in {X_{k - 1}}} {\Theta _k}({x_{k - 1}})]\text{。}$

(5)

其中： ${\Theta _k}({x_{k - 1}})$ 为源于真实目标的量测集合； ${K_k}$ 为由于杂波或虚警引起的量测集合。

2.2 PHD滤波原理

针对多目标跟踪中Bayes公式难求解的问题，Mahler提出了概率假设密度（PHD）滤波近似实现多目标Bayes滤波器，PHD滤波通过递推多目标后验密度的一阶统计矩来降低计算复杂度。

在 $k > 0$ 时， ${D_k}$ 表示联合多目标后验分布的强度函数。PHD滤波通过预测和更新操作对强度函数 ${D_k}$ 进行传播。

PHD滤波预测方程：

$\begin{split} & {D_{k\left| {k - 1} \right.}}(x\left| {{Z_{k - 1}}} \right.) = \\ & \int {{\phi _{k\left| {k - 1} \right.}}} (x,{x_{k - 1}}){D_{k - 1\left| {k - 1} \right.}}({x_{k - 1}}\left| {{Z_{k - 1}}} \right.){\rm d}{x_{k - 1}} + {\gamma _k}(x)\end{split}\text{，}$

(6)

${\phi _{k\left| {k - 1} \right.}}(x,\xi ) = {b_{k\left| {k - 1} \right.}}(x\left| \xi \right.) + {e_{k\left| {k - 1} \right.}}(\xi ){f_{k\left| {k - 1} \right.}}(x\left| \xi \right.)\text{。}$

(7)

其中： ${\gamma _k}(x)$ 和 ${b_{k\left| {k - 1} \right.}}( \cdot \left| \xi \right.)$ 分别为新生目标随机集 ${\Gamma _k}$ 的PHD和前一时刻状态为 $\xi $ 的目标衍生随机集 ${B_{k\left| {k - 1} \right.}}$ 的PHD； ${e_{k\left| {k - 1} \right.}}(\xi )$ 为 $k - 1$ 刻状态为 $\xi $ 的目标存活概率； ${f_{k\left| {k - 1} \right.}}( \cdot \left| \cdot \right.)$ 为单个目标的状态转移概率密度。

PHD滤波更新方程：

$\begin{split} & {D_k}(x\left| {{Z_k}} \right.) = \Bigg[1 - {P_{D,k}} +\Bigg.\\ & \Bigg.\sum\limits_{z \in {Z_k}} {\frac{{{\psi _{k,z}}(x)}}{{{K_k}(z) + \langle {\psi _{k,z}},{D_{k\left| {k - 1} \right.}}\rangle }}}\Bigg] \cdot {D_{k\left| {k - 1} \right.}}(x\left| {{Z_{k - 1}}} \right.)\end{split}\text{。}$

(8)

其中： ${\psi _{k,z}}(x) = {P_{D,k}}{g_k}(z\left| x \right.)$ ， ${P_{D,k}}$ 为 $k$ 时刻的检测概率， ${g_k}( \cdot \left| \cdot \right.)$ 为单个目标的观测似然函数； ${K_k}( \cdot )$ 为 $k$ 时刻杂波随机集的PHD，有 ${K_k}(z) = {\lambda _k}{c_k}(z)$ ， ${\lambda _k}$ 为观测范围内杂波的平均数； ${c_k}(z)$ 为杂波概率密度。

对式（8）中的PHD函数 ${D_k}$ 的积分，即为 ${N_k} = \int {{D_k}} (x\left| {{Z_k}} \right.){\rm d}x$ ，然后对 ${N_k}$ 取整得到目标数的估计值 $\mathop {{N_k}}\limits^ \wedge = {[{N_k}]_{\operatorname{int} }}$ 。进一步由 ${D_k}(x\left| {{Z_k}} \right.)$ 的 $\mathop {{N_k}}\limits^ \wedge $ 个峰值点所在位置，得到 $k$ 时刻 $\mathop {{N_k}}\limits^ \wedge $ 个目标的状态向量。

2.3 多目标跟踪的PF-PHD滤波

PF-PHD滤波是利用一系列离散的带权重的样本近似相应的PHD函数，通过聚类算法提取目标的状态估计，当粒子数目达到一定程度，PF-PHD滤波逼近Bayes最优估计，文献[2]详细介绍了粒子滤波实现PHD滤波的原理及具体步骤。

传统的PHD滤波器的应用前提是单传感器，对于多传感器的情况，也有相应的算法提出，比如一种利用所有传感器的量测更新PHD的迭代更新近似算法^[8]，以及Mahler提出的一种乘积形式的多传感器PHD滤波算法^[9]，这些算法的共同缺点是计算复杂度大。本文根据文献[10]，采用集中式融合策略，在第k时刻,将各个传感器的观测统一送至数据融合中心，假定量测数据已完成配准及关联，这样就可以将多传感器问题转化为单传感器问题，进一步简化多目标跟踪算法。

3 仿真实验

考虑二维监测区域范围为[–3000 m，3000 m] × [–3000 m，3000 m]，进入监测区域的目标都作匀速直线运动且航速约40 kn，即20 ${\rm m}/{\rm s}$ ，利用双传感器进行纯方位角跟踪。仿真中采样间隔 $T = 5{\rm s}$ ，总步长为60，目标检测概率 ${p_D} = 1$ ，目标存活概率 $e = 0.99$ ，杂波数服从泊松分布且杂波在跟踪场景内在 $[ - {\text{π}} ,{\text{π}} ]$ 上服从均匀分布，新生目标强度函数的高斯分量权值均为0.1，不考虑目标的衍生。采用OSPA距离作为多目标跟踪的评价指标，侧重评价目标跟踪的距离误差，因此OSPA距离参数设为： ${{c}} = 1$ ， $p = 150$ 。初始粒子数为500，每个实验均进行50次Monte Carlo仿真。

1）实验1　单目标跟踪

假设监测区域内只有一个目标，沿横坐标 $x$ 方向航行，航速为40 kn。为研究目标和两传感器连线初始链距（简称目标初始链距）、双站距离对算法跟踪性能的影响，将目标初始链距记为 $D$ ，两传感器间距记为 $S$ ，单位为 ${\rm m}$ ，令 $S,D = 500:500:6\;000$ ，即初始距离为 $500\;{\rm m}$ 且以 $500\;{\rm m}$ 的步长变化。

仿真中单个目标在 $k = 1$ 时刻出现， $k = 60$ 时刻消失，位置单位为 ${\rm m}$ ，速度单位为 ${\rm m}/{\rm s}$ ，假设目标初始状态为 ${x_{1,1}} = {[ - 3\;000,20,(D - 3\;000),0]^{\rm T}}$ ，两传感器位置坐标分别为 $\left( { - 3\;000, - 3\;000} \right)$ ， $(S - 3\;000, - 3\;000)$ 。目标状态噪声方差矩阵 ${ Q} = {\rm diag}[0.01,0.01]$ ，测量角度方差 $\sigma _v^2{\rm{ = (1}} \times \pi /180{)^2}$ ，不考虑杂波的干扰。新生目标强度函数为 ${\gamma _k}(x\left| r \right.) = 0.1 \times N(x;{x_{1,1}},{p_\gamma })$ ，协方差矩阵 ${{ p}_\gamma } = {\rm diag}[100,4,100,4]$ 。

二维搜索得到的结果是当两传感器间距和目标初始链距以 $500\;{\rm m}$ 步长同时从 $500\;{\rm m}$ 变化到 $6\;000\;{\rm m}$ 时，单目标位置估计的均方根误差变化情况。从得到数据初步判断，在某一传感器间距大小一定时，随目标初始链距增大，则目标估计的均方根误差总体呈上升趋势，算法的估计性能下降；同理，固定某一初始链距大小，目标估计的均方根误差随着传感器间距增大而总体呈下降趋势。这一直观结果与文献[7]吻合。假设单目标位置估计的均方根误差在 $10\;{\rm m}$ 内是可以接受的，那么，当 $D \in ({\rm{500}} \sim 3\;500)$ 且 $S \in (3\;000 \sim 6\;000)$ 时，算法的性能更优越，这就表明二维搜索能够有效确定参数的最佳取值范围。

进一步分析数据标准差，判断目标估计均方根误差对参数变化的敏感程度。固定每一个传感器间距 $S$ 的值，对变量目标初始链距 $D$ 求其标准差，同样固定每一个目标初始链距 $D$ 的值，计算变量传感器间距 $S$ 的标准差，数据标准差变化显示在图1。根据图1，当传感器间距变化时，目标初始链距对应的均方根误差的标准差一直在 $10\;{\rm m}$ 以上，而且不断增大，相比之下，当目标初始链距变化，传感器间距对应的均方误差的标准差有超过 $80\% $ 的值在 $10\;{\rm m}$ 以下。这就表明目标估计的均方根误差对于目标初始链距的改变更敏感，因此，利用2个传感器进行目标跟踪时，可以首先确定出合适的传感器间距。

图2表示目标估计均方根误差均值随传感器间距变化情况，其中均值指的是各初始链距对应均方根误差的均值。将实验2和实验3的传感器间距设置为 $3\;000\;{\rm m}$ ，此时图2中曲线趋向平稳。

2）实验2　多目标编队航行的场景

假设目标编队航行，各目标的间距为 $DIS$ ，单位为 ${\rm m}$ ，目标速度设为40 kn，沿横坐标 $x$ 方向航行，为保证在正确估计目标数情况下研究目标间距对算法的影响，故将目标数目设置成随时间变化。

假设监测区域内存在3个目标，位置单位为 ${\rm m}$ ，速度单位为 ${\rm m}/{\rm s}$ 。目标1初始状态为 ${x_{1,1}} = {[ - {\rm{300}}0,20,{\rm{0}},0]^{\rm T}}$ ，目标2初始状态为 ${x_{1,2}} = {[ - 3\;000,20,DIS,0]^{\rm T}}$ ；目标3初始状态为 ${x_{1,3}} = {[ - 3\;000,20, - DIS,0]^{\rm T}}$ ，取各目标间距 $DIS = 2\;500,1500,500,250,125,62.5$ ，分别进行实验仿真。其中各目标出现和消失的时刻与实验一相同。2个传感器位置分别为 $( - 3\;000, - {\rm{300}}0)$ ， $(0, - 3\;000)$ 。状态噪声方差矩阵 ${ Q} = {\rm diag}[0.{\rm{01}},0.{\rm{01}}]$ ，测量角度方差为 $\sigma _v^2{\rm{ = (1}} \times \pi /180{)^2}$ 。杂波数服从 $r = 5$ 的泊松分布，则虚警PHD为 ${\kappa _k}(y) = r/{(2\pi )^2}$ 。新生目标强度函数为： ${\gamma _k}(x\left| r \right.) \!=\! 0.1 \!\times \![N(x;{x_{1,1}},{p_\gamma }) \!+\! N(x;{x_{1,2}},{p_\gamma }) \!+\! N(x;{x_{1,3}},{P_\gamma })]$ ，协方差矩阵 ${p_\gamma }$ 与实验1相同。

图3描述了不同 $DIS$ 取值时，利用本文算法对3个编队航行的目标进行跟踪，评价指标OSPA距离的变化情况。当 $DIS = 2\;500,1\;500,500,250$ 时，整个跟踪时段内，OSPA距离始终维持在 $40\;{\rm m}$ 以下，没有明显的波动。对比之下，当 $DIS = 125,62.5$ 时，OSPA距离在不同时刻出现突变，而且与目标间距的取值大小相当，这是由于目标之间的相关性增强导致算法无法准确分辨哪一个估计值属于哪一个目标。因此需要合理选取目标间距，避免算法失效。

3）实验3　航迹交叉的场景

假设监测区域内存在3个目标，位置单位为 ${\rm m}$ ，速度单位为 ${\rm m}/{\rm s}$ 。目标1初始状态 ${x_{1,1}} = {[ - {\rm{250}}0,18,100{\rm{0}},0]^{\rm T}}$ ，目标2初始状态 ${x_{1,2}} = {[2\;000, - 15,2\;000, - 12]^{\rm T}}$ ，目标3初始状态为 ${x_{1,3}} = {[ - 2\;000,18,2\;000, - 14]^{\rm T}}$ ，其中各目标出现和消失的时刻与实验1相同。两传感器位置、虚警PHD和新生目标强度函数与实验2相同。考虑实际情况，假设目标状态噪声方差矩阵 ${ Q} = {\rm diag}[0.{\rm{04}},0.{\rm{04}}]$ ，测量角度误差从1°变化到3°。

由图4~图6知，当目标航迹发生交叉时，量测角度误差保持在2°以下时，采用PF-PHD算法能够对监测区域内每时刻存在的目标进行跟踪，且目标跟踪精度很高。对于角度量测误差超过一定范围时，在航迹交叉点和交叉区域内目标跟踪偏离程度高，效果很差，因此算法对于声呐基阵的测向精度有一定要求。

4 结　语

本文研究了基于概率假设密度滤波的水下多目标跟踪技术，针对算法中两声呐距离和目标初始链距2个重要参数取值的问题，提出的二维搜索法能够有效指导实际目标跟踪中算法的参数选取，并且通过数据分析得出单目标跟踪的精度对于目标初始链距的变化比两声呐间距变化更敏感这一重要结论。同时，对于多目标编队航行和航迹交叉情况，仿真结果表明算法的高效性和稳定性。