A Robust Interacting Multi-Model Multi-Bernoulli Mixture Filter for Maneuvering Multitarget Tracking under Glint Noise

In practical radar systems, changes in the target aspect toward the radar will result in glint noise disturbances in the measurement data. The glint noise has heavy-tailed characteristics and cannot be perfectly modeled by the Gaussian distribution widely used in conventional tracking algorithms. In this article, we investigate the challenging problem of tracking a time-varying number of maneuvering targets in the context of glint noise with unknown statistics. By assuming a set of models for the possible motion modes of each single-target hypothesis and leveraging the multivariate Laplace distribution to model measurement noise, we propose a robust interacting multi-model multi-Bernoulli mixture filter based on the variational Bayesian method. Within this filter, the unknown noise statistics are adaptively learned while filtering and the predictive likelihood is approximately calculated by means of the variational lower bound. The effectiveness and superiority of our proposed filter is verified via computer simulations.


Introduction
The aim of multitarget tracking (MTT) lies in estimating the kinematic state (e.g., position, velocity, acceleration) of each target within the region under surveillance from measurement data provided by sensing devices (e.g., radar, sonar, microphone) [1].Often, the number of targets may change over time as a result of the stochastic births, deaths, and spawns of targets.Moreover, merely part of the available measurements originates from targets, and the association maps between targets and measurements are not clear.The random finite set (RFS) [2] offers a natural formulation of the multitarget states and multitarget measurements.Using RFS modeling, the probability hypothesis density (PHD) filter [3], cardinalized PHD (CPHD) filter [4], and multi-Bernoulli filters [2,5] were developed as tractable MTT algorithms.Recently, driven by the development of multitarget conjugate priors, the advanced generalized labeled multi-Bernoulli (GLMB) filter [6], multi-Bernoulli mixture (MBM) filter [7], and Poisson MBM (PMBM) filter [8] were reported by researchers, which admit closed-form Bayesian recursion without certain approximations necessary in the PHD, CPHD, and multi-Bernoulli filters.By approximating a GLMB with a single term, an efficient alternative to the GLMB filter named the labeled multi-Bernoulli (LMB) filter [9] was proposed to make compromise between complexity and accuracy.Performance evaluations of the GLMB, PMBM, MBM, and LMB filters can be found in [10].
Most MTT algorithms assume that all targets follow the same dynamic model throughout.However, this assumption is too ideal in practice.For example, in a battlefield environment, a fighter jet has to carry out a series of tactical maneuvers to avoid being locked on by enemy weapon systems.In this case, these algorithms demonstrate poor performance.The jump Markov system or multi-model (MM) approach [11] has proven to be effective for maneuvering target tracking, in which the target can switch among a set of dynamic models in a Markovian fashion.By applying the MM approach in conjunction with RFS-based multitarget filters, numerous maneuvering MTT (MMTT) algorithms have been reported in [12][13][14][15][16][17][18][19][20].A fundamental premise of these algorithms is that the measurement noise is Gaussian distributed with known statistics.However, this is not necessarily the case in practical radar systems, where changes in the target aspect with respect to the radar can cause the apparent center of radar reflections to wander significantly.The random wandering of the apparent radar reflecting center dramatically increases the radar cross-section fluctuations [21][22][23][24][25], resulting in significant glint noise.It was found that the glint noise has a heavy-tailed probability density function (PDF) and cannot be perfectly modeled by the Gaussian distribution.In general, a priori knowledge of the glint noise statistics is not available.
The common approach for modeling glint noise is to exploit a type of distribution or the combination of multiple distributions.In [26], the Student's t (ST) distribution was used to model glint noise, whereas the mixture of Gaussian distributions was used in [27].In [28], the glint noise was modeled by the mixture of a Gaussian distribution and a Laplacian distribution.In particular, the ST distribution is immune to measurement outliers and has been widely used in RFS-based MTT [29][30][31], which can accurately characterize the tailed behavior by carefully selecting its degree of freedom (DOF) parameter.Based on the ST distribution and variational Bayesian (VB) method [32], a robust MMTT algorithm was proposed under the marginal distribution Bayes (MDB) filtering framework [33].The resulting ST-MM-MDB filter can adaptively learn the unknown scale matrix and DOF parameter of the ST distribution while filtering.In addition, the ST-MM-LMB filter was also reported in [34], which adopts a similar idea to estimate unknown noise statistics.However, as indicated in [35], these ST-based filters cannot estimate the DOF parameter accurately by means of limited measurement samples.To mitigate this drawback, attempts have been made to model glint noise using the multivariate Laplace (ML) distribution [36], which avoids the selection of the DOF parameter.A related Kalman filter was presented in [37], which shows higher estimation accuracy than the existing ST-based counterpart [38].Extensions and improvements of this filter were reported in [39,40].These works, however, are limited to the case where there is only a single non-maneuvering target to be tracked using one noisy measurement per scan.When multiple maneuvering targets are involved, a series of complicated factors, e.g., the stochastic appearances and disappearances of targets, unknown target-measurement association, misdetection (a target is detected by radar with a certain probability), and false alarms (spurious measurements not originating from any target) need to be systematically considered.It is infeasible to extend the single non-maneuvering case to the maneuvering multitarget case straightforwardly.
Motivated by the above discussions, in this article, we propose a robust MMTT algorithm in the context of glint measurement noise based on the MBM filter and ML distribution.Specifically, we assume a set of models for each Bernoulli component (BC) in the MBM representing a single-target hypothesis so as to capture the real-time maneuvers of targets.Moreover, we employed the ML distribution to model measurement noise and formulate the corresponding likelihood function as a Gaussian scale mixture form, leading to a hierarchical measurement model.Based on the MM assumption and such a hierarchical model, we derive a robust MBM filter for maneuvering targets using strategies from the interacting MM (IMM) algorithm (a celebrated MM approach with high cost effectiveness) [41], in which the unknown noise statistics are adaptively estimated by means of the VB method accompanied by IMM-MBM filtering.Since the involved predictive likelihood cannot be exactly calculated, we make use of the variational lower bound to obtain an approximate alternative.To summarize, the main contributions of this article are as follows: (1) Based on ML modeling and the VB method, a robust IMM-MBM filter is proposed to adaptively learn unknown glint noise statistics while filtering.
(2) A series of numerical simulations is performed to test the robustness of the proposed algorithm and compare its performance with the existing solutions.
The remainder of this article is organized as follows.Section 2 provides the necessary background knowledge.Section 3 reviews the standard (single-model) MBM filter and its Gaussian implementation.Section 4 details the robust MMTT algorithm proposed in this article.Section 5 presents the simulation results.Section 6 gives closing remarks.

Background
This section gives the notations used throughout this article and reviews several kinds of RFSs relevant to the development of our key results, as well as the ST and ML distributions.

Notation
In this article, tr(•) and det(•) denote the trace and determinant of a matrix, respectively; E[•] denotes the expectation operator; N(•; m, P) denotes the Gaussian PDF with mean vector m and covariance matrix P; IW(•; Ψ, ν) denotes the inverse Wishart (IW) PDF with scale matrix Ψ and DOF parameter ν; GIG(•; a, b, c) denotes the generalized inverse Gaussian (GIG) PDF with shape parameters a, b, and c; E(•; λ) denotes the exponential PDF with rate parameter λ; the superscripts −1 and ⊤ denote the inverse operation and transpose operation, respectively; exp denotes the natural exponential; log denotes the natural logarithm.

RFS Statistics
An RFS is essentially a set-valued random variable of which both the cardinality (number of elements) and the elements are random.Akin to conventional random variables, the randomness of an RFS is entirely described by its probability density.Several kinds of RFSs of interest are given next.
A Poisson RFS X is an RFS with the cardinality |X| being Poisson distributed with mean N = D(x)dx, where D(x) denotes the intensity, also known as the PHD, and the elements x ∈ X are independent and identically distributed in light of the spatial density D(x)/ N for any finite cardinality.The probability density of the Poisson RFS X is given by [2] (pp.366) Clearly, the Poisson RFS density f p (X) is entirely determined by the intensity D(x).
A Bernoulli RFS X can either be empty, i.e., X = ∅, with probability 1 − r, or be composed of a single element x, i.e., X = {x}, with probability r.Conditional on being nonempty, x is distributed according to the spatial density p.The probability density of the Bernoulli RFS X takes the form [2] (pp.368) A multi-Bernoulli (MB) RFS X is a union of a fixed number of independent Bernoulli RFSs, i.e., X = ⊎ i∈I X i , where ⊎ denotes the disjoint union, I is an index set for BCs in the MB RFS, and X i denotes the i-th BC.The probability density of the MB RFS X is given by [7] An MBM RFS X is a normalized and weighted sum of MB RFSs, whose probability density is given by [7] where ∝ stands for proportionality, J is the index set for MB components in the MBM, I j is the index set for BCs in the j-th MB RFS, and w i,j and f i,j (X i ) are the respective weight and probability density of the i-th BC in the j-th MB RFS.

ST and ML Distributions
The ST distribution is widely used to model glint noise.For a d-dimensional ST distributed random variable x, its PDF is given by [42] where Γ(•) denotes the Gamma function, ν denotes the DOF parameter, ∆ 2 = (x − m) ⊤ P −1 (x − m), m denotes the mean vector, and P denotes the scale matrix.Note that the covariance matrix of x is ν ν−2 P (ν > 2) instead of P. The tailed behavior of the ST PDF f st (x) is dominated by the DOF parameter ν.More specifically, the smaller the DOF parameter ν, the heavier the tail, and vice versa.When the DOF parameter ν tends to infinity, the ST PDF f st (x) becomes a Gaussian PDF.
When exact statistics of glint noise are not available, robust ST-based Kalman filters [43,44] have been proposed to jointly estimate the state vector and unknown noise parameters.However, as indicated in [35], these filters cannot estimate the DOF parameter accurately by means of limited measurement samples.Such a drawback motivates the employment of the ML distribution to model glint noise, which gets rid of the intractable DOF parameter.For a d-dimensional ML distributed random variable x, its PDF is given by [45] where ϱ(x) = (x − m) ⊤ P −1 (x − m), m is the mean vector, P is the covariance matrix, and B ν denotes the modified Bessel function of the second kind with order ν.As shown in [37], the ML PDF f ml (x) can be reformulated as a Gaussian-scale mixture form: where ζ is an auxiliary variable.

Gaussian MBM Filter
The MBM filter is closed to the Bayes prediction and update steps, in which the multitarget density at time t, t ∈ {k, k + 1}, conditional on the measurements up to time k is an MBM density of the form [46] An explanation of expression ( 9) is given as follows.To begin with, X t is the multitarget state RFS at time t, n t|k is the number of BCs, and a i = τ i , ℓ i , χ i τ i :k denotes the single-target hypothesis in regard to the i-th BC with τ i , ℓ i , and χ i t i :k = (χ i τ i , . . ., χ i k ) being the birth time, birth index, and data associations up to time k, respectively.At a specific time j, BC is associated with the p-th measurement, where m k is the number of measurements available at time k.Moreover, a = a 1 , • • • , a n t|k denotes a global hypothesis, which is composed of all single-target hypotheses, and A t|k denotes the collection of all global hypotheses.The i-th BC with single-target hypothesis a i has an associated weight w i,a i t|k and probability density f i,a i t|k of the form (2) characterized by the existence probability r i,a i t|k and spatial density p i,a i t|k .One recursion of the MBM filter is summarized below; see [7,8] for the mathematical proofs.
Proposition 1 (prediction).Assume that the filtering density at time k is of the form of (9) with t = k.Then, the predicted density involves n k+1|k BCs, where n k+1|k = n k|k + n b k+1 , and n k|k and n b k+1 are the number of BCs characterizing surviving targets and newborn targets, respectively.For each surviving BC i ∈ 1, • • • , n k|k , the predicted parameters are calculated by where p S,k (x) and φ k+1|k (x|x ′ ) denote the probability of survival and single-target state transition density, respectively.For each newborn BC i ∈ n k|k + 1, • • • , n k+1|k , we have Here, the single-target hypothesis reduces to a i = (τ i , ℓ i ) since there has not been a data association event for newborn targets yet.
Proposition 2 (update).The number of BCs does not change in the update step, so n k|k = n k|k−1 .
For each BC i ∈ 1, • • • , n k|k with single-target hypothesis a i , it generates a misdetection hypothesis and one association hypothesis for each measurement z . Note that Z k is composed of target-originated measurements and false alarms not originating from any target.To distinguish these two hypotheses, we make use of an ordered pair of integers (a i , p) with p ∈ {0, • • • , m k }.Under the misdetection hypothesis, we have w i,(a i ,0) k|k where p D,k (x) denotes the probability of detection.Under the measurement association hypothesis, the updated parameters are given by w i,(a i ,j) k|k where l k (z|x) denotes the single-target likelihood function and κ k (z) denotes the intensity of a Poisson RFS modeling false alarms.
It has been shown in [46] that the MBM recursion stated in the above propositions admits a closed-form solution under the following assumptions: Assumption 1.The single-target state transition density φ k+1|k (x|x ′ ) satisfies where F k and Q k are the state transition matrix and process noise covariance matrix, respectively.
Assumption 2. The single-target measurement likelihood function l k (z|x) satisfies where H k is the measurement matrix and R k is the measurement noise covariance matrix.
Assumption 3. The probabilities of survival and detection are constants, i.e., Assumption 4. Each BC in the predicted and filtering MBM densities has a Gaussian spatial density: The following propositions show how the parameters characterizing MBM densities are analytically propagated as time progresses.

Proposition 3 (prediction). The prediction of parameters characterizing the filtering density is
given in (10)-( 16).For each surviving BC i ∈ 1, • • • , n k|k , we have where For each newborn BC i ∈ n k|k + 1, • • • , n k+1|k , we have where the quantities r b,i−n k|k k+1 , m b,i−n k|k k+1 , and P b,i−n k|k k+1 are given birth model parameters.
Proposition 4 (update).For each BC i ∈ 1, • • • , n k|k with single-target hypothesis a i , the updated parameters for the misdetection case are given in (18) and (19), which simplify as w i,(a i ,0) k|k The updated parameters for the measurement association case are given in (20)-( 22), which simplify as w i,(a i ,j) k|k where q i,(a i ,j) k (z

P
i,(a i ,j) k|k After the prediction and update steps, the pruning of global hypotheses and BCs is further required to guarantee the computing efficiency.These procedures are detailed in [46] and, thus, will not be shown here for clarity.Moreover, the estimators given in [7] can be straightforwardly employed to extract multitarget state estimates.

Robust MMTT Under Glint Noise
The Gaussian MBM filter stated in the previous section relies on both the single-model and Gaussian noise assumptions.In some practical applications, these assumptions are no longer valid.On the one hand, target maneuvers are inevitable, and it is sufficient to describe the target motion mode by a single model.On the other hand, if a target of interest is observed by a radar system, non-Gaussian glint noise may occur as a result of changes in the target aspect toward the radar [47], and exact noise statistics are generally not available.In this section, we will show how the Gaussian MBM filter is tailored to accommodate maneuvering targets in the context of glint measurement noise with unknown statistics.Specifically, in order to capture the real-time maneuvers of targets, a set of models is assumed for possible motion modes of each BC representing a single-target hypothesis.The switching among models follows a homogeneous Markov process with transition probability t k+1|k (ξ ′ |ξ), where ξ ′ , ξ ∈ M and M denotes the discrete set of model labels.In addition, the ML distribution is employed to model the measurement noise.As a result, the likelihood function l k (z|x) in (24) now becomes For the sake of convenience, we make use of an equivalent form of the likelihood function where ML(•; m, P) denotes the ML PDF of the form (6). According to ( 7) and ( 8), l(z k |x k ) can be reformulated as The prior distribution of R k is selected as The rationality behind such a selection is that the IW distribution is the conjugate prior of a positive definite matrix [48].
The following propositions show how the unknown noise statistics are adaptively learned by means of the VB approximation accompanied by the IMM-MBM filtering.

Proposition 5 (prediction).
The predicted parameters for each surviving BC i ∈ 1, • • • , n k|k and each model with label ξ ∈ M are given by w i,a i k+1|k (ξ) = w i,a i k|k (ξ), where Here, µ i,a i k (ξ ′ |ξ) is the mixing weight and µ i,a i k+1|k (ξ) is the predicted model probability.For each newborn BC i ∈ n k|k + 1, • • • , n k+1|k , the corresponding parameters are given by Note that the birth parameters r b,i−n k|k k+1 (ξ), m b,i−n k|k k+1 (ξ), and P b,i−n k|k k+1 (ξ) are known a priori.

Proposition 6 (update).
For the i-th BC and model with label ξ ∈ M, under the misdetection hypothesis, the updated parameters are given by w i,(a i ,0) k|k The relative parameters for the measurement association hypothesis are given by w i,(a i ,j) k|k i,(a i ,j) k|k where m i,(a i ,j) k|k (ξ), P i,(a i ,j) k|k (ξ), and q i,(a i ,j) k (z Here, [•] = VB(•) denotes the VB approximation and N denotes the number of VB iterations.
The updated model probability for the misdetection hypothesis case denoted by µ i,(a i ,0) k (ξ) equals µ i,a i k|k−1 (ξ), while the updated model probability for the measurement association hypothesis case denoted by µ i,(a i ,j) k Moreover, the model-dependent quantities are then combined to obtain the final estimates of the i-th BC according to Note that the pruning and estimate-extraction procedures are needed after the prediction and update steps.They are the same as in the standard MBM filter.In addition, the above IMM-MBM filter is not restricted to the linear dynamic and measurement models.It can accommodate nonlinear models using approximate strategies such as the unscented transform [49] and spherical-radial cubature rule [50].Now, we specify the VB approximation adopted in (72).Without loss of generality, we consider a notationally convenient form: For the hierarchical measurement model given in ( 49)-( 51), an analytic solution of the joint posterior PDF p(Λ k |z 1:k ) is not available, where Λ k ≜ {x k , R k , ζ k } and z 1:k denotes the measurement sequence accumulated to time k.Consequently, the VB method is exploited to obtain an approximate solution: Each term of the approximated posterior PDF satisfies where Λ −ϑ k ∪ {ϑ} ≜ Λ k and c ϑ is a ϑ-independent constant.To address the mutual coupling of the variational parameter, a fixed-point iteration scheme is adopted to obtain an approximation of q(ϑ) by iteratively solving (81).This means that q(ϑ) is updated as q (d+1) (ϑ) by exploiting q (d) (Λ −ϑ k ) to compute the expectation in (81) at the (d + 1)-th iteration.Using Bayes' theorem, the joint PDF p(Λ k , z 1:k ) can be formed as Proposition 7. Let ϑ = x k , and plugging (82) into (81), q (d+1) (x k ) is updated as a Gaussian PDF, i.e., where Here, n z denotes the dimensions of measurement vector z k and I n denotes the n × n identity matrix, and the Kalman gain G (d+1) k|k is calculated by where the modified measurement noise covariance matrix Proof.See Appendix A.
To implement the fixed-point iteration, the required expectations need to be calculated as follows.Using (83), the auxiliary parameter U (d+1) k|k in (92) is calculated by Employing (88), the expectation Utilizing (93), the expectation In addition, abiding by [29,31], the predictive likelihood q k|k (z Since the first term on the right-hand side of (99) is minimized by the VB method, L(q(x k ) q(R k )q(ζ k )) can be treated as the variational lower bound of log q k|k (z k ).As a consequence, q k|k (z k ) can be approximately computed as where The explicit expression for L(q(x k )q(R k )q(ζ k )) can be found in Appendix D. A summary of the VB function of the form (79) is given in Algorithm 1, where Rk is the nominal measurement covariance matrix.
Algorithm 1 A summary of the VB function.

Numerical Studies
A 2D tracking scenario involving four maneuvering targets is considered to demonstrate the performance of the proposed MMTT algorithm (hereafter refereed to as the ML-IMM-MBM filter).The target kinematic state x k is composed of position [p x,k , p y,k ] ⊤ and velocity [ ṗx,k , ṗy,k ] ⊤ , i.e., x k = [p x,k , ṗx,k , p y,k , ṗy,k ] ⊤ .All targets follow the state dynamics: with different turn rates ω ∈ {−10 • /s, 0, 10 • /s}, where T = 1s is the scan period and w k is the zero-mean Gaussian process noise with covariance matrix More specifically, ω = 0 corresponds to a constant velocity model indexed by M1 with σ w = 10m/s 2 , while ω = −10 The simulation lasted 80 s, corresponding to 80 scans.The dynamics of each target is given as follows.Target 1 appears at time k = 1 s and disappears at time k = 60 s.It follows model M2 during [1,25] s, model M1 during [26,40] s, and model M3 during [41,60]  A radar observes the targets of interest and provides both the range and angle according to the measurement equation: where [p s,x , p s,y ] ⊤ = [1000 m, 2000 m] ⊤ denotes the location of the radar and v k is the glint measurement noise satisfying [31] Here, 0 n denotes the n-dimensional zero vector, Rk = diag [10m, 5 • /s] ⊤ 2 is the nominal measurement covariance matrix, τ = 100 is a scale factor, and γ = 0.1 is the glint probability.The probability of detection is p D = 0.9.Apart from target-originated measurements, at each scan, an average of 10 false alarms are also received by the sensor.
We first tested the ML-IMM-MBM filter with 100 Monte Carlo (MC) trials.In each trial, the true target tracks remained unchanged, whereas the sensor measurements were randomly generated.In the ML-IMM-MBM filter, the prior parameters were configured as λ 0 = 1, ν 0 = β, and Ψ k = β Rk , where β is a user-defined parameter satisfying β > n z + 1.The thresholds for pruning, merging, and state extraction were set as 10 −5 , 4, and 0.4, respectively.In addition, the ML-IMM-MBM filter uses a maximum number of global hypotheses N h = 100 and adopts the unscented transform [49] to tackle measurement model nonlinearity.We evaluated the filter performance using the generalized optimal subpattern assignment (GOSPA) metric [51] with α = 2, p = 2, and c = 10m.The GOSPA error (GOSPAE) along with its decomposed localization error (LE), missed target error (MTE), and false target error (FTE) for different settings of the VB iteration number N are presented in Figure 2. It can be seen that the ML-IMM-MBM filter has a similar performance when N ≥ 4.This can be attributed to the fast convergence of the VB approach [32].The relevant results for different settings of the prior parameter β are shown in Figure 3.It can be seen that the ML-IMM-MBM filter behaves almost identically for varying parameter β, indicating excellent robustness to the prior parameter β.We then compared the ML-IMM-MBM filter with the ST-MM-MDB filter, ST-MM-LMB filter, and the IMM-MBM filter with the true measurement noise covariance matrix (hereafter referred to as the IMM-MBM-T filter) with 100 MC trials.In the ML-IMM-MBM filter, the prior parameter β was set as β = 8.Necessary parameters for the ST-MM-MDB filter and ST-MM-LMB filter were configured consistent with [33,34], respectively.In all filters under study, the VB iteration number N was set as N = 6.Figures 4 and 5 show the performance of these filters in terms of cardinality estimates and the GOSPA metric, respectively.It is shown that the ML-IMM-MBM filter is comparable to the MM-MBM-T filter, but outperforms the ST-MM-MDB and ST-MM-LMB filters.This can be attributed to two reasons.First, the MBM filter is an exact solution to the Bayes multitarget recursion (centerpiece of RFS-based MTT), while the MDB filter [52] and LMB filter [9] are approximate alternatives.Second, the ML-IMM-MBM filter models glint noise by the ML distribution, which gets rid of the selection of the DOF parameter, which is necessary in the ST-MM-MDB filter.The GOSPAE and its decomposition of different filters for varying scale factors and glint probability are shown in Figure 6 and Figure 7, respectively.It can be observed that the glint probability has a more significant influence on performance, and all the results indicate the superiority of the proposed ML-MM-MBM filters over the existing ST-MM-MDB and ST-MM-LMB filters.

Conclusions
In this article, we proposed a robust algorithm to address the problem of tracking multiple maneuvering targets with random appearances and disappearances from radar measurements corrupted by glint noise.Specifically, in order to capture the real-time maneuvers of targets, a set of models was assumed for the possible motion modes of each single-target hypothesis.Furthermore, the ML distribution was employed to model measurement noise, and the relevant likelihood function was formulated by a Gaussian scale form, leading to a hierarchical Gaussian model.Using this model, a robust IMM-MBM filter was derived based on the VB method, which can adaptively learn unknown noise statistics while filtering.In particular, the involved predicted likelihood was approximately calculated by means of the variational lower bound.The simulation results showed that our proposed algorithm not only possesses strong robustness to its freedom parameters, but also outperforms the existing ST-MM-MDB and ST-MM-LMB filters.For future work, we will extend the proposed algorithm to accommodate amplitude information under different SNRs and a network of cooperative sensors with limited sensing ranges.