CN104794735A

CN104794735A - Extended target tracking method based on variational Bayesian expectation maximization

Info

Publication number: CN104794735A
Application number: CN201510152626.9A
Authority: CN
Inventors: 李翠芸; 王晋斌; 姬红兵; 王荣
Original assignee: Xidian University
Current assignee: Xidian University
Priority date: 2015-04-02
Filing date: 2015-04-02
Publication date: 2015-07-22
Anticipated expiration: 2035-04-02
Also published as: CN104794735B

Abstract

The invention discloses an extended target tracking method based on Variational Bayesian Expectation Maximization (VBEM), which mainly solves the problem that in the traditional extended target tracking field, when the covariance of measurement noise is unknown, the tracking performance of the target will decrease. problem of a sharp decline. This method first predicts the relevant parameters of the Gaussian inverse gamma component in the joint probability hypothesis density of the target state and measurement noise covariance; then updates the parameters of the Gaussian inverse gamma component; finally obtains the extended target state and number by pruning and merging . Simulation experiments show that the present invention can well track multi-spread targets under unknown number and unknown measurement noise covariance, and has high tracking precision, and can be used for tracking aircraft and submarine targets.

Description

Extended Object Tracking Method Based on Variational Bayesian Expectation Maximization

技术领域technical field

本发明属于信息处理技术领域，特别涉及一种目标跟踪方法,可用于跟踪多扩展目标。The invention belongs to the technical field of information processing, and in particular relates to a target tracking method, which can be used to track multiple extended targets.

背景技术Background technique

在传统的目标跟踪领域里，由于雷达的分辨率有限，因此通常情况下将目标看作是点目标，即一个时刻每个目标只能产生一个量测。近年来，随着雷达探测技术的发展和实际应用的需要，更多的将目标视为扩展目标，即每个目标在每一个时刻都可以产生多个量测。In the traditional target tracking field, due to the limited resolution of the radar, the target is usually regarded as a point target, that is, each target can only produce one measurement at a time. In recent years, with the development of radar detection technology and the needs of practical applications, more targets are regarded as extended targets, that is, each target can generate multiple measurements at each moment.

实际的目标跟踪场景中，目标的数目是无法提前预知的，因此随机集理论的提出极大的满足了目标跟踪理论的需要。在诸多对目标的模型假设中，尤其属扩展目标理论的提出更加贴近目前跟踪理论的需要，并且在实际生活中得到了广泛的应用，近十年来，成为了目标跟踪领域中的研究热点。2003年，Mahler将随机有限集理论应用于多目标跟踪问题，提出了概率假设密度PHD滤波。2005年，Gilholm和Salmond提出一种空间分布服从泊松分布的扩展目标模型。2009年，Mahler推导出了扩展目标PHD滤波，即用每一时刻的量测随机集对目标随机集进行预测、更新，继而可以准确提取目标的运动状态和估计目标的数目。2010年，等给出了扩展目标PHD的高斯混合实现形式。2011年，Orguner等又提出了带势分布的扩展目标PHD(ET-CPHD)滤波，很好的解决了ET-PHD估计目标数目时的缺陷。然而，传统的扩展目标跟踪算法处理的都是量测噪声协方差是已知的情形，实际中当量测噪声协方差未知时，扩展目标的跟踪性能将会急剧下降。In the actual target tracking scene, the number of targets cannot be predicted in advance, so the proposal of random set theory greatly meets the needs of target tracking theory. Among many model assumptions on the target, especially the proposed extended target theory is closer to the needs of the current tracking theory, and has been widely used in real life. In the past ten years, it has become a research hotspot in the field of target tracking. In 2003, Mahler applied stochastic finite set theory to the problem of multi-target tracking, and proposed the probability hypothesis density PHD filter. In 2005, Gilholm and Salmond proposed an extended objective model in which the spatial distribution obeys the Poisson distribution. In 2009, Mahler derived the extended target PHD filter, that is, the random set of measurements at each moment is used to predict and update the target random set, and then the motion state of the target and the number of estimated targets can be accurately extracted. year 2010, etc. gave the Gaussian mixture realization form of extended objective PHD. In 2011, Orguner et al. proposed the extended target PHD (ET-CPHD) filter with potential distribution, which solved the defects of ET-PHD in estimating the number of targets. However, the traditional extended target tracking algorithms deal with the situation where the measurement noise covariance is known. In practice, when the measurement noise covariance is unknown, the tracking performance of the extended target will drop sharply.

发明内容Contents of the invention

本发明的目的在于针对上述问题，提出一种基于变分贝叶斯期望最大化的扩展目标跟踪方法，以提高在量测噪声协方差未知条件下的跟踪性能。The object of the present invention is to address the above problems, and propose an extended target tracking method based on variational Bayesian expectation maximization, so as to improve the tracking performance under the condition of unknown measurement noise covariance.

实现本发明的技术关键是：在势概率假设密度滤波框架下，引入变分贝叶斯期望最大化VBEM技术，估计目标状态和量测噪声协方差的联合概率假设密度，实现未知量测噪声协方差下的目标跟踪问题。The technical key to realize the present invention is: under the potential probability hypothesis density filtering framework, introduce the variational Bayesian expectation maximization VBEM technology, estimate the joint probability hypothesis density of the target state and the measurement noise covariance, and realize the unknown measurement noise covariance The problem of object tracking under variance.

VB是一类用于贝叶斯估计和机器学习领域中近似计算复杂积分的方法，本文中VB被用来近似线性高斯系统中，扩展目标状态和量测噪声协方差的联合概率假设密度，该算法的主要思想是对扩展目标状态和量测噪声协方差的联合概率假设密度进行参数化近似,并给出其参数化表达形式。在用VB理论对扩展目标状态和量测噪声协方差的联合概率假设密度进行近似的过程中，为了判断所估计的高斯逆伽玛分量相关参数的性能，在VB的基础上又引入了期望最大化EM算法。在期望E步，估计未知参数的期望值，给出当前的参数估计；在最大化M步，重新估计分布参数，以使该似然函数达到最大。M步上得到的参数估计值被用于下一个E步计算，这个过程不断交替进行。这种方法可以广泛的应用于数据有缺损的情况，且具有简单性和稳定性等优点。VB is a class of methods used in Bayesian estimation and machine learning to approximate complex integrals. In this paper, VB is used to approximate the joint probability hypothesis density of the extended target state and measurement noise covariance in a linear Gaussian system. The main idea of the algorithm is to parametrically approximate the joint probability hypothesis density of the extended target state and measurement noise covariance, and give its parametric expression. In the process of using the VB theory to approximate the joint probability hypothesis density of the expanded target state and the measurement noise covariance, in order to judge the performance of the parameters related to the estimated Gaussian inverse gamma component, the expected maximum Chem algorithm. In the expectation E step, estimate the expected value of the unknown parameter, and give the current parameter estimate; in the maximization M step, re-estimate the distribution parameters to maximize the likelihood function. The parameter estimates obtained in the M step are used for the calculation of the next E step, and this process is carried out alternately. This method can be widely used in the case of data defects, and has the advantages of simplicity and stability.

本发明利用上述VBEM技术进行扩展目标跟踪的技术步骤包括如下：The technical steps that the present invention utilizes above-mentioned VBEM technology to carry out extended target tracking include as follows:

(1)在时刻k＝0时，初始化扩展目标状态和量测噪声协方差的联合概率假设密度为v₀(x,R)；(1) At time k=0, initialize the joint probability assumption density of the extended target state and the measurement noise covariance as 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)；(2) When k ≥ 1, the joint probability assumption density v _k-1|k-1 (x, R) for the extended target state and measurement noise covariance and the potential distribution P _k- for calculating the number of extended targets _1|k-1 (num) is predicted, and the predicted extended target joint probability hypothesis density v _k|k-1 (x, R) and the predicted potential distribution P _k|k-1 (num) are obtained;

(3)对预测的扩展目标状态和量测噪声协方差的联合概率假设密度v_k|k-1(x,R)及用于计算扩展目标数目的势分布P_k|k-1(num)进行更新：(3) The joint probability assumption density v _k|k-1 (x, R) for the predicted extended target state and the measurement noise covariance and the potential distribution P _k|k-1 (num) used to calculate the number of extended targets Make an update:

(3a)利用变分贝叶斯VB方法对联合概率假设密度v_k|k-1(x,R)进行近似，得到用高斯分布的求和形式表示的扩展目标状态的概率假设密度Q_x,k|k-1(x)和用逆伽玛分布的求和形式表示的量测噪声协方差的概率假设密度Q_R,k|k-1(R)；(3a) Use the variational Bayesian VB method to approximate the joint probability hypothesis density v _k|k-1 (x, R), and obtain the probability hypothesis density Q _{x of the extended target state represented by the summation form of the Gaussian distribution, k|k-1} (x) and the probability hypothesis density Q _R,k|k-1 (R) of the measurement noise covariance expressed in the sum form of the inverse gamma distribution;

(3b)利用变分贝叶斯期望最大化VBEM方法对扩展目标状态的概率假设密度Q_x,k|k-1(x)中的高斯分量和量测噪声协方差的概率假设密度Q_R,k|k-1(R)中的逆伽玛分量进行迭代更新，得到表示扩展目标运动状态x的高斯分量和表示量测噪声协方差R的逆伽玛分量；(3b) Using the Variational Bayesian Expectation Maximization VBEM method to expand the probability hypothesis density Q _x,k|k-1 (x) of the target state and the probability hypothesis density Q _{R of the measurement noise covariance in the Gaussian component,} The inverse gamma component in _k|k-1 (R) is iteratively updated to obtain the Gaussian component representing the extended target motion state x and the inverse gamma component representing the measurement noise covariance R;

(3c)对步骤(2)预测得到的势分布P_k|k-1(num)进行更新，得到更新后的势分布P_k|k(num)；(3c) Update the potential distribution P _k|k-1 (num) predicted in step (2) to obtain the updated potential distribution P _k|k (num);

(4)对更新后的高斯分量和逆伽玛分量进行修剪与合并，并提取合并后的高斯分量和逆伽玛分量的位置和速度作为扩展目标的状态；(4) pruning and merging the updated Gaussian component and the inverse gamma component, and extracting the position and velocity of the merged Gaussian component and the inverse gamma component as the state of the expanded target;

(5)对步骤(3)更新得到的势分布P_k|k(num)进行加权平均，得到扩展目标的数目： ${num}_{k | k} = Σ_{num = 1}^{\infty} num \times p_{k | k} (num);$ (5) Perform a weighted average on the potential distribution P _k|k (num) updated in step (3) to obtain the number of expanded targets: ${num}_{k | k} = Σ_{num = 1}^{\infty} num \times p_{k | k} (num);$

(6)重复步骤(2)-(5)，继续跟踪扩展目标。(6) Repeat steps (2)-(5) to continue tracking the extended target.

本发明具有以下优点：The present invention has the following advantages:

第一，本发明引入了变分贝叶斯EM技术，通过估计扩展目标状态和量测噪声协方差的联合概率假设密度，有效估计了各个目标在不同时刻的真实量测噪声，为复杂环境下多扩展目标跟踪场景的分析提供了帮助，保证了CPHD滤波算法可以有效地实现对未知量测噪声协方差环境下扩展目标跟踪。First, the present invention introduces the variational Bayesian EM technology, by estimating the joint probability hypothesis density of the expanded target state and the covariance of the measurement noise, effectively estimating the real measurement noise of each target at different times, providing a solution for complex environments The analysis of multi-extended target tracking scene provides help, which ensures that the CPHD filtering algorithm can effectively realize the extended target tracking in the environment of unknown measurement noise covariance.

第二，本发明由于采用先预测高斯逆伽玛分量的相关参数，再对高斯逆伽玛分量参数进行更新，最后通过修剪与合并得到扩展目标状态和数目的跟踪扩展目标的过程，与传统的GM-CPHD滤波算法相比，提高了跟踪精度。The second, because the present invention adopts the relevant parameter that predicts Gaussian inverse gamma component earlier, then Gaussian inverse gamma component parameter is updated, obtains the process of tracking extended target state and number by pruning and merging at last, and traditional Compared with the GM-CPHD filter algorithm, the tracking accuracy is improved.

附图说明Description of drawings

图1是本发明的总流程图；Fig. 1 is a general flowchart of the present invention;

图2是单次实验条件下，用本发明方法跟踪扩展目标的仿真结果图；Fig. 2 is under single experiment condition, uses the simulation result figure of the inventive method to track extended target;

图3是在100次实验条件下，用本发明方法与传统GM-CPHD方法估计目标数目的仿真结果对比图；Fig. 3 is under 100 experimental conditions, uses the simulation result contrast figure of the method of the present invention and traditional GM-CPHD method to estimate target number;

图4是在100次实验条件下，用本发明方法与传统GM-CPHD方法借助OSPA距离来判断目标跟踪精度的仿真结果对比图。Fig. 4 is a comparison diagram of the simulation results of judging the target tracking accuracy by using the method of the present invention and the traditional GM-CPHD method by means of OSPA distance under the condition of 100 experiments.

具体实施方式detailed description

下面结合附图对本发明做进一步的描述。The present invention will be further described below in conjunction with the accompanying drawings.

参照图1，本发明的具体实现步骤包括如下：With reference to Fig. 1, concrete implementation steps of the present invention include as follows:

步骤1，在时刻k＝0时，初始化目标状态和量测噪声协方差的联合概率假设密度：Step 1, at time k=0, initialize the joint probability hypothesis density of target state and measurement noise covariance:

${v v}_{00} ((x x,, R R)) = = {Σ Σ}_{i i = = 11}^{{J J}_{00}} [[{w w}_{00}^{((i i))} N N ((x x;; {m m}_{00}^{((i i))},, {P P}_{00}^{((i i))})) {Π Π}_{l l = = 11}^{d d} IG IG (({(({σ σ}_{00,, l l}^{((i i))}))}^{22};; {α α}_{00,, l l}^{((i i))},, {β β}_{00,, l l}^{((i i))}))]]$

其中，J₀表示高斯分量的个数,表示第i个高斯分量的权值，N(·)表示高斯分布，表示第i个高斯分量的均值，表示第i个高斯分量的协方差；IG(·)表示逆伽玛分布，表示第i个逆伽玛分量的协方差,为第i个逆伽玛分量的常量因子，为第i个逆伽玛分量的迭代因子，l表示量测噪声协方差的第l维，d表示量测噪声协方差的维数。Among them, J ₀ represents the number of Gaussian components, Represents the weight of the i-th Gaussian component, N( ) represents the Gaussian distribution, Indicates the mean value of the i-th Gaussian component, Represents the covariance of the i-th Gaussian component; IG( ) represents the inverse gamma distribution, Indicates the covariance of the i-th inverse gamma component, is a constant factor of the ith inverse gamma component, is the iteration factor of the i-th inverse gamma component, l represents the l-th dimension of the measurement noise covariance, and d represents the dimension of the measurement noise covariance.

步骤2，在k≥1时，对扩展目标状态和量测噪声协方差的联合概率假设密度v_k-1|k-1(x,R)进行预测，得到预测的扩展目标联合概率假设密度v_k|k-1(x,R)。Step 2. When k≥1, predict the joint probability hypothesis density v _k-1|k-1 (x, R) of the extended target state and measurement noise covariance, and obtain the predicted joint probability hypothesis density of the extended target v _k|k-1 (x,R).

2a)对存活目标状态和量测噪声协方差的联合概率假设密度v_S,k|k-1(x,R)中高斯分量的均值和协方差进行预测，得到预测的存活目标高斯分量的均值和协方差 2a) The mean of the Gaussian component in the joint probability hypothesis density vS _,k|k-1 (x,R) for the surviving target state and the measurement noise covariance and covariance Make a prediction and get the mean of the Gaussian component of the predicted survival target and covariance

${m m}_{S S,, k k | | k k - - 11}^{((i i))} = = {F f}_{S S,, k k - - 11}^{((i i))} {m m}_{S S,, k k - - 11}^{((i i))}$

${P P}_{S S,, k k | | k k - - 11}^{((i i))} = = {Q Q}_{k k - - 11} + + {F f}_{S S,, k k - - 11}^{((i i))} {P P}_{S S,, k k - - 11}^{((i i))} {(({F f}_{S S,, k k - - 11}^{((i i))}))}^{T T}$

其中，表示存活目标状态转移矩阵，Q_k-1表示存活目标过程噪声协方差，(·)^T表示转置；in, Represents the survival target state transition matrix, Q _k-1 represents the noise covariance of the survival target process, ( ) ^T represents the transposition;

2b)对存活目标状态和量测噪声协方差的联合概率假设密度v_S,k|k-1(x,R)中逆伽玛分量的常量因子和迭代因子进行预测，得到预测的存活目标逆伽玛分量的常量因子和逆伽玛分量的迭代因子 2b) Constant factor of the inverse gamma component in the joint probability assumption density vS _,k|k-1 (x,R) for the surviving target state and measurement noise covariance and iteration factor Make predictions to get a constant factor of the inverse gamma component of the predicted survival target and the iteration factor of the inverse gamma component

${α α}_{S S,, k k | | k k - - 11,, l l}^{((i i))} = = {ρ ρ}_{l l} {α α}_{S S,, k k - - 11,, l l}^{((i i))}$

${β β}_{S S,, k k | | k k - - 11,, l l}^{((i i))} = = {ρ ρ}_{l l} {β β}_{S S,, k k - - 11,, l l}^{((i i))}$

其中ρ_l表示遗忘因子，且ρ_l∈(0,1]；Where ρ _l represents the forgetting factor, and ρ _l ∈ (0,1];

2c)在高斯混合框架下，利用预测的存活目标高斯分量的均值协方差和逆伽玛分量的常量因子迭代因子计算存活目标状态和量测噪声协方差的联合概率假设密度v_S,k|k-1(x,R)：2c) Under the Gaussian mixture framework, use the mean of the Gaussian components of the predicted survival targets Covariance and a constant factor of the inverse gamma component iteration factor Compute the joint probability hypothesis density v _S,k|k-1 (x,R) of the surviving target state and the measurement noise covariance:

${v v}_{S S,, k k | | k k - - 11} ((x x,, R R)) = = {P P}_{S S,, k k} {Σ Σ}_{i i = = 11}^{{J J}_{k k - - 11}} [[{w w}_{k k - - 11}^{((i i))} N N ((x x;; {m m}_{S S,, k k | | k k - - 11}^{((i i))},, {P P}_{S S,, k k | | k k - - 11}^{((i i))})) {Π Π}_{l l = = 11}^{d d} IG IG (({(({σ σ}_{S S,, k k | | k k - - 11}^{((i i))}))}^{22};; {α α}_{S S,, k k | | k k - - 11,, l l}^{((i i))},, {β β}_{S S,, k k | | k k - - 11 . . l l}^{((i i))}))]]$

其中，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个逆伽玛分量的迭代因子；Among them, P _S,k represents the survival probability of the extended target, J _k-1 represents the number of Gaussian components at time k-1, Indicates the weight of the i-th Gaussian component at time k-1, N(·) indicates the Gaussian distribution, represents the mean of the i-th Gaussian component in the joint probability hypothesis density v _S,k|k-1 (x,R) of the predicted survival target state and the measurement noise covariance, Represents the covariance of the i-th Gaussian component in the joint probability hypothesis density v _S,k|k-1 (x,R) of the predicted survival target state and the measurement noise covariance; IG(·) represents the inverse gamma distribution, d represents the dimensionality of the measurement noise covariance, represents the constant factor of the ith inverse gamma component in the joint probability hypothesis density vS _,k|k-1 (x,R) of the predicted survival target state and the measurement noise covariance, Represents the iteration factor of the i-th inverse gamma component in the joint probability hypothesis density v _S,k|k-1 (x,R) of the predicted survival target state and the measurement noise covariance;

2d)对衍生目标状态和量测噪声协方差的联合概率假设密度b_k|k-1(x,R)中高斯分量的均值和协方差进行预测，得到预测的衍生目标高斯分量的均值和协方差 2d) The mean of the Gaussian component in the joint probability hypothesis density b _k|k-1 (x,R) for the derived target state and measurement noise covariance and covariance Make a prediction and get the mean of the predicted Gaussian component of the derived target and covariance

${m m}_{b b,, k k | | k k - - 11}^{((i i,, j j))} = = {F f}_{b b,, k k - - 11}^{((j j))} {m m}_{b b,, k k - - 11}^{((i i,, j j))} + + {d d}_{b b,, k k - - 11}^{((j j))}$

${P P}_{b b,, k k | | k k - - 11}^{((i i,, j j))} = = {Q Q}_{b b,, k k - - 11}^{((j j))} + + {F f}_{b b,, k k - - 11}^{((j j))} {P P}_{b b,, k k - - 11}^{((i i))} {(({F f}_{b b,, k k - - 11}^{((j j))}))}^{T T}$

其中，i表示k-1时刻高斯逆伽玛分量的第i个,j表示由k-1时刻到k时刻衍生出的高斯逆伽玛分量的第j个，表示衍生目标的状态转移矩阵，表示衍生目标的状态修正量，表示衍生目标过程噪声协方差；Among them, i represents the i-th Gaussian inverse gamma component at time k-1, and j represents the j-th Gaussian inverse gamma component derived from time k-1 to time k, represents the state transition matrix of the derived target, Indicates the status modifier of the derived target, Denotes the noise covariance of the derived target process;

2e)对衍生目标状态和量测噪声协方差的联合概率假设密度b_k|k-1(x,R)中逆伽玛分量的常量因子和迭代因子进行预测，得到预测的衍生目标逆伽玛分量的常量因子和迭代因子 2e) Constant factor of the inverse gamma component in the joint probability assumption density b _k|k-1 (x,R) for the derived target state and measurement noise covariance and iteration factor Make predictions to get a constant factor of the predicted inverse gamma component of the derived target and iteration factor

${α α}_{b b,, k k | | k k - - 11,, l l}^{((i i,, j j))} = = {ρ ρ}_{l l} {α α}_{b b,, k k - - 11,, l l}^{((i i))}$

${β β}_{b b,, k k | | k k - - 11,, l l}^{((i i,, j j))} = = {ρ ρ}_{l l} {β β}_{b b,, k k - - 11,, l l}^{((i i))};;$

2f)在高斯混合框架下,利用预测的衍生目标高斯分量的均值协方差和逆伽玛分量的常量因子迭代因子计算衍生目标状态和量测噪声协方差的联合概率假设密度b_k|k-1(x,R)：2f) Under the Gaussian mixture framework, use the mean of the predicted Gaussian components of the derived target Covariance and a constant factor of the inverse gamma component iteration factor Compute the joint probability hypothesis density b _k|k-1 (x,R) of the derived target state and measurement noise covariance:

${b b}_{k k | | k k - - 11} ((x x,, R R)) = = {Σ Σ}_{i i = = 11}^{{J J}_{k k - - 11}} {Σ Σ}_{j j = = 11}^{{J J}_{b b,, k k}} [[{w w}_{k k - - 11}^{((i i))} {w w}_{b b,, k k}^{((j j))} N N ((x x;; {m m}_{b b,, k k | | k k - - 11}^{((i i,, j j))},, {P P}_{b b,, k k | | k k - - 11}^{((i i,, j j))})) {Π Π}_{l l = = 11}^{d d} IG IG (({(({σ σ}_{b b,, k k | | k k - - 11,, l l}^{((i i,, j j))}))}^{22};; {α α}_{b b,, k k | | k k - - 11,, l l}^{((i i,, j j))},, {β β}_{b b,, k k | | k k - - 11,, l l}^{((i i,, j j))}))]]$

其中，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个逆伽玛分量的迭代因子；Among them, J _b,k represents the number of Gaussian components of the derived target from time k-1 to time k, Indicates the weight of the jth derived target Gaussian component at time k, Indicates the mean value of the j-th Gaussian component derived from the i-th Gaussian component from time k-1 to k-time, Indicates the covariance of the j-th Gaussian component derived from the i-th Gaussian component from time k-1 to k-time; Indicates the constant factor of the j-th inverse gamma component derived from the i-th inverse gamma component from time k-1 to k time, Indicates the iteration factor of the j-th inverse gamma component derived from the i-th inverse gamma component from time k-1 to k time;

2g)计算新生目标状态和量测噪声协方差的联合概率密度γ_k(x,R)：2g) Compute the joint probability density γ _k (x,R) of the nascent target state and measurement noise covariance:

${γ γ}_{k k} ((x x,, R R)) = = {Σ Σ}_{i i = = 11}^{{J J}_{γ γ,, k k}} [[{w w}_{γ γ . . k k}^{((i i))} N N ((x x;; {m m}_{γ γ,, k k}^{((i i))},, {P P}_{γ γ,, k k}^{((i i))})) {Π Π}_{l l = = 11}^{d d} IG IG (({(({σ σ}_{γ γ,, k k,, l l}^{((i i))}))}^{22};; {α α}_{γ γ,, k k,, l l}^{((i i))},, {β β}_{γ γ,, k k,, l l}^{((i i))}))]]$

其中，J_γ,k为k时刻新生目标高斯分量数，为为k时刻新生目标第i个高斯分量的权值，为新生目标第i个高斯分量的状态均值，为新生目标第i个高斯分量的运动状态协方差；为新生目标第i个逆伽玛分量的常量因子，为新生目标第i个逆伽玛分量的迭代因子。Among them, J _γ,k is the number of Gaussian components of the newborn target at time k, is the weight of the i-th Gaussian component of the newborn target at time k, is the state mean of the i-th Gaussian component of the newborn target, is the motion state covariance of the i-th Gaussian component of the newborn target; is the constant factor of the ith inverse gamma component of the nascent target, It is the iteration factor of the ith inverse gamma component of the nascent target.

2h)利用步骤2a)到步骤2g)得到的参数，计算扩展目标状态和量测噪声协方差的联合概率假设密度v_k|k-1(x,R)：2h) Using the parameters obtained from steps 2a) to 2g), calculate the joint probability hypothesis density v _k|k-1 (x,R) of the extended target state and measurement noise covariance:

v_k|k-1(x,R)＝v_S,k|k-1(x,R)+b_k|k-1(x,R)+γ_k(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)：Step 3, when k≥1, predict the potential distribution P _k-1|k-1 (num), and obtain the predicted potential distribution P _k|k-1 (num):

${P P}_{k k | | k k - - 11} ((num num)) = = {Σ Σ}_{j j = = 00}^{num num} {P P}_{birth birth,, k k} ((num num - - j j)) {Σ Σ}_{h h = = j j}^{\infty \infty} \frac{h h!!}{j j!! ((h h - - j j))!!} {P P}_{k k - - 11 | | k k - - 11} ((h h)) {p p}_{S S,, k k}^{j j} {(({11 - - p p}_{S S,, k k}))}^{h h - - j 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个扩展目标消亡的概率。Among them, P _birth,k (num-j) represents the probability of newborn num-j extended targets from time k-1 to time k, p _S,k represents the survival probability of extended targets at time k, P _k-1|k-1 (h) represents the probability that there are h extended targets at k-1 time,! means factorial, Indicates the probability that j extended targets survive from time k-1 to time k, and (1-p _S,k ) ^hj indicates the probability that hj extended targets die from time k-1 to time k.

步骤4.对扩展目标状态和量测噪声协方差的联合概率假设密度v_k|k-1(x,R)及预测势分布P_k|k-1(num)进行更新：Step 4. Update the joint probability hypothesis density v _k|k-1 (x, R) and the predicted potential distribution P _k|k-1 (num) of the extended target state and measurement noise covariance:

4a)利用变分贝叶斯VB方法将联合概率假设密度v_k|k-1(x,R)近似为：4a) Approximate the joint probability hypothesis density v _k|k-1 (x, R) using the variational Bayesian VB method as:

式中，Q_x,k|k-1(x)为高斯分布的求和形式， In the formula, Q _x,k|k-1 (x) is the sum form of Gaussian distribution,

Q_R,k|k-1(R)为逆伽玛分布的求和形式,Q _R,k|k-1 (R) is the summation form of inverse gamma distribution,

${Q Q}_{R R,, k k | | k k - - 11} ((R R)) = = {Σ Σ}_{i i = = 11}^{{J J}_{k k}} [[{Π Π}_{l l = = 11}^{d d} IG IG (({(({σ σ}_{k k | | k k - - 11,, l l}^{((i i))}))}^{22};; {α α}_{k k | | k k - - 11,, l l}^{((i i))},, {β β}_{k k | | k k - - 11,, l l}^{((i i))}))]];;$

其中，表示第i个高斯分量在第k个时刻的预测权值，i＝1,...,J_k，J_k表示第k个时刻扩展目标高斯分量的个数，N(·)表示高斯分布，为第k个时刻预测得到的第i个高斯分量的均值，为第k个时刻预测得到的第i个高斯分量的协方差；IG(·)表示逆伽玛分布，为第k个时刻预测得到的第i个逆伽玛分量的常量因子，为第k个时刻预测得到的第i个逆伽玛分量的迭代因子，l＝1,...,d，d表示量测噪声协方差的维数；in, Represents the prediction weight of the i-th Gaussian component at the k-th moment, i=1,...,J _k , J _k represents the number of extended target Gaussian components at the k-th moment, N( ) represents the Gaussian distribution, is the mean value of the i-th Gaussian component predicted at the k-th moment, is the covariance of the i-th Gaussian component predicted at the k-th moment; IG(·) represents the inverse gamma distribution, is the constant factor of the i-th inverse gamma component predicted at the k-th moment, It is the iteration factor of the i-th inverse gamma component predicted at the k-th moment, l=1,...,d, where d represents the dimension of the measurement noise covariance;

4b)利用变分贝叶斯期望最大化VBEM方法对扩展目标状态的概率假设密度Q_x,k|k-1(x)中的高斯分量和量测噪声协方差的概率假设密度Q_R,k|k-1(R)中的逆伽玛分量进行迭代更新：4b) Using the Variational Bayesian Expectation Maximization VBEM method to expand the probability hypothesis density Q _x,k|k-1 (x) of the target state in the Gaussian component and the probability hypothesis density Q _R,k of the measurement noise covariance The inverse gamma component in _|k-1 (R) is iteratively updated:

(4b1)设定逆伽玛分量的常量因子和迭代因子其中l＝1,...,d，d为量测噪声协方差R的维数；(4b1) Set the constant factor of the inverse gamma component and iteration factor Where l=1,...,d, d is the dimension of the measurement noise covariance R;

(4b2)根据设定的逆伽玛分量的两个因子，计算量测噪声协方差：(4b2) Calculate the measurement noise covariance according to the two factors of the set inverse gamma component:

其中n＝1,...,N，N为最大迭代次数，diag[...]表示对角化其中的元素； Among them, n=1,...,N, N is the maximum number of iterations, diag[...] means to diagonalize the elements;

(4b3)利用量测噪声协方差计算更新因子 (4b3) Using the measurement noise covariance Calculate update factor

${S S}_{W W}^{((n no))} = = {H h}_{W W} {P P}_{k k | | k k - - 11}^{((i i))} {H h}_{W W}^{T T} + + {R R}_{W W}^{((n no))}$

其中，表示对量测噪声协方差矩阵进行对角线连接当前单元W的量测数目后的矩阵，blkdiag(·)表示对其中的元素进行对角线连接，|W|表示当前单元W的量测数目；H_W表示对观测矩阵H_k垂直连接当前单元W的量测数目后的矩阵，表示k时刻的观测矩阵H_k的转置；表示k-1到k时刻预测的扩展目标高斯分量运动状态协方差，表示矩阵H_W的转置；in, Denotes the measurement noise covariance matrix The matrix after diagonally connecting the number of measurements of the current unit W, blkdiag(·) indicates that the elements in it are connected diagonally, |W| indicates the number of measurements of the current unit W; H _W indicates the observation matrix H _k vertically connects the matrix after the measurement number of the current unit W, Represents the transpose of the observation matrix H _k at time k; Represents the covariance of the extended target Gaussian component motion state predicted from k-1 to k time, Represents the transpose of matrix H _W ;

(4b4)利用更新因子计算增益矩阵 (4b4) Utilize update factor Calculate the gain matrix

${K K}_{k k}^{((i i)) ((n no))} = = {P P}_{k k | | k k - - 11}^{((i i))} {H h}_{W W}^{T T} {[[{S S}_{W W}^{((n no))}]]}^{- - 11}$

其中[·]^-1表示对矩阵求逆；Where [ ] ^-1 means to invert the matrix;

(4b5)利用增益矩阵计算扩展目标高斯分量运动状态和扩展目标高斯分量运动状态协方差 (4b5) Using the gain matrix Compute the extended target Gaussian component motion state and the extended target Gaussian component motion state covariance

${m m}_{k k | | k k}^{((i i)) ((n no))} = = {m m}_{k k | | k k - - 11}^{((i i))} + + {K K}_{k k}^{((i i)) ((n no))} (({z z}_{W W} - - {H h}_{W W} {m m}_{k k | | k k - - 11}^{((i i))}))$

${P P}_{k k | | k k}^{((i i)) ((n no))} = = [[I I - - {K K}_{k k}^{((i i)) ((n no))} {H h}_{W W}]] {P P}_{k k | | k k - - 11}^{((i i))}$

其中，I表示一个单位矩阵，z_W表示某个划分单元W中的所有量测；Among them, I represents an identity matrix, and z _W represents all measurements in a certain division unit W;

(4b6)提取扩展目标高斯分量运动状态的位置信息，将该位置信息用表示；(4b6) Extract the motion state of the extended target Gaussian component location information of the , use the location information express;

(4b7)利用该高斯分量位置信息和量测噪声协方差计算量测Y_n′由位置处的高斯分量生成的概率γ_n′i：(4b7) Use the Gaussian component position information and the measurement noise covariance Calculate the measurement Y _n′ by the position The probability γ _n′i generated by the Gaussian component at :

${γ γ}_{{n no}^{' '} i i} = = \frac{{π π}_{i i} N N (({Y Y}_{{n no}^{' '}} | | {m m}_{k k | | k k}^{' ' ((i i)) ((n no))},, {R R}_{k k}^{((i i)) ((n no))}))}{{Σ Σ}_{i i = = 11}^{{J J}_{k k}} {π π}_{i i} N N (({Y Y}_{{n no}^{' '}} | | {m m}_{k k | | k k}^{' ' ((i i)) ((n no))},, {R R}_{k k}^{((i i)) ((n no))}))},,$

其中，J_k表示扩展目标高斯分量数目，Y_n′表示当前单元W的第n′个量测，n′＝1,...,|W|；N(·)表示高斯分布；π_i代表混合系数，N_i表示由位置处的高斯分量产生的有效量测数目， Among them, J _k represents the number of Gaussian components of the extended target, Y _n' represents the n'th measurement of the current unit W, n'=1,...,|W|; N(·) represents the Gaussian distribution; π _i represents mixing coefficient, N _i is represented by the position The number of valid measurements produced by the Gaussian component at ,

(4b8)利用量测Y_n′由位置处的高斯分量生成的概率γ_n′i，迭代更新得到扩展目标高斯分量运动状态的位置信息 (4b8) Using the measurement Y _n' from the position The probability γ _n′i generated by the Gaussian component at , iteratively updated to obtain the position information of the motion state of the extended target Gaussian component

${m m}_{k k | | k k}^{' '' ' ((i i)) ((n no))} = = \frac{11}{{N N}_{i i}} {Σ Σ}_{{n no}^{' '} = = 11}^{| | W W | |} {γ γ}_{{n no}^{' '} i i} {Y Y}_{{n no}^{' '}};;$

(4b9)利用扩展目标高斯分量运动状态的位置信息混合系数π_i，量测噪声协方差计算最大似然函数L⁽ⁱ⁾⁽ⁿ⁾：(4b9) Use the position information of the extended target Gaussian component motion state Mixing coefficient π _i , measurement noise covariance Compute the maximum likelihood function L ⁽ⁱ⁾⁽ⁿ⁾ :

$L^{(i) (n)} = Σ_{n^{'} = 1}^{| W |} \ln Σ_{i = 1}^{J_{k}} π_{i} N (Y_{n^{'}} | m_{k | k}^{'' (i) (n)}, R_{k}^{(i) (n)}),$ 其中J_k表示扩展目标高斯分量数目； $L^{(i) (no)} = Σ_{{no}^{'} = 1}^{| W |} \ln Σ_{i = 1}^{J_{k}} π_{i} N (Y_{{no}^{'}} | m_{k | k}^{'' (i) (no)}, R_{k}^{(i) (no)}),$ where J _k represents the number of Gaussian components of the extended target;

(4b10)判断|L⁽ⁱ⁾⁽ⁿ⁾-L^(i)(n-1)|是否小于常量ε＝0.01，同时判断当前迭代次数n是否小于最大迭代次数N＝100，如果是，则停止迭代；否则，返回步骤(4b2)，更新逆伽玛分量迭代因子： $β_{k, 1}^{(i) (n + 1)} = \frac{Σ β_{W}^{(n + 1)}}{d | W |},$ (4b10) Judging whether |L ⁽ⁱ⁾⁽ⁿ⁾ -L ^(i)(n-1) | is less than the constant ε=0.01, and judging whether the current number of iterations n is less than the maximum number of iterations N=100, if yes, stop iteration; otherwise, return to step (4b2), and update the iteration factor of the inverse gamma component: $β_{k, 1}^{(i) (no + 1)} = \frac{Σ β_{W}^{(no + 1)}}{d | W |},$

其中，表示对向量中所有元素进行相加，in, Represents a pair of vectors All elements in are added together,

$β_{W}^{(n + 1)} = β_{W}^{(n)} + \frac{1}{2} {(z_{W} - H_{W} m_{k | k}^{(i) (n)})}_{j}^{2} + \frac{1}{2} {(H_{W} P_{k | k}^{(i) (n)} H_{W}^{T})}_{jj}, {(\cdot)}_{j}^{2}$ 表示对向量的第j维元素的平方，(·)_jj表示取矩阵的对角线元素，表示对迭代因子垂直连接当前单元W的量测数目后的向量，z_W表示当前单元的量测； $β_{W}^{(no + 1)} = β_{W}^{(no)} + \frac{1}{2} {(z_{W} - h_{W} m_{k | k}^{(i) (no)})}_{j}^{2} + \frac{1}{2} {(h_{W} P_{k | k}^{(i) (no)} h_{W}^{T})}_{jj}, {(&Center Dot;)}_{j}^{2}$ Indicates the square of the j-th dimension element of the vector, (·) _jj represents the diagonal element of the matrix, represents the iteration factor The vector after the number of measurements of the current unit W is vertically connected, z _W represents the measurement of the current unit;

(4b11)提取扩展目标状态分量扩展目标运动状态协方差迭代因子即， $m_{k | k}^{(i)} = m_{k | k}^{(i) (n)}, P_{k | k}^{(i)} = P_{k | k}^{(i) (n)}, β_{k, l}^{(i)} = β_{k, l}^{(i) (n)},$ 其中扩展目标状态分量中的位置信息用的是步骤(4b8)中迭代更新得到的扩展目标分量运动状态的位置信息 (4b11) Extract extended target state components Extended target motion state covariance iteration factor Right now, $m_{k | k}^{(i)} = m_{k | k}^{(i) (no)}, P_{k | k}^{(i)} = P_{k | k}^{(i) (no)}, β_{k, l}^{(i)} = β_{k, l}^{(i) (no)},$ where the extended target state component The position information in is the position information of the motion state of the extended target component obtained by iterative update in step (4b8)

4c)对势分布P_k|k-1(num)进行更新，得到更新后的势分布P_k|k(num)表示如下：4c) The potential distribution P _k|k-1 (num) is updated, and the updated potential distribution P _k|k (num) is expressed as follows:

${P P}_{k k | | k k} ((num num)) = = \{\begin{matrix} \frac{{Σ Σ}_{p p &angle; &angle; Z Z} {Σ Σ}_{W W &Element; &Element; p p} {ψ ψ}_{p p,, W W} {G G}_{k k | | k k - - 11}^{((num num))} ((00)) (\begin{matrix} {G G}_{FA FA} ((00)) \frac{ηW ηW}{| | p p | |} \frac{{ρ ρ}^{num num - - | | P P | |}}{((num num - - | | p p | |))!!} {δ δ}_{num num &GreaterEqual; &Greater Equal; | | p p | |} \\ + + {G G}_{FA FA}^{((| | W W | |))} ((00)) \frac{{ρ ρ}^{num num - - p p + + 11}}{((num num - - | | p p | | + + 11))!!} {δ δ}_{num num &GreaterEqual; &Greater Equal; | | p p | | - - 11} \end{matrix})}{{Σ Σ}_{p p &angle; &angle; Z Z} {Σ Σ}_{W W &Element; &Element; p p} {Ψ Ψ}_{p p,, W W} {l l}_{p p,, W W}},, & | | Z Z | | &NotEqual; &NotEqual; 00 \\ \frac{{ρ ρ}^{num num} {G G}_{k k | | k k - - 11}^{((num num))} ((00))}{{G G}_{k k | | k k - - 11} ((ρ ρ))} & | | Z Z | | = = 00 \end{matrix}$

其中，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表示目标产生量测概率的乘积，ρ表示扩展目标分量没有被检测到的概率：Among them, p∠Z indicates that the measurement set Z is divided into p non-empty subsets, W∈p indicates a certain unit under the p-th non-empty subset, and G _k|k-1 (ρ) indicates the state prediction probability generation function , Represents the num order partial derivative of the state prediction probability generation function, G _FA (0) represents the false alarm probability generation function when there is no measurement, ηW represents the measurement probability of the expanded target, |p| represents all non- the number of empty cells, Indicates the |W| order partial derivative of the false alarm probability generation function, δ _num≥|p| indicates that when the target number num is greater than the division unit |p|, the value is 1, otherwise it is 0, and |Z|=0 indicates that the extended target has no Generate measurement, |W| indicates the number of measurements in each non-empty unit W, l _{p, W} indicates the false alarm constant coefficient when the measurement division unit is |p|-1, ψ _{p, W} indicates the target generation amount The product of the detection probability, ρ represents the probability that the extended target component is not detected:

$ρ ρ = = Σ Σ {\overset{&OverBar; &OverBar;}{ω ω}}_{k k | | k k - - 11}^{j j} [[11 - - {P P}_{D D.} ((\cdot \cdot)) + + {P P}_{D D.} ((\cdot \cdot)) {G G}_{z z} ((00 | | \cdot &Center Dot;))]]$

$ηW ηW = = {p p}_{k k | | k k - - 11} [[{P P}_{D D.} ((\cdot \cdot)) {G G}_{z z}^{((| | W W | |))} ((00 | | \cdot &Center Dot;)) \underset{{z z}^{' '} &Element; &Element; W W}{Π Π} \frac{{p p}_{z z} (({z z}^{' '} | | \cdot &Center Dot;))}{{p p}_{FA FA} (({z z}^{' '}))}]]$

${l l}_{p p,, W W} = = {G G}_{FA FA} ((00)) {G G}_{k k | | k k - - 11}^{((| | p p | |))} ((ρ ρ)) \frac{ηW ηW}{| | p p | |} + + {G G}_{FA FA}^{((| | W W | |))} ((00)) {G G}_{k k | | k k - - 11}^{((| | p p | | - - 11))} ((ρ ρ))$

${ψ ψ}_{p p,, W W} = = \underset{{W W}^{' '} &Element; &Element; p p - - W W}{Π Π} {ηW η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阶偏导。Among them, Π is the multiplication symbol, Indicates the proportion of the j-th Gaussian inverse gamma component in all current Gaussian inverse gamma components, p _k|k-1 represents the state transition probability density function of a single extended target, and p _z (z′|·) represents the extended target Measurement likelihood, p _FA (z′) represents false alarm measurement likelihood, G _z (0|·) represents measurement probability generating function, Indicates the |W| order partial derivative of the measurement probability generation function, z′∈W indicates that the measurement z′ belongs to the unit W, _PD ( ) indicates the detection probability, W′∈pW indicates that the unit W is removed from all units divided by p The remaining unit, ηW' represents the probability of false alarm measurement of the extended target, Represents the |p| order partial derivative of the state prediction probability generating function, Indicates the |p|-1 order partial derivative of the state prediction probability generating function.

步骤5.对更新后的高斯分量和逆伽玛分量进行修剪与合并，其步骤如下：Step 5. pruning and merging the updated Gaussian component and inverse gamma component, the steps are as follows:

(5a)设置两个修剪门限T1和T2，一个合并门限U：T1＝10^-5，T2＝120，U＝10；设置最大高斯逆伽玛分量数目：J_max＝100；(5a) Set two pruning thresholds T1 and T2, one merging threshold U: T1=10 ⁻⁵ , T2=120, U=10; set the maximum number of Gaussian inverse gamma components: J _max =100;

(5b)计算每个扩展目标分量对应的量测噪声协方差: (5b) Calculate the measurement noise covariance corresponding to each extended target component:

(5c)设变量l′＝0，对更新后的扩展目标分量进行修剪，得到修剪后的扩展目标分量对应的序号集合I为： $I = {i = 1, . . ., J_{k} | w_{k}^{(i)} > T 1, {| | R_{k}^{(i)} | |}_{2} < T 2};$ (5c) Setting variable l'=0, pruning the extended target component after the update, obtaining the serial number set I corresponding to the extended target component after pruning is: $I = {i = 1, . . ., J_{k} | w_{k}^{(i)} > T 1, {| | R_{k}^{(i)} | |}_{2} < T 2};$

(5d)令l′＝l′+1，取表示取最大权值所对应的集合I中元素i′，对修剪后的扩展目标分量中满足合并门限U的分量进行提取，得到适合合并的扩展目标分量对应的序号集合为： $\overset{&OverBar;}{Q} = {i &Element; I | {(m_{k}^{(i)} - m_{k}^{(i^{'})})}^{T} {(P_{k}^{(i)})}^{- 1} (m_{k}^{(i)} - m_{k}^{(i^{'})}) \leq U};$ (5d) Let l'=l'+1, take means to take the maximum weight The element i' in the corresponding set I extracts the components that meet the merging threshold U in the pruned extended target components, and obtains the sequence number set corresponding to the extended target components suitable for merging for: $\overset{&OverBar;}{Q} = {i &Element; I | {(m_{k}^{(i)} - m_{k}^{(i^{'})})}^{T} {(P_{k}^{(i)})}^{- 1} (m_{k}^{(i)} - m_{k}^{(i^{'})}) \leq u};$

(5e)分别对序号集合中对应的扩展目标分量的权值运动状态常量因子迭代因子协方差进行合并，得到合并后的扩展目标分量的权值运动状态常量因子迭代因子协方差如下：(5e) Separately set the serial numbers The weight of the corresponding extended target component in exercise state constant factor iteration factor Covariance Merge to obtain the weight of the merged extended target component exercise state constant factor iteration factor Covariance as follows:

${\overset{~ ~}{w w}}_{k k}^{(({l l}^{' '}))} = = \underset{i i &Element; &Element; \overset{&OverBar; &OverBar;}{Q Q}}{Σ Σ} {w w}_{k k}^{((i i))},,$

${\overset{~ ~}{m m}}_{k k}^{(({l l}^{' '}))} = = \frac{11}{{\overset{~ ~}{w w}}_{k k}^{(({l l}^{' '}))}} \underset{i i &Element; &Element; \overset{&OverBar; &OverBar;}{Q Q}}{Σ Σ} {w w}_{k k}^{((i i))} {m m}_{k k}^{((i i))},,$

${\overset{~ ~}{α α}}_{k k}^{(({l l}^{' '}))} = = \frac{11}{{\overset{~ ~}{w w}}_{k k}^{(({l l}^{' '}))}} \underset{i i &Element; &Element; \overset{&OverBar; &OverBar;}{Q Q}}{Σ Σ} {w w}_{k k}^{((i i))} {α α}_{k k}^{((i i))},,$

${\overset{~ ~}{β β}}_{k k}^{(({l l}^{' '}))} = = \frac{11}{{\overset{~ ~}{w w}}_{k k}^{(({l l}^{' '}))}} \underset{i i &Element; &Element; \overset{&OverBar; &OverBar;}{Q Q}}{Σ Σ} {w w}_{k k}^{((i i))} {β β}_{k k}^{((i i))},,$

${\overset{~ ~}{P P}}_{k k}^{(({l l}^{' '}))} = = \frac{11}{{\overset{~ ~}{w w}}_{k k}^{(({l l}^{' '}))}} \underset{i i &Element; &Element; \overset{&OverBar; &OverBar;}{Q Q}}{Σ Σ} {w w}_{k k}^{((i i))} (({P P}_{k k}^{((i i))} + + (({\overset{~ ~}{m m}}_{k k}^{(({l l}^{' '}))} - - {m m}_{k k}^{((i i))})) {(({\overset{~ ~}{m m}}_{k k}^{(({l l}^{' '}))} - - {m m}_{k k}^{((i i))}))}^{T T}));;$

(5f)将步骤(5c)得到的修剪后的扩展目标分量对应的序号集合I中与步骤(5d)得到的适合合并的扩展目标分量对应的序号集合中相同的元素去除掉，然后判断修剪后的扩展目标分量对应的序号集合I是否为空集，如果不为空集，则返回步骤(5d)，否则执行(5g)；(5f) in the sequence number set I corresponding to the expanded target component after pruning that step (5c) obtains and the sequence number set corresponding to the extended target component that step (5d) obtains is suitable for merging Identical elements in are removed, then judge whether the sequence number set I corresponding to the extended target component after pruning is an empty set, if not an empty set, then return to step (5d), otherwise perform (5g);

(5g)判断变量l′是否大于最大高斯逆伽玛分量数目J_max，如果l′＞J_max，则将权值对应的高斯逆伽玛分量按从大到小排列，并取前J_max个权值大于0.5的高斯逆伽玛分量的位置和速度作为扩展目标的状态；如果l′＜J_max，则将所有权值大于0.5对应的高斯逆伽玛分量的位置和速度作为扩展目标的状态。(5g) Determine whether the variable l' is greater than the maximum number of Gaussian inverse gamma components J _max , if l'>J _max , then the weight The corresponding Gaussian inverse gamma components are arranged from large to small, and take the first J _max weights The position and velocity of the Gaussian inverse gamma component greater than 0.5 are used as the state of the extended target; if l′<J _max , all values The position and velocity of the Gaussian inverse gamma component corresponding to greater than 0.5 are used as the state of the extended target.

步骤6.根据步骤4c)更新得到的势分布P_k|k(num)进行加权平均，得到扩展目标的数目： ${num}_{k | k} = Σ_{num = 1}^{\infty} {num \times p}_{k | k} (num);$ Step 6. According to step 4c), the potential distribution P _k|k (num) obtained by updating the weighted average is obtained to obtain the number of expanded targets: ${num}_{k | k} = Σ_{num = 1}^{\infty} {num \times p}_{k | k} (num);$

得到扩展目标的状态和数目以后即可完成对扩展目标的跟踪。After obtaining the state and number of the extended object, the tracking of the extended object can be completed.

本发明用于跟踪扩展目标的效果可通过以下仿真实验进一步说明：The effect of the present invention for tracking extended target can be further illustrated by the following simulation experiments:

1.仿真条件1. Simulation conditions

考虑二维平面中做匀速直线运动4个目标，并且在两个位置出现运动轨迹交叉的情况，采样周期为整个观测过程持续40个时刻，目标的运动方程和量测方程分别为：Consider the case of 4 targets moving in a straight line at a uniform speed in a two-dimensional plane, and the intersection of motion trajectories occurs at two positions, and the sampling period is The whole observation process lasts for 40 moments, and the motion equation and measurement equation of the target are respectively:

x_k＝Fx_k-1+Gw_k x _k ＝Fx _k-1 +Gw _k

y_k＝Hx_k+v_k y _k ＝Hx _k +v _k

其中，高斯噪声w_k的标准差σ_w＝2。扩展目标状态转移矩阵F，输入矩阵G，高斯噪声w_k，观测矩阵H设置为：Wherein, the standard deviation σ _w =2 of the Gaussian noise w _k . The extended target state transition matrix F, input matrix G, Gaussian noise w _k , and observation matrix H are set as:

$F f = = [\begin{matrix} 11 & 00 & T T & 00 \\ 00 & 11 & 00 & T T \\ 00 & 00 & 11 & 00 \\ 00 & 00 & 00 & 11 \end{matrix}],, G G = = [\begin{matrix} 11 / / 22 & 00 \\ 00 & 11 / / 22 \\ 11 & 00 \\ 00 & 11 \end{matrix}],, {w w}_{k k} ~ ~ N N (\begin{matrix} 00,, & [\begin{matrix} {σ σ}_{w w}^{22} & 00 \\ 00 & {σ σ}_{w w}^{22} \end{matrix}] \end{matrix}),, H h = = [\begin{matrix} 11 & 00 & 00 & 00 \\ 00 & 11 & 00 & 00 \end{matrix}],,$

实际目标跟踪场景中量测噪声v_k的协方差是未知的，本次实验中设真实的量测噪声标准差σ_v＝1；The covariance of the measurement noise v _k in the actual target tracking scene is unknown. In this experiment, the real measurement noise standard deviation σ _v =1;

设传统GM-CPHD算法中采用的三个量测噪声标准差分别为σ＝0.05,1,50；Suppose the standard deviations of the three measurement noises used in the traditional GM-CPHD algorithm are σ=0.05, 1, 50;

设目标产生量测数目服从泊松分布，泊松分布参数β＝20，目标产生量测位置服从高斯分布；Assume that the number of target generation measurements obeys the Poisson distribution, the Poisson distribution parameter β=20, and the target generation measurement positions obey the Gaussian distribution;

设目标存活概率p_S,k＝0.99，检测概率p_D,k＝0.98；Set the target survival probability p _S,k =0.99, detection probability p _D,k =0.98;

设杂波数服从均值为5的泊松分布，杂波在整个观测区域内服从均匀分布。Assume that the number of clutter obeys the Poisson distribution with a mean of 5, and the clutter obeys a uniform distribution in the entire observation area.

设置新生目标的状态为 $m_{γ}^{(1)} = {[- 100,200,0,0]}^{T}, m_{γ}^{(2)} = {[0, - 100,0,0]}^{T},$ 新生目标的协方差均为P_γ＝diag[10,10,10,10]；Set the status of the newborn target to $m_{γ}^{(1)} = {[- 100,200,0,0]}^{T}, m_{γ}^{(2)} = {[0, - 100,0,0]}^{T},$ The covariance of the newborn target is P _γ =diag[10,10,10,10];

设置修剪门限T1＝10^-5，修剪门限T2＝120，合并门限U＝10，最大高斯逆伽玛分量数目J_max＝100，遗忘因子ρ＝0.9，实验仿真次数为100次。Set pruning threshold T1=10 ⁻⁵ , pruning threshold T2=120, merging threshold U=10, maximum number of Gaussian inverse gamma components J _max =100, forgetting factor ρ=0.9, and the number of experimental simulations is 100.

2.仿真内容与结果2. Simulation content and results

仿真1，采用本发明方法对未知量测噪声协方差下扩展目标运动轨迹的单次实验进行跟踪，结果如图2所示。从图2可以看出，本发明方法可以很好的跟踪扩展目标运动轨迹。In simulation 1, the method of the present invention is used to track a single experiment of the extended target trajectory under the unknown measurement noise covariance, and the result is shown in FIG. 2 . It can be seen from FIG. 2 that the method of the present invention can track the trajectory of the extended target very well.

仿真2，采用本发明方法与传统的分别采用不同量测噪声协方差的三个GM-CPHD方法进行100次实验的目标数目估计对比，结果如图3所示。从图3可以看出，本发明方法与传统的给定真实量测噪声协方差下的GM-CPHD方法对目标数目的估计效果一样良好。但是，传统的GM-CPHD方法，当量测噪声协方差未知时，且使用的量测噪声协方差与真实值偏差比较大时，估计的目标数目会发生很大的偏差，同时在k＝20,k＝35这两个目标运动轨迹出现交叉时刻，目标数目出现漏检；Simulation 2, using the method of the present invention and the traditional three GM-CPHD methods respectively using different measurement noise covariances for 100 times of target number estimation comparison, the results are shown in Figure 3 . It can be seen from FIG. 3 that the method of the present invention is as good as the traditional GM-CPHD method in estimating the number of targets given the covariance of real measurement noise. However, in the traditional GM-CPHD method, when the measurement noise covariance is unknown and the deviation between the measurement noise covariance and the true value is relatively large, the estimated number of targets will have a large deviation. At the same time, when k=20 ,k=35 When the two target trajectories intersect, the number of targets is missed;

仿真3，采用本发明方法与传统的分别采用不同量测噪声协方差的三个GM-CPHD方法借助OSPA距离进行100次实验的跟踪精度对比，结果如图4所示。从图4可以看出，本发明方法与传统的给定真实量测噪声协方差下的GM-CPHD算法的跟踪精度一样良好。但是，传统的GM-CPHD方法，当使用的量测噪声协方差与真实的量测噪声协方差偏差比较大时，跟踪精度会变得很差。Simulation 3, using the method of the present invention and the traditional three GM-CPHD methods respectively using different measurement noise covariances to compare the tracking accuracy of 100 experiments with the help of OSPA distance, the results are shown in Figure 4. It can be seen from Fig. 4 that the tracking accuracy of the method of the present invention is as good as that of the traditional GM-CPHD algorithm given the real measurement noise covariance. However, in the traditional GM-CPHD method, when the deviation between the used measurement noise covariance and the real measurement noise covariance is relatively large, the tracking accuracy will become very poor.

实验表明，本发明方法在处理量测噪声协方差是未知条件下的扩展目标跟踪场景时，其跟踪效果优于传统的GM-CPHD扩展目标跟踪方法。Experiments show that the method of the present invention is better than the traditional GM-CPHD extended target tracking method when dealing with the extended target tracking scene under the condition that the measurement noise covariance is unknown.

Claims

1. An extended target tracking method based on variational Bayesian expectation maximization, comprising the steps of:

(1) At time k=0, initialize the joint probability assumption density of the extended target state and the measurement noise covariance as v ₀ (x,R);

(2) When k ≥ 1, the joint probability assumption density v _k-1|k-1 (x, R) for the extended target state and measurement noise covariance and the potential distribution P _k- for calculating the number of extended targets _1|k-1 (num) is predicted, and the predicted extended target joint probability hypothesis density v _k|k-1 (x, R) and the predicted potential distribution P _k|k-1 (num) are obtained;

(3) The joint probability assumption density v _k|k-1 (x, R) for the predicted extended target state and the measurement noise covariance and the potential distribution P _k|k-1 (num) used to calculate the number of extended targets Make an update:

(3a) Use the variational Bayesian VB method to approximate the joint probability hypothesis density v _k|k-1 (x, R), and obtain the probability hypothesis density Q _{x of the extended target state represented by the summation form of the Gaussian distribution, k|k-1} (x) and the probability hypothesis density Q _R,k|k-1 (R) of the measurement noise covariance expressed in the sum form of the inverse gamma distribution;

(3b) Using the Variational Bayesian Expectation Maximization VBEM method to expand the probability hypothesis density Q _x,k|k-1 (x) of the target state and the probability hypothesis density Q _{R of the measurement noise covariance in the Gaussian component,} The inverse gamma component in _k|k-1 (R) is iteratively updated to obtain the Gaussian component representing the extended target motion state x and the inverse gamma component representing the measurement noise covariance R;

(3c) Update the potential distribution P _k|k-1 (num) predicted in step (2) to obtain the updated potential distribution P _k|k (num);

(4) pruning and merging the updated Gaussian component and the inverse gamma component, and extracting the position and velocity of the merged Gaussian component and the inverse gamma component as the state of the expanded target;

(5) Perform a weighted average on the potential distribution P _k|k (num) updated in step (3) to obtain the number of expanded targets:

{num}_{k | k} = Σ_{num = 1}^{\infty} num \times p_{k | k} (num);

(6) Repeat steps (2)-(5) to continue tracking the extended target.

2. The extended target tracking method based on variational Bayesian expectation maximization according to claim 1, wherein, the utilization of variational Bayesian VB method described in step (3a) to the joint probability hypothesis density v _{k|k -1} (x,R) for approximation, according to the following formula:

{Q Q}_{x x,, k k | | k k - - 11} ((x x)) = = {Σ Σ}_{i i = = 11}^{{J J}_{k k}} [[{w w}_{k k | | k k - - 11}^{((i i))} N N ((x x;; {m m}_{k k | | k k - - 11}^{((i i))},, {P P}_{k k | | k k - - 11}^{((i i))}))]]

{Q Q}_{R R,, k k | | k k - - 11} ((R R)) = = {Σ Σ}_{l l = = 11}^{{J J}_{k k}} [[{Π Π}_{l l = = 11}^{d d} IG IG (({(({σ σ}_{k k | | k k - - 11,, l l}^{((i i))}))}^{22};; {α α}_{k k | | k k - - 11,, l l}^{((i i))},, {β β}_{k k | | k k - - 11,, l l}^{((i i))}))]]

Among them, Q _x,k|k-1 (x) is the sum form of Gaussian distribution, expressed as

Q _{R, k|k-1} (R) is the summation form of the inverse gamma distribution, expressed as

{Q Q}_{R R,, k k | | k k - - 11} ((R R)) = = {Σ Σ}_{l l = = 11}^{{J J}_{k k}} [[{Π Π}_{l l = = 11}^{d d} IG IG (({(({σ σ}_{k k | | k k - - 11,, l l}^{((i i))}))}^{22};; {α α}_{k k | | k k - - 11,, l l}^{((i i))},, {β β}_{k k | | k k - - 11,, l l}^{((i i))}))]];;

Represents the weight of the i-th Gaussian component at the k-th moment, i=1,...,J _k , J _k represents the number of expanded target Gaussian components at the k-th moment, N( ) represents the Gaussian distribution, is the mean value of the i-th Gaussian component predicted at the k-th moment, is the covariance of the i-th Gaussian component predicted at the k-th moment; IG(·) represents the inverse gamma distribution, is the constant factor of the i-th inverse gamma component predicted at the k-th moment, is the iteration factor of the i-th inverse gamma component predicted at the k-th moment, l=1,...,d, where d represents the dimension of the measurement noise covariance.

3. the extended target tracking method based on variational Bayesian expectation maximization according to claim 1, wherein, the described utilization variational Bayesian expectation maximization VBEM method of step (3b) is to the probability of extended target state The Gaussian component in the assumption density Q _x,k|k-1 (x) and the inverse gamma component in the probability assumption density Q _R,k|k-1 (R) of the measurement noise covariance are updated iteratively as follows Steps to proceed:

(3b1) Set the constant factor of the inverse gamma component and iteration factor Where l=1,...,d, d is the dimension of the measurement noise covariance R;

(3b2) Calculate the measurement noise covariance according to the two factors of the set inverse gamma component: Among them, n=1,...,N, N is the maximum number of iterations, diag[...] means to diagonalize the elements;

(3b3) Using measurement noise covariance Calculate update factor

{S S}_{W W}^{((n no))} = = {H h}_{W W} {P P}_{k k | | k k - - 11}^{((i i))} {H h}_{W W}^{T T} + + {R R}_{W W}^{((n no))}

in, Denotes the measurement noise covariance matrix The matrix after diagonally connecting the number of measurements of the current unit W, blkdiag( ) indicates that the elements in it are connected diagonally, |W| indicates the number of measurements of the current unit W; H _W indicates the observation matrix H _k The matrix after the number of measurements of the current unit W is vertically connected, Represents the transpose of the observation matrix H _k at time k; Represents the covariance of the extended target Gaussian component motion state predicted from k-1 to k time, Represents the transpose of matrix H _W ;

(3b4) Utilize update factor Calculate the gain matrix

{K K}_{k k}^{((i i)) ((n no))} = = {P P}_{k k | | k k - - 11}^{((i i))} {H h}_{W W}^{T T} {[[{S S}_{W W}^{((n no))}]]}^{- - 11},,

Where [ ] ^-1 means to invert the matrix;

(3b5) Using the gain matrix Compute the extended target Gaussian component motion state and the extended target Gaussian component motion state covariance

{m m}_{k k | | k k}^{((i i)) ((n no))} = = {m m}_{k k | | k k - - 11}^{((i i))} + + {K K}_{k k}^{((i i)) ((n no))} (({z z}_{w w} - - {H h}_{w w} {m m}_{k k | | k k - - 11}^{((i i))}))

{P P}_{k k | | k k}^{((i i)) ((n no))} = = [[I I - - {K K}_{k k}^{((i i)) ((n no))} {H h}_{W W}]] {P P}_{k k | | k k - - 11}^{((i i))}

Among them, I represents an identity matrix, and z _W represents all measurements in a certain division unit W;

(3b6) Extract the motion state of the extended target Gaussian component location information of the , use the location information express;

(3b7) Use the Gaussian component position information and the measurement noise covariance The calculated measure Y _n′ is determined by the position The probability γ _n′i generated by the Gaussian component at :

{γ γ}_{{n no}^{' '} i i} = = \frac{{π π}_{i i} N N (({Y Y}_{{n no}^{' '}} | | {m m}_{k k | | k k}^{' ' ((i i)) ((n no))},, {R R}_{k k}^{((i i)) ((n no))}))}{{Σ Σ}_{i i = = 11}^{{J J}_{k k}} {π π}_{i i} N N (({Y Y}_{{n no}^{' '}} | | {m m}_{k k | | k k}^{' ' ((i i)) ((n no))},, {R R}_{k k}^{((i i)) ((n no))}))}

Among them, J _k represents the number of Gaussian components of the extended target, Y _n' represents the n'th measurement of the current unit W, n'=1,...,|W|; N(·) represents the Gaussian distribution; π _i represents mixing coefficient, N _i is represented by the position The number of valid measurements produced by the Gaussian component at ,

(3b8) Using the measurement Y _n' from the position The probability γ _n′i generated by the Gaussian component at , iteratively updated to obtain the position information of the motion state of the extended target Gaussian component

{m m}_{k k | | k k}^{' '' ' ((i i)) ((n no))} = = \frac{11}{{N N}_{i i}} {Σ Σ}_{{n no}^{' '} = = 11}^{| | W W | |} {γ γ}_{{n no}^{' '} i i} {Y Y}_{{n no}^{' '}};;

(3b9) Use the position information of the extended target Gaussian component motion state Mixing coefficient π _i , measurement noise covariance Compute the maximum likelihood function L ⁽ⁱ⁾⁽ⁿ⁾ :

L^{(i) (no)} = Σ_{{no}^{'} = 1}^{| W |} 1 no Σ_{i = 1}^{J_{k}} π_{i} N (Y_{{no}^{'}} | m_{k | k}^{'' (i) (no)}, R_{k}^{(i) (no)}),

where J _k represents the number of Gaussian components of the extended target;

(3b10) Judging whether |L ⁽ⁱ⁾⁽ⁿ⁾ -L ^(i)(n-1) | is less than the constant ε=0.01, and judging whether the current number of iterations n is less than the maximum number of iterations N=100, if yes, stop Iterate, otherwise return to step (3b2), and update the iteration factor of the inverse gamma component:

in, Represents a pair of vectors All elements in are added together,

β_{W}^{(no + 1)} = β_{W}^{(no)} + \frac{1}{2} {(z_{W} - h_{W} m_{k | k}^{(i) (no)})}_{j}^{2} + \frac{1}{2} {(h_{W} P_{k | k}^{(i) (no)} h_{W}^{T})}_{jj},

Indicates the square of the j-th dimension element of the vector, (·) _jj represents the diagonal element of the matrix, represents the iteration factor The vector after vertically connecting the current unit W to measure the number, z _W represents the measurement of the current unit;

(3b11) Extract extended target state components Extended target motion state covariance iteration factor Right now,

m_{k | k}^{(i)} = m_{k | k}^{(i) (no)}, P_{k | k}^{(i)} = P_{k | k}^{(i) (no)}, β_{k . l}^{(i)} = β_{k, l}^{(i) (no)},

where the extended target state component The position information in is the position information of the motion state of the extended target component obtained by iterative update in step (3b8)

4. The extended target tracking method based on variational Bayesian expectation maximization according to claim 1, wherein said step (3c) updates the potential distribution P _k|k-1 (num), and obtains the updated The potential distribution P _k|k (num) of is expressed as follows:

{P P}_{k k | | k k} ((num num)) = = \{\begin{matrix} \frac{{Σ Σ}_{p p &angle; &angle; Z Z} {Σ Σ}_{W W &Element; &Element; p p} {ψ ψ}_{p p,, W W} {G G}_{k k | | k k - - 11}^{((num num))} ((00)) (\begin{matrix} {G G}_{FA FA} ((00)) \frac{ηW ηW}{| | p p | |} \frac{{ρ ρ}^{num num - - | | P P | |}}{((num num - - | | p p | |))!!} & {δ δ}_{num num &GreaterEqual; &Greater Equal; | | p p | |} \\ + + {G G}_{FA FA}^{((| | W W | |))} ((00)) \frac{{ρ ρ}^{num num - - p p + + 11}}{((num num - - | | p p | | + + 11))!!} & {δ δ}_{num num &GreaterEqual; &Greater Equal; | | p p | | - - 11} \end{matrix})}{{Σ Σ}_{p p &angle; &angle; Z Z} {Σ Σ}_{W W &Element; &Element; p p} {ψ ψ}_{p p,, W W} {l l}_{p p,, W W}},, & | | Z Z | | &NotEqual; &NotEqual; 00 \\ \frac{{ρ ρ}^{num num} {G G}_{k k | | k k - - 11}^{((num num))} ((00))}{{G G}_{k k | | k k - - 11} ((ρ ρ))} & | | Z Z | | = = 00 \end{matrix}

Among them, p∠Z indicates that the measurement set Z is divided into p non-empty subsets, W∈p indicates a certain unit under the p-th non-empty subset, and G _k|k-1 (ρ) indicates the state prediction probability generation function , Represents the num order partial derivative of the state prediction probability generation function, G _FA (0) represents the false alarm probability generation function when there is no measurement, ηW represents the measurement probability of the expanded target, |p| represents all non- the number of empty cells, Indicates the |W| order partial derivative of the false alarm probability generation function, δ _num≥p means that when the target number num is greater than the division unit |p|, the value is 1, otherwise it is 0, |Z|=0 means that the extended target has no generation |W| represents the number of measurements in each non-empty unit W, l _{p, W} represents the false alarm constant coefficient when the measurement division unit is |p|-1, ψ _{p, W} represents the measurement probability of target generation The product of , ρ represents the probability that the extended target component is not detected:

ρ ρ = = Σ Σ {\overset{&OverBar; &OverBar;}{ω ω}}_{k k | | k k - - 11}^{j j} [[11 - - {P P}_{D D.} ((\cdot \cdot)) + + {P P}_{D D.} ((\cdot &Center Dot;)) {G G}_{z z} ((00 | | \cdot &Center Dot;))]]

ηW ηW = = {p p}_{k k | | k k - - 11} [[{P P}_{D D.} ((\cdot &Center Dot;)) {G G}_{z z}^{((| | W W | |))} ((00 | | \cdot &Center Dot;)) \underset{{z z}^{' '} &Element; &Element; W W}{Π Π} \frac{{p p}_{z z} (({z z}^{' '} | | \cdot &Center Dot;))}{{p p}_{FA FA} (({z z}^{' '}))}]]

{l l}_{p p,, W W} = = {G G}_{FA FA} ((00)) {G G}_{k k | | k k - - 11}^{((| | p p | |))} ((ρ ρ)) \frac{ηW ηW}{| | p p | |} + + {G G}_{FA FA}^{((| | W W | |))} ((00)) {G G}_{k k | | k k - - 11}^{((| | p p | | - - 11))} ((ρ ρ))

{ψ ψ}_{p p,, W W} = = \underset{{W W}^{' '} &Element; &Element; p p - - W W}{Π Π} η η {W W}^{' '}

Among them, Π is the multiplication symbol, Indicates the proportion of the j-th Gaussian inverse gamma component in all current Gaussian inverse gamma components, p _k|k-1 represents the state transition probability density function of a single extended target, and p _z (z′|·) represents the extended target Measurement likelihood, p _FA (z′) represents false alarm measurement likelihood, G _z (0|·) represents measurement probability generating function, Indicates the |W| order partial derivative of the measurement probability generation function, z′∈W indicates that the measurement z′ belongs to the unit W, _PD ( ) indicates the detection probability, W′∈pW indicates that the unit W is removed from all units divided by p The remaining unit, ηW' represents the probability of false alarm measurement of the extended target, Represents the |p| order partial derivative of the state prediction probability generating function, Indicates the |p|-1 order partial derivative of the state prediction probability generating function.

5. the extended target tracking method based on variational Bayesian expectation maximization according to claim 1, wherein, the described Gaussian component and inverse gamma component after step (4) are carried out pruning and merging, press Follow the steps below:

(5a) Set two pruning thresholds T1 and T2, one merging threshold U: T1=10 ⁻⁵ , T2=120, U=10; set the maximum number of Gaussian inverse gamma components: J _max =100;

(5b) Calculate the measurement noise covariance corresponding to each extended target component:

(5c) Setting variable l'=0, pruning the extended target component after the update, obtaining the serial number set I corresponding to the extended target component after pruning is:

I = {i = 1, . . ., J_{k} | w_{k}^{(i)} > T 1, {| | R_{k}^{(i)} | |}_{2} < T 2};

(5d) Let l'=l'+1, take means to take the maximum weight The element i' in the corresponding set I extracts the components that meet the merging threshold U in the pruned extended target components, and obtains the sequence number set corresponding to the extended target components suitable for merging for:

\overset{&OverBar; &OverBar;}{Q Q} = = {{i i &Element; &Element; I I | | {(({m m}_{k k}^{((i i))} - - {m m}_{k k}^{(({i i}^{' '}))}))}^{T T} {(({P P}_{k k}^{((i i))}))}^{- - 11} (({m m}_{k k}^{((i i))} - - {m m}_{k k}^{(({i i}^{' '}))})) \leq \leq U u}};;

(5e) Separately set the serial numbers The weight of the corresponding extended target component in exercise state constant factor iteration factor Covariance Merge to obtain the weight of the merged extended target component exercise state constant factor iteration factor Covariance as follows:

{\overset{~ ~}{w w}}_{k k}^{(({l l}^{' '}))} = = \underset{i i &Element; &Element; \overset{&OverBar; &OverBar;}{Q Q}}{Σ Σ} {w w}_{k k}^{((i i)),,}

{\overset{~ ~}{m m}}_{k k}^{(({l l}^{' '}))} = = \frac{11}{{\overset{~ ~}{w w}}_{k k}^{(({l l}^{' '}))}} \underset{i i &Element; &Element; \overset{&OverBar; &OverBar;}{Q Q}}{Σ Σ} {w w}_{k k}^{((i i))} {m m}_{k k}^{((i i))},,

{\overset{~ ~}{α α}}_{k k}^{(({l l}^{' '}))} = = \frac{11}{{\overset{~ ~}{w w}}_{k k}^{(({l l}^{' '}))}} \underset{i i &Element; &Element; \overset{&OverBar; &OverBar;}{Q Q}}{Σ Σ} {w w}_{k k}^{((i i))} {α α}_{k k}^{((i i))},,

{\overset{~ ~}{β β}}_{k k}^{(({l l}^{' '}))} = = \frac{11}{{\overset{~ ~}{w w}}_{k k}^{(({l l}^{' '}))}} \underset{i i &Element; &Element; \overset{&OverBar; &OverBar;}{Q Q}}{Σ Σ} {w w}_{k k}^{((i i))} {β β}_{k k}^{((i i))},,

{\overset{~ ~}{P P}}_{k k}^{(({l l}^{' '}))} = = \frac{11}{{\overset{~ ~}{w w}}_{k k}^{(({l l}^{' '}))}} \underset{i i &Element; &Element; \overset{&OverBar; &OverBar;}{Q Q}}{Σ Σ} {w w}_{k k}^{((i i))} (({P P}_{k k}^{((i i))} + + (({\overset{~ ~}{m m}}_{k k}^{(({l l}^{' '}))} - - {m m}_{k k}^{((i i))})) {(({\overset{~ ~}{m m}}_{k k}^{(({l l}^{' '}))} - - {m m}_{k k}^{((i i))}))}^{T T}));;

(5f) in the sequence number set I corresponding to the expanded target component after pruning that step (5c) obtains and the sequence number set corresponding to the extended target component that step (5d) obtains is suitable for merging Identical elements in are removed, then judge whether the sequence number set I corresponding to the extended target component after pruning is an empty set, if not an empty set, then return to step (5d), otherwise execute (5g);

(5g) Determine whether the variable l' is greater than the maximum number of Gaussian inverse gamma components J _max , if l'>J _max , then the weight The corresponding Gaussian inverse gamma components are arranged from large to small, and take the first J _max weights The position and velocity of the Gaussian inverse gamma component greater than 0.5 are used as the state of the extended target; if l′<J _max , all values The position and velocity of the Gaussian inverse gamma component corresponding to greater than 0.5 are used as the state of the extended target.