Multisensor Multiobject Tracking With Improved Sampling Efficiency

Passive monitoring of acoustic or radio sources has important applications in modern convenience, public safety, and surveillance. A key task in passive monitoring is multiobject tracking (MOT). This paper presents a Bayesian method for multisensor MOT for challenging tracking problems where the object states are high-dimensional, and the measurements follow a nonlinear model. Our method is developed in the framework of factor graphs and the sum-product algorithm (SPA) and implemented using random samples or “particles”. The multimodal probability density functions provided by the SPA are effectively represented by a Gaussian mixture model (GMM). To perform the operations of the SPA with improved sample efficiency, we make use of particle flow (PFL). Here, particles are migrated towards regions of high likelihood based on the solution of a partial differential equation. This makes it possible to obtain good object detection and tracking performance even in challenging multisensor MOT scenarios with single sensor measurements that have a lower dimension than the object positions. We perform a numerical evaluation in a passive acoustic monitoring scenario where multiple sources are tracked in 3-D from 1-D time-difference-of-arrival (TDOA) measurements provided by pairs of hydrophones. Our numerical results, obtained by processing synthetic and real data, demonstrate favorable detection and estimation accuracy compared to state-of-the-art reference techniques.


I. INTRODUCTION
MOT is an important capability for a variety of applications, including surveillance, autonomy, and marine mammal research.MOT is a high-dimensional nonlinear filtering problem complicated by measurement-origin uncertainty (MOU), i.e., the associations between measurements and objects, and an unknown number of objects to be tracked.In this paper, we develop a sequential Bayesian MOT framework for the particularly challenging scenarios where object states are highdimensional and measurement models are nonlinear.We expect that our approach is particularly useful for the passive monitoring of acoustic [1] or radio [2] sources in 3-D.

A. State-of-the-Art
Traditional methods for MOT include probabilistic data association (PDA) [3], multi-hypothesis tracking (MHT) [4], and methods based on random finite sets (RFS) [5]- [8].Most of these traditional approaches suffer from a computational Wenyu Zhang is with the Department of Electrical and Computer Engineering, University of California San Diego, La Jolla, CA, USA (e-mail: wez078@eng.ucsd.edu).
Florian Meyer is with the Scripps Institution of Oceanography and the Department of Electrical and Computer Engineering, University of California San Diego, La Jolla, CA, USA (e-mail: flmeyer@ucsd.edu).
complexity that is exponential in important system parameters, including the number of measurements, objects, and sensors.MOT methods that are scalable with respect to these parameters have been recently developed in the framework of factor graphs and the SPA [9]- [13].Factor graphs represent statistical independencies of random variables.The SPA is known to provide accurate solutions to high-dimensional Bayesian estimation problems efficiently.In particular, by performing local operations ("messages") on the factor graph, accurate approximations ("beliefs") of the marginal posterior pdfs of unknown states [14] are computed.SPA-based methods are versatile and have been successfully applied to a variety of applications, including cooperative localization [15]- [18], simultaneous localization and mapping (SLAM) [19]- [21], and focalization for underwater localization [22].
To calculate messages that, due to nonlinearities in the system model, cannot be evaluated in closed form, SPAbased methods for MOT typically rely on particle-based computations that closely follow the bootstrap particle filter (BPF) [23]- [25] and rely on importance sampling.A known drawback of this approach is that it typically fails in tracking problems where (i) the states of individual objects have dimensions higher than four, (ii) measurements are very informative compared to the predicted/prior pdfs.In particular, tracking of objects in 3-D Cartesian coordinates or employing sensors that yield low measurement variance often leads to a failure of particle-based computations due to particle degeneracy [26].The particle degeneracy problem is related to the fact that predicted pdfs are used as proposal pdfs for sampling.Since predicted pdfs can have completely different shapes than the posterior pdfs, this sampling strategy is highly inefficient, i.e., few or none of the generated particles are suitable to represent the posterior pdfs.Fig. 1 shows an example of particle degeneracy in a 3-D tracking scenario with a single object and a single TDOA measurement.Particle degeneracy is exacerbated in high dimensional problems and in problems with low measurement variance.In particular, it can lead to the unwanted behavior that filter performance degrades as measurement variance is reduced.As the dimension of the problem increases, or measurement variance is reduced, the likelihood function becomes "peakier," and it becomes more unlikely that a particle, sampled from the prior distribution, is located in a region of high likelihood.
Sometimes particle degeneracy can be avoided by using vast numbers of particles or by implementing regularization strategies [16], [25], [27].A straightforward approach to improve sampling efficiency and avoid particle degeneracy is to design proposal pdfs that are similar to the posterior pdfs [23]- [25].However, finding a distribution that is easy to sample from and simultaneously similar to posterior pdfs The object, shown in black, is located on the hyperboloid.Note that any other location on the hyperboloid will lead to the same measurement in the case without noise.(a): The prior pdf, f (x), is Gaussian, with the mean depicted as a big blue dot.2000 particles, shown as small light blue dots, are drawn from the prior distribution.On the right, the prior and posterior pdfs for the case with measurement noise are shown in three separate 2-D plots.Each of these plots is obtained by depicting the prior and posterior pdfs along the three axes of the coordinate system.(b): After importance sampling, as performed by the conventional "bootstrap" particle filter, only a single particle has a nonzero weight.This single particle does not accurately represent the posterior pdf p(x|z1) for future processing, e.g., of a measurement, z2, provided by a second sensor.[23]- [25] is often challenging.To improve samping efficiency, adaptive importance sampling can be employed [28], [29].In particular, auxiliary particle filters use a delayed resampling strategy to increase the number of particles with significant weights after importance sampling [30].This approach can improve sampling efficiency but can only be applied in combination with a prediction step, which may be unavailable for newly introduced object states in MOT scenarios.Furthermore, multiple particle filtering [31], similar to a particle-based implementation of the SPA, aims to increase sample efficiency by exploiting factorization of the underlying statistical model [28].Since it relies on a suitable factorization of the conditional posterior, its applicability is restricted.Incorporating sequential Markov chain Monte Carlo (SMCMC) methods into particle filters [32], [33] is another general approach for nonlinear sequential Bayesian estimation, which is known to be very computationally expensive in high-dimensional state spaces.An alternative approach to improve sampling efficiency in sequential estimation is to perform the update step of an unscented Kalman filter [34] and use the resulting Gaussian pdf as a proposal pdf for particle filtering.The unscented particle filter [35], [36] combines this idea with a Gaussian mixture representation of predicted and posterior pdfs.To the best of our knowledge, the unscented particle filtering approach has not yet been extended to problems with MOU and an unknown number of states to be estimated.
PFL [37]- [41] is a promising strategy for challenging nonlinear estimation problems that has recently received significant attention [42].It has the potential to avoid particle degeneracy due to its ability to actively move particles representing a prior or predicted pdf to locations of high likelihood. 1   1 Conventional strategies that rely on resampling also move particles but do so in a more passive way.In particular, in conventional strategies, after particles are randomly drawn, only those that correspond to locations of high likelihood remain after resampling.In challenging problems, this is prone to particle degeneracy, i.e., the number of remaining particles can become too low to be a representative description of the underlying posterior pdf.
This active motion is illustrated Fig. 2. For PFL a homotopy function is defined to formulate a pdf that can be smoothly deformed from the predicted pdf (or prior pdf) to the posterior pdf.PFL then makes use of the homotopy function to incrementally move a set of particles sampled from the predicted pdf.In particular, a partial differential equation (PDE) for particle velocity is obtained by combining the homotopy function with the Fokker-Planck equation.The particle velocity solution to the PDE can be discretized and used as a transport equation for particle migration.After migration, the set of particles represents the posterior pdf.There are two different types of PFL resulting in the exact Daum and Huang (EDH) filter and the localized exact Daum and Huang (LEDH) filter.In the EDH filter, the PFL equations are computed once for the mean of all particles.In contrast, in the computationally more demanding LEDH filter, the PFL equations are computed for each particle individually.PFL has been demonstrated to achieve a superior performance complexity tradeoff compared to existing approaches that aim at improving sampling efficiency [42].As other particle filtering approaches, PFL is highly parallelizable [38]- [40] and thus ideal for real-time processing on graphical processing units (GPUs).
Traditional PFL methods avoid importance sampling and can only provide an approximate representation of posterior pdfs in general nonlinear systems [37]- [41].Nevertheless, these "proposal-free" methods often lead to accurate estimation results at a significantly reduced computational complexity compared to BPF [42].Recently, it has been shown that PFL can be described by an invertible mapping and can thus be used as a measurement-driven proposal pdf for importance sampling [42].The resulting invertible PFL filter [42] is an asymptotically optimal approach to nonlinear filtering that avoids particle degeneracy and can provide accurate estimation results in high-dimensional and nonlinear problems.
A significant limitation of the PFL filter presented in [42] is that it assumes that the prior or predicted pdfs follow Gaussian distributions.It is thus unsuitable for problems that involve multimodal pdfs.For problems where the measurement noise follows a Gaussian mixture pdf, [43] introduces the Gaussian sum PFL filter.Here, the means of the Gaussian mixture components are updated by performing an update step similar to the LEDH.On the other hand, the covariance matrices of the components are updated by extended Kalman filters that also run in parallel.An extension of [43] to the case where both driving noise and measurement noise are distributed by a Gaussian mixture pdf is presented in [44].Here, invertible flow is used for particle weight update in an importance sampling step.For problems where both driving noise and measurement noise can be multimodal, [45] combines the invertible PFL with a SMCMC method that relies on the Metropolis-Hastings approach, i.e., a Metropolis-Hastings kernel is constructed using a PFL algorithm based on a GMM.However, aforementioned PFL approaches that can represent multimodal pdfs [43]- [45] are unsuitable for MOT since neither model MOU nor an unknown number of states to be estimated.For the cooperative localization problem, a method that relies on invertible PFL is presented in [46], [47].This method is not suitable for the more challenging multiobject tracking problems since it can only be applied to problems without measurement origin uncertainty, known number of states to be estimated, and posterior pdfs with simple, unimodal shapes.A variant of the PFL filter has been proposed for MOT [8].In particular, EDH and LEDH variants of the singlesensor δ-Generalized Labeled Multi-Bernoulli filter [7] with invertible flow are presented.These approaches are unsuitable for multiobject tracking problems where measurements are provided by multiple sensors.

B. Contributions, Paper Organization, and Notation
We develop a method for multisensor MOT with improved sample efficiency that can be used in scenarios with highdimensional object states and informative measurements.Of particular interest are multisensor MOT problems, where inexpensive sensors are used and the tracking of objects in Cartesian coordinates is impossible based on the measurements provided by a single sensor.In this type of tracking problems, the measurement of a single sensor typically has a lower dimension than the positions of objects.Consider a scenario where object positions are 3-D, but sensors only provide 1-D measurements, e.g., times of arrival (TOAs), time differences of arrival (TDOAs), or directions of arrival (DOAs).In this type of MOT problem, prior or predicted pdfs can have complicated multimodal shapes, e.g., spheres, hyperboloids, or cones at the initial step after the appearance of a new object.As an example, Figs. 1 and 2 show the hyperboloid-shaped pdfs resulting from a TDOA measurement model in a 3-D tracking scenario.
Our approach performs SPA-based message passing on the factor graph for scalable multisensor MOT developed in [10].The messages of the SPA are computed sequentially across sensors.To improve sampling efficiency, we embed invertible particle flow into SPA computations.For the evaluations of particle weights, invertible particle flow relies on a Gaussian representation of the prior or predicted pdf at the onset of the flow.To represent beliefs of object states with complicated non-Gaussian shapes, such as, e.g., hyperboloids, as in the example in Figs. 1 and 2, we make use of a GMM representation that is known to be asymptotically optimal [48].Combining a GMM with an efficient sampling approach to represent pdfs with complicated shapes in high dimensions is inspired by unscented particle filtering [35], [36].A general proposal pdf that takes MOU into account and consists of a mixture of pdfs related to different particle flows is developed.The resulting computations are asymptotically optimal.In particular, since particles are migrated towards regions of high likelihood, an accurate approximation of SPA messages with a relatively small number of particles is obtained.
The technical novelty of the proposed method lies in a new method for multitarget tracking that can achieve a superior runtime-estimation accuracy tradeoff in nonlinear and highdimensional problems by improving sampling efficiency.The improved tradeoff is obtained by carefully embedding invertible PFL.In particular, to address MOU, association probabilities are computed by performing parallel flows, one for each component of the GMM and each possible measurement-to-object association.The particles of the parallel flows are weighted based on association probabilities and combined into a mixture of flows.The mixture of flows provides samples of the proposal for importance sampling.Since all flows are invertible, it is possible to evaluate the proposal pdf represented by the mixture of flows at each particle.Thus, the resulting SPA-based computation of beliefs is asymptotically optimal in the sense that the resulting particle representation of the beliefs provided by the SPA is arbitrarily accurate for an increasingly large number of Gaussian components and a number of particles.Our method, for the first time, performs MOT with probabilistic data association based on PFL.
We further demonstrate that the proposed multisensor MOT can outperform reference methods based on conventional ("bootstrap") and unscented particle filtering in a 3-D passive source tracking scenario.In particular, in the considered realistic source tracking scenario, graph-based MOT based on conventional particle filtering [9] cannot provide acceptable estimation accuracy.The also considered, yet unpublished, implementation of graph-based MOT based on unscented particle filtering, has a lower estimation accuracy but a higher runtime compared to the proposed method that embeds invertible PFL.Key contributions of this paper are as follows.
• We develop a graph-based MOT method based on a GMM and invertible PFL for challenging scenarios with high-dimensional object states and arbitrarily shaped posterior pdfs.
• We demonstrate that the proposed method can significantly outperform reference techniques in a challenging 3-D passive source multisensor MOT scenario and show tracking results using real passive acoustic data.
This paper advances over the preliminary account of our method provided in the conference publication [49] by (i) introducing a GMM for multimodal state distribution with dynamic kernel resampling; (ii) considering the multisensor MOT problem; (iii) presenting an improved proposal distribution based on PFL; (iv) performing a comprehensive numerical evaluation in a 3-D passive source tracking scenario; and (v) applying the proposed method to an underwater acoustic dataset 2 .Contrary to the approach presented in [8], the proposed method is suitable for multisensor scenarios.In addition, SPA-based processing makes our approach scalable with respect to relevant system parameters.
Notation: Random variables are displayed in sans serif, upright fonts and their realizations in serif, italic fonts.Vectors and matrices are denoted by bold lowercase and uppercase letters, respectively.For example, a random variable and its realization are denoted by x and x, respectively, and a random vector and its realization by x and x, respectively.Furthermore, x and x T denote the Euclidean norm and the transpose of vector x, respectively; and ∝ indicates equality up to a normalization factor.N (x; x * , P ) denotes the Gaussian pdf (of random vector x) with mean x * and covariance matrix 2 More details on the application of the proposed method to the problem of tracking multiple whales underwater by performing TDOA measurements of their echolocation clicks, is presented in the companion paper [50] P .The trace of matrix M is denoted as Tr{M }.Finally, 1(a) denotes the indicator function of the event a = 0, i.e., 1(a) = 1 if a = 0 and 0 otherwise.

II. REVIEW OF INVERTIBLE PFL
We consider the general setting of calculating the posterior pdf based on Bayes' rule f (x|z) ∝ f (x)f (z|x) with the state of interest x and the observed (fixed) measurement z.If the prior pdf f (x) follows a Gaussian distribution and the likelihood function f (z|x) represents a linear measurement model z = Hx + v with Gaussian measurement noise v, the posterior pdf f (x|z) also follows a Gaussian distribution.In this special case, the mean and covariance of the Gaussian posterior pdf f (x|z) can be calculated in closed form by the Kalman update step [51].
If the measurement model is nonlinear, e.g., z = h(x) + v, a popular approach is to approximate the posterior pdf f (x|z) by a set of N p weighted particles {(x (i) , w (i) )} Np i=1 .Note that the weights are normalized to one, i.e., Np i=1 w (i) = 1 and can be computed based on the importance sampling principle [25] as follows Here, the proposal pdf q(x|z) is used to sample the particles It is an arbitrary pdf that has the same support as f (x|z).Importance sampling is asymptotically optimal if q(x|z) is "heavier tailed", i.e., less informative, than f (x|z) [24].In particular, importance sampling can provide an approximation of f (x|z) that can be made arbitrarily good by choosing N p sufficiently large [25].For N p fixed, if the proposal q(x|z) is "more similar" to the posterior f (x|z) [24], importance sampling is "more accurate".
A simple choice for the proposal pdf used in the update step of the conventional "bootstrap" particle filter [23], [25] is the prior pdf f (x).However, for a feasible number of particles N p and most choices of the proposal pdf, importance sampling can suffer from particle degeneracy [26].Particle degeneracy is especially severe if the state x is high-dimensional and the measurement z is informative (i.e., the likelihood function has narrow and sharp peaks).

A. Particle Flow (PFL)
PFL [37]- [40] is an approach that aims at avoiding particle degeneracy.Here, particles are smoothly migrated in the state space from a representation of the prior pdf to a representation of the posterior pdf by solving a PDE.Let us introduce the homotopy function π λ (x) = f (x)l λ (x) where λ ∈ [0, 1] is the pseudo time of the flow process and l(x) = f (z|x) is the likelihood function.Note that for λ = 1, the homotopy function is equal to the unnormalized posterior pdf, i.e., π(x) π 1 (x) = f (x)l(x).The log-homotopy function is then given by [37], [38], The log-homotopy function is a pseudo posterior pdf in the log domain that defines a smooth and continuous deformation from φ(x, 0) = log f (x) to φ(x, 1) = log π(x).This deformation describes the PFL process.
It can be shown that the stochastic process defined by homotopy function π λ (x) satisfies the Fokker-Planck equation [39]- [41].Combining the Fokker-Planck equation for the zerodiffusion case with (2) results in the following PDE [39], [40] where ζ(x, λ) = dx dλ describes particle velocity (samples of x) as the pseudo time λ increases from 0 to 1, i.e., as the homotopy function is deformed from the prior pdf to the posterior pdf.This migration is referred to as the PFL.

B. PFL Update Step
If f (x) and l(x) are Gaussians or in another exponential family, then an exact and closed form solution for (3) is available.The EDH filter [39], [52] makes use of this closed-form solution in its update step.More precisely, let f (x) = N (x; x * 0 , P ) and z = Hx+v be a linear measurement model with measurement noise v ∼ N (v; 0, R).The exact flow solution [39], [52] where we introduce Note that in (5), z is the observed and thus fixed measurement.This solution is extended to the nonlinear measurement model z = h(x) + v by performing a suboptimal linearization step.In particular, in a first-order approximation, a Jacobian matrix is computed, i.e.H(λ) = ∂h(x) ∂x x=x * λ where x * λ is the approximated mean of x at pseudo time λ.
In a practical implementation, we calculate ζ(x, λ) at N λ discrete values of λ, i.e., 0 = λ 0 < λ 1 < ... < λ N λ = 1, to perform the PFL.Here, we first sample N p particles x . Next, at each discrete pseudo time step l ∈ {1, . . ., N λ }, particles are migrated according to for all i ∈ {1, . . ., N p }. Here, the linearized flow solution ζ(x (i) λ l−1 +b l is computed based on A l and b l given by (cf.( 4) and ( 5)) Note that e l = h(x * λ l−1 , 0) − H l x * λ l−1 is the error of linearization and that the linearized measurement model, H l , is computed based on the mean of the last step l − 1, i.e., The mean at which the measurement model is linearized is typically propagated in parallel to the particles, i.e., that approximately represent the unnormalized posterior pdf π(x) are finally obtained.Pseudocode for the PFL update step is provided in Algorithm 1. PFL based on a linearized model has no optimality guarantees.However, it has been demonstrated numerically to typically provide an accurate representation of the posterior pdf f (x|z) [37]- [42].In what follows, the PFL related to the measurement z as defined by ( 6)-( 8), is denoted as x 0 − → z − → x 1 or for notational convenience in future derivations as x 0 − → z − → x.
Calculate the linearized measurement model H l according to (9); 5 Compute A l and b l according to ( 7) and ( 8); As an alternative to EDH PFL as discussed above, LEDH PFL has been introduced in [42], [53].Here, a linearization is performed for each particle location individually instead of only at the mean particle location, i.e., individual flow parameters are computed to migrate the particles along their individual flows.Although the LEDH flow usually outperforms the EDH flow, it suffers from considerable computational complexity since the main computational burden is related to calculating the flow parameters.
In Fig. 2, it is depicted how PFL actively move particles representing a prior or predicted pdfs to locations of high likelihood.Active motion of particles leads to a significantly improved approximation of the posterior pdf compared to conventional importance sampling shown in Fig. 1.

C. Importance Sampling with Invertible Flow
PFL can be used to compute a measurement-driven proposal pdf q(x|z) for importance sampling (cf.( 1)) to perform asymptotically optimal estimation [42].Here, the mapping as performed by PFL x 0 − → z − → x is invertible, i.e., there exists an invertible mapping of the particles after the flow x (i) Np i=1 to the particles x (i) 0 Np i=1 if certain constraints on the differences of consecutive discrete pseudo times λ l −λ l−1 , l ∈ {1, . . ., N λ } are satisfied [42].
By exploiting the invertible mapping, the proposal pdf resulting from PFL can be evaluated at the particles as [42] Here, the "mapping factor" θ is defined as By plugging (10) into (1) the weight of the particle x (i) is obtained as The resulting particle set {x (i) , w (i) } Np i=1 is an asymptotically optimal sample representation of the posterior pdf f (x|z) that can often provide accurate estimation results in nonlinear and high-dimensional estimation problems even if the number of particles is moderate [42].Pseudocode for importance sampling with invertible PFL is provided in Algorithm 2. Note that since the PFL used for the measurement-driven proposal pdf is typically based on the EDH filter update step, a Gaussian prior pdf is assumed.
Algorithm 2: Importance Sampling with Invertible Flow Perform weight update according to (12), i.e., Np i=1 An approximate Gaussian representation of this posterior distribution can be subsequently obtained by applying Algorithm 3 which calculates a mean x * 1 and a covariance matrix P 1 from the unnormalized weighted particles x Note that in Algorithm 2, the same mapping factor θ is used to calculate all particle weights.If Algorithm 3 is applied after Algorithm 2, this factor is irrelevant since all weights are normalized in Algorithm 3.However, making use of θ is important if multiple flows are performed in parallel, as will be discussed in Section III. ; 5 Compute mean and covariance matrix from particles, i.e., Note that instead of a particle-based covariance matrix computation as performed in Line 7 of Algorithm 3, an extended or unscented Kalman update step can be used [42].For a large number of particles N p , a particle-based computation is more accurate than a computation based on the extended or unscented Kalman update step.

III. GAUSSIAN MIXTURE REPRESENTATION FOR NONLINEAR ESTIMATION IN HIGH-DIMENSIONS
In this section, we use GMM representations [54]- [56] for nonlinear estimation in high dimensions.In particular, we developed methods for updating the parameters of GMMs based on PFL.This approach is suitable for high-dimensional pdfs that are multimodal and thus relevant for estimation problems in multiobject tracking and SLAM [19].

A. GMM Importance Sampling with Invertible Flow
As discussed in the previous Section II, for the evaluations of particle weights, invertible particle flow relies on a Gaussian representation of the prior pdf at the onset of the flow.In challenging multisensor MOT problems, the complicated multimodal shapes of prior and posterior pdfs (see, e.g., the hyperboloid-shaped posterior pdf in Fig. 2(c)) can often not be approximated accurately by a single Gaussian.A GMM aims at representing multimodal distributions based on an additively weighted combination of multiple Gaussian components.Each Gaussian component is typically referred to as a "kernel".Let N k be the total number of kernels and let h ∈ {1, . . ., N k } be the kernel index.A multimodal prior pdf that follows a GMM representation can then be written as (h) .The corresponding multimodal posterior pdf f (x|z) can be computed by performing Algorithm 2 and Algorithm 3 N k times in parallel, i.e., one instance of both algorithms is performed for each kernel N x; x * (h) , P (h) , h ∈ {1, . . ., N k }.To obtain a GMM representation composed of an arbitrary number N ′ k of kernels, a resampling step is then performed, i.e., N ′ k particles are drawn from the overall N k N p particles based on their weights w (h,i) 1 . The resampled particles represent the mean of N ′ k new kernels.The covariance of the new kernels is inherited from the original kernel the mean was sampled from.Pseudocode for GMM importance sampling with invertible PFL is provided in Algorithm 4. Importance sampling with invertible PFL makes use of resampling as presented in Algorithm 5. Algorithm 4: GMM Importance Samp.with Invertible Flow // see Alg. 5 9 Output: x * (h) For N k = 1 and N p > 1, this importance sampling approach is equivalent to invertible PFL based on the EDH update step.Furthermore, if N p = 1, N k > 1, this importance sampling approach is equivalent to invertible PFL based on the LEDH update step.Note that for N p = 1, as performed by the LEDH, an additional extended or unscented Kalman update step needs to be used to calculate an approximate covariance matrix [42]. ; and

B. Measurement-Origin Uncertainty (MOU)
In a variety of estimation problems, the measurement model suffers from a deficiency beyond measurement noise referred to as MOU [3].Here, there is a single object but multiple measurements and it is not known which measurement was generated by the object.Consider a single object with state x and measurements z (m) , m ∈ {1, . . ., M }.If m ′ is the measurement that was generated by the object, the corresponding measurement model is given by z (m ′ ) = h(x) + v. Based on this model, the conditional pdf of the object-originated measurement z (m ′ ) reads f o z (m ′ ) |x .All the other measurements are false positives (FPs) that follow the FP pdf f fp z (m ′ ) .It is assumed that at most one measurement originates from the object.The probability that the object generates a measurement is p d , and the mean number of FPs is Poisson distributed with mean µ fp .Since it is unknown which measurement was generated by the object, a discrete and random association variable a ∈ {0, 1, . . ., M } is introduced.Here, a = 0 describes the event where no measurement originated from the object and a = m, m ∈ {1, . . ., M } describes the event where measurement z (m) was originated from the object.Let z = [z (1)T , . . ., z (M)T ] T be the joint measurement vector.Following common assumptions [3], conditioned on x, the joint pdf of z and a is given by For z fixed, one can use this conditional pdf to directly compute the MOU likelihood function and the marginal probability mass function (pmf) The values of the pmf p(a|z) are also referred to as marginal association probabilities [3], i.e., they represent the probability of a particular association event a ∈ {0, . . ., M } conditioned on an observed z.

C. Importance Sampling with Invertible PFL for Problems with MOU
In principle, the MOU likelihood function in ( 14) can be directly used for importance sampling as in (1).However, in problems with MOU, importance sampling based on invertible PFL is complicated by the fact that there are multiple measurements, and it is thus not clear which measurement should be used to compute PFL parameters ( 7) and ( 8), i.e., the PFL proposal pdf q PFL (x 10) cannot be directly used.To address this problem, we propose the combined proposal pdf where we used the marginal association probabilities p(a|z) to weight the proposal pdfs q PFL x|z (m) related to PFL based on measurements z (m) , m ∈ {1, . . ., M } (cf.(10)) and the Gaussian prior pdf f (x) = N (x; x * 0 , P ).In particular, recall that the proposal distribution q PFL x|z (m) ) related to the PFL x 0 − → z (m) − → x (m) , can be evaluated as where θ (m) is the mapping factor.A total of (M +1)N p particles representing the proposal pdf in ( 16) is obtained by drawing N p particles for each of the M + 1 components in (16) and calculating corresponding marginal association probabilities and weights.For the first component related to association event a = 0, N p particles x (i,0) Np i=1 are directly drawn from f (x), the corresponding marginal association probability is obtained by using ( 13) in (15), i.e., p(a = 0|z) ∝ 1−p d , and the corresponding combined proposal weights are set to ω ′(i,0) = p(a = 0|z) N (x i=1 are obtained for each m ∈ {1, . . ., M }.An approximation of each marginal association probability p(a = m|z), m ∈ {1, . . ., M } is finally calculated from these particles by using ( 13) in (15) and performing Monte Carlo integration [24] based on the proposal pdf q PFL x (i,m) |z (i,m) in (17), i.e., The corresponding combined proposal weights are set according to (cf. ( 16) and (17)) Finally, we reindex the resulting particles and weights ω ′(i,a) , x (i,a) Np i=1 , a ∈ {0, . . ., M } to obtain ω (l) , x (l) L l=1 where l = i(a + 1) and L = N p (M + 1).Following the importance sampling principle, we next aim to compute particles w (l) , x (l) L l=1 that represent the posterior pdf f (x|z).In particular, by plugging ( 14) into (1) and by using ω (l) , x (l) (cf.( 19)) to represent q(x|z) in (1), we obtain ω (l)  . ( The resulting set of particles w (l) , x (l) L l=1 is an asymptotically optimal representation of f (x|z) for scenarios with MOU.

D. GMM Importance Sampling with Invertible Flow for Problems with MOU
For problems where the prior distribution is non-Gaussian and potentially multimodal, GMM PFL with invertible flow discussed in Section III-A can be directly applied to problems with MOU.Pseudocode for GMM importance sampling with the invertible flow for MOU problems is provided in Algorithm 6. Algorithm 6 relies on the measurement evaluation presented in Algorithm 7 and the measurement update presented in Algorithm 8.Note that in Algorithm 8, for future reference, we have also introduced extrinsic data association information denoted as κ (a) , a ∈ {0, . . ., M }.In the single object tracking considered here, we have κ (a) = 1.Note that for later use in Section V-D, Algorithm 7 performs measurements evaluation by also taking a probability of existence, p, into account.Similarly, Algorithm 8 also updates p.By using p = 1 as input for Algorithm 7 (see Algorithm 6, line 2), Algorithm 6 is equivalent to the estimation method discussed in Sections III-B and III-C.Furthermore, note that for consistency with Section V-D, we introduced the notation β (a) ∝ p(a|z), a ∈ {0, . . ., M } in Algorithm 6, Algorithm 7, and Algorithm 8.

A. PO States
As in [9], [10], we consider MOT for an unknown, timevarying number of objects by introducing PO states.The number of POs J k−1 at discrete time k−1 0 is the maximum possible number of objects that have generated a measurement up to time k − 1.At time k, a new PO is introduced for each of the M k observed measurements, and the total number of POs is updated as J k = J k−1 + M k .All POs that have been introduced at previous time steps are referred to as legacy POs, i.e., at time k, there are J k−1 legacy POs and M k new POs.
The augmented state of PO j ∈ {1, . . ., J k } is given by , where the state x (j) k of PO j consists of the position and possibly further parameters of the object represented by the PO.Furthermore, the existence variable r (j) k ∈ {0, 1} models the existence/nonexistence of PO j in the sense that PO j exists at time k if and only if r For nonexistent POs, i.e., r k is obviously irrelevant.Thus, all pdfs of augmented PO states f y is an arbitrary "dummy pdf" and f is a constant.To distinguish between legacy and new POs, we denote by y The concept of legacy and new POs can be extended to scenarios with S sensors as follows.Let M k,s be the number of measurements at time k and sensor s ∈ (1, . . ., S), where (1, . . ., S) is an arbitrary processing order of the sensors.The maximum possible number of objects that generated a measurement up to time k and sensor s is J k,s = J k,s−1 + M k,s , with J k,0 J k−1 .

B. Problem Formulation and Selected Messages of the SPA
At each time step k 1, we consider the tracking of an unknown number of objects based on measurements z 1:k .Object detection is performed by comparing the existence probability p r (j) k = 1 z 1:k with a threshold P th , i.e., PO j ∈ {1, . . ., J k } is declared to exist if p r For existent POs, state estimation is performed by calculating the minimum mean-square error (MMSE) estimate [58] as Both object detection and estimation require the marginal posterior pdfs f x k z 1:k by direct marginalization is infeasible due to the large number parameters in the joint posterior distribution in [57, Eq. ( 1)].
As in [9], [10], we consider approximate calculation by performing the loopy SPA on the factor graph in Fig. 1 of [57] and passing messages only forward in time.This makes it possible to efficiently calculate so-called beliefs f x , j ∈ {1, . . ., J k } which accurately approximate the marginal posterior pdfs f x k z 1:k , j ∈ {1, . . ., J k } needed for object detection and estimation.To keep computational complexity feasible, at the end of each time k with all sensors processed, a suboptimal pruning step has to be performed.Here, POs with probability of existence p (j) k p(r (j) k = 1|z 1:n ) below a threshold P pr are removed from the state space.
Next, we review the SPA messages that will later be calculated based on PFL.We will limit our discussion to messages and beliefs related to legacy PO states.Messages and beliefs related to new PO states are obtained by performing similar steps.(A complete description of message passing for MOT is provided in [10, Section IX-A].) We consider sequential sensor processing, where the SPA incorporates sensor measurements sequentially and in an arbitrary processing order at each time step.The part of the factor graph that represents the processing of the measurements of one sensor at one time step is shown in Fig. 1 of [57].At the processing step related to time k and sensor s, the "prior messages" of legacy PO states are denoted as α k,s .At time k and sensor s = 1, these message are computed by a prediction step [10, Section IX-A1], i.e., α that makes use of the state-transition function f x At the time k and sensor s > 1, the "prior message" is the same as the belief after the previous sensor update, i.e., α k,s−1 .For future reference, we also introduce p k,s as the predicted probability of existence for each legacy PO.Note that p After the prior messages have been computed, a "measurement evaluation" step for each legacy and each new PO is performed.Here, we denote by a k,s the association variable related to measurement m at the update step related to sensor s at time k (see [57,Sec. 1.2] for details).The SPA messages that are passed from the factor nodes q(x in [57, Eq. ( 3)] to the adjacent variables nodes a k,s , respectively, are computed.For legacy POs, these messages are given by (see [10, Section IX]) For new POs, the corresponding messages are denoted as ξ k,s , m ∈ {1, . . ., M k,s } and are calculated similarly (see [10,Section IX]).
Next, probabilistic data association (DA) is performed by means of iterative SPA message passing with input messages β [10, Section IX-A3] for details).After convergence, corresponding output messages κ k,s , m ∈ {1, . . ., M k,s } are available for legacy POs and new POs, respectively.Probabilistic DA is followed by a "measurement update" step.Here, for legacy POs, messages γ and as γ k,s 0 .Measurement update for new POs is performed by following similar steps [10, Section IX-A].
Finally, beliefs are calculated to approximate the posterior pdfs of POs.For legacy POs, beliefs f x and as f x k,s , m ∈ {1, . . ., M k,s } also shown in Fig. 1 of [57].The resulting beliefs for legacy and new POs are used as the prior messages for measurement update of sensor s + 1 as discussed above, i.e., α k,s , j ∈ {1, . . ., J k,s }.When the measurements of the last sensor in the sequence have been processed, i.e., s = S, the resulting beliefs are used in the prediction steps (21) of the next time step k + 1.

V. GRAPH-BASED MULTISENSOR MOT
WITH INVERTIBLE PFL In nonlinear MOT scenarios, calculation of β k,s in (24) related to legacy PO states as well as their counterparts ξ k,s , related to new PO states cannot be performed in closed form.We propose a particle-based implementation where a proposal pdf is established using invertible PFL as introduced in Section II-B.This makes it possible to implement multisensor MOT with high dimensional states and nonlinear measurement models.A single time step of the proposed particle-based implementation is discussed next.At first, we assume a single Gaussian kernel as the prior knowledge for each PO.An extension to GMM is also presented.In what follows, we consider a single time step, remove the time index k, and use the index − short for k − 1.

A. Prediction
It is assumed that the beliefs of legacy POs at time k − 1 are represented by a single Gaussian distribution, i.e., f x In MOT problems, the state transition function underlying the state transition model f x discussed in [10, Section VIII-C] is typically linear with additive Gaussian noise, i.e., j) where G is the state transition matrix and u (j) is an additive Gaussian noise vector with mean u * and covariance matrix P u .Consequently, the messages computed in the prediction step are also represented by a Gaussian distribution, i.e., α , P 1 ), j ∈ {1, . . ., J − } with mean, covariance matrix, and existence probability given by x * (j) 1 = Gx * (j) − +u * , P − G T + P u , and p (j) 1 = p su p (j) − , respectively.Here, p su is the survival probability, i.e., the probability that an object that exists at time step k − 1, still exists at time step k.Here If the state transition function is not linear with additive Gaussian noise, for each j ∈ {1, . . ., J − }, N p particles are drawn from N (x − ), the prediction step of a conventional particle filter is performed [25], and a predicted Gaussian representation N (x , P ) is computed from the resulting particles using Algorithm 3.

B. Measurement Evaluation
The following steps are performed sequentially for each sensor s = 1, . . ., S. First, a particle representation Np i=1 is obtained by drawing particles x s ) and setting the corresponding weights to ω (i,j) 0,s = p (j)  s /N p .Next, we compute an extended set x (i,a, Ms a=0 , that consists of N p particles and weights for each value of a (j) s ∈ {0, . . ., M s } and j ∈ {1, . . ., J s−1 }. (Note that measurement gating [3] can be employed to reduce the number of measurements used for PFL.)For a (j) s = 0, we perform no flow, i.e., we set . By making use of the invertible PFL principle (cf.( 10)) [42], the weights ω (i,m,j) s corresponding to the migrated particles x (i,m,j) s are obtained as with mapping factor θ (j) m,s (cf.(11)).Note that the sets of weighted particles x s , 1 .The result of this particle migration along the flow defined by a measurement z (m) s , m ∈ {1, . . ., M s }, is that the particles are now at locations where the evaluation of the corresponding likelihood function f z (m) s x (j) s will produce a significant particle weight.Particle degeneracy and thus approximate message-passing operations accurately even if the dimension of the state is high [9].
The measurement evaluation step can now be performed on weighted particles x (i,a, Ms a=0 by calculating an approximation β(j) s (a) of the messages β (j) s (a) in ( 22) for all j ∈ {1, . . ., J s−1 }, a ∈ {0, . . ., M s } as For the computation of a particle representation of new PO states, we first draw particles x .These messages are used as an input for the iterative SPA for data association [9], [10], [59] performed next (see [10, Sec.VI] for details).

C. Measurement Update and Belief Calculation
After the iterative loopy SPA for data association has been converged, the messages κ (see [10, Section IX] for details).Consequently, we obtain a new set of reindex particles and weights x Ls l=1 from x (i,a,j) , ω (i,a,j) Np i=1 , a ∈ {0, . . ., M s } by using l = i(a + 1), ω (l,j) s = ω (i,a,j) /p(a (j) s = a|z s ), and L s = N p (M s + 1).Next, based on (24), we update the particle weights of the legacy POs j ∈ {1, . . ., J s−1 } by first computing and then calculating normalized weights as (25) Note that denominator of ( 25) is a particle-based approximation of C (j) s in (24).The resulting particles and weights x x (i,h,a,j) , w (i,h,a,j) Np i=1 , θ (h,a,j) N k h=1 , β (a,j) Ms a=0 = Evaluation x * (h,j) , P (h,j) N k h=1 , p (j) , zs , m ∈ {1, . . ., Ms} w (i,h,a,j) Np i=1 , θ (h,a,j) N k h=1 , β (a,j) , κ (a,j) M a=0 , p (j)    For new POs j ∈ {J s−1 +1, . . ., J s }, approximate existence probabilities p(j) s and a Gaussian representation of the beliefs f x , for m = j −J s−1 are calculated by performing similar steps as discussed above for legacy POs.Existence probabilities and Gaussian representations of the beliefs related to PO that have not been pruned are then used as input for processing measurements of the next sensor s + 1 or, in case s = S, for processing at the next time step.
For POs that have been declared to exist after the last sensor update, i.e., p(j) S > P th , an approximate MMSE state estimate is directly given by the mean of the Gaussian representation, i.e., x(j) ≈ x * (j) S .

D. The Proposed Multisensor MOT Method
In multisensor MOT problems with nonlinear measurement models, object beliefs are non-Gaussian and potentially multimodal.Here, graph-based multisensor MOT with invertible PFL has to be combined with a GMM discussed in Section III-A.Pseudocode for a single time step of the resulting multisensor MOT method is provided in Algorithm 11.Algorithm 11 relies on the single sensor update step provided x * (h,j) , P (h,j) N k h=1 , p (j) J j=1 = xS * (h,j) , P where J = JS 13 Output: x * (h,j) , P (h,j) N k h=1 , p (j) J j=1 Note that for the execution of the proposed method, for each time step and each object, we need to update each particle per kernel, pseudo-time, measurement, and sensor.The asymptotic complexity with respect to these parameters, per object and time step, thus reads O N λ N k N p S s=1 M s .The complexity of a conventional bootstrap implementation per object and time step is O N b S s=1 M s , where N b is the number of particles.The improved runtime-complexity tradeoff of the proposed method is due to the fact that, in challenging problems as the one considered in Sections VI, for N λ N k N p = N b , the proposed method strongly outperforms graph-based MOT that relies on conventional particle filtering.Since memory requirements per time step and object are N k N p for the proposed method, and N b for bootstrap particle filtering, and since N λ ≫ 1, for N λ N k N p = N b , the memory requirements of the proposed method are significantly lower compared to conventional particle filtering.

VI. NUMERICAL RESULTS
Next, we report simulation results assessing the performance of our method and comparing it with that of two reference methods for multisensor-multiobject tracking .

A. Tracking Scenario and Reference Methods
We consider an underwater 3-D surveillance scenario where eight objects are tracked by two static sonar hydrophone arrays.200 time steps are considered.The hydrophones are deployed about 1300 m below sea level.The object states at time k consist of 3-D position and velocity, i.e., x 3,0 = −1000 m.The velocity is obtained by setting the vertical speed to zero and the horizontal velocity vector pointing to the circle center with ẋ(1)2 1,0 + ẋ(1)2 2,0 = 0.3m/s.For the appearance of each further object j = 2, . . ., 8, an initial position is randomly generated near the initial position of the previously appeared object, i.e., around a circle with radius 50 m centered at [x 3,k ].The initial velocity vector is set with respect to the center circle, as discussed above.As a result of this initialization procedure, tracks start in close proximity in time and space.This makes it challenging to perform data association and declare the existence of newborn objects.The times when objects appear and disappear, as well as the states of appearing objects, are unknown to all simulated tracking methods.All methods aim to detect the presence of a new object and sequentially estimate its state across time merely based on TDOA measurements and the statistical model.The prior intensity of object birth is modeled at each timestep by a Poisson point process with mean µ b = 0.05.The prior pdf for newborn object, f b x k , is uniform on the ROI for the 3-D position and Gaussian with zero mean and covariance matrix 5m 2 /s 2 I 3 for the 3-D velocity.The survival probability is p su = 0.95.The object declaration threshold is set to P th = 0.5 and the pruning threshold to P pr = 10 −4 .
The two hydrophone arrays have the geometry of the array described in [61] and are located at [519m 137m −1300m] T and [−519m − 137m − 1300m] T , respectively.Each hydrophone array consists of 4 receivers.Hence, there are six receiver pairs at each array.Each receiver pair generates TDOA measurements and is considered a sensor for MOT.This means that multisensor measurements z where p sL and p sR are the paired receiver positions of sensor s, c = 1500m/s is the propagation speed, and v (m) k,s is additive zero-mean Gaussian noise with standard deviation σ v that is assumed statistically independent across s, k, and m.The pdf of FP measurements, f fp z We compare the proposed SPA-based MOT with the embedded particle flow sampling strategy ("SPA-PF") with two reference sampling strategies.The first ("SPA-PM") follows the sampling strategy of the bootstrap particle filter [23], [25] and uses predicted beliefs as proposal pdf [9].The second (yet unpublished) ("SPA-UT") follows the sampling strategy of the unscented particle filter [35], i.e., it uses a Gaussian mixture representation and the unscented transformation to calculate an informative proposal pdfs.For SPA-PM, we use N b = 10 6 particles for newborn POs and N b = 6 • 10 4 particles for legacy POs.For all other simulated methods, we use N k = 100 kernels.For each kernel representing a newborn PO, we set N p = 500 for SPA-UT and SPA-PF and for each kernel representing a legacy PO, we se N p = 30 for SPA-UT and SPA-PF.Since N b /N p N k = 20 the memory requirements of SPA-PM are 20 time higher compared to SPA-PF and SPA-UT.
We also simulate two variants of the second reference method to obtain a similar runtime for SPA-PF.In particular, for "SPA-UT-1", we set N p = 4000 for kernels representing newborn POs and N p = 250 for kernels representing legacy POs.Furthermore, for "SPA-UT-2", we set N p = 6000 for kernels representing newborn PO and N p = 30 for kernels representing legacy PO.Note that the memory requirements of SPA-UT-1 and SPA-UT-2 are higher than the ones of SPA-UT and SPA-PF due to their higher values of N p .Finally, we also simulate a method ("SPA-PF-H") that uses the sampling strategy of SPA-PF for kernels representing newborn POs and the sampling strategy of SPA-PM for kernels representing legacy POs.Using fewer samples for the kernels representing legacy PO, or even using the strategy of SPA-PF, is motivated by the fact that the beliefs of legacy POs are typically unimodal and quite informative.We set N λ = 20 for SPA-PF.For all considered methods, 100 simulation runs are performed.All methods are implemented in MATLAB, and each simulation run is processed on a single core of a 2.6GHz Intel Xeon Gold 6240 processor.
The performance of the six MOT methods is evaluated w.r.t. to changes in four system parameters (i) object driving noise standard deviation σ w , (ii) mean number of FPs µ fp , (iii) detection probability p d and (iv) measurement noise standard deviation σ v .Note that for the setting σ v = 5 × 10 −7 s, the number of samples N p was doubled for all methods to yield high tracking performance.The tracking accuracy of the various methods is measured by the Euclidean distance-based optimal sub-pattern assignment (OSPA) metric with cutoff parameter C = 50 [62].

B. Performance Comparison
In what follows, we discuss two scenarios where particle degeneracy is particularly pronounced.Particle degeneracy can be caused by uninformative prior information or a very informative likelihood function.To obtain a scenario with uninformative prior information, we consider the case σ w = 1 m/s 2 .In addition, to obtain a scenario with an informative likelihood function, we consider the case σ v = 5×10 −7 s.Fig. 3 and Fig. 4 show the mean OSPA (MOSPA) error-averaged over 100 simulation runs-of all methods versus time k for these two scenarios.It can be seen that all methods yield error peaks at time steps where the objects appear.This is because the MOT methods do not know when a new object appears and often need few time steps after the appearance of an object to declare its existence based on TDOA measurements and the statistical model.However, the proposed methods, i.e., SPA-PF and SPA-PF-H, have lower error peaks at the time steps when objects appear, i.e., SPA-PF and SPA-PF-H can often declare the existence of an object faster.In addition, SPA-PF and SPA-PF-H can outperform the other reference methods at almost all time steps.SPA-PM, in particular performs very poorly due to particle degeneracy.It can be noted that SPA-UT-1 and SPA-UT-2 have improved performance compared to SPA-UT but are still outperformed by SPA-PF and SPA-PF-H despite their more extensive memory requirements.Furthermore, the MOSPA error and the runtime per time step of the six methods for different system parameter values are shown in Table I.The default value of the four parameters are σ w = 0.1 m/s 2 , σ v = 1×10 −6 s, µ fp = 5 and p d = 0.9.For each row corresponding to a particular system parameter value, the value of the other three parameters is set to the default value.For each system parameter value corresponding to one row, the MOSPA and runtime are averaged over 100 simulation runs and 200 time steps.The best and the second best MOSPA value corresponding to each system parameter value is marked by an underline and dashed underline, respectively.As can be noted, the performance of SPA-PM significantly degrades when measurement noise standard deviation is reduced to σ v = 5 × 10 −7 s.This unwanted and counterintuitive behavior is a clear indicator of particle degeneracy in SPA-PM.Only with the proposed method is it possible to yield improved tracking performance as measurement noise variance is reduced.It can be seen that for almost all system parameter values, the proposed SPA-PF and SPA-PF-H outperform the reference methods.At the same time, their runtime is comparable with SPA-UT-1 and SPA-UT-2, which yield higher memory requirements.The only scenario where SPA-PF and SPA-PF-H do not result in the lowest MOSPA is when σ v = 2×10 −6 s.In this case, since measurements are not very informative, SPA-PM does not suffer from particle degeneracy.It can also be seen that SPA-PF-H outperforms SPA-PF both in terms of MOSPA and runtime, while both SPA-PF and SPA-PF-H typically outperform SPA-UT.We can thus conclude that the main challenge for the sampling method is the initial time step after a new object appears in the scene.Here, beliefs of new POs are highly uninformative and have complicated shapes.At later time steps, beliefs become informative and unimodal and can thus be computed accurately with fewer samples and the sampling strategy of the bootstrap particle filter.For real-time processing, it is expected that an adaptation of the proposed method for execution on graphical processing units (GPUs), can strongly reduce runtimes by exploiting the highly parallelizable nature of PFL.
In what follows, we numerically investigate the effect of the number of kernels N k on system performance.In particular, we set the number of kernels as N k ∈ {1, 5, 10, 50, 100} in SPA-PF-H and compare SPA-PF-H for these different values of N k with SPA-PM.The parameters σ v , µ fp , σ w , and p d are set as dicussed above.For comparable runtimes of all SPA-PF-H variants, we set N p ∈ {100000, 20000, 8000, 1300, 500} for newborn object and N p ∈ {15000, 1000, 400, 70, 30} for legacy objects.As expected, it can be seen in Fig. 5 that the accuracy of SPA-PF-H improves with increasing N k .While SPA-PF-H performs very poorly for N k = 1, notably, it can already outperform SPA-PM for N k = 5.The average runtime per time step of SPA-PF-H is 27.14s, 15.06s, 10.85s, 11.35s, and 12.72s for N k = 1,N k = 5, N k = 10, N k = 50, and N k = 100, respectively, as well as 22.66s for SPA-PM.

C. Real Data Experiment: Echolocation in Oceanography
We further validate our proposed MOT method in an underwater acoustic tracking scenario.The acoustic signals "clicks" emitted by two Cuvier's beaked whales are recorded by two high-frequency acoustic recording package (HARP) [63], each of which is equipped with four hydrophones.HARPs are deployed at a depth of 1330 m and approximately 1 km apart.Preprocessing of acoustic data is described in [50].Each pair of hydrophones on each HARP acts as a sensor that provides TDOA measurements every 7s.Since there are 6 pairs of hydrophones on each HARP, there are a total of S = 12 sensors providing TDOA measurements.We use a dataset that consists of 172 time steps and has a total duration of roughly 20 minutes.It was recorded on July 1st, 2018, in Southern California.Tracking results are shown in Fig. 6.The red solid lines show the estimated tracks of two whales provided by SPA-PF-H.The two whales were initially detected at depths of about 450 m and then kept diving until a depth of 1300 m.For obvious reasons, no ground truth information exists for this scenario.However, we have added reference tracks of the two whales that are the result of a trained operator handannotating preprocessed acoustic data.
In Figure 7, we show the estimate tracks of SPA-PF-H compared with SPA-PM in 2-D.SPA-PM does break and merge tracks, which makes it very difficult to determine how many whales are actually there.Most importantly SPA-PF-H can potentially replace the human operator, while SPA-PM cannot.The overall number of particles of SPA-PM is 20 times that of SPA-PF-H.The runtime per time step of SPA-PM is 4.65 s while 2.35 s of SPA-PF-H, i.e., SPA-PF-H is faster than   SPA-PM.Since the measurement interval is 7 s, both SPA-PF-H and SPA-PM can be used in real-time.For larger scenarios with more than two whales, a GPU implementation is required for real-time processing.

VII. CONCLUSION
We presented a graph-based Bayesian method for multisensor MOT with high-dimensional object states.Particle degeneracy is avoided by performing operations on the graph using PFL.Our numerical results indicate that the main challenge for sampling is representing the posterior distribution at the initial time step after a new object appears in the scene.Compared to state-of-the-art reference methods, we show favorable tracking performance in a 3-D MOT scenario.The introduced approach is expected to be particularly appealing for passive surveillance problems [64].Future research avenues include graph-based processing with embedded stochastic PFL [65], [66] and applications including extended object tracking [13], [67], simultaneous localization and object tracking [16], and information-seeking control [68].We also aim to demonstrate real-time processing capabilities of the proposed approach by execution on graphical processing units (GPUs), exploiting the highly parallelizable nature of PFL.
This material is based upon work supported by the Under Secretary of Defense for Research and Engineering under Air Force Contract No. FA8702-15-D-0001 and by the National Science Foundation (NSF) under CAREER Award No. 2146261.Parts of this work have been presented at the IEEE RADAR-21, Atlanta, Georgia, May 2021.

Fig. 1 .
Fig. 1.Particle degeneracy in a tracking scenario with 3-D object state, x = [x1 x2 x3] T , and a single 1-D TDOA measurement z1.A single time step in considered.The 1-D TDOA measurement is generated by the sensor shown as gray circle.Assuming no measurement noise, the 1-D TDOA measurement describes potential 3-D object locations on the hyperboloid shown in red.The object, shown in black, is located on the hyperboloid.Note that any other location on the hyperboloid will lead to the same measurement in the case without noise.(a): The prior pdf, f (x), is Gaussian, with the mean depicted as a big blue dot.2000 particles, shown as small light blue dots, are drawn from the prior distribution.On the right, the prior and posterior pdfs for the case with measurement noise are shown in three separate 2-D plots.Each of these plots is obtained by depicting the prior and posterior pdfs along the three axes of the coordinate system.(b): After importance sampling, as performed by the conventional "bootstrap" particle filter, only a single particle has a nonzero weight.This single particle does not accurately represent the posterior pdf p(x|z1) for future processing, e.g., of a measurement, z2, provided by a second sensor.

Fig. 2 .
Fig. 2. Example of PFL in the tracking scenario with 3-D object state and a single 1-D TDOA measurement discussed in Fig. 1. (a): 2000 particles represent the prior pdf at the onset of the flow, i.e., at pseudo time λ = 0, are depicted.(b): An intermediate flow state corresponding to λ = 2 × 10 −9 is shown.The tracks of 8 selected particles are indicated as red dashed line with arrows.(c): At pseudo time λ = 1, particle migration is completed and the resulting particles represent the hyperboloid-shaped posterior pdf.(d):The histogram of the flowed particles together with 1-D prior and posterior pdfs is drawn.The representation of the posterior pdf provided by the particles after the flow is much more accurate than the single "degenerated" particle resulting from conventional particle filtering discussed in Fig.1.Due to approximations performed in PFL, there can be a small mismatch of particles after the flow and the true posterior pdf.Such a mismatch can also be seen in (d) by comparing the posterior pdf with the histogram of particles at λ = 1.Invertible PFL can eliminate such mismatch and provide an asymptotical optimal representation of the posterior pdf by making it possible to compute particle weights for importance sampling.
)).For each other component related to association event a = m, m ∈ {1, . . ., M }, first N p particles x (i,m) 0 Np i=1 are drawn from f (x).Next, the PFL x 0 − →z (m) − → x (m) is applied to the particles x (i,m) 0 Np i=1 and new particles x (i,m) Np

M a=0 Algorithm 8 :
GMM Importance Sampling with Invertible Flow and DA -Measurement Update 1

k
and by y (m) k the augmented state of a legacy PO and a new PO states, respectively.
the association variable related to PO j at the update step related to sensor s at time k and by b (m) , m ∈ {1, . . ., M k,s } for new POs is performed by following similar steps [10, Section IX-A].Note that this calculation of new POs involves the messages ι , i ∈ {1, . . ., N p } from f b x as introduced in [Sec.1.1] [57].Next, for each new PO j = J s−1 + m, m ∈ {1, . . ., M s }, new particles and corresponding weights x (m,i) Np i=1 by performing the invertible PFL x 0,s − → z (m) s − → x (m) s .Note that this flow relies on the mean x b and covariance matrix P b of f b x .Finally, for each m ∈ {1, . . ., M s } approximate messages ξ(m) s b (m) s are calculated from x (m,i) s , w (m,i) s Np i=1 by performing the same steps as described above for the calculation of β(j)

s , 1 ,
{1, . . ., J s−1 } and ι(m) s b (m) s , m ∈ {1, . . ., M s } are available.These messages are used to obtain an approximation γ(j) j ∈ {1, . . ., J s−1 } in (23) as well as an approximation ς(m) s x (m) s , 1 of the messages ς (m) s x (m) s , 1 , m ∈ {1, . . ., M s } in [10, Section IX].Beliefs approximating the posterior pdf of POs are now computed by means of importance sampling.As in Section III-C, we use the marginal association probabilities p(a (j) s |z s ) to weight the proposal pdfs q PFL x (j) s |z (m) s related to PFL based on measurements z (m) s , m ∈ {1, . . ., M }.Note that in MOT scenarios, accurate approximations of p(a (j) s |z s ) can be obtained as p(a (j)

17
j ∈ {1, . . ., J s−1 }.These particles can be used to calculate an approximation of the existence probability as p(j) ] T , j = 1, . . ., 8 and evolve according to a constant-velocity model[60, Sec.6.3.2],where the dynamic noise has the physical interpretation as an acceleration with variance σ 2 w .The region of interest (ROI) is[−1000m, 1000m]×[−1000m, 1000m]×[−1500m, −500m].Objects appear at k ∈ {1, 10, 20, 30, 40, 50, 60, 70} and disappear at k ∈ {130, 140, 150, 160, 170, 180, 190, 200}.To simulate a tracking scenario with challenging data association, we generate the initial state of the objects as follows.The initial state of the first object is randomly generated by setting its position at a circle centered at the origin with a radius of 50 m and a depth of 1000 m, i.e., for its position, we have x , m ∈ 1, . . ., M k,s and s ∈ 1, . . ., 12 are obtained by the two arrays at time k.At each sensor s, a random number of M k,s measurements are generated.In particular, the TDOA measurement z (m) k,s of a detected object with state x (j) k is modeled as z

Fig. 5 .
Fig. 5. OSPA performance of SPA-PM and SPA-PF-H for different number number of kernels, N k .

Fig. 7 .
Fig. 7.A comparison between the estimated tracks provided by SPA-PF-H and SPA-PM.The hand-annotated tracks of the two whales are marked for reference (gray dashed line).Each axis of the 3-D domain is shown individually.

TABLE I MOSPA
AND RUNTIME PER TIME STEP OF DIFFERENT ALGORITHMS FOR THE CONSIDERED TRACKING SCENARIO W.R.T RELATED PARAMETERS.