CN104794735A - 基于变分贝叶斯期望最大化的扩展目标跟踪方法 - Google Patents

Abstract

本发明公开了一种基于变分贝叶斯期望最大化VBEM的扩展目标跟踪方法，主要解决传统的扩展目标跟踪领域中，在量测噪声协方差是未知的情况下，目标的跟踪性能将会急剧下降的问题。该方法首先预测了目标状态和量测噪声协方差的联合概率假设密度中高斯逆伽玛分量的相关参数；再对高斯逆伽玛分量参数进行更新；最后通过修建与合并得到扩展目标状态和数目。仿真实验表明，本发明能很好地跟踪未知数目和未知量测噪声协方差下的多扩展目标，且具有较高的跟踪精度，可用于对飞行器和潜艇目标的跟踪。

Description

基于变分贝叶斯期望最大化的扩展目标跟踪方法

技术领域

本发明属于信息处理技术领域，特别涉及一种目标跟踪方法,可用于跟踪多扩展目标。

背景技术

在传统的目标跟踪领域里，由于雷达的分辨率有限，因此通常情况下将目标看作是点目标，即一个时刻每个目标只能产生一个量测。近年来，随着雷达探测技术的发展和实际应用的需要，更多的将目标视为扩展目标，即每个目标在每一个时刻都可以产生多个量测。

实际的目标跟踪场景中，目标的数目是无法提前预知的，因此随机集理论的提出极大的满足了目标跟踪理论的需要。在诸多对目标的模型假设中，尤其属扩展目标理论的提出更加贴近目前跟踪理论的需要，并且在实际生活中得到了广泛的应用，近十年来，成为了目标跟踪领域中的研究热点。2003年，Mahler将随机有限集理论应用于多目标跟踪问题，提出了概率假设密度PHD滤波。2005年，Gilholm和Salmond提出一种空间分布服从泊松分布的扩展目标模型。2009年，Mahler推导出了扩展目标PHD滤波，即用每一时刻的量测随机集对目标随机集进行预测、更新，继而可以准确提取目标的运动状态和估计目标的数目。2010年，等给出了扩展目标PHD的高斯混合实现形式。2011年，Orguner等又提出了带势分布的扩展目标PHD(ET-CPHD)滤波，很好的解决了ET-PHD估计目标数目时的缺陷。然而，传统的扩展目标跟踪算法处理的都是量测噪声协方差是已知的情形，实际中当量测噪声协方差未知时，扩展目标的跟踪性能将会急剧下降。

发明内容

本发明的目的在于针对上述问题，提出一种基于变分贝叶斯期望最大化的扩展目标跟踪方法，以提高在量测噪声协方差未知条件下的跟踪性能。

实现本发明的技术关键是：在势概率假设密度滤波框架下，引入变分贝叶斯期望最大化VBEM技术，估计目标状态和量测噪声协方差的联合概率假设密度，实现未知量测噪声协方差下的目标跟踪问题。

VB是一类用于贝叶斯估计和机器学习领域中近似计算复杂积分的方法，本文中VB被用来近似线性高斯***中，扩展目标状态和量测噪声协方差的联合概率假设密度，该算法的主要思想是对扩展目标状态和量测噪声协方差的联合概率假设密度进行参数化近似,并给出其参数化表达形式。在用VB理论对扩展目标状态和量测噪声协方差的联合概率假设密度进行近似的过程中，为了判断所估计的高斯逆伽玛分量相关参数的性能，在VB的基础上又引入了期望最大化EM算法。在期望E步，估计未知参数的期望值，给出当前的参数估计；在最大化M步，重新估计分布参数，以使该似然函数达到最大。M步上得到的参数估计值被用于下一个E步计算，这个过程不断交替进行。这种方法可以广泛的应用于数据有缺损的情况，且具有简单性和稳定性等优点。

本发明利用上述VBEM技术进行扩展目标跟踪的技术步骤包括如下：

(1)在时刻k＝0时，初始化扩展目标状态和量测噪声协方差的联合概率假设密度为v₀(x,R)；

(2)在k≥1时，对扩展目标状态和量测噪声协方差的联合概率假设密度v_k-1|k-1(x,R)及用于计算扩展目标数目的势分布P_k-1|k-1(num)进行预测，得到预测的扩展目标联合概率假设密度v_k|k-1(x,R)和预测势分布P_k|k-1(num)；

(3)对预测的扩展目标状态和量测噪声协方差的联合概率假设密度v_k|k-1(x,R)及用于计算扩展目标数目的势分布P_k|k-1(num)进行更新：

(3a)利用变分贝叶斯VB方法对联合概率假设密度v_k|k-1(x,R)进行近似，得到用高斯分布的求和形式表示的扩展目标状态的概率假设密度Q_x,k|k-1(x)和用逆伽玛分布的求和形式表示的量测噪声协方差的概率假设密度Q_R,k|k-1(R)；

(3b)利用变分贝叶斯期望最大化VBEM方法对扩展目标状态的概率假设密度Q_x,k|k-1(x)中的高斯分量和量测噪声协方差的概率假设密度Q_R,k|k-1(R)中的逆伽玛分量进行迭代更新，得到表示扩展目标运动状态x的高斯分量和表示量测噪声协方差R的逆伽玛分量；

(3c)对步骤(2)预测得到的势分布P_k|k-1(num)进行更新，得到更新后的势分布P_k|k(num)；

(4)对更新后的高斯分量和逆伽玛分量进行修剪与合并，并提取合并后的高斯分量和逆伽玛分量的位置和速度作为扩展目标的状态；

(5)对步骤(3)更新得到的势分布P_k|k(num)进行加权平均，得到扩展目标的数目： ${\mathrm{num}}_{k|k}=\underset{\mathrm{num}=1}{\overset{\∞}{\mathrm{\Σ}}}\mathrm{num}\×p_{k|k}\left(\mathrm{num}\right);$

(6)重复步骤(2)-(5)，继续跟踪扩展目标。

本发明具有以下优点：

第一，本发明引入了变分贝叶斯EM技术，通过估计扩展目标状态和量测噪声协方差的联合概率假设密度，有效估计了各个目标在不同时刻的真实量测噪声，为复杂环境下多扩展目标跟踪场景的分析提供了帮助，保证了CPHD滤波算法可以有效地实现对未知量测噪声协方差环境下扩展目标跟踪。

第二，本发明由于采用先预测高斯逆伽玛分量的相关参数，再对高斯逆伽玛分量参数进行更新，最后通过修剪与合并得到扩展目标状态和数目的跟踪扩展目标的过程，与传统的GM-CPHD滤波算法相比，提高了跟踪精度。

附图说明

图1是本发明的总流程图；

图2是单次实验条件下，用本发明方法跟踪扩展目标的仿真结果图；

图3是在100次实验条件下，用本发明方法与传统GM-CPHD方法估计目标数目的仿真结果对比图；

图4是在100次实验条件下，用本发明方法与传统GM-CPHD方法借助OSPA距离来判断目标跟踪精度的仿真结果对比图。

具体实施方式

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

参照图1，本发明的具体实现步骤包括如下：

步骤1，在时刻k＝0时，初始化目标状态和量测噪声协方差的联合概率假设密度：

$v_0(x,R)=\underset{i=1}{\overset{J_0}{\mathrm{\Σ}}}\left[w_0^{\left(i\right)}N(x;m_0^{\left(i\right)},P_0^{\left(i\right)})\underset{l=1}{\overset{d}{\mathrm{\Π}}}\mathrm{IG}({\left({\mathrm{\σ}}_{0,l}^{\left(i\right)}\right)}^2;{\mathrm{\α}}_{0,l}^{\left(i\right)},{\mathrm{\β}}_{0,l}^{\left(i\right)})\right]$

其中，J₀表示高斯分量的个数,表示第i个高斯分量的权值，N(·)表示高斯分布，表示第i个高斯分量的均值，表示第i个高斯分量的协方差；IG(·)表示逆伽玛分布，表示第i个逆伽玛分量的协方差,为第i个逆伽玛分量的常量因子，为第i个逆伽玛分量的迭代因子，l表示量测噪声协方差的第l维，d表示量测噪声协方差的维数。

步骤2，在k≥1时，对扩展目标状态和量测噪声协方差的联合概率假设密度v_k-1|k-1(x,R)进行预测，得到预测的扩展目标联合概率假设密度v_k|k-1(x,R)。

2a)对存活目标状态和量测噪声协方差的联合概率假设密度v_S,k|k-1(x,R)中高斯分量的均值和协方差进行预测，得到预测的存活目标高斯分量的均值和协方差

$m_{S,k|k-1}^{\left(i\right)}=F_{S,k-1}^{\left(i\right)}{m}_{S,k-1}^{\left(i\right)}$

$P_{S,k|k-1}^{\left(i\right)}=Q_{k-1}+F_{S,k-1}^{\left(i\right)}{P}_{S,k-1}^{\left(i\right)}{\left(F_{S,k-1}^{\left(i\right)}\right)}^T$

其中，表示存活目标状态转移矩阵，Q_k-1表示存活目标过程噪声协方差，(·)^T表示转置；

2b)对存活目标状态和量测噪声协方差的联合概率假设密度v_S,k|k-1(x,R)中逆伽玛分量的常量因子和迭代因子进行预测，得到预测的存活目标逆伽玛分量的常量因子和逆伽玛分量的迭代因子

${\mathrm{\α}}_{S,k|k-1,l}^{\left(i\right)}={\mathrm{\ρ}}_l{\mathrm{\α}}_{S,k-1,l}^{\left(i\right)}$

${\mathrm{\β}}_{S,k|k-1,l}^{\left(i\right)}={\mathrm{\ρ}}_l{\mathrm{\β}}_{S,k-1,l}^{\left(i\right)}$

其中ρ_l表示遗忘因子，且ρ_l∈(0,1]；

2c)在高斯混合框架下，利用预测的存活目标高斯分量的均值协方差和逆伽玛分量的常量因子迭代因子计算存活目标状态和量测噪声协方差的联合概率假设密度v_S,k|k-1(x,R)：

$v_{S,k|k-1}(x,R)=P_{S,k}\underset{i=1}{\overset{J_{k-1}}{\mathrm{\Σ}}}\left[w_{k-1}^{\left(i\right)}N(x;m_{S,k|k-1}^{\left(i\right)},P_{S,k|k-1}^{\left(i\right)})\underset{l=1}{\overset{d}{\mathrm{\Π}}}\mathrm{IG}({\left({\mathrm{\σ}}_{S,k|k-1}^{\left(i\right)}\right)}^2;{\mathrm{\α}}_{S,k|k-1,l}^{\left(i\right)},{\mathrm{\β}}_{S,k|k-1.l}^{\left(i\right)})\right]$

其中，P_S,k表示扩展目标存活概率，J_k-1表示k-1时刻的高斯分量数，表示k-1时刻第i个高斯分量的权值，N(·)表示高斯分布，表示预测的存活目标状态和量测噪声协方差的联合概率假设密度v_S,k|k-1(x,R)中第i个高斯分量的均值，表示预测的存活目标状态和量测噪声协方差的联合概率假设密度v_S,k|k-1(x,R)中第i个高斯分量的协方差；IG(·)表示逆伽玛分布，d表示量测噪声协方差的维数，表示预测的存活目标状态和量测噪声协方差的联合概率假设密度v_S,k|k-1(x,R)中第i个逆伽玛分量的常量因子，表示预测的存活目标状态和量测噪声协方差的联合概率假设密度v_S,k|k-1(x,R)中第i个逆伽玛分量的迭代因子；

2d)对衍生目标状态和量测噪声协方差的联合概率假设密度b_k|k-1(x,R)中高斯分量的均值和协方差进行预测，得到预测的衍生目标高斯分量的均值和协方差

$m_{b,k|k-1}^{(i,j)}=F_{b,k-1}^{\left(j\right)}{m}_{b,k-1}^{(i,j)}+d_{b,k-1}^{\left(j\right)}$

$P_{b,k|k-1}^{(i,j)}=Q_{b,k-1}^{\left(j\right)}+F_{b,k-1}^{\left(j\right)}{P}_{b,k-1}^{\left(i\right)}{\left(F_{b,k-1}^{\left(j\right)}\right)}^T$

其中，i表示k-1时刻高斯逆伽玛分量的第i个,j表示由k-1时刻到k时刻衍生出的高斯逆伽玛分量的第j个，表示衍生目标的状态转移矩阵，表示衍生目标的状态修正量，表示衍生目标过程噪声协方差；

2e)对衍生目标状态和量测噪声协方差的联合概率假设密度b_k|k-1(x,R)中逆伽玛分量的常量因子和迭代因子进行预测，得到预测的衍生目标逆伽玛分量的常量因子和迭代因子

${\mathrm{\α}}_{b,k|k-1,l}^{(i,j)}={\mathrm{\ρ}}_l{\mathrm{\α}}_{b,k-1,l}^{\left(i\right)}$

${\mathrm{\β}}_{b,k|k-1,l}^{(i,j)}={\mathrm{\ρ}}_l{\mathrm{\β}}_{b,k-1,l}^{\left(i\right)};$

2f)在高斯混合框架下,利用预测的衍生目标高斯分量的均值协方差和逆伽玛分量的常量因子迭代因子计算衍生目标状态和量测噪声协方差的联合概率假设密度b_k|k-1(x,R)：

$b_{k|k-1}(x,R)=\underset{i=1}{\overset{J_{k-1}}{\mathrm{\Σ}}}\underset{j=1}{\overset{J_{b,k}}{\mathrm{\Σ}}}\left[w_{k-1}^{\left(i\right)}{w}_{b,k}^{\left(j\right)}N(x;m_{b,k|k-1}^{(i,j)},P_{b,k|k-1}^{(i,j)})\underset{l=1}{\overset{d}{\mathrm{\Π}}}\mathrm{IG}({\left({\mathrm{\σ}}_{b,k|k-1,l}^{(i,j)}\right)}^2;{\mathrm{\α}}_{b,k|k-1,l}^{(i,j)},{\mathrm{\β}}_{b,k|k-1,l}^{(i,j)})\right]$

其中，J_b,k表示k-1时刻到k时刻的衍生目标高斯分量数，表示k时刻第j个衍生目标高斯分量的权值，表示k-1时刻到k时刻由第i个高斯分量衍生得到第j个高斯分量的均值，表示k-1时刻到k时刻由第i个高斯分量衍生得到第j个高斯分量的协方差；表示k-1时刻到k时刻由第i个逆伽玛分量衍生得到第j个逆伽玛分量的常量因子，表示k-1时刻到k时刻由第i个逆伽玛分量衍生得到第j个逆伽玛分量的迭代因子；

2g)计算新生目标状态和量测噪声协方差的联合概率密度γ_k(x,R)：

${\mathrm{\γ}}_k(x,R)=\underset{i=1}{\overset{J_{\mathrm{\γ},k}}{\mathrm{\Σ}}}\left[w_{\mathrm{\γ}.k}^{\left(i\right)}N(x;m_{\mathrm{\γ},k}^{\left(i\right)},P_{\mathrm{\γ},k}^{\left(i\right)})\underset{l=1}{\overset{d}{\mathrm{\Π}}}\mathrm{IG}({\left({\mathrm{\σ}}_{\mathrm{\γ},k,l}^{\left(i\right)}\right)}^2;{\mathrm{\α}}_{\mathrm{\γ},k,l}^{\left(i\right)},{\mathrm{\β}}_{\mathrm{\γ},k,l}^{\left(i\right)})\right]$

其中，J_γ,k为k时刻新生目标高斯分量数，为为k时刻新生目标第i个高斯分量的权值，为新生目标第i个高斯分量的状态均值，为新生目标第i个高斯分量的运动状态协方差；为新生目标第i个逆伽玛分量的常量因子，为新生目标第i个逆伽玛分量的迭代因子。

2h)利用步骤2a)到步骤2g)得到的参数，计算扩展目标状态和量测噪声协方差的联合概率假设密度v_k|k-1(x,R)：

v_k|k-1(x,R)＝v_S,k|k-1(x,R)+b_k|k-1(x,R)+γ_k(x,R)。

步骤3，在k≥1时，对势分布P_k-1|k-1(num)进行预测，得到预测势分布P_k|k-1(num)：

$P_{k|k-1}\left(\mathrm{num}\right)=\underset{j=0}{\overset{\mathrm{num}}{\mathrm{\Σ}}}{P}_{\mathrm{birth},k}(\mathrm{num}-j)\underset{h=j}{\overset{\∞}{\mathrm{\Σ}}}\frac{h!}{j!(h-j)!}{P}_{k-1|k-1}\left(h\right)p_{S,k}^j{\left({1-p}_{S,k}\right)}^{h-j}$

其中，P_birth,k(num-j)表示k-1时刻到k时刻新生num-j个扩展目标的概率，p_S,k表示k时刻扩展目标的存活概率，P_k-1|k-1(h)表示k-1时刻有h个扩展目标的概率，！表示阶乘，表示k-1时刻到k时刻有j个扩展目标存活的概率，(1-p_S,k)^h-j表示k-1时刻到k时刻有h-j个扩展目标消亡的概率。

步骤4.对扩展目标状态和量测噪声协方差的联合概率假设密度v_k|k-1(x,R)及预测势分布P_k|k-1(num)进行更新：

4a)利用变分贝叶斯VB方法将联合概率假设密度v_k|k-1(x,R)近似为：

式中，Q_x,k|k-1(x)为高斯分布的求和形式，

Q_R,k|k-1(R)为逆伽玛分布的求和形式,

$Q_{R,k|k-1}\left(R\right)=\underset{i=1}{\overset{J_k}{\mathrm{\Σ}}}\left[\underset{l=1}{\overset{d}{\mathrm{\Π}}}\mathrm{IG}({\left({\mathrm{\σ}}_{k|k-1,l}^{\left(i\right)}\right)}^2;{\mathrm{\α}}_{k|k-1,l}^{\left(i\right)},{\mathrm{\β}}_{k|k-1,l}^{\left(i\right)})\right];$

其中，表示第i个高斯分量在第k个时刻的预测权值，i＝1,...,J_k，J_k表示第k个时刻扩展目标高斯分量的个数，N(·)表示高斯分布，为第k个时刻预测得到的第i个高斯分量的均值，为第k个时刻预测得到的第i个高斯分量的协方差；IG(·)表示逆伽玛分布，为第k个时刻预测得到的第i个逆伽玛分量的常量因子，为第k个时刻预测得到的第i个逆伽玛分量的迭代因子，l＝1,...,d，d表示量测噪声协方差的维数；

4b)利用变分贝叶斯期望最大化VBEM方法对扩展目标状态的概率假设密度Q_x,k|k-1(x)中的高斯分量和量测噪声协方差的概率假设密度Q_R,k|k-1(R)中的逆伽玛分量进行迭代更新：

(4b1)设定逆伽玛分量的常量因子和迭代因子 其中l＝1,...,d，d为量测噪声协方差R的维数；

(4b2)根据设定的逆伽玛分量的两个因子，计算量测噪声协方差：

其中n＝1,...,N，N为最大迭代次数，diag[...]表示对角化其中的元素；

(4b3)利用量测噪声协方差计算更新因子

$S_W^{\left(n\right)}=H_W{P}_{k|k-1}^{\left(i\right)}{H_W}^T+R_W^{\left(n\right)}$

其中，表示对量测噪声协方差矩阵进行对角线连接当前单元W的量测数目后的矩阵，blkdiag(·)表示对其中的元素进行对角线连接，|W|表示当前单元W的量测数目；H_W表示对观测矩阵H_k垂直连接当前单元W的量测数目后的矩阵，表示k时刻的观测矩阵H_k的转置；表示k-1到k时刻预测的扩展目标高斯分量运动状态协方差，表示矩阵H_W的转置；

(4b4)利用更新因子计算增益矩阵

$K_k^{\left(i\right)\left(n\right)}=P_{k|k-1}^{\left(i\right)}{H}_W^T{\left[S_W^{\left(n\right)}\right]}^{-1}$

其中[·]^-1表示对矩阵求逆；

(4b5)利用增益矩阵计算扩展目标高斯分量运动状态和扩展目标高斯分量运动状态协方差

$m_{k|k}^{\left(i\right)\left(n\right)}=m_{k|k-1}^{\left(i\right)}+K_k^{\left(i\right)\left(n\right)}(z_W-H_W{m}_{k|k-1}^{\left(i\right)})$

$P_{k|k}^{\left(i\right)\left(n\right)}=[I-K_k^{\left(i\right)\left(n\right)}{H}_W]P_{k|k-1}^{\left(i\right)}$

其中，I表示一个单位矩阵，z_W表示某个划分单元W中的所有量测；

(4b6)提取扩展目标高斯分量运动状态的位置信息，将该位置信息用表示；

(4b7)利用该高斯分量位置信息和量测噪声协方差计算量测Y_n′由位置处的高斯分量生成的概率γ_n′i：

${\mathrm{\γ}}_{n^{\′}i}=\frac{{\mathrm{\π}}_{i}N\left(Y_{n^{\′}}\right|m_{k|k}^{\′\left(i\right)\left(n\right)},R_k^{\left(i\right)\left(n\right)})}{\underset{i=1}{\overset{J_k}{\mathrm{\Σ}}}{\mathrm{\π}}_{i}N\left(Y_{n^{\′}}\right|m_{k|k}^{\′\left(i\right)\left(n\right)},R_k^{\left(i\right)\left(n\right)})},$

其中，J_k表示扩展目标高斯分量数目，Y_n′表示当前单元W的第n′个量测，n′＝1,...,|W|；N(·)表示高斯分布；π_i代表混合系数，N_i表示由位置处的高斯分量产生的有效量测数目，

(4b8)利用量测Y_n′由位置处的高斯分量生成的概率γ_n′i，迭代更新得到扩展目标高斯分量运动状态的位置信息

$m_{k|k}^{\′\′\left(i\right)\left(n\right)}=\frac{1}{N_i}\underset{n^{\′}=1}{\overset{\left|W\right|}{\mathrm{\Σ}}}{\mathrm{\γ}}_{n^{\′}i}{Y}_{n^{\′}};$

(4b9)利用扩展目标高斯分量运动状态的位置信息混合系数π_i，量测噪声协方差计算最大似然函数L⁽ⁱ⁾⁽ⁿ⁾：

$L^{\left(i\right)\left(n\right)}=\underset{n^{\′}=1}{\overset{\left|W\right|}{\mathrm{\Σ}}}\ln \underset{i=1}{\overset{J_k}{\mathrm{\Σ}}}{\mathrm{\π}}_{i}N\left(Y_{n^{\′}}\right|m_{k|k}^{\′\′\left(i\right)\left(n\right)},R_k^{\left(i\right)\left(n\right)}),$ 其中J_k表示扩展目标高斯分量数目；

(4b10)判断|L⁽ⁱ⁾⁽ⁿ⁾-L^(i)(n-1)|是否小于常量ε＝0.01，同时判断当前迭代次数n是否小于最大迭代次数N＝100，如果是，则停止迭代；否则，返回步骤(4b2)，更新逆伽玛分量迭代因子： ${\mathrm{\β}}_{k,1}^{\left(i\right)(n+1)}=\frac{\mathrm{\Σ}{\mathrm{\β}}_W^{(n+1)}}{d\left|W\right|},$

其中，表示对向量中所有元素进行相加，

${\mathrm{\β}}_W^{(n+1)}={\mathrm{\β}}_W^{\left(n\right)}+\frac{1}{2}{(z_W-H_W{m}_{k|k}^{\left(i\right)\left(n\right)})}_j^2+\frac{1}{2}{\left(H_W{P}_{k|k}^{\left(i\right)\left(n\right)}{H}_W^T\right)}_{\mathrm{jj}},{(\·)}_j^2$ 表示对向量的第j维元素的平方，(·)_jj表示取矩阵的对角线元素，表示对迭代因子垂直连接当前单元W的量测数目后的向量，z_W表示当前单元的量测；

(4b11)提取扩展目标状态分量扩展目标运动状态协方差迭代因子即， $m_{k|k}^{\left(i\right)}=m_{k|k}^{\left(i\right)\left(n\right)},P_{k|k}^{\left(i\right)}=P_{k|k}^{\left(i\right)\left(n\right)},{\mathrm{\β}}_{k,l}^{\left(i\right)}={\mathrm{\β}}_{k,l}^{\left(i\right)\left(n\right)},$ 其中扩展目标状态分量中的位置信息用的是步骤(4b8)中迭代更新得到的扩展目标分量运动状态的位置信息

4c)对势分布P_k|k-1(num)进行更新，得到更新后的势分布P_k|k(num)表示如下：

$P_{k|k}\left(\mathrm{num}\right)=\left\{\begin{array}{cc}\frac{{\mathrm{\Σ}}_{p\∠Z}{\mathrm{\Σ}}_{W\∈p}{\mathrm{\ψ}}_{p,W}{G}_{k|k-1}^{\left(\mathrm{num}\right)}\left(0\right)\left(\begin{array}{c}{G}_{\mathrm{FA}}\left(0\right)\frac{\mathrm{\ηW}}{\left|p\right|}\frac{{\mathrm{\ρ}}^{\mathrm{num}-\left|P\right|}}{(\mathrm{num}-|p\left|\right)!}{\mathrm{\δ}}_{\mathrm{num}\≥\left|p\right|}\\ +G_{\mathrm{FA}}^{\left(\right|W\left|\right)}\left(0\right)\frac{{\mathrm{\ρ}}^{\mathrm{num}-p+1}}{(\mathrm{num}-|p|+1)!}{\mathrm{\δ}}_{\mathrm{num}\≥\left|p\right|-1}\end{array}\right)}{{\mathrm{\Σ}}_{p\∠Z}{\mathrm{\Σ}}_{W\∈p}{\mathrm{\Ψ}}_{p,W}{l}_{p,W}},& \left|Z\right|\≠0\\ \frac{{\mathrm{\ρ}}^{\mathrm{num}}{G}_{k|k-1}^{\left(\mathrm{num}\right)}\left(0\right)}{G_{k|k-1}\left(\mathrm{\ρ}\right)}& \left|Z\right|=0\end{array}\right.$

其中，p∠Z表示把量测集Z划分成p个非空子集，W∈p表示第p个非空子集下的某一个单元，G_k|k-1(ρ)表示状态预测概率生成函数，表示状态预测概率生成函数的num阶偏导，G_FA(0)表示没有量测时的虚警概率生成函数，ηW表示扩展目标产生量测概率，|p|表示第p个划分下的所有非空单元数目，表示虚警概率生成函数的|W|阶偏导，δ_num≥|p|表示当目标数目num大于划分单元|p|时取值为1，否则为0，|Z|＝0表示扩展目标没有产生量测，|W|表示各非空单元W中的量测个数，l_p,W表示量测划分单元为|p|-1时的虚警常系数，ψ_p,W表示目标产生量测概率的乘积，ρ表示扩展目标分量没有被检测到的概率：

$\mathrm{\ρ}=\mathrm{\Σ}{\stackrel{\‾}{\mathrm{\ω}}}_{k|k-1}^j[1-P_D(\·)+P_D(\·)G_z\left(0\right|\·)]$

$\mathrm{\ηW}=p_{k|k-1}\left[P_D(\·)G_z^{\left(\right|W\left|\right)}\left(0\right|\·)\underset{z^{\′}\∈W}{\mathrm{\Π}}\frac{p_z\left(z^{\′}\right|\·)}{p_{\mathrm{FA}}\left(z^{\′}\right)}\right]$

$l_{p,W}=G_{\mathrm{FA}}\left(0\right)G_{k|k-1}^{\left(\right|p\left|\right)}\left(\mathrm{\ρ}\right)\frac{\mathrm{\ηW}}{\left|p\right|}+G_{\mathrm{FA}}^{\left(\right|W\left|\right)}\left(0\right)G_{k|k-1}^{\left(\right|p|-1)}\left(\mathrm{\ρ}\right)$

${\mathrm{\ψ}}_{p,W}=\underset{W^{\′}\∈p-W}{\mathrm{\Π}}{\mathrm{\ηW}}^{\′}$

其中，Π为连乘符号，表示第j个高斯逆伽玛分量在当前所有高斯逆伽玛分量中所占的比重，p_k|k-1表示单扩展目标状态转移概率密度函数，p_z(z′|·)表示扩展目标量测似然，p_FA(z′)表示虚警量测似然，G_z(0|·)表示量测概率生成函数，表示量测概率生成函数的|W|阶偏导，z′∈W表示量测z′属于单元W，P_D(·)表示检测概率，W′∈p-W表示p划分下所有单元中除去单元W后剩下的单元，ηW′表示扩展目标虚警量测产生概率，表示状态预测概率生成函数的|p|阶偏导，表示状态预测概率生成函数的|p|-1阶偏导。

步骤5.对更新后的高斯分量和逆伽玛分量进行修剪与合并，其步骤如下：

(5a)设置两个修剪门限T1和T2，一个合并门限U：T1＝10^-5，T2＝120，U＝10；设置最大高斯逆伽玛分量数目：J_max＝100；

(5b)计算每个扩展目标分量对应的量测噪声协方差:

(5c)设变量l′＝0，对更新后的扩展目标分量进行修剪，得到修剪后的扩展目标分量对应的序号集合I为： $I=\{i=1,...,J_k|w_k^{\left(i\right)}>T1,{\left|\right|R_k^{\left(i\right)}\left|\right|}_2<T2\};$

(5d)令l′＝l′+1，取表示取最大权值所对应的集合I中元素i′，对修剪后的扩展目标分量中满足合并门限U的分量进行提取，得到适合合并的扩展目标分量对应的序号集合为： $\stackrel{\&OverBar;}{Q}=\{i\&Element;I|{(m_k^{\left(i\right)}-m_k^{\left(i^{\&prime;}\right)})}^T{\left(P_k^{\left(i\right)}\right)}^{-1}(m_k^{\left(i\right)}-m_k^{\left(i^{\&prime;}\right)})\&le;U\};$

(5e)分别对序号集合中对应的扩展目标分量的权值运动状态常量因子迭代因子协方差进行合并，得到合并后的扩展目标分量的权值运动状态常量因子迭代因子协方差如下：

${\tilde{w}}_k^{\left(l^{\&prime;}\right)}=\underset{i\&Element;\stackrel{\&OverBar;}{Q}}{\mathrm{\&Sigma;}}{w}_k^{\left(i\right)},$

${\tilde{m}}_k^{\left(l^{\&prime;}\right)}=\frac{1}{{\tilde{w}}_k^{\left(l^{\&prime;}\right)}}\underset{i\&Element;\stackrel{\&OverBar;}{Q}}{\mathrm{\&Sigma;}}{w}_k^{\left(i\right)}{m}_k^{\left(i\right)},$

${\widetilde{\mathrm{\&alpha;}}}_k^{\left(l^{\&prime;}\right)}=\frac{1}{{\tilde{w}}_k^{\left(l^{\&prime;}\right)}}\underset{i\&Element;\stackrel{\&OverBar;}{Q}}{\mathrm{\&Sigma;}}{w}_k^{\left(i\right)}{\mathrm{\&alpha;}}_k^{\left(i\right)},$

${\widetilde{\mathrm{\&beta;}}}_k^{\left(l^{\&prime;}\right)}=\frac{1}{{\tilde{w}}_k^{\left(l^{\&prime;}\right)}}\underset{i\&Element;\stackrel{\&OverBar;}{Q}}{\mathrm{\&Sigma;}}{w}_k^{\left(i\right)}{\mathrm{\&beta;}}_k^{\left(i\right)},$

${\tilde{P}}_k^{\left(l^{\&prime;}\right)}=\frac{1}{{\tilde{w}}_k^{\left(l^{\&prime;}\right)}}\underset{i\&Element;\stackrel{\&OverBar;}{Q}}{\mathrm{\&Sigma;}}{w}_k^{\left(i\right)}(P_k^{\left(i\right)}+({\tilde{m}}_k^{\left(l^{\&prime;}\right)}-m_k^{\left(i\right)}){({\tilde{m}}_k^{\left(l^{\&prime;}\right)}-m_k^{\left(i\right)})}^T);$

(5f)将步骤(5c)得到的修剪后的扩展目标分量对应的序号集合I中与步骤(5d)得到的适合合并的扩展目标分量对应的序号集合中相同的元素去除掉，然后判断修剪后的扩展目标分量对应的序号集合I是否为空集，如果不为空集，则返回步骤(5d)，否则执行(5g)；

(5g)判断变量l′是否大于最大高斯逆伽玛分量数目J_max，如果l′＞J_max，则将权值对应的高斯逆伽玛分量按从大到小排列，并取前J_max个权值大于0.5的高斯逆伽玛分量的位置和速度作为扩展目标的状态；如果l′＜J_max，则将所有权值大于0.5对应的高斯逆伽玛分量的位置和速度作为扩展目标的状态。

步骤6.根据步骤4c)更新得到的势分布P_k|k(num)进行加权平均，得到扩展目标的数目： ${\mathrm{num}}_{k|k}=\underset{\mathrm{num}=1}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{\mathrm{num}\&times;p}_{k|k}\left(\mathrm{num}\right);$

得到扩展目标的状态和数目以后即可完成对扩展目标的跟踪。

本发明用于跟踪扩展目标的效果可通过以下仿真实验进一步说明：

1.仿真条件

考虑二维平面中做匀速直线运动4个目标，并且在两个位置出现运动轨迹交叉的情况，采样周期为整个观测过程持续40个时刻，目标的运动方程和量测方程分别为：

x_k＝Fx_k-1+Gw_k

y_k＝Hx_k+v_k

其中，高斯噪声w_k的标准差σ_w＝2。扩展目标状态转移矩阵F，输入矩阵G，高斯噪声w_k，观测矩阵H设置为：

$F=\left[\begin{array}{cccc}1& 0& T& 0\\ 0& 1& 0& T\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right],G=\left[\begin{array}{cc}1/2& 0\\ 0& 1/2\\ 1& 0\\ 0& 1\end{array}\right],w_k~N\left(\begin{array}{cc}0,& \left[\begin{array}{cc}{\mathrm{\&sigma;}}_w^2& 0\\ 0& {\mathrm{\&sigma;}}_w^2\end{array}\right]\end{array}\right),H=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\end{array}\right],$

实际目标跟踪场景中量测噪声v_k的协方差是未知的，本次实验中设真实的量测噪声标准差σ_v＝1；

设传统GM-CPHD算法中采用的三个量测噪声标准差分别为σ＝0.05,1,50；

设目标产生量测数目服从泊松分布，泊松分布参数β＝20，目标产生量测位置服从高斯分布；

设目标存活概率p_S,k＝0.99，检测概率p_D,k＝0.98；

设杂波数服从均值为5的泊松分布，杂波在整个观测区域内服从均匀分布。

设置新生目标的状态为 $m_{\mathrm{\&gamma;}}^{\left(1\right)}={[-\mathrm{100,200,0,0}]}^T,m_{\mathrm{\&gamma;}}^{\left(2\right)}={[0,-\mathrm{100,0,0}]}^T,$ 新生目标的协方差均为P_γ＝diag[10,10,10,10]；

设置修剪门限T1＝10^-5，修剪门限T2＝120，合并门限U＝10，最大高斯逆伽玛分量数目J_max＝100，遗忘因子ρ＝0.9，实验仿真次数为100次。

2.仿真内容与结果

仿真1，采用本发明方法对未知量测噪声协方差下扩展目标运动轨迹的单次实验进行跟踪，结果如图2所示。从图2可以看出，本发明方法可以很好的跟踪扩展目标运动轨迹。

仿真2，采用本发明方法与传统的分别采用不同量测噪声协方差的三个GM-CPHD方法进行100次实验的目标数目估计对比，结果如图3所示。从图3可以看出，本发明方法与传统的给定真实量测噪声协方差下的GM-CPHD方法对目标数目的估计效果一样良好。但是，传统的GM-CPHD方法，当量测噪声协方差未知时，且使用的量测噪声协方差与真实值偏差比较大时，估计的目标数目会发生很大的偏差，同时在k＝20,k＝35这两个目标运动轨迹出现交叉时刻，目标数目出现漏检；

仿真3，采用本发明方法与传统的分别采用不同量测噪声协方差的三个GM-CPHD方法借助OSPA距离进行100次实验的跟踪精度对比，结果如图4所示。从图4可以看出，本发明方法与传统的给定真实量测噪声协方差下的GM-CPHD算法的跟踪精度一样良好。但是，传统的GM-CPHD方法，当使用的量测噪声协方差与真实的量测噪声协方差偏差比较大时，跟踪精度会变得很差。

实验表明，本发明方法在处理量测噪声协方差是未知条件下的扩展目标跟踪场景时，其跟踪效果优于传统的GM-CPHD扩展目标跟踪方法。

Claims (5)

1.一种基于变分贝叶斯期望最大化的扩展目标跟踪方法，包括如下步骤：

(1)在时刻k＝0时，初始化扩展目标状态和量测噪声协方差的联合概率假设密度为v₀(x,R)；

(2)在k≥1时，对扩展目标状态和量测噪声协方差的联合概率假设密度v_k-1|k-1(x,R)及用于计算扩展目标数目的势分布P_k-1|k-1(num)进行预测，得到预测的扩展目标联合概率假设密度v_k|k-1(x,R)和预测势分布P_k|k-1(num)；

(3)对预测的扩展目标状态和量测噪声协方差的联合概率假设密度v_k|k-1(x,R)及用于计算扩展目标数目的势分布P_k|k-1(num)进行更新：

(3a)利用变分贝叶斯VB方法对联合概率假设密度v_k|k-1(x,R)进行近似，得到用高斯分布的求和形式表示的扩展目标状态的概率假设密度Q_x,k|k-1(x)和用逆伽玛分布的求和形式表示的量测噪声协方差的概率假设密度Q_R,k|k-1(R)；

(3b)利用变分贝叶斯期望最大化VBEM方法对扩展目标状态的概率假设密度Q_x,k|k-1(x)中的高斯分量和量测噪声协方差的概率假设密度Q_R,k|k-1(R)中的逆伽玛分量进行迭代更新，得到表示扩展目标运动状态x的高斯分量和表示量测噪声协方差R的逆伽玛分量；

(3c)对步骤(2)预测得到的势分布P_k|k-1(num)进行更新，得到更新后的势分布P_k|k(num)；

(4)对更新后的高斯分量和逆伽玛分量进行修剪与合并，并提取合并后的高斯分量和逆伽玛分量的位置和速度作为扩展目标的状态；

(5)对步骤(3)更新得到的势分布P_k|k(num)进行加权平均，得到扩展目标的数目： ${\mathrm{num}}_{k|k}=\underset{\mathrm{num}=1}{\overset{\&infin;}{\mathrm{\&Sigma;}}}\mathrm{num}\&times;p_{k|k}\left(\mathrm{num}\right);$

(6)重复步骤(2)-(5)，继续跟踪扩展目标。

2.根据权利要求1所述的基于变分贝叶斯期望最大化的扩展目标跟踪方法，其中，步骤(3a)所述的利用变分贝叶斯VB方法对联合概率假设密度v_k|k-1(x,R)进行近似，按如下公式进行：

$Q_{x,k|k-1}\left(x\right)=\underset{i=1}{\overset{J_k}{\mathrm{\&Sigma;}}}\left[w_{k|k-1}^{\left(i\right)}N(x;m_{k|k-1}^{\left(i\right)},P_{k|k-1}^{\left(i\right)})\right]$

$Q_{R,k|k-1}\left(R\right)=\underset{l=1}{\overset{J_k}{\mathrm{\&Sigma;}}}\left[\underset{l=1}{\overset{d}{\mathrm{\&Pi;}}}\mathrm{IG}({\left({\mathrm{\&sigma;}}_{k|k-1,l}^{\left(i\right)}\right)}^2;{\mathrm{\&alpha;}}_{k|k-1,l}^{\left(i\right)},{\mathrm{\&beta;}}_{k|k-1,l}^{\left(i\right)})\right]$

其中，Q_x,k|k-1(x)为高斯分布的求和形式，表示为

Q_R,k|k-1(R)为逆伽玛分布的求和形式,表示为

$Q_{R,k|k-1}\left(R\right)=\underset{l=1}{\overset{J_k}{\mathrm{\&Sigma;}}}\left[\underset{l=1}{\overset{d}{\mathrm{\&Pi;}}}\mathrm{IG}({\left({\mathrm{\&sigma;}}_{k|k-1,l}^{\left(i\right)}\right)}^2;{\mathrm{\&alpha;}}_{k|k-1,l}^{\left(i\right)},{\mathrm{\&beta;}}_{k|k-1,l}^{\left(i\right)})\right];$

表示第i个高斯分量在第k个时刻的权值，i＝1,...,J_k，J_k表示第k个时刻扩展目标高斯分量的个数，N(·)表示高斯分布，为第k个时刻预测得到的第i个高斯分量的均值，为第k个时刻预测得到的第i个高斯分量的协方差；IG(·)表示逆伽玛分布，为第k个时刻预测得到的第i个逆伽玛分量的常量因子，为第k个时刻预测得到的第i个逆伽玛分量的迭代因子，l＝1,...,d，d表示量测噪声协方差的维数。

3.根据权利要求1所述的基于变分贝叶斯期望最大化的扩展目标跟踪方法，其中，步骤(3b)所述的利用变分贝叶斯期望最大化VBEM方法对扩展目标状态的概率假设密度Q_x,k|k-1(x)中的高斯分量和量测噪声协方差的概率假设密度Q_R,k|k-1(R)中的逆伽玛分量进行迭代更新，按如下步骤进行：

(3b1)设定逆伽玛分量的常量因子和迭代因子 其中l＝1,...,d，d为量测噪声协方差R的维数；

(3b2)根据设定的逆伽玛分量两个因子，计算量测噪声协方差：其中n＝1,...,N，N为最大迭代次数，diag[...]表示对角化其中的元素；

(3b3)利用量测噪声协方差计算更新因子

$S_W^{\left(n\right)}=H_W{P}_{k|k-1}^{\left(i\right)}{H_W}^T+R_W^{\left(n\right)}$

其中，表示对量测噪声协方差矩阵进行对角线连接当前单元W量测数目后的矩阵，blkdiag(·)表示对其中的元素进行对角线连接，|W|表示当前单元W的量测数目；H_W表示对观测矩阵H_k垂直连接当前单元W的量测数目后的矩阵， 表示k时刻的观测矩阵H_k的转置；表示k-1到k时刻预测的扩展目标高斯分量运动状态协方差，表示矩阵H_W的转置；

(3b4)利用更新因子计算增益矩阵

$K_k^{\left(i\right)\left(n\right)}=P_{k|k-1}^{\left(i\right)}{H}_W^T{\left[S_W^{\left(n\right)}\right]}^{-1},$

其中[·]^-1表示对矩阵求逆；

(3b5)利用增益矩阵计算扩展目标高斯分量运动状态和扩展目标高斯分量运动状态协方差

$m_{k|k}^{\left(i\right)\left(n\right)}=m_{k|k-1}^{\left(i\right)}+K_k^{\left(i\right)\left(n\right)}(z_w-H_w{m}_{k|k-1}^{\left(i\right)})$

$P_{k|k}^{\left(i\right)\left(n\right)}=[I-K_k^{\left(i\right)\left(n\right)}{H}_W]P_{k|k-1}^{\left(i\right)}$

其中，I表示一个单位矩阵，z_W表示某个划分单元W中的所有量测；

(3b6)提取扩展目标高斯分量运动状态的位置信息，将该位置信息用表示；

(3b7)利用该高斯分量位置信息和量测噪声协方差计算量测Y_n′是由位置处的高斯分量生成的概率γ_n′i：

${\mathrm{\&gamma;}}_{n^{\&prime;}i}=\frac{{\mathrm{\&pi;}}_{i}N\left(Y_{n^{\&prime;}}\right|m_{k|k}^{\&prime;\left(i\right)\left(n\right)},R_k^{\left(i\right)\left(n\right)})}{\underset{i=1}{\overset{J_k}{\mathrm{\&Sigma;}}}{\mathrm{\&pi;}}_{i}N\left(Y_{n^{\&prime;}}\right|m_{k|k}^{\&prime;\left(i\right)\left(n\right)},R_k^{\left(i\right)\left(n\right)})}$

其中，J_k表示扩展目标高斯分量数目，Y_n′表示当前单元W的第n′个量测，n′＝1,...,|W|；N(·)表示高斯分布；π_i代表混合系数，N_i表示由位置处的高斯分量产生的有效量测数目，

(3b8)利用量测Y_n′由位置处的高斯分量生成的概率γ_n′i，迭代更新得到扩展目标高斯分量运动状态的位置信息

$m_{k|k}^{\&prime;\&prime;\left(i\right)\left(n\right)}=\frac{1}{N_i}\underset{n^{\&prime;}=1}{\overset{\left|W\right|}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_{n^{\&prime;}i}{Y}_{n^{\&prime;}};$

(3b9)利用扩展目标高斯分量运动状态的位置信息混合系数π_i，量测噪声协方差计算最大似然函数L⁽ⁱ⁾⁽ⁿ⁾：

$L^{\left(i\right)\left(n\right)}=\underset{n^{\&prime;}=1}{\overset{\left|W\right|}{\mathrm{\&Sigma;}}}1n\underset{i=1}{\overset{J_k}{\mathrm{\&Sigma;}}}{\mathrm{\&pi;}}_{i}N\left(Y_{n^{\&prime;}}\right|m_{k|k}^{\&prime;\&prime;\left(i\right)\left(n\right)},R_k^{\left(i\right)\left(n\right)}),$ 其中J_k表示扩展目标高斯分量数目；

(3b10)判断|L⁽ⁱ⁾⁽ⁿ⁾-L^(i)(n-1)|是否小于常量ε＝0.01，同时判断当前迭代次数n是否小于最大迭代次数N＝100，如果是，则停止迭代，否则返回步骤(3b2)，更新逆伽玛分量迭代因子：

其中，表示对向量中所有元素进行相加，

${\mathrm{\&beta;}}_W^{(n+1)}={\mathrm{\&beta;}}_W^{\left(n\right)}+\frac{1}{2}{(z_W-H_W{m}_{k|k}^{\left(i\right)\left(n\right)})}_j^2+\frac{1}{2}{\left(H_W{P}_{k|k}^{\left(i\right)\left(n\right)}{H}_W^T\right)}_{\mathrm{jj}},$ 表示对向量的第j维元素的平方，(·)_jj表示取矩阵的对角线元素，表示对迭代因子垂直连接当前单元W量测数目后的向量，z_W表示当前单元的量测；

(3b11)提取扩展目标状态分量扩展目标运动状态协方差迭代因子即， $m_{k|k}^{\left(i\right)}=m_{k|k}^{\left(i\right)\left(n\right)},P_{k|k}^{\left(i\right)}=P_{k|k}^{\left(i\right)\left(n\right)},{\mathrm{\&beta;}}_{k.l}^{\left(i\right)}={\mathrm{\&beta;}}_{k,l}^{\left(i\right)\left(n\right)},$ 其中扩展目标状态分量中的位置信息用的是步骤(3b8)中迭代更新得到的扩展目标分量运动状态的位置信息

4.根据权利要求1所述的基于变分贝叶斯期望最大化的扩展目标跟踪方法，其中，所述步骤(3c)对势分布P_k|k-1(num)进行更新，得到更新后的势分布P_k|k(num)表示如下：

$P_{k|k}\left(\mathrm{num}\right)=\left\{\begin{array}{cc}\frac{{\mathrm{\&Sigma;}}_{p\&angle;Z}{\mathrm{\&Sigma;}}_{W\&Element;p}{\mathrm{\&psi;}}_{p,W}{G}_{k|k-1}^{\left(\mathrm{num}\right)}\left(0\right)\left(\begin{array}{cc}{G}_{\mathrm{FA}}\left(0\right)\frac{\mathrm{\&eta;W}}{\left|p\right|}\frac{{\mathrm{\&rho;}}^{\mathrm{num}-\left|P\right|}}{(\mathrm{num}-|p\left|\right)!}& {\mathrm{\&delta;}}_{\mathrm{num}\&GreaterEqual;\left|p\right|}\\ +G_{\mathrm{FA}}^{\left(\right|W\left|\right)}\left(0\right)\frac{{\mathrm{\&rho;}}^{\mathrm{num}-p+1}}{(\mathrm{num}-|p|+1)!}& {\mathrm{\&delta;}}_{\mathrm{num}\&GreaterEqual;\left|p\right|-1}\end{array}\right)}{{\mathrm{\&Sigma;}}_{p\&angle;Z}{\mathrm{\&Sigma;}}_{W\&Element;p}{\mathrm{\&psi;}}_{p,W}{l}_{p,W}},& \left|Z\right|\&NotEqual;0\\ \frac{{\mathrm{\&rho;}}^{\mathrm{num}}{G}_{k|k-1}^{\left(\mathrm{num}\right)}\left(0\right)}{G_{k|k-1}\left(\mathrm{\&rho;}\right)}& \left|Z\right|=0\end{array}\right.$

其中，p∠Z表示把量测集Z划分成p个非空子集，W∈p表示第p个非空子集下的某一个单元，G_k|k-1(ρ)表示状态预测概率生成函数，表示状态预测概率生成函数的num阶偏导，G_FA(0)表示没有量测时的虚警概率生成函数，ηW表示扩展目标产生量测概率，|p|表示第p个划分下的所有非空单元数目，表示虚警概率生成函数的|W|阶偏导，δ_num≥p表示当目标数目num大于划分单元|p|时取值为1，否则为0，|Z|＝0表示扩展目标没有产生量测，|W|表示各非空单元W中的量测个数，l_p,W表示量测划分单元为|p|-1时的虚警常系数，ψ_p,W表示目标产生量测概率的乘积，ρ表示扩展目标分量没有被检测到的概率：

$\mathrm{\&rho;}=\mathrm{\&Sigma;}{\stackrel{\&OverBar;}{\mathrm{\&omega;}}}_{k|k-1}^j[1-P_D(\&CenterDot;)+P_D(\&CenterDot;)G_z\left(0\right|\&CenterDot;)]$

$\mathrm{\&eta;W}=p_{k|k-1}\left[P_D(\&CenterDot;)G_z^{\left(\right|W\left|\right)}\left(0\right|\&CenterDot;)\underset{z^{\&prime;}\&Element;W}{\mathrm{\&Pi;}}\frac{p_z\left(z^{\&prime;}\right|\&CenterDot;)}{p_{\mathrm{FA}}\left(z^{\&prime;}\right)}\right]$

$l_{p,W}=G_{\mathrm{FA}}\left(0\right)G_{k|k-1}^{\left(\right|p\left|\right)}\left(\mathrm{\&rho;}\right)\frac{\mathrm{\&eta;W}}{\left|p\right|}+G_{\mathrm{FA}}^{\left(\right|W\left|\right)}\left(0\right)G_{k|k-1}^{\left(\right|p|-1)}\left(\mathrm{\&rho;}\right)$

${\mathrm{\&psi;}}_{p,W}=\underset{W^{\&prime;}\&Element;p-W}{\mathrm{\&Pi;}}\mathrm{\&eta;}{W}^{\&prime;}$

其中，Π为连乘符号，表示第j个高斯逆伽玛分量在当前所有高斯逆伽玛分量中所占的比重，p_k|k-1表示单扩展目标状态转移概率密度函数，p_z(z′|·)表示扩展目标量测似然，p_FA(z′)表示虚警量测似然，G_z(0|·)表示量测概率生成函数，表示量测概率生成函数的|W|阶偏导，z′∈W表示量测z′属于单元W，P_D(·)表示检测概率，W′∈p-W表示p划分下所有单元中除去单元W后剩下的单元，ηW′表示扩展目标虚警量测产生概率，表示状态预测概率生成函数的|p|阶偏导，表示状态预测概率生成函数的|p|-1阶偏导。

5.根据权利要求1所述的基于变分贝叶斯期望最大化的扩展目标跟踪方法，其中，步骤(4)所述的对更新后的高斯分量和逆伽玛分量进行修剪与合并，按如下步骤进行：

(5a)设置两个修剪门限T1和T2，一个合并门限U：T1＝10^-5，T2＝120，U＝10；设置最大高斯逆伽玛分量数目：J_max＝100；

(5b)计算每个扩展目标分量对应的量测噪声协方差:

(5c)设变量l′＝0，对更新后的扩展目标分量进行修剪，得到修剪后的扩展目标分量对应的序号集合I为： $I=\{i=1,...,J_k|w_k^{\left(i\right)}>T1,{\left|\right|R_k^{\left(i\right)}\left|\right|}_2<T2\};$

(5d)令l′＝l′+1，取表示取最大权值所对应的集合I中元素i′，对修剪后的扩展目标分量中满足合并门限U的分量进行提取，得到适合合并的扩展目标分量对应的序号集合为：

$\stackrel{\&OverBar;}{Q}=\{i\&Element;I|{(m_k^{\left(i\right)}-m_k^{\left(i^{\&prime;}\right)})}^T{\left(P_k^{\left(i\right)}\right)}^{-1}(m_k^{\left(i\right)}-m_k^{\left(i^{\&prime;}\right)})\&le;U\};$

(5e)分别对序号集合中对应的扩展目标分量的权值运动状态常量因子迭代因子协方差进行合并，得到合并后的扩展目标分量的权值运动状态常量因子迭代因子协方差如下：

${\tilde{w}}_k^{\left(l^{\&prime;}\right)}=\underset{i\&Element;\stackrel{\&OverBar;}{Q}}{\mathrm{\&Sigma;}}{w}_k^{\left(i\right),}$

${\tilde{m}}_k^{\left(l^{\&prime;}\right)}=\frac{1}{{\tilde{w}}_k^{\left(l^{\&prime;}\right)}}\underset{i\&Element;\stackrel{\&OverBar;}{Q}}{\mathrm{\&Sigma;}}{w}_k^{\left(i\right)}{m}_k^{\left(i\right)},$

${\widetilde{\mathrm{\&alpha;}}}_k^{\left(l^{\&prime;}\right)}=\frac{1}{{\tilde{w}}_k^{\left(l^{\&prime;}\right)}}\underset{i\&Element;\stackrel{\&OverBar;}{Q}}{\mathrm{\&Sigma;}}{w}_k^{\left(i\right)}{\mathrm{\&alpha;}}_k^{\left(i\right)},$

${\widetilde{\mathrm{\&beta;}}}_k^{\left(l^{\&prime;}\right)}=\frac{1}{{\tilde{w}}_k^{\left(l^{\&prime;}\right)}}\underset{i\&Element;\stackrel{\&OverBar;}{Q}}{\mathrm{\&Sigma;}}{w}_k^{\left(i\right)}{\mathrm{\&beta;}}_k^{\left(i\right)},$

${\tilde{P}}_k^{\left(l^{\&prime;}\right)}=\frac{1}{{\tilde{w}}_k^{\left(l^{\&prime;}\right)}}\underset{i\&Element;\stackrel{\&OverBar;}{Q}}{\mathrm{\&Sigma;}}{w}_k^{\left(i\right)}(P_k^{\left(i\right)}+({\tilde{m}}_k^{\left(l^{\&prime;}\right)}-m_k^{\left(i\right)}){({\tilde{m}}_k^{\left(l^{\&prime;}\right)}-m_k^{\left(i\right)})}^T);$

(5f)将步骤(5c)得到的修剪后的扩展目标分量对应的序号集合I中与步骤(5d)得到的适合合并的扩展目标分量对应的序号集合中相同的元素去除掉，然后判断修剪后的扩展目标分量对应的序号集合I是否为空集，如果不为空集则返回步骤(5d)，否则执行(5g)；

(5g)判断变量l′是否大于最大高斯逆伽玛分量数目J_max，如果l′＞J_max，则将权值对应的高斯逆伽玛分量按从大到小排列，并取前J_max个权值大于0.5的高斯逆伽玛分量的位置和速度作为扩展目标的状态；如果l′＜J_max，则将所有权值大于0.5对应的高斯逆伽玛分量的位置和速度作为扩展目标的状态。