CN105182291A - 自适应目标新生强度的phd平滑器的多目标跟踪方法 - Google Patents

Abstract

本发明公开了一种自适应目标新生强度的PHD平滑器的多目标跟踪方法，主要解决在杂波环境下，自适应目标新生强度的PHD滤波器在目标新生时刻存在目标确认的滞后现象，并给出其线性高斯条件下的实现形式。通过前向滤波和后向平滑，能够准确估计目标数目和状态，减小目标新生时确认滞后对航迹形成带来的影响。其步骤包括目标新生率估计、前向滤波和后向平滑，即首先根据先验杂波数均值来估计k时刻的新生率；其次，采用k时刻的量测对目标进行预测和更新来完成前向滤波；然后，用滞后的L时刻量测对滤波结果进行后向平滑；最后，通过修剪合并和状态提取完成跟踪结果的输出。

Description

自适应目标新生强度的PHD平滑器的多目标跟踪方法

技术领域

本发明属于雷达信号处理技术领域，涉及多目标跟踪。具体地说是一种自适应目标新生强度的概率假设密度平滑器(λ-ATBI-PHDSmoother)的多目标跟踪方法，可用于杂波环境下的火控、航管等探测***中。

背景技术

无论现代的防御探测还是空中交通管制ATC(AirTrafficControl)***中，多目标跟踪MTT(MutipleTargetTracking)都是其中的关键技术，也是历来受到最多关注的方向之一。多目标环境下的跟踪问题，存在以下难点：(1)每一时刻都可能存在目标的出现、衍生以及消失，使得目标数目处在一个不断变化的过程中；(2)量测信息不确定，若对漏检、虚警等问题处理不当，将会极大地影响跟踪效果；(3)进行跟踪滤波的同时，还需要完成数据关联。

因此，在复杂环境下，存在的虚警、漏警等量测信息的不确定性和目标生灭带来的数目变化给多目标的稳定跟踪造成很大困难，一直是该领域的研究热点和难点。

传统的多目标跟踪需要进行数据关联，在目标和量测数目较多时计算量会急剧增大。近年来，Mahler等学者提出基于随机有限集(RandomFiniteSet,RFS)理论的多目标跟踪方法，并推导了PHD(ProbabilityHypothesisDensity，概率假设密度)、CPHD(CardinalisedProbabilityHypothesisDensity，势概率假设密度)、多伯努利等实现形式。该理论框架下，多目标的量测和状态分别被视为一个随机有限集，避免了复杂的数据关联，因此受到广泛的关注。

PHD滤波是通过递推目标状态的后验概率密度的一阶矩，可从中来估计多目标的状态和个数。PHD的实现方式有两种：高斯混合(GaussianMixture,GM)和粒子。在线性高斯条件下，高斯混合PHD(GM-PHD)利用加权的高斯混合分布来拟合目标的后验概率密度函数，状态提取不再需要进行聚类，可以用较小的计算量完成滤波估计。

在传统的PHD和CPHD滤波中，认为新生目标多出现在感兴趣的特定位置(例如：机场，视场的边界处等)，新生目标的初始强度是作为先验已知量。而实际应用中，新生目标的初始强度是很难获得的，这给工程应用造成了困难。

B.Ristic等人于2012年提出了基于量测驱动的自适应目标新生强度概率假设密度滤波器(AdaptiveTargetBirthIntensityProbabilityHypothesisDensityFilter，ATBI-PHDFilter)，即ATBI-PHD滤波器，能够从量测中估计目标新生强度并同时对新生目标和存活目标进行递推滤波，摆脱对新生强度先验的依赖，并将ATBI方法和多个固定点检测新生的方法进行了对比。LiangMa等人于2014年针对杂波环境下新生目标的快速生成，提出对PHD滤波器的目标新生率的估计方法，得到λ-ATBI-PHDFilter，即λ-ATBI-PHD滤波器。

但该方法在杂波环境下，自适应目标新生强度的PHD滤波器在目标新生时刻存在目标确认的滞后现象，即对新生目标出现时刻的估计存在延迟，对后续完整的目标航迹形成和关联造成影响。

发明内容

本发明为解决上述现有方法的问题，提出一种自适应目标新生强度的PHD平滑器的多目标跟踪方法，即λ-ATBI-PHD平滑器的多目标跟踪方法。本发明的主要方法：首先，根据先验杂波数均值来估计目标的新生率，能够在得到特定杂波数均值的情况下，尽快检测到新生目标；其次，通过后向平滑能够更加准确的估计目标出现时刻，减小目标新生时刻的确认滞后情况对后续航迹关联的影响，使得目标数目的估计更加准确。

本发明实现上述目的的技术方法包括如下步骤：

1)对k-1时刻跟踪滤波器中存活目标PHD的高斯混合形式 $D_{k-1/k-1}(y,0)=\underset{i}{\Σ}{w_{s,k-1}}^{\left(i\right)}N({m_{s,k-1}}^{\left(i\right)};{P_{s,k-1/k-1}}^{\left(i\right)})$ 进行初始化，设定第i个高斯项均值m_s,k-1 ⁽ⁱ⁾，高斯项滤波协方差P_s,k-1/k-1 ⁽ⁱ⁾和高斯项权值w_s,k-1 ⁽ⁱ⁾的初始值，其中i为高斯项的标号，D_k-1/k-1(y,β)表示全体目标PHD，D_k-1/k-1(y,0)表示存活目标PHD，D_k-1/k-1(y,1)表示新生目标PHD，y为目标状态，N(m_s,k-1 ⁽ⁱ⁾；P_s,k-1/k-1 ⁽ⁱ⁾)表示均值为m_s,k-1 ⁽ⁱ⁾，方差为P_s,k-1/k-1 ⁽ⁱ⁾的高斯项；

2)根据先验杂波信息对新生目标强度估计：

2a)根据k时刻传感器的测量方差Σ_k/k-1、杂波检测概率p_D ^(c)、目标检测概率p_D ^(t)、目标存活概率p_S ^(t)、先验杂波数均值N^(c)、目标新生权值门限T_e和监控区域体积V_S信息来估计k时刻的目标新生率λ_b,k；

2b)根据估计的目标新生率λ_b,k和量测数据z_k，把新生目标的PHD中的高斯项加入到k-1时刻跟踪结果的全体目标的高斯项中，完成对监测区域内新生目标强度的检测，式中w_b,k ⁽ⁱ⁾＝λ_b,k，w_b,k ⁽ⁱ⁾，m_b,k ⁽ⁱ⁾，P_b,k ⁽ⁱ⁾分别对应新生目标高斯项的权值，均值，方差；

3)前向滤波，包括PHD预测和PHD更新：

3a)PHD预测：对k-1时刻的全体目标PHD向k时刻进行预测，其中全体目标包括存活目标和新生目标；

3b)PHD更新：利用k时刻的量测数据z_k，对预测PHD中的新生目标D_k/k-1(y,1)和存活目标D_k/k-1(y,0)分别进行更新，得到前向滤波结果其中w_f,k，m_f,k，P_f,k分别对应前向滤波后的高斯项权值、均值、方差；

4)后向平滑：利用滞后的L时刻的量测信息z_L来平滑前向滤波后的高斯项，得到后向平滑结果其中w_k/L，m_k/L，P_k/L分别对应后向平滑后的高斯项权值，均值，方差，当L＝k+1时为一步后向平滑；

5)删剪合并高斯项：对平滑后权值小于经验门限T_prun的高斯项进行删剪，对均值m_k ⁽ⁱ⁾之间的距离小于门限U_merg的高斯项进行合并，得到删剪合并后的状态估计结果其中w_k ⁽ⁱ⁾，m_k ⁽ⁱ⁾，P_k ⁽ⁱ⁾分别对应删剪合并后的高斯项权值，均值，方差；

6)估计全体目标数目：对修剪合并后的PHD进行权值求和得到全体目标数目估计 $\stackrel{\‾}{N_k}=\underset{i}{\mathrm{\Σ}}{w_k}^{\left(i\right)};$

7)输出最终的状态估计和目标数估计结果得到多目标跟踪的结果。

本发明与现有技术相比具有以下优点：

1.针对目标新生强度未知条件下，ATBI-PHD滤波器和λ-ATBI-PHD滤波器在目标新生时刻，对目标数的估计和目标确认存在滞后的这一问题提出了一种自适应目标新生强度的PHD平滑器，即λ-ATBI-PHD平滑器，分别对新生目标和存活目标进行前向滤波和后向平滑，给出了高斯混合条件下的实现形式。对于杂波环境下的多目标跟踪问题，由于PHD平滑器能够利用更多的滞后时间测量，因此相比于PHD滤波器，它在杂波密度较大以及传感器检测概率较低时可以有效提高多目标的个数和状态估计精度。

2.引入目标新生率估计方法，根据先验杂波强度来估计目标新生率，能够在得到特定杂波数均值的情况下，尽快检测到新生目标。

仿真结果表明，相较于ATBI-PHD滤波器、λ-ATBI-PHD滤波器，本发明的λ-ATBI-PHD平滑器能够更好的估计目标状态，更加准确的估计目标出现时刻，可以减小目标新生时确认的滞后对航迹生成的影响并准确估计目标数目。火控和航管等探测***中常常需要获得准确的航迹信息，本发明的λ-ATBI-PHD平滑器在杂波环境下对于保证多目标跟踪精度和形成正确航迹方面具有积极的意义。

附图说明

图1是本发明方法的流程图

图2是传感器获得的量测数据图

图3是多目标跟踪的目标真实航迹

图4是ATBI-PHD滤波器经过50次蒙特卡洛仿真的跟踪结果

图5是λ-ATBI-PHD滤波器经过50次蒙特卡洛仿真的跟踪结果

图6是λ-ATBI-PHD平滑器经过50次蒙特卡洛仿真的跟踪结果

图7是多目标跟踪的目标数估计图

图8是多目标跟踪的OSPA距离误差图

具体实施方式

依据附图，对本发明的技术方案作具体说明。

本发明使用的量测和状态模型如下：

PHD滤波是将所有量测和目标状态看作两个随机有限集，分别为

Z_k＝{z_k,1,........,z_k,m}∈F(Ζ)(1)

Ξ_k＝S_k(X_k-1)∪Γ_k(X_k-1)∈F(χ)(2)

其中Z_k为k时刻量测集合，Ξ_k为全体目标，Sk为存活目标，Γk为新生目标，F(Z)为量测随机有限集，F(χ)为目标状态的随机有限集，X_k-1为k-1时刻的目标状态。

在传感器每一时刻获得的量测是来自于目标或杂波，而全体目标的量测来源又可分为新生目标和存活目标。自适应目标新生强度方法是把目标的状态空间分为存活目标和新生目标，根据每一时刻的量测信息反演出可能的新生目标位置，同时对全体目标进行递推滤波，使得跟踪方法不再依赖已知目标新生强度的先验信息。

本发明中使用的平滑器分为前向滤波、后向平滑两个步骤。在前向滤波时，后验密度是贝叶斯递推的前向传播；在后向平滑时，利用滞后的L时刻的信息来更新k时刻的状态，完成多目标的后向平滑递推，其中k<L。

如附图1的流程图所示，本发明实现步骤如下：

1)对k-1时刻跟踪滤波器中存活目标PHD的高斯混合形式 $D_{k-1/k-1}(y,0)=\underset{i}{\&Sigma;}{w_{s,k-1}}^{\left(i\right)}N({m_{s,k-1}}^{\left(i\right)};{P_{s,k-1/k-1}}^{\left(i\right)})$ 进行初始化，设定第i个高斯项均值m_s,k-1 ⁽ⁱ⁾，高斯项滤波协方差P_s,k-1/k-1 ⁽ⁱ⁾和高斯项权值w_s,k-1 ⁽ⁱ⁾的初始值，其中i为高斯项的标号，D_k-1/k-1(y,β)表示全体目标PHD，D_k-1/k-1(y,0)表示存活目标PHD，D_k-1/k-1(y,1)表示新生目标PHD，y为目标状态，N(m_s,k-1 ⁽ⁱ⁾；P_s,k-1/k-1 ⁽ⁱ⁾)表示均值为m_s,k-1 ⁽ⁱ⁾，方差为P_s,k-1/k-1 ⁽ⁱ⁾的高斯项；

用β＝0,1来分别表示目标状态x中区分新生目标和存活目标的标志，即

用y表示目标状态x中可以观测到的部分，如目标位置、速度等；

2)根据先验杂波信息对新生目标强度估计：

2a)根据k时刻传感器的测量方差Σ_k/k-1、杂波检测概率p_D ^(c)、目标检测概率p_D ^(t)、目标存活概率p_S ^(t)、先验杂波数均值N^(c)、目标新生权值门限T_e和监控区域体积V_S信息来估计k时刻的目标新生率λ_b,k；

在杂波条件下，第j个量测的来源是杂波或目标，故PHD中第i个高斯项的预测值为：

${D_{k/k-1}}^{\left(i\right)}=\frac{{p_{D,k}}^{\left(c\right)}{N}^{\left(c\right)}+{\mathrm{\&lambda;}}_{b,k}}{V_S}+{p_{D,k}}^{\left(t\right)}{w_{k/k-1}}^i{q_k}^{i,j}---\left(4\right)$

其中p_D,k ^(c)，p_D,k ^(t)分别为杂波和目标的检测概率，λ_b,k为目标新生率；

可以看出，预测的PHD中分为杂波的PHD、新生目标的PHD、存活目标的PHD三部分，p_D,k ^(c)N^(c)/V_S为杂波部分，λ_b,k/V_S为新生目标部分，p_D,k ^(t)w_k/k-1 ⁱq_k ^i,j为存活目标部分，q_k ^(i,j)＝N(z_k ^j；Hm_k/k-1 ⁽ⁱ⁾,R+HP_k/k-1 ⁽ⁱ⁾H^T)为似然函数；

由于B.Ristic提出的自适应目标新生强度模型是假设新生目标强度在监测视场内满足均匀分布，并没有考虑杂波对目标新生的影响。由(4)式可以看出，当杂波数目均值N^(c)增加时，新生部分所占的比例就会相应的下降，即新生目标初始强度降低，会造成目标新生时刻估计的滞后；当新生目标强度部分增加时，又会低估杂波数均值；

LiangMa提出了目标新生率的估计方法，其目标新生率估计式为：

$\stackrel{\&OverBar;}{{\mathrm{\&lambda;}}_{b,k}}=\frac{2\mathrm{\&pi;}\sqrt{\left|{\mathrm{\&Sigma;}}_{k/k-1}\right|}{\left({p_{D,k}}^{\left(c\right)}{N}^{\left(c\right)}\right)}^2}{{p_{D,k}}^{\left(t\right)}{p_{S,k}}^{\left(t\right)}{V}_S({T_e}^{-1}-1)}---\left(5\right)$

将引入目标新生率估计的自适应目标新生强度的PHD滤波称为λ-ATBI-PHD滤波，通过对目标新生率的估计，可以根据杂波环境来设置新生目标的初始权值，能够在得到特定杂波数均值的情况下，使得新生目标在连续两帧数据出现后就可以被检测到；

2b)将估计的新生率赋值给对应的新生高斯项权值w_b ⁽ⁱ⁾，并根据量测数据把新生高斯项的参数w_b,k ⁽ⁱ⁾，m_b,k ⁽ⁱ⁾，P_b,k ⁽ⁱ⁾赋值给并把D_k-1/k-1(y,1)加入到原有的中，完成对监测区域内新生目标的检测；

3)前向滤波：

3a)PHD预测：对全体目标的PHD进行预测，其中全体目标分为新生和存活两部分：

$D_{k/k-1}(y,\mathrm{\&beta;})={\mathrm{\&gamma;}}_{k/k-1}(y,\mathrm{\&beta;})+\underset{{\mathrm{\&beta;}}^{,}=0}{\overset{1}{\&Sigma;}}\&Integral;D_{k-1/k-1}(y^{\&prime;},{\mathrm{\&beta;}}^{\&prime;})\&times;p_s(y^{\&prime;},{\mathrm{\&beta;}}^{\&prime;})f_{k/k-1}(y,\mathrm{\&beta;}|y^{\&prime;},{\mathrm{\&beta;}}^{\&prime;}){\mathrm{dy}}^{\&prime;}---\left(6\right)$

D_k/k-1(y,β)为预测的全体目标PHD强度，y和y’分别表示k和k-1时刻目标的可观测状态，p_s(y,β)为目标的存活概率，f_k/k-1为目标从k-1到k时刻的状态转移函数；

新生目标：D_k/k-1(y,1)＝γ_k(y,1)＝Σw_b,kN(m_b,k；P_b,k)(7)

式中γ_k(y,1)表示根据k时刻的量测反演的新生目标的强度；

存活目标:

由于γ_k(y,0)＝0，故有

D_k/k-1(y,0)＝＜D_k-1/k-1(y,1)+D_k-1/k-1(y,0),ps(y',β')f_k/k-1(y,β|y',β')>＝<p_s(y',1)Σw_b,k-1N(m_b,k-1；P_b,k-1)+p_s(y',0)Σw_s,k-1N(m_s,k-1；P_s,k-1),f_k/k-1(y,β|y',β')＞＝Σw_b,k/k-1N(m_b,k/k-1；P_b,k/k-1)+Σw_s,k/k-1N(m_s,k/k-1；P_s,k/k-1)

(8)

(8)式中

w_b,k/k-1＝w_b,k-1p_s(y',1)(8-1)

m_b,k/k-1＝F_km_b,k-1(8-2)

P_b,k/k-1＝F_kP_b,k-1F_k ^T(8-3)

w_s,k/k-1＝w_s,k-1p_s(y',0)(8-4)

m_s,k/k-1＝F_km_s,k-1(8-5)

P_s,k/k-1＝F_kP_s,k-1F_k ^T+Q(8-6)

w_b,k-1，m_b,k-1，P_b,k-1和w_s,k-1，m_s,k-1，P_s,k-1分别为k-1时刻滤波结果中PHD的新生和存活部分的高斯项权值，均值，方差；

w_b,k/k-1，m_b,k/k-1，P_b,k/k-1和w_s,k/k-1，m_s,k/k-1，P_s,k/k-1分别为k-1时刻预测PHD的新生和存活部分的高斯项权值，均值，方差，Q为状态噪声协方差，D_k-1/k-1(y,0)＝Σw_s,k-1N(m_s,k-1；P_s,k-1)为k-1时刻的存活目标PHD；

3b)PHD更新：根据k时刻的量测数据z_k，对预测PHD中的新生目标D_k/k-1(y,1)和存活目标D_k/k-1(y,0)分别进行更新：

$\begin{array}{c}{D}_{k/k}(y,\mathrm{\&beta;})=\&lsqb;1-p_{D,k}(y,\mathrm{\&beta;})\&rsqb;D_{k/k-1}(y,\mathrm{\&beta;})\\ +\underset{z\&Element;Z_k}{\mathrm{\&Sigma;}}\frac{p_{D,k}(y,\mathrm{\&beta;})g_k\left(z\right|y,\mathrm{\&beta;})D_{k/k-1}(y,\mathrm{\&beta;})}{{\mathrm{\&kappa;}}_k\left(z\right)+\underset{\mathrm{\&beta;}=0}{\overset{1}{\mathrm{\&Sigma;}}}<p_{D,k}(y,\mathrm{\&beta;})g_k\left(z\right|y,\mathrm{\&beta;})D_{k/k-1}(y,\mathrm{\&beta;})>}\\ ={\mathrm{\&Sigma;w}}_{m,k}N(m_{m,k};P_{m,k})+{\mathrm{\&Sigma;w}}_{s,k}N(m_{s,k};P_{s,k})+{\mathrm{\&Sigma;w}}_{b,k}N(m_{b,k};P_{b,k})\\ =\underset{i}{\mathrm{\&Sigma;}}{w_{k/k}}^{\left(i\right)}N({m_{k/k}}^{\left(i\right)};{P_{k/k}}^{\left(i\right)})\end{array}---\left(9\right)$

其中D_k/k(y,β)为更新后的全体目标的PHD强度，g_k(z|y,β)为目标的似然函数，κ_k(z)＝N⁽ⁱ⁾/V_S为监视区域内的杂波强度，为得到的全体目标的PHD前向滤波结果，p_D,k(y,β)为目标的检测概率，有

$p_{D,k}(y,\mathrm{\&beta;})=\left\{\begin{array}{cc}1,& \mathrm{\&beta;}=0\\ p_{D,k}\left(y\right),& \mathrm{\&beta;}=1\end{array}\right\};$

(9)式中

w_m,k ⁽ⁱ⁾＝(1-P_D,k)w_s,k/k-1 ⁽ⁱ⁾(9-1)

${w_{s,k}}^{\left(i\right)}=\frac{P_{D,k}{g}_k{\left(z\right|y,0)}^{\left(i\right)}{w_{s,k/k-1}}^{\left(i\right)}}{{\mathrm{\&kappa;}}_k\left(z\right)+P_{D,k}\underset{i}{\&Sigma;}{w_{k/k-1}}^{\left(i\right)}{g}_k{\left(z\right|y,0)}^{\left(i\right)}+{w_{b,k}}^{\left(i\right)}/V_S}---(9-2)$

${w_{b,k}}^{\left(i\right)}=\frac{w_{b,k}/V_S}{{\mathrm{\&kappa;}}_k\left(z\right)+P_{D,k}\underset{i}{\&Sigma;}{w_{k/k-1}}^{\left(i\right)}{g}_k{\left(z\right|y,0)}^{\left(i\right)}+w_{b,k}/V_S}---(9-3)$

均值m_m,k ⁽ⁱ⁾，m_s,k ⁽ⁱ⁾，m_b,k ⁽ⁱ⁾的更新式都形如其中m_k ⁽ⁱ⁾＝m_m,k ⁽ⁱ⁾，m_s,k ⁽ⁱ⁾，m_b,k ⁽ⁱ⁾；

协方差P_m,k ⁽ⁱ⁾，P_s,k ⁽ⁱ⁾，P_b,k ⁽ⁱ⁾的更新式都形如P_k ⁽ⁱ⁾＝[I-K_k ⁽ⁱ⁾H]P_k/k-1 ⁽ⁱ⁾，其中P_k ⁽ⁱ⁾＝P_m,k ⁽ⁱ⁾，P_s,k ⁽ⁱ⁾，P_b,k ⁽ⁱ⁾；

$\stackrel{\&OverBar;}{{z_{k/k-1}}^{\left(i\right)}}={{\mathrm{Hm}}_{k/k-1}}^{\left(i\right)}---\left(10\right)$

K_k ⁽ⁱ⁾＝P_k/k-1 ⁽ⁱ⁾H^T[S_k/k-1 ⁽ⁱ⁾]^-1(11)

S_k/k-1 ⁽ⁱ⁾＝HP_k/k-1 ⁽ⁱ⁾H^T+R(12)

其中w_m,k为更新后漏检部分权值，w_s,k为更新后存活部分权值，w_b,k为更新后新生部分权值，H为观测矩阵，R为观测噪声协方差，K_k ⁽ⁱ⁾为卡尔曼滤波增益；

4)后向平滑：利用滞后的L时刻的量测信息z_L来更新k时刻的状态，完成多目标的后向平滑递推，其中当L＝k+1时为一步后向平滑；

平滑器的后向更新方程为：

$D_{k/L}\left(y_k\right|Z_{1:L})=D_{k/k}\left(y_k\right|Z_{1:k})\&lsqb;1-p_{S,k+1/k}+\&Integral;D_{k+1/L}\left(y_{k+1}\right|Z_{1:L})\frac{F_{k+1/k}\left(y_{k+1}\right|y_k)}{D_{k+1/k}\left(y_{k+1}\right|Z_{1:k})}{\mathrm{\&delta;y}}_{k+1}\&rsqb;---\left(13\right)$

其中y_k为k时刻目标的可观测状态，Z_1:k为1到k时刻的量测集合，p_S,k+1/k为存活概率，D_k/k(y_k|Z_1:k)为k时刻存活目标和新生目标的PHD滤波结果，F_k+1/k(y_k+1|y_k)为k-1到k时刻的状态转移函数，D_k+1/k(y_k+1|Z_1:k)为k到k+1时刻存活目标和新生目标的PHD预测，D_k+1/L(y_k+1|Z_1:k)为L向k+1时刻存活目标和新生目标的后向递推结果；

上式可分解为：

D_k/L(y_k|Z_1:L)＝D_k/k(y_k|Z_1:k)[1-p_S,k+1/k+B_k/L(y_k|Z_1:L)](14)

其中B_k/L(y_k|Z_1:L)为L时刻对k时刻的后向平滑算子：

B_k/L(y_k|Z_1:L)＝p_s,k+1/k(y_k)＜B_k+1/kL_k+1(Z_k+1|y_k'),f_k+1/k(y_k'|y_k)＞+q_s,k+1/k(y_k)(15)

自适应目标新生强度的PHD平滑器的高斯混合实现如下：

$\begin{array}{c}{D}_{k/L}\left(y_k\right|Z_{1:L})=\underset{I\&SubsetEqual;k+1:L}{\&Sigma;}\underset{z_L\&Element;Z_L}{\&Sigma;}\underset{i}{\&Sigma;}{w_{k/k}}^{\left(i\right)}N({m_{k/k}}^{\left(i\right)},{P_{k/k}}^{\left(i\right)})\&lsqb;1-p_{S,k+1/k}+\frac{F_{k+1/k}{q}_{k+1/I}\left(z_I\right)}{{w_{k+1/k}}^{\left(i\right)}N({m_{k+1/k}}^{\left(i\right)},{P_{k+1/k}}^{\left(i\right)})}\&rsqb;\\ =\underset{i}{\&Sigma;}{w_{k/L}}^{\left(i\right)}N({m_{k/L}}^{\left(i\right)},{P_{k/L}}^{\left(i\right)})\end{array}---\left(16\right)$

(16)式中

q_k+1/I(z_I)＝N(z_I|m_k+1/k+1 ⁽ⁱ⁾(z_I),R_I/k+H_I/kP_k/kH_I/k ^(T))(16-1)

${w_{k/L}}^{\left(i\right)}\left(z_I\right)={w_{k/k}}^{\left(i\right)}\&times;\&lsqb;1-p_{S,k+1/k}+\underset{i}{\&Sigma;}\frac{q_{k+1/I}\left(z_I\right)}{{w_{k+1/k}}^{\left(i\right)}}\&rsqb;---(16-2)$

w_k+1/k ⁽ⁱ⁾＝w_k/k ⁽ⁱ⁾p_S,k+1/k(16-3)

m_k+1/k ⁽ⁱ⁾＝F_k+1/km_k/k ⁽ⁱ⁾(16-4)

P_k+1/k ⁽ⁱ⁾＝F_k+1/kP_k/k ⁽ⁱ⁾F_k+1/k ^T(16-5)

其中q_k+1/I(z_I)为量测z_I的似然函数， $D_{k/k}\left(y_k\right|Z_{1:k})=\underset{i}{\&Sigma;}{w_{k/k}}^{\left(i\right)}N({m_{k/k}}^{\left(i\right)},{P_{k/k}}^{\left(i\right)}),$ D_k+1/k(y_k+1|Z_1:k)＝w_k+1/k ⁽ⁱ⁾N(m_k+1/k ⁽ⁱ⁾,P_k+1/k ⁽ⁱ⁾)，I表示平滑使用的滞后量测的所在时刻，z_I为I时刻的量测，F_k+1/k为k到k+1时刻的状态转移矩阵，w_k+1/k ⁽ⁱ⁾，m_k+1/k ⁽ⁱ⁾，P_k+1/k ⁽ⁱ⁾分别对应k到k+1时刻的预测高斯项权值，预测均值，预测协方差矩阵，w_k/L ⁽ⁱ⁾，m_k/L ⁽ⁱ⁾，P_k/L ⁽ⁱ⁾分别对应L时刻对k时刻平滑结果中的高斯项权值，预测均值，预测协方差矩阵；

5)删剪合并高斯项：

根据经验门限T_prun，对平滑后PHD中权值小于门限的高斯项进行删剪，即若w_k/L ⁽ⁱ⁾＜T_prun，则将w_k/L ⁽ⁱ⁾,m_k/L ⁽ⁱ⁾,P_k/L ⁽ⁱ⁾对应的高斯项从中剔除，得到删剪合并后的状态估计结果其中w_k ⁽ⁱ⁾，m_k ⁽ⁱ⁾，P_k ⁽ⁱ⁾分别对应删剪合并后的高斯项权值，均值，方差；

将高斯项均值m_k/L ⁽ⁱ⁾之间的距离小于门限U_merg的高斯项进行合并，合并后的权值w_k ⁽ⁱ⁾为合并前各高斯项权值w_k/L ⁽ⁱ⁾之和，即若有|m_k ⁽ⁱ⁾-m_k ⁽ⁱ⁺¹⁾|＜U_merg,则

w_k ⁽ⁱ⁾＝w_k/L ⁽ⁱ⁾+w_k/L ⁽ⁱ⁺¹⁾(17)

m_k ⁽ⁱ⁾＝w_k/L ⁽ⁱ⁾m_k/L ⁽ⁱ⁾+w_k/L ⁽ⁱ⁺¹⁾m_k/L ⁽ⁱ⁺¹⁾(18)

${P_k}^{\left(i\right)}=\frac{1}{w_{k/L}}\underset{i}{\&Sigma;}{w_{k/L}}^{\left(i\right)}({P_{k/L}}^{\left(i\right)}+(m_{k/L}-{m_{k/L}}^{\left(i\right)}\left){\left(m_{k/L}-{m_{k/L}}^{\left(i\right)}\right)}^T\right)---\left(19\right)$

6)估计目标数目：对修剪合并后的PHD进行权值求和得到目标数目估计 $\stackrel{\&OverBar;}{N_k}=\underset{i}{\&Sigma;}{w_k}^{\left(i\right)};$

7)输出最终的状态估计和目标数估计结果得到多目标跟踪的结果。

下面结合附图对本发明的仿真效果做进一步的描述。

1.仿真条件：

本发明的仿真是在主频3.0GHZ的Intel(R)Pentium(R)CPUG2030、内存4.00GB的硬件环境和MATLABR2009b的软件环境下进行的。

本实验为对比ATBI-PHD滤波器、λ-ATBI-PHD滤波器和λ-ATBI-PHD平滑器三种方法的跟踪效果，取L＝k+1，即一步后向平滑，仿真场景设置如下：

整个观测区域为[-100,100]×[-100,100]㎡，采样周期为1s，观测时刻k＝1～40,观测过程持续40帧，相继出现4个目标，不考虑衍生目标情况，杂波数均值λ_c,k＝2；

假设在线性高斯条件下，目标在k时刻的状态向量为状态分别为目标的x轴坐标，x方向速度，y轴坐标，y方向速度，其运动方程为：

X_k＝F_kX_k-1+G_kw_k

$\begin{array}{ccc}{F}_k=\left[\begin{array}{cccc}1& T& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& T\\ 0& 0& 0& 1\end{array}\right]& G_k=\left[\begin{array}{cc}{T}^2/2& 0\\ T& 0\\ 0& T^2/2\\ 0& T\end{array}\right]& w_k=\left[\begin{array}{c}{\mathrm{\&sigma;}}_{1,w}\\ {\mathrm{\&sigma;}}_{2,w}\end{array}\right]\end{array}$

其中，符号Fk代表目标匀速(CV)运动模型，G_k代表扰动矩阵，w_k为状态噪声，w_k服从 $N(0,\left[\begin{array}{cc}{\mathrm{\&sigma;}}_{1,w}& 0\\ 0& {\mathrm{\&sigma;}}_{2,w}\end{array}\right])$ 高斯分布，σ_w为状态噪声根方差，设置σ_1,w＝0.01，σ_2,w＝0.01；

观测方程为：Z_k＝H_kX_k+n_k

其中，观测矩阵 $H_k=\left[\begin{array}{cc}1& 0\\ 0& 0\\ 0& 1\\ 0& 0\end{array}\right],$ n_k为观测噪声，n_k服从 $N\left(\begin{array}{cc}0,& \left[\begin{array}{cc}{\mathrm{\&sigma;}}_n& 0\\ 0& {\mathrm{\&sigma;}}_n\end{array}\right]\end{array}\right)$ 高斯分布，设置σ_n＝0；

假设杂波强度在场景内服从均匀分布，杂波数目服从参数λ_c,k的泊松分布。目标存活概率Ps＝0.9，新生目标的检测概率Pd_birth＝1，存活目标的检测概率Pd＝0.99；

4个目标在视场内相继出现后运动至仿真结束，新生时刻分别为：目标1在时刻k＝1，目标2在时刻k＝9，目标3在时刻k＝13，目标4在时刻k＝20；

4个目标初始状态设定：目标1为[0,0.5,10,1]，目标2为[10,0.3,0,-1.1]，目标3为[20,0.6,0,0.9]，目标4为[30,0.7,0,-1.1]；

2.仿真结果及分析

传感器量测如图2所示，目标的真实航迹如图3所示，图中“□”表示目标的起始点。

采用目标状态估计、目标个数均值和最优次模式分配(optimalsubpatternassignment，OSPA)距离对各方法的跟踪性能进行评价。

$d_{p,c}(X,Z)={\{\frac{1}{n}\&lsqb;\underset{\mathrm{\&pi;}\&Element;\&Pi;n}{min}\underset{i=1}{\overset{m}{\&Sigma;}}{d}_c{(x_i,z_{\mathrm{\&pi;}\left(i\right)})}^p+c^p(n-m)\&rsqb;\}}^{1/p}$

其中，设定参数c＝70，p＝2，p为距离敏感性参数，c是一个水平调节数，用于调节集合势误差的影响，X和Z为任意两个有限集合，m，n分别对应X，Z两个集合的元素个数，且m≤n，d(x_i,z_π(i))表示x_i和z_π(i)之间的一阶距离，。

图4，图5，图6分别为ATBI-PHD滤波器，λ-ATBI-PHD滤波器和λ-ATBI-PHD平滑器经过50次蒙特卡洛仿真的跟踪结果。由于λ-ATBI-PHD平滑器采用了一步平滑，所以有效数据处理的时刻为k＝2～39。可以看出，在有效时刻上，三者在未知新生目标强度的情况下，都能完成对多目标的跟踪。

图7为三种方法的目标数的估计与真实目标数目。通过对比可以看出，估计目标新生率的λ-ATBI-PHD滤波器比ATBI-PHD滤波器的目标数估计效果更加平稳，这是由于估计目标的新生率，可以得到在得到特定杂波数均值的情况下合适的新生目标初始权值，但是对新生目标的确认滞后也更大。ATBI-PHD滤波器和λ-ATBI-PHD滤波器在k＝9,13,20的目标新生的时刻，对目标数目都会产生低估，其原因如(4)式的分析。相比之下，λ-ATBI-PHD平滑器在杂波环境下对目标数的估计更加准确，在目标新生时刻能够及时检测到目标，有利于后续的航迹起始处理。

图8为三种方法的OSPA距离对比，同样，在k＝9,13,20的目标新生的时刻，ATBI-PHD滤波器和λ-ATBI-PHD滤波器的OSPA距离明显增大，但本发明的λ-ATBI-PHD平滑器的OSPA距离要明显小于ATBI-PHD滤波器和λ-ATBI-PHD滤波器，表明本发明的λ-ATBI-PHD平滑器有更好的跟踪精度。

综上所述，在杂波环境下，本发明的λ-ATBI-PHD平滑器通过后向平滑能够更准确地估计目标新生时刻，减小目标新生时刻的确认滞后情况对后续航迹形成的影响，获得更好的目标跟踪精度。火控和航管等探测***中常常需要获得准确的航迹起始信息，本发明的λ-ATBI-PHD平滑器在杂波环境下对于保证多目标跟踪精度和形成正确航迹方面具有积极的意义。

Claims (3)

1.一种自适应目标新生强度的PHD平滑器的多目标跟踪方法，包括如下步骤：

(1)对k-1时刻跟踪滤波器中存活目标的PHD的高斯混合形式 $D_{k-1/k-1}(y,0)=\underset{i}{\&Sigma;}{w_{s,k-1}}^{\left(i\right)}N({m_{s,k-1}}^{\left(i\right)};{P_{s,k-1/k-1}}^{\left(i\right)})$ 进行初始化，设定第i个高斯项均值m_s,k-1 ⁽ⁱ⁾，高斯项滤波协方差P_s,k-1/k-1 ⁽ⁱ⁾和高斯项权值w_s,k-1 ⁽ⁱ⁾的初始值，其中i为高斯项的标号，D_k-1/k-1(y,β)表示全体目标的PHD，β为区分新生目标和存活目标的标志，D_k-1/k-1(y,0)表示存活目标的PHD，D_k-1/k-1(y,1)表示新生目标的PHD，y为目标的可观测状态，N(m_s,k-1 ⁽ⁱ⁾；P_s,k-1/k-1 ⁽ⁱ⁾)表示均值为m_s,k-1 ⁽ⁱ⁾，方差为P_s,k-1/k-1 ⁽ⁱ⁾的高斯项；

(2)根据先验杂波信息对新生目标强度估计：

2a)根据k时刻传感器的测量方差Σ_k/k-1、杂波检测概率p_D ^(c)、目标检测概率p_D ^(t)、目标存活概率p_S ^(t)、先验杂波数均值N^(c)、目标新生权值门限T_e和监控区域体积V_S信息来估计k时刻的目标新生率λ_b,k；

2b)根据估计的目标新生率λ_b,k和量测数据z_k，把新生目标的PHD中的高斯项加入到k-1时刻跟踪结果的全体目标的高斯项中，完成对监测区域内新生目标强度的检测，式中w_b,k ⁽ⁱ⁾＝λ_b,k，w_b,k，m_b,k，P_b,k分别对应新生目标高斯项的权值，均值，方差；

(3)前向滤波，包括PHD预测和PHD更新：

3a)PHD预测：对k-1时刻的全体目标PHD向k时刻进行预测，其中全体目标包括存活目标和新生目标；

3b)PHD更新：利用k时刻的量测数据z_k，对预测PHD中的新生目标D_k/k-1(y,1)和存活目标D_k/k-1(y,0)分别进行更新，得到前向滤波结果其中w_f,k，m_f,k，P_f,k分别对应前向滤波后的高斯项权值，均值，方差；

(4)后向平滑：利用滞后的L时刻的量测信息z_L来平滑前向滤波后的高斯项，得到后向平滑结果其中w_k/L，m_k/L，P_k/L分别对应后向平滑后的高斯项的权值，均值，方差；当L＝k+1时为一步后向平滑；

(5)删剪合并高斯项：对平滑后权值小于经验门限T_prun的高斯项进行删剪，对均值m_k ⁽ⁱ⁾之间的距离小于门限U_merg的高斯项进行合并，得到删剪合并后的状态估计结果其中w_k ⁽ⁱ⁾，m_k ⁽ⁱ⁾，P_k ⁽ⁱ⁾分别对应删剪合并后的高斯项权值，均值，方差；

(6)估计全体目标数目：对修剪合并后的PHD进行权值求和得到全体目标数目估计 $\stackrel{\&OverBar;}{N_k}=\underset{i}{\&Sigma;}{w_k}^{\left(i\right)};$

(7)输出最终的状态估计和目标数估计结果得到多目标跟踪的结果。

2.根据权利要求1所述的自适应目标新生强度的PHD平滑器的多目标跟踪方法，其特征在于：步骤(2)所述的目标新生率的估计计算式为：

$\stackrel{\&OverBar;}{{\mathrm{\&lambda;}}_{b,k}}=\frac{2\mathrm{\&pi;}\sqrt{\left|{\mathrm{\&Sigma;}}_{k/k-1}\right|}{\left({p_{D,k}}^{\left(c\right)}{N_{k/k-1}}^{\left(c\right)}\right)}^2}{{p_{D,k}}^{\left(t\right)}{p_{S,k}}^{\left(t\right)}{V}_S({T_e}^{-1}-1)}$

式中Σ_k/k-1为传感器的测量方差，p_D,k ^(c)，p_D,k ^(t)分别为杂波和目标的检测概率，p_S,k ^(t)为目标的存活概率，N^(c)为先验杂波数均值，V_S为监控区域体积，T_e为新生一个目标的新生权值门限，表示目标新生率λ_b,k的估计值。

3.根据权利要求1所述的自适应目标新生强度的PHD平滑器的多目标跟踪方法，其特征在于：步骤(4)所述的根据L时刻的量测信息对k时刻全体目标的前向滤波后的高斯项进行平滑，得到k时刻的状态估计，全体目标包括存活目标和新生目标，其中PHD平滑的计算式为：

$D_{k/L}\left(y_k\right|Z_{1:L})=D_{k/k}\left(y_k\right|Z_{1:k})\&lsqb;1-p_{S,k+1/k}+\&Integral;D_{k+1/L}\left(y_{k+1}\right|Z_{1:L})\frac{F_{k+1/k}\left(y_{k+1}\right|y_k)}{D_{k+1/k}\left(y_{k+1}\right|Z_{1:k})}{\mathrm{\&delta;y}}_{k+1}\&rsqb;$

其中y_k为k时刻目标的可观测状态，Z_1:k为1到k时刻的量测集合，p_S,k+1/k为存活概率，D_k/k(y_k|Z_1:k)为k时刻存活目标和新生目标的PHD滤波结果，F_k+1/k(y_k+1|y_k)为k-1到k时刻的状态转移函数，D_k+1/k(y_k+1|Z_1:k)为k到k+1时刻存活目标和新生目标的PHD预测，D_k+1/L(y_k+1|Z_1:k)为L向k+1时刻存活目标和新生目标的后向递推结果。