New Box Particle Filter with Improved Resampling Method and Extended Inclusion Volume Criteria for Multi-Target Tracking

In the resampling procedure of traditional box particle filtering, selected box particles are divided in a randomly chosen dimension. This resampling procedure may fail when some elements in the target state vector are unmeasured. To deal with this problem, an improved resampling method for box particle filtering is proposed, where limits on the sizes of box particles are imposed to restrain the box particles from growing too large. In addition, we extend the inclusion and the volume criteria from single-target tracking to multi-target tracking. Instead of indicating whether the true target state is included in the support of the posterior track probability density in single target tracking, the inclusion value in multi-target tracking indicates how many true targets are included in the supports of the posterior probability densities. And the volume value in multi-target tracking is redefined as the mean volume of the supports of the posterior probability densities. Simulation results are provided to illustrate the effectiveness of the proposed approach.


Introduction
The problem of multi-target tracking is the on-line recursive estimation of an unknown and time-varying number of targets, their states and trajectories.Recently, the random finite set (RFS) approach [1] attracts much attention because it casts multi-target tracking into a Bayesian framework.The propagation of the multi-target state cannot be implemented directly due to prohibitive numerical complexity.Therefore, feasible approximations of the Bayes multi-target filter have been developed, including the probability hypothesis density (PHD) [2], [3], cardinalized PHD (CPHD) [4], multi-target multi-Bernoulli (MeMBer) [1], cardinalitybalanced MeMBer (CBMeMBer) [5], [6], δ-generalized la-beled multi-Bernoulli (δ-GLMB) [8], [7], and labeled multi-Bernoulli (LMB) [9], [10] filters.In this paper, the LMB filter is adopted as the basic framework for tracking.
Though the aforementioned filters are numerically less complex, they still involve multiple integrals without closedform expressions.To address this problem, two main implementation methods are developed: Gaussian mixture (GM) [11], and sequential Monte Carlo (SMC) or particle filter (PF) [12].The GM method is suitable for linear Gaussian dynamic and measurement models, whereas the SMC\PF method is suitable for nonlinear non-Gaussian scenarios.The GM method assumes precise point measurements affected by additive measurement noises of a Gaussian probability density function (pdf).However, many practical applications involve not only nonlinear non-Gaussian models, but also imprecise stochastic measurements with unknown biases and error distributions.Though the SMC\PF method can be used to handle such situations, it requires a large number of particles to achieve satisfactory performance.As a result, it can be so time-consuming that real-time implementation is hardly achievable.
In recent years, box particle filter (BPF) [13][14][15][16], an emerging implementation method for RFS-based filters, is gaining in popularity because it accommodates imprecise stochastic measurements with unknown biases and error distributions.A box particle is a vector comprised of intervals.The average computation cost per box particle is higher than that per point particle.Nevertheless, to achieve a similar level of performance, the BPF requires less computation power than the PF does, because the number of particles used by the BPF is much smaller than that used by the PF.
For the PF resampling, particles with large weights are replicated [12], while for the traditional BPF resampling, boxes with large weights are divided into sub-boxes as many times as they are selected [13].The dimension in which the selected boxes are divided is chosen randomly.This method is only suitable for the situations where all elements in the target state vector are measured.Because in such situations intervals in all dimensions of a box particle can be contracted by the measurements.So the choice of dimension is insignificant and thus can be random.However, this method may fail in situations where some elements in the target state vector are unmeasured.To deal with this problem, an improved resampling method for the BPF is proposed.A new parameter called box resolution vector is defined to set an upper limit to the interval width in every dimension of every box particle after resampling.Another contribution of this paper is that the inclusion and volume criteria [13] initially defined for single-target tracking are extended for multi-target tracking.We propose that the inclusion value for multi-target tracking indicate how many true targets are included in the supports of the posterior track pdfs, and that the volume value for multi-target tracking be redefined as the mean volume of the supports of the posterior track pdfs.
The paper is organized as follows.Background on LMB filter and interval analysis is presented in Sec. 2. Section 3 presents the proposed approach.Section 4 outlines BPF-LMB and PF-LMB implementations using the proposed approach.Section 5 presents a numerical example to examine the effect of the proposed resampling procedure.Conclusions are drawn in Sec. 6.

LMB Filter
A labeled RFS [8] is a finite-set-valued random variable, where each state x ∈ X is coupled with a unique label ∈ L.Here X denotes the target state space and L denotes the finite label space.

1) LMB Filter Prediction:
Suppose that the multitarget posterior density and birth model are both LMB RFSs with parameter sets π = {r ( ) , p ( ) : ∈ L} and π γ = {r ( )  γ , p ( ) γ : ∈ B}, respectively, where B denotes the birth label space.Then the multi-target predicted density is characterized by the parameter set where f , g ∫ f (x)g(x)dx denotes the standard inner product, P s (x, ) is the state dependent survival probability, f (x|x , ) is the single-target transition density for track .
2) LMB Filter Update: Suppose that the multitarget predicted density is an LMB RFS with parameter set π + = {r ( ) + , p ( ) + : ∈ L + }.Then the updated parameter set is π Z = {r ( ) , p ( ) : ∈ L + }, where F (L + ) denotes the class of finite subsets of L + , Θ I + is the space of mappings θ : is the state dependent detection probability, g(z|x, ) is the single-target measurement likelihood, and κ(z) is the clutter intensity.

Interval Analysis
where IR is the class of closed and connected subsets of the set R of real numbers, x and x denote its lower and upper bounds respectively.The width of x + x /2.An n-dimensional real interval vector or a box ). Elementary arithmetic and set operations for real intervals and boxes are summarized in [17], [18].

1) Inclusion Function
Inclusion functions may result in bounding boxes larger than necessary [16], as shown in Fig. 1.Consider a function h from R n to R expressed as a finite composition of the operators +, −, * , / and elementary functions (sin, cos, exp, sqrt, etc.).Its natural inclusion function [h] from IR n to IR is obtained by replacing each real variable x i , i = 1, • • • , n, with an interval variable [x i ] and each operator or function with its interval counterpart.If each of the variables x i , i = 1, • • • , n, occurs at most once in the formal expression of h, then [h] is minimal.

2) Contraction:
) .And the solution set of H is defined as Contraction is used in the definition of likelihood functions and the update step.Existing methods for building contractors, such as the Gauss elimination, Gauss-Seidel algorithms, linear programming and constraint propagation (CP), are elaborated in [17].In this paper, the CP technique is adopted.

Improved Resampling
For situations where not all elements in the target state vector can be measured, assume that a prior for the unmeasured elements is available.We propose to contract the measured elements of a box particle by the measurements and the unmeasured elements by the support of the prior, which may result in different levels of uncertainty in different dimensions of the box particle.So instead of being random, the choice of the dimension for the division of boxes to take place is proposed to be based on calculation and a new step is added to the traditional box particle resampling procedure, as is explained as follows.
We define a box resolution vector are user specified and application dependent and stand for an upper limit to the interval widths in the jth dimension of all box particles after resampling.has the same dimensionality and units as a single-target state vector does.
Step 1: Suppose that we are to resample where d i is calculated as: After dividing all boxes [ xi ], i = 1, • • • , J, and renormalizing their weights, we have a new set of box parti- is added to the traditional box particle resampling procedure to reduce the uncertainty each box represents.And the new step is explained below.
Step 2: To redivide and x denotes the ceiling of x.
Next, divide [x k ] in the first dimension into a set of a k1 boxes and denote the resulting set of boxes as Similarly, divide each box in Ξ k in the second dimension into a k2 boxes, then we have a renewed set of a k1 a k2 boxes as ] .Continue this operation in each of the remaining dimensions sequentially.And finally, have been redivided in this way, the set of resampled box particles is obtained as Hereafter, the proposed box particle resampling procedure is refered to as planned division resampling (PDR), on the other hand the resampling procedure in [13] is refered to as random division resampling (RDR).

Extended Inclusion and Volume Criteira
As is proposed in [13], the inclusion value for singletarget tracking measures if the true target state is contained in the support of the posterior pdf p(x).And the volume value measures the spread of particles that approximate p(x).We propose that the inclusion value for multi-target tracking indicate how many true targets are included in the supports of the posterior pdfs, and that the volume value for multi-target tracking be redefined as the mean volume of the supports of the posterior pdfs.
To compute the inclusion value, first we define τ ( ) t as an indicator function showing whether or not the true target state x * t is contained in the support of posterior pdf p ( ) (x): where C ( ) (α) is a credible set [13] with α 1 satisfying The computation of τ ( ) t is given in Appendix 1.
Note that in multi-target tracking one true target state may be contained in the supports of multiple posterior pdfs and that one posterior pdf may have multiple true target states contained in its support.We propose that the inclusion value for multi-target tracking be the number of one-to-one assignments between the true target states x * t , t = 1, • • • , T, and posterior pdfs p ( ) (x), ∈ L. The meaning of a one-toone assignment is threefold: a) x * t is assigned to p ( ) (x) if x * t is included in the support of p ( ) (x); b) Once x * t is assigned to p ( ) (x), it cannot be assigned to another posterior pdf p ( ) (x), , even though it may be included in the supports of both pdfs; c) Once p ( ) (x) is assigned to x * t , it cannot be assigned to another true target state x * t , t t, even though it may contain both true target states in its support.It is unnecessary to know exactly which target is assigned to which pdf.Only the number of assignments matters here.Therefore, a simple way to compute the inclusion value ρ is where n α is the number of nonzero t , ∈ L.Moreover, α t can be interpreted as the number of posterior pdfs that have x * t in their supports and β ( ) as the number of true target states that are in the support of p ( ) (x).
Suppose that the volume of p ( ) (x) is given as ν ( ) , then volume value for multi-target tracking is computed by where |L| is the cardinality of the label space L. The computation of ν ( ) is given in Appendix 1

Implementations
The implementations are based on the following assumptions.

A.1. The sensor provides biased1 box measurements
[z] ∈ IZ, where IZ denotes the set of closed and connected subsets of the point measurement space Z.
A.2.The survival probability P s and detection probability P d are state independent.

A.3.
In PF, the measurement noise v is modeled as zero mean white Gaussian.Assume that the stochastic uncertainty due to v is small so that p(v) can be approximated by p(v) ≈ U [ε] (v) in BPF, where [ε] is the measurement noise support.
2) Update: Suppose that the predicted parameter set is given as π + = {r ( ) + , p ( ) i=1 .Then p (θ) (x, |Z) in ( 7) can be approximated by {w ( ,θ)  U,i , [x ( ,θ) , where ) is a function that returns a contracted version of [x] under the constraints given by the measurement func- is the measurement likelihood for box measurements conditioned on box particles [13].The rest of the computations for (4)-( 8) are the same as in [9].
3) Resampling, Pruning and Estimation: Suppose the updated parameter set is given as π = {r ( ) , p ( ) : ∈ L} with p ( ) (x) approximated by {w ( )  i , [x ( ) i ]} J ( ) i=1 .The proposed resampling method has been explained in Sec.3.1.In PDR, we resample J ( ) 1 = r ( ) J max1 boxes for each updated track in Step 1, where J max1 is the maximum number of boxes per track after Step 1.In RDR, we resample J ( ) = r ( ) J max2 boxes for each updated track, where J max2 is the maximum number of boxes per track after the whole resampling.Note that if J max1 = J max2 , then the number of boxes per track after PDR is larger than that after RDR due to the redivision of boxes in Step 2 of PDR.After resampling, tracks with an existence probability r ( ) lower than a specific threshold are pruned.Box estimates of the target states are extracted in the same way as in [9], then we take the centers of the box estimates as the output point estimates.
3) Resampling, Pruning and Estimation: Suppose the updated parameter set is given as π = {r ( ) , p ( ) : ∈ L} with p ( ) (x) approximated by {w ( ) i , x ( ) i } J ( ) i=1 .For simplicity, the multinomial resampling method is used in this paper.We resample J ( ) = r ( ) J max particles for each updated track, where J max is the maximum number of particles per track after resampling.Tracks with r ( ) lower than a specific threshold are pruned.Estimates of the target states are obtained in the same way as in [9].

Numerical Examples
In this section, we first present the simulated scenario.Then we examine the performances of PDR-BPF-LMB, RDR-BPF-LMB and PF-LMB.We show that in the scenario considered PDR-BPF-LMB is better than the RDR-BPF-LMB, and that to achieve a similar level of performance, PDR-BPF-LMB is more time efficient than the PF-LMB.In this paper, we use cardinality statistics, the proposed extended inclusion and volume, and OSPA distance [20] to measure the performance of a multi-target tracking filter.
Consider a non-linear multi-target tracking scenario with a maximum of six targets appearing on the scene using range and bearings measurements.The number of targets varies in time due to births and deaths, and observations are subject to missed detections and clutter.The true trajectories are shown by different colors in Fig. 2 along with the start and stop positions of each track.The initial states, birth times and death times of the true trajectories are given in Tab. 1.The sensor is placed at the origin.The target state  A nearly constant turn model is considered.The state transition model for the PF is where The state transition model for the BPF is where The measurement function is given by The measurement noise v k is zero mean white Gaussian with a covariance R k = diag((σ 2 θ , σ 2 r ) T ).The parameters σ θ and σ r are given in Tab. 2. The sensor returns interval measurements with an interval length ∆ = (∆θ, ∆r) T , where ∆θ = π/90 rad is the length of bearing interval and ∆r = 60 m is the length of range interval.Thus a biased interval measurement at time k is defined as: The birth process is modeled by a multi-Bernoulli RFS with density parameters π γ = {r (i)  γ , p (i) γ } 3 i=1 , where p (i) γ (x) = N (x; m (i) γ , P γ ).Parameters of the birth process are given in Tab. 3. Clutter is modeled as a Poisson RFS with intensity κ k (z) = λ c U(z)/V, where λ c = 1.6 For the scenario considered, the resolution vector for the position and velocity vector xT k is = (200, 120, 200, 120) T .The MATLAB code for the PF-LMB filter is from Ba-Ngu Vo's website2.Note that no dynamic grouping or adaptive birth is implemented in this code, only the standard filter with static birth is considered.The implementation of the BPF-LMB involves the using of INTLAB toolbox [19], which provides users with built-in functions for interval calculations.
As is suggested by the scenario setting, the sensor returns no range rate measurements.Therefore, the velocities of a target state cannot be measured by the sensor, nor can they be directly contracted by the measurements.Assume that the target velocities vary within the range [-60,60] m/s.We use the measurements to contract the target positions and the assumed velocity range to contract target velocities.And the constraint propagation algorithm for target position and velocity vector [ xk ] is given in Tab. 4.
Figure 3 and 4 present the performance of the PDR-BPF-LMB for a single sample run.Different tracks are plotted in different colors.It can be seen that the estimates of individual target states are accurate in general, though a few outliers exist because of the high degree of ambiguity involved in interval calculations.No track label switching is observed in this sample run, but dropped or false tracks occur as expected.However, label switching does happen occasionally during target crossings due to process noise.In addition, it can also be seen in Fig. 3 that the farther the box measurements are from the sensor, the larger they are.This is due to the fact that corresponding to the same length of angle interval, the farther a circular arc is from the center, the longer it is.Less deviation from the ground truth is observed for the estimates of the trajectories that are closer to the origin, where the sensor is located, because the box particles representing those trajectories are contracted by smaller boxes and smaller boxes embody less uncertainty.Input:   In Fig. 5, a comparison of the cardinality statistics for the PDR-BPF-LMB with that for the RDR-BPF-LMB is de-picted3 .To have a fair comparison, we set J max1 = 40 for the PDR-BPF-LMB and J max2 = 130 for the RDR-BPF-LMB so that their average computation times for one single Monte Carlo run are very close.The average computation times are presented in Tab. 5.The cardinality statistics for the two filters are averaged over 100 Monte Carlo runs.As is shown in Fig. 5, the PDR-BPF-LMB produces reliable cardinality estimates throughout the simulation period, whereas the RDR-BPF-LMB only produces reliable cardinality estimates during the period from 1 s to 20 s and underestimates the cardinality for the rest of the simulation period.The reason for the performance gap is presented as follows.At every time step, after the addition of interval process noises to the predicted box particles in (15), the widths of intervals in every dimension of the box particles grow wider.During the update, the unmeasured target velocities are not contracted by measurements but contracted by the range of velocities.After a few recursions, because no redivision is applied in the RDR-BPF-LMB, the widths of velocity intervals in box particles will grow wider than the range of velocities.Consequently the results of the contraction of velocities in RDR-BPF-LMB will always be the range of velocities itself, namely [−60, 60] in this scenario.Because we use the center of a selected box particle to produce a point estimate, so the velocities will become zero in this case.This leads to wrong cardinality estimates for the RDR-BPF-LMB.It is obvious that the RDR-BPF-LMB is not suitable for this scenario.Furthermore, the cardinality statistics of the PDR-BPF-LMB is compared against that of the PF-LMB.As is shown in Fig. 6, the PDR-BPF-LMB performs better than the PF-LMB with J max = 1000, as well as the PF-LMB with J max = 2000, and slightly worse than the PF-LMB with J max = 5000 in terms of cardinality statistics.
Next, the inclusion and volume values for the PDR-BPF-LMB, RDR-BPF-LMB, and PF-LMB with J max = 1000, 2000 and 5000 are depicted in Fig. 7 and 8.We see that the inclusion values of PDR-BPF-LMB are generally lower than the inclusion values of RDR-BPF-LMB, and that the volume values of the former are smaller than the volume values of the latter.This is because PDR produces smaller boxes than RDR does due to the PDR redivision step, and smaller boxes (ultimately point particles) tend to result in lower inclusion and volume values, although PDR-BPF-LMB with J max1 = 40 and RDR-BPF-LMB with J max2 = 130 have very close computation times.We can also see that the PDR-BPF-LMB has higher inclusion values than the PF-LMB with J max = 1000 does, has a similar level of inclusion values than the PF-LMB with J max = 2000 does, and has generally lower 3Hereafter, different filters are shown in different colors and the choice of colors is irrelevant with the choice of colors for tracks in Figs.2-4.
inclusion values than the PF-LMB with J max = 5000 does.In terms of volume, the two BPF-LMB filters have higher values than the PF-LMB does.
Figure 9 compares the PDR-BPF-LMB, RDR-BPF-LMB, and PF-LMB with J max = 1000, 2000 and 5000 in terms of OSPA miss distance (p = 1, c = 100 m).The OSPA distance of the PDR-BPF-LMB is smaller than that of the RDR-BPF-LMB and is very close to those of the PF-LMB with different J max during most of the simulation time.But the BPF-LMB filters have significantly larger OSPA distance than the PF-LMB filter does during the periods from 1 s to 20 s when there are new births of targets, and during 50 s to 60 s when there are targets crossing.This is mainly due to the imprecision of the measurements.Based on our observations of the cardinality statistics, inclusion and volume values, and OSPA distance, it is reasonable to say that the overall performance of the PDR-BPF-LMB is better than that of RDR-BPF-LMB, and is comparable to that of the PF-LMB with J max = 2000.But the computation time for one single Monte Carlo run of the PDR-BPF-LMB is 49.16 s, about half of the time for the PF-LMB with J max = 2000, namely 106.77 s.So it is proved again that for a comparable level of performance, the BPF is more time efficient than the PF.

Conclusion
In this paper, we proposed an improved resampling method for the box particle filter by introducing the box resolution vector to the resampling procedure of the standard box particle filter.This vector is in fact a user specified and application dependent parameter that can be used to control the size of box particles after resampling.It can be applied to the situations where not all elements of the target state can be measured.For example, it can be applied to the situation described in this paper where the velocities of a target state cannot be measured due to the fact that the sensor returns only azimuth and range measurements but no range rate measurements.Additionally, we extend the inclusion and volume criteria from single-target tracking to multi-target tracking.In the multi-target tracking scenarios where the measurements are imprecise and biased, it is necessary to use inclusion and volume to measure the tracker's performance.It has been shown with a numerical example that the PDR-BPF-LMB outperforms the RDR-BPF-LMB and that the PDR-BPF-LMB reaches a comparable level of performance as the PF-LMB does with much lower computation time.

Fig. 1 .
Fig. 1.Image of a 2-dimensional box by a function f and two of its inclusion functions [ f ] and [ f ] * .[ f ] * is minimal.

Fig. 3 .Fig. 4 .Fig. 5 .
Fig. 3.The tracking result of a single sample run.Box measurements of target positions and clutter are shown as green rectangles, while estimates of targets are circles of other colors.
Average computation times for the considered filters for one single Monte Carlo run.