Extended Target Marginal Distribution Poisson Multi-Bernoulli Mixture Filter

The existence of clutter, unknown measurement sources, unknown number of targets, and undetected probability are problems for multi-extended target tracking, to address these problems; this paper proposes a gamma-Gaussian-inverse Wishart (GGIW) implementation of a marginal distribution Poisson multi-Bernoulli mixture (MD-PMBM) filter. Unlike existing multiple extended target tracking filters, the GGIW-MD-PMBM filter computes the marginal distribution (MD) and the existence probability of each target, which can shorten the computing time while maintaining good tracking results. The simulation results confirm the validity and reliability of the GGIW-MD-PMBM filter.


Introduction
Multiple target tracking (MTT) involves estimating the state of an unknown number of targets in the presence of clutter [1][2][3]. The traditional MTT algorithm assumes each target is a point as one target generates at most one measurement in the sensor at each time step [3]. However, the rapid development of the modern high-resolution sensors makes the "point assumption" impractical because with such sensors, one target generates at least one measurement per time scan. The MTT problem with such sensors becomes an extended MTT problem [4,5]. Compared with the point target, the extended target can not only provide accurate target movement information, but also the target' shape and size information, making it useful in the fields of robot recognition and positioning, moving crowd tracking, and tracking of close cars or large ships using the high-resolution sensors or automotive radars. Due to its wide application, the research on extended target tracking has received considerable attention from many scholars, making it a growing research hotspot.
As one of the most used extended target measurement models, The Poisson Point Process (PPP) model assumes a Poisson distributed random number of measurements are spatially distributed around the target at each time step [4,5]. To efficiently and correctly represent the spatially distributed measurements, many shape models have been developed, such as the random matrix model (RMM) [6][7][8][9][10][11][12][13][14][15], random hypersurface model (RHM) [16][17][18][19][20][21] and Gaussian process model (GPM) [22]. RMM assumes that the shape of the measurements can be approximated by an ellipse, and the measurements surrounding the centroid of the target obey a Gaussian distribution. RHM uses a more general star-convex shape to approximate the target contour. The GPM automatically learns the shape of the target through a Gaussian process and can give an analytic expression of its contour for an arbitrary shape target. Among the above three measurement models, the RMM has the smallest amount of calculation and is easy to implement, and can meet the requirements of many practical situations, so we use the RMM method in this paper.
In the above-listed RFS-based filters, many kinds of conjugate priors are widely used, to give a closed-form solution of the posterior density. The prior conjugate refers to a time series of random finite set that satisfies the conjugate prior property; if the prediction is a Gaussian (Multi-Bernoulli) distribution, the update is also a Gaussian (Multi-Bernoulli) distribution. Using the conjugate prior, given enough parameters, the posterior distribution can be approximated arbitrarily.
In recent years, developing conjugate prior-based MTT filters has been a significant trend; the two kinds of most used conjugate priors in various MTT filters are the PMBM conjugate prior and the δ-Generalized Labelled Multi-Bernoulli (δ-GLMB) conjugate prior. In the δ-GLMB conjugate prior, the state of each target is added with a unique label, which is helpful to find each target's trajectory. The PMBM prior conjugate divides the target set into two disjoint subsets, targets that have been detected and targets that have not yet been detected, and then processes the two subsets separately. This strategy can greatly reduce the target missed detection rate and improve tracking accuracy. In reference [40], it has been shown that the point target tracking filter based on PMBM conjugate prior outperforms the filter based on the GLMB conjugate prior in terms of tracking accuracy and computing time; thus, in this paper we used the PMBM conjugate prior.
In this paper, to solve the extended MTT problem and combine the PPP model and PMBM conjugate prior, we propose a GGIW implementation of a marginal distribution Poisson multi-Bernoulli mixture (MD-PMBM) filter on the basis of the GGIW-PMBM filter. The main difference between the GGIW-MD-PMBM filter and the GGIW-PMBM filter is that in the prediction step and the update step, instead of recursively propagating the joint state distribution of targets (GGIW-PMBM), the GGIW-MD-PMBM filter recursively propagates the marginal distribution and the existence probability of each target; such a strategy can effectively reduce the calculation cost while maintaining good tracking results.
The rest of the paper is organized as follows. Several model assumptions, such as the Bayesian model, PPP model, MBM model, motion model, measurement model, and PMBM conjugate prior, are discussed in Section 2. Section 3 presents the prediction step and the update step of the GGIW-MD-PMBM filter. Section 4 gives the four extended target filters' experimental simulation results, which are summarized in Section 5.

Background
In this section, we outline various model assumptions, such as the Bayesian model [45][46][47], PPP model, Multi-Bernoulli (MB) process model, motion and measurement model, and PMBM conjugate prior model, which were used in the prediction and update steps of the GGIW-MD-PMBM filter.

Bayesian Model
We define the state of the i-th extended target at time step k as ξ (i) k , and the state of the target set at time step k can be written as where N ξ,k = X k is the number of the target, which is unknown and time-varying. The measurements (the union of target generated measurements and clutter) obtained from the sensor at time k are where N z,k = Z k is the number of measurements. Let Z k denote the union of measurements set from time step 1 to k, that is, Z k = Z 1 ∪ Z 2 ∪ · · · ∪ Z k . In many point target-tracking filters and extended target-tracking filters, the main purpose of using the Bayesian model is to recursively propagate the multi-target density distribution. In a Bayesian model, the multi-target set density, multi-target transition density, and the multi-target measurement likelihood are represented by f k|k (X k Z k ) , f k|k−1 (X k X k−1 ) , and f k|k (Z k X k ) , respectively. Using the Chapman-Kolmogorov equation, the multi-target set density can be defined as and the Bayes update as

PPP Model and Multi-Bernoulli (MB) Process Model
The PPP model, proposed by Gilholm et al. [4] and Granstrom et al. [5], is widely used in point target tracking filters and extended target tracking filters. In a PPP model, each target generates a Poisson distributed number of measurements, and each measurement is independent and identically distributed (i.i.d.). The Poisson rate µ and the spatial distribution f (·) determine the PPP intensity.
We define the spatial distribution of the target state ξ The existence of a single target can be modeled by a Bernoulli RFS [14,[31][32][33][34][40][41][42]. The cardinality of a Bernoulli RFS X An MB RFS X k is the union of a limited number of independent Bernoulli RFS X In this paper, we divide the targets into the undetected target set and the detected target set. The PPP represents the distribution of undetected target measurements, and the MB mixture (MBM) represents the distribution of detected target measurements.

The Motion Model and Measurement Model
To simplify the calculation, at each time step, we assume each target evolves independently from the other targets, and the measurements generated by each target are independent of each other and independent of those generated by other targets.
We assume each birth target can be described by a PPP model, and its intensity is D b (x k ). All the existing targets either survive with a probability P D (x k ) or die with a probability 1 − P D (x k ) from time step k to time step k + 1. Each extended target state consists of three variables, the measurement rate γ k , kinematic state x k , and extent state E k , that is, ξ k . Given γ k , x k and E k , the Markov motion model can be defined as where w k is Gaussian process noise with zero mean and covariance Q, f (·) and M(·) are the transformation matrix of kinematic state and extent state, respectively. The clutter measurements at time step k can be modeled as a PPP model with measurement rate λ and spatial Poisson distribution c(z) . The clutter PPP density function is k(z) = λ c(z) . The measurements generated by each target at time step k are a Poisson distribution with a measurement rate γ k . The measurement model is where H k is the known observation model function, and υ k is the Gaussian observation noise with zero mean and covariance E k . Assuming the measurements generated by a single target follow a Gaussian distribution, that is The distribution of each target state can be expressed as where G(γ k ; α k|k , β k|k ) is the Gamma distribution, α k|k > 0 and β k|k > 0 are its shape parameter and inverse scale parameter; N(x k ; m k|k , P k|k ⊗ E k ) is the Gaussian distribution, m k|k and P k|k ⊗ E k are its mean and covariance; IW d (E k ; v k|k , V k|k ) is the inverse Wishart distribution, v k|k and V k|k are its degrees of freedom and shape matrix. Note that ζ k|k = α k|k , β k|k , m k|k , P k|k , v k|k , V k|k is a shorthand for the set of GGIW density parameters. The detailed implementation of GGIW prediction and update is given in [21,[23][24][25][26].

PMBM Conjugate Prior
Existing studies [40][41][42] show that developing PMBM conjugate priors-based filters is a significant trend, due to their excellent tracking performance. In a PMBM model, the extended target set X k at time k is divided into two disjoint sets, undetected targets X u k and detected targets X d k , then the PPP model is used to describe the distribution of undetected targets X u k , and the MBM model to describe the distribution of targets X d k that have been detected at least once.

of 15
The PMBM set density can be defined as where f d,j,i k|k (·) is the Bernoulli set density. The MBM has |J| components, and J is the index set (or MBM component) of the MB in the MBM. I j is the Bernoulli index set of the j-th MB, that is, the j-th component has I j Bernoulli components. ω d,j is the probability of the j-th MB component. The conjugacy property determines the intensities of the birth PPP and the initial undetected PPP, which are expressed by both GGIW mixture density.
Given the PMBM conjugate prior assumption, the extended target tracking is equivalent to recursively propagate the PMBM density parameters in the prediction step and the update step, which is shown in Section 3.

The GGIW-MD-PMBM Filter
In this section, the GGIW-MD-PMBM filter is presented. Instead of recursively propagating the joint state distribution of targets (GGIW-PMBM filter), the proposed GGIW-MD-PMBM filter recursively propagates the marginal distributions and the existence probabilities of each extended target. The GGIW-MD-PMBM filter can be modeled by the following assumptions: (1) Clutter is uniformly distributed in the surveillance area with a given Poisson rate λ and is independent of targets' distribution; (2) Each extended target generates at least one measurement per time scan and evolves independently of the other targets; (3) The birth PPP intensity is a GGIW mixture with a given Poisson rate µ b , weight ω b , and various density parameters ς b of GGIW; (4) The initial undetected PPP intensity is also a GGIW mixture with a known Poisson rate µ u 0 , weight ω u 0 , and various density parameters ς u 0 of GGIW.
We assume the birth PPP and the initial undetected targets are all GGIW mixture intensities; according to the PMBM conjugate property, each target density in the PMBM filter is also a GGIW density. where The existence probability of undetected targets ω ub,j k and ω uu,j k|k−1 can be normalized, as Through the above analysis, in the prediction step, all extended target distributions and existence probabilities can be described as Sensors 2020, 20, 5387 where

Update
If each target predicted density is a PMBM at time k, according to the conjugate prior property, the updated density is also a PMBM. The updated MBM is given by the formulas (31) and (32), which contain the MB components predicted by each target in the previous processing step and their related data associations.
Here, the predicted likelihood L P of partition and the predicted likelihood L w of each MB are as follows: In formulas (33) and (34), for each MB component in each partition cell, the data association can be divided into the following three types: (1) The measurement set Z can be divided into the union of two disjoint sets, Z = Y ∪ (Z\Y), where Y corresponds to the clutter and the previously undetected targets' measurements, and Z\Ycorresponds to the previously detected targets' measurements.
(2) For set Y, form a partition of P of non-empty cells C; each cell contains clutter or target-generated measurements.
(3) For set Z\Y, the union of a set of index subsets I j can be used to represent the measurements related to the j-th MB component of the detected targets.
In the update step, the Bernoulli parameters of each detected target and the PPP parameters of each undetected target are determined, as given below.

Detected Targets
In this paper, there are two types of detected targets at time k: The targets detected for the first time in the measurement set Y, and the previously detected targets are detected again in the measurement set Z\Y. The first detected targets are processed by the PPP model, and the second detected target is estimated with the existing MB. According to Equations (9) and (10), the spatial distribution f Sensors 2020, 20, 5387 For each partition cell, its predicted likelihood L C is related to the predicted likelihood of each GGIW component in the spatial density, such as The predicted likelihood L C is determined by C ∩ I j and the measurement cell C C , as follows: (1) When C ∩ I j = ∅ , C C = 1 , the predicted likelihood L C is an approximation of the predicted likelihood of the cell composed of the undetected target-generated measurements and the clutter measurements in the first detection process. Since this cell contains the clutter and undetected targets' measurements, L C is related to κ C C and l C , where κ C C is the predicted likelihood of clutter, and l C is the predicted likelihood of undetected target-generated measurements.
(2) When C ∩ I j = ∅ , C C > 1 , it means the targets detected for the first time are divided into multiple cells. If there are undetected targets, they must be contained in the above cells, thus L C is only related to l C .
(3) When C ∩ I j ∅ , C C = ∅, it means L C is the predicted likelihood when the j-th MB component of the previously detected target is an empty set.
(4) When C ∩ I j ∅ , C C ∅, the predicted likelihood L C is the approximation when the j-th MB component of the previously detected target is not empty. The calculation L u,j,C k is presented in [39], Table 2.
The predicted likelihood l C of undetected targets can be defined as According to the different values of C ∩ I j and C C , the likelihood probability r j,C can be defined as The effective probability of the missed detection can be defined as Sensors 2020, 20, 5387 9 of 15 The spatial distribution f d,j,i k|k (ξ k ) of each measurement cell at time step k is as follows: The second equation in (41) uses a gamma-mixture reduction to reduce its bi-modal GGIW distribution to a uni-modal GGIW distribution [44]. Detailed calculations of parameter ζ d,j,C k|k and prediction likelihood L u,j,C k in the GGIW update are provided in [41].

Undetected Targets
The detection probability p D is always less than 1, and the Poisson parameter γ k may be zero; these two cases may cause missed detection of targets. Therefore, we need to consider the above two situations, and the updated spatial density is where q u,j D is the effective probability of the undetected target, and µ u k|k is the PPP Poisson rate of the undetected targets, The existence probability is Note: Most of the mathematical details of formulas (42) and (45) are given in the Appendix A. The Gaussian and inverse Wishart parameters in formula (42) are the same in both cases, but the gamma parameters are different. Therefore, formula (42) can be reduced by gamma mixing [44] to a single-peak GGIW distribution and the number of GGIW components remains unchanged.
The main contribution of the GGIW-MD-PMBM proposed in this paper is to recursively propagate the marginal distribution and existence probability of each target through formulas (21), (22), (29), (34), (35), (42), and (45), which is different from the GGIW-PMBM in [40] that it recursively propagates the joint target state distribution and existence probability.

Complexity Reduction and Data Association
In the tracking process, as the number of targets increases, the number of unknown data associations, the number of components in the MBM, and the number of parameters in the PMBM greatly increase, which brings the problem of "combination explosion". It is necessary to approximate the target density through reduction methods. Possible reduction methods include gating, clustering, pruning, merging, and recycling. In this paper, we focused on developing the approximate method to recursively propagate the marginal distribution of targets in the prediction and update steps. The complexity reduction and the data association process used in this paper follow the same methods used in [40], Section V; hence the complexity reduction and the data association are briefly discussed.

Simulation
In this simulation, we compared the proposed GGIW-MD-PMBM filter with the GGIW-PMBM, GGIW-LMB, and GGIW-PHD filters for multi-target trajectories in the four classic scenarios, which are provided in the excellent sample MATLAB code [40][41][42][43]. Each extended target can be defined as where p k ∈ R 2 is each target's position and v k ∈ R 2 is the velocity. The motion model's parameters used in Formulas (8)-(10) are where σ a is the acceleration standard deviation and T s = 1s is the sampling time.
To evaluate the four algorithms' performance, we used the computing time and the generalized optimal sub-pattern assignment (GOSPA) metric [47], the cutoff parameter c = 10, and the ordered parameter p = 1, among which, for the distance measurement, we use the Gaussian Wasserstein Distance metric [48]. We divided the GOSPA metric into three categories: localization error, false detection error, and miss detection error. For GGIW-LMB filter and GGIW-PHD filter, we extracted the target states by taking the mean vector of all Bernoullis whose existence probability is larger than 0.5. For the GGIW-PMBM filter and GGIW-MD-PMBM filter, target state extraction was performed similarly, but only from the MB component with the highest weight.
In the four scenarios, the target detection probabilityp D , target survival probability p S , the clutter Poisson rate λ, and the measurement Poisson rate γ used in different scenarios are given in Table 1. In the first scenario, there are 100 time steps, and 27 highly time-varying targets are randomly generated in four positions; this scenario aims to compare the four algorithms' tracking performances with a high clutter density and high target number scenario. In the second scenario, there are 100 time steps, and two targets are born well separated, move close to each other, and then split; this scenario aims to compare the four algorithms' tracking performance when highly close to each other. In the third scenario, there are 10 time steps, and five targets are born closely at the same time; this scenario aims to compare the four algorithms' tracking performances of handling dense birth. In the fourth scenario, there are 300 time steps, and two targets first get close and then they maneuver closely before splitting; this scenario aims to compare the four algorithms' tracking performance of seriously handling the data association problem. The target trajectories of different scenarios are given in Figure 1. time; this scenario aims to compare the four algorithms' tracking performances of handling dense birth. In the fourth scenario, there are 300 time steps, and two targets first get close and then they maneuver closely before splitting; this scenario aims to compare the four algorithms' tracking performance of seriously handling the data association problem. The target trajectories of different scenarios are given in Figure 1.    Table 2 shows the estimation errors and cycling time of the four algorithms in the four   Table 2 shows the estimation errors and cycling time of the four algorithms in the four scenarios. From the simulation results in Figure 2 and Table 2, we can see that the proposed GGIW-MD-PMBM obtains the general minimum average GOSPA error, which outperforms the other three algorithms; GGIW-PMBM is the second, and GGIW-PHD obtains the highest average GOSPA error. As for the time consumption of each MC run, the GGIW-PHD has the lowest computational cost, GGIW-LMB is the second, and GGIW-MD-PMBM is the third. Compared with other results of normalized location error, the number of missed targets, and the number of false detection, the proposed GGIW-MD-PMBM algorithm is generally better than the other three algorithms. In the third scenario, the GGIW-PMBM filter outperformed the GGIW-MD-PMBM filter in terms of the location error and the number of false detections; this is because the targets in this scenario are intensively generated at the same location and at the same time. Compared with the joint state distribution (GGIW-PMBM), the marginal distribution (GGIW-MD-PMBM) has a larger deviation, but as time increases, the targets gradually move away, and the average GOSPA error of the GGIW-MD-PMBM filter is still lower than the GGIW-PMBM filter. The tracking performance of each algorithm in the fourth scenario is worse than in the other three scenarios; this is because the data association in the fourth scenario is worse than in the other three scenarios.
Sensors 2020, 20, x FOR PEER REVIEW 11 of 15 scenarios. From the simulation results in Figure 2 and Table 2, we can see that the proposed GGIW-MD-PMBM obtains the general minimum average GOSPA error, which outperforms the other three algorithms; GGIW-PMBM is the second, and GGIW-PHD obtains the highest average GOSPA error. As for the time consumption of each MC run, the GGIW-PHD has the lowest computational cost, GGIW-LMB is the second, and GGIW-MD-PMBM is the third. Compared with other results of normalized location error, the number of missed targets, and the number of false detection, the proposed GGIW-MD-PMBM algorithm is generally better than the other three algorithms. In the third scenario, the GGIW-PMBM filter outperformed the GGIW-MD-PMBM filter in terms of the location error and the number of false detections; this is because the targets in this scenario are intensively generated at the same location and at the same time. Compared with the joint state distribution (GGIW-PMBM), the marginal distribution (GGIW-MD-PMBM) has a larger deviation, but as time increases, the targets gradually move away, and the average GOSPA error of the GGIW-MD-PMBM filter is still lower than the GGIW-PMBM filter. The tracking performance of each algorithm in the fourth scenario is worse than in the other three scenarios; this is because the data association in the fourth scenario is worse than in the other three scenarios.

Conclusions
In this paper, we propose an efficient filter to solve the extended target tracking problem. Unlike the existing GGIW-PMBM filter that recursively propagates the joint state distribution of targets, the proposed GGIW-MD-PMBM filter recursively propagates the marginal distributions and the existence probabilities of each extended target. Comparing the other three extended target filters in four classic scenarios, the experimental results show that the proposed filter in this paper improves tracking accuracy and reduces computing time.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Let p e k,D be the effective probability of target detection [35] (the set of target detections is nonempty with probability 1 − e −γ( j) ), then the updated PPP intensity corresponding to previously existing targets that are not detected is At this time, all existence probabilities of the undetected targets can be normalized, and the existence probability and spatial density of each target can be as follows: where f u,j 1,k|k (ξ k ) represents the spatial density in the case of a missed detection caused by p D < 1, and f u,j 2,k|k (ξ k ) represents the spatial density in the case of a missed detection caused by γ k = 0.