Multi-Target Localization and Tracking Using TDOA and AOA Measurements Based on Gibbs-GLMB Filtering

This paper deals with mobile multi-target detection and tracking. In the traditional method, there are uncertainties such as misdetection and false alarm in the measurement data, and it will be inevitable having to deal with the data association. To solve the target trajectory and state estimation problem under a cluttered environment, this paper proposes a non-concurrent multi-target acoustic localization tracking method based on the Gibbs-generalized labelled multi-Bernoulli (Gibbs-GLMB) filter and considers an acoustic array of a fixed arrangement for the tracking of targets by joint time difference of arrival (TDOA) and angle of arrival (AOA) measurements. Firstly, the TDOAs are calculated by using the generalized cross-correlation algorithm (GCC) and the AOAs are derived from the received signal directions. Secondly, we assume the independence of the targets and fuse the measurements which are used to track the multiple targets via the Gibbs-GLMB filter. Finally, the effectiveness of the method is verified by Monte Carlo simulation experiments.


Introduction
Passive detection, such as multi-sensor array localization of acoustic sources, plays an important role in the field of target tracking [1]. Localization through acoustic signals has broad applications in both civil and military fields, for example, the detection of unknown objects in the airport, detection of illegal traffic whistles/horns, localization of submarines or marine animals and the localization of explosion sources. Existing passive localization and tracking techniques include Time Difference of Arrival (TDOA) [2], and Angle of arrival (AOA) [3]. TDOA requires multiple observation devices to be used at the same time to ensure that the clock time of the sensor is consistent among the multiple groups of sensors [4].
In multi-target scenarios, it is difficult to determine the number and states of targets due to clutter effects, miss detection and data association uncertainty. Reference [5] proposed the route-based dynamic modeling to improve data association performance. For the complex and non-linear acoustic signals, traditional signal processing techniques such as the Fang algorithm, the music algorithm and the Taylor series expansion method are used [2]. To address the nonlinearity of the measurement process, the extended Kalman filter, one of the common methods in target tracking, can be used to provide better information for multi-sensor information fusion. In Reference [6], they performed data fusion that combines the active detection and the passive interception using Multi-acoustic array localization is typical of multi-sensor, passive localization and nonlinear problems. In the case of traditional algorithms dealing with multi-sensor measurement uncertainty and target uncertainty, the target tracking is based on the method of associating data and there is no effective way to estimate the number of targets. The Gibbs-GLMB filtering can effectively solve these problems. We change from the original active tracking to passive tracking based on the original Gibbs-GLMB filter [36]. TDOA and AOA measurements are generated by calculating the true signal correlation, which is a nonlinear estimation problem. The TDOAs and the AOAs are computed by the generalized cross-correlation (GCC) [4] and receiving direction of acoustic signal, respectively.
For the reasons above, we use the Gibbs-GLMB filter. First, we assume that the target obeys a cv motion model and the observation information includes TDOA and AOA. Secondly, target tracking in a multi-sensor array is done using the Gibbs sampling implementation of the GLMB filter (Gibbs-GLMB filter), to reduce the computational complexity of the algorithm without sacrificing accuracy. The effectiveness of the algorithm is verified by three pairs of acoustic array sensors are deployed to track three targets.

•
Single-target state is expressed by a small letter, (e.g., x).

•
Multi-target states are represented by an italic capital letter, (e.g., X).

•
The labeled states and distribution are bolded, (e.g., x, X, π). • The spaces are represented by blackboard bold (e.g., the state space X and measurement space Z).
The inner product symbol is abbreviated as: The following multi-target exponential notation h X ∏ x∈X h (x), where h is a real-valued function, The generalization of the Kronecker delta for sets, vectors and integers: where inclusion function is denoted as: • X m:n is shorthand for the list of variables X m , X m+1 , · · · , X n .

Random Finite Set
In the multi-target environment based on RFS framework in time k, the states of multiple targets can be denoted by a set as X k = x k,1 , · · · , x k,N(k) ∈ F (X ) [13], where F (X ) is the all of finite subsets of state space and N (k) is the number of surviving targets. In the similar way, the observation about TDOAs and AOAs of the qth sensor pair can be described as Z is the space of finite subsets of observation space Z and M q (k) is the number of observed measurement.
In the location tracking area, there are many uncertainties in the number and measurement of the detection process, such as birth, death, derivation, false alarm and missed detection. Consequently, The multiple targets state in time k can be defined as [59,60]: where S k|k−1 (x), B k|k−1 (x) and Γ k are the RFS of survival target at time k − 1, the RFS of the target spawn at time k from the survival target at time k − 1 and the RFS of target new-born, respectively. There are clutter or false alarms in the tracking area, which can be expressed as: where Θ k (x) is the measurements with the RFS which produced by targets in the tracking area: here, theH miss and H miss are the hypotheses of detection and miss detection, more generally, whether the sensor has received the signal generated by targets. Moreover, K [q] k is the measurement set of alarms or clutter false which follows a poisson distribution with a uniform density U (z) on the observation area and is given by: In RFS framework, the probability density function that the state of multi-target makes a transition from state X k−1 to X k can be described as: where π T,k|k−1 (· |· ) is the probability density of spontaneous target birth and π Γ,k (·) is the probability density of target new-born.

Multi-Bernoulli RFS
The state of the target and the measurements are random variables. The Bernoulli distribution can be used to describe a single target X ∈ X. Hence, a singleton target probability X is r which satisfies the spatial distribution of a probability density p (x) and the probability that the target does not exist is 1 − r. The probability density distribution of the Bernoulli RFS is written as follows: The Multi-Bernoulli RFS is given by combining of the M independent Bernoulli RFS X (i) ∈ X, i = 1, · · · , M(satisfying X = M i=1 X (i) ) with existence probability r (i) ∈ (0, 1), which is described by . ∑ M i = 1 r (i) is the mean cardinality of the Multi-Bernoulli RFS. Therefore, the probability density distribution of multi-Bernoulli is expressed by [12,18]:

Labeled Multi-Bernoulli RFS
A labeled-random finite set (L-RFS) [34,35] means that each state of the RFS has a unique tag. This means we attach a unique label l ∈ L = {α i : i ∈ N} to each state x ∈ X where L is discrete countable space and N is the positive integer set space. The single target state is expressed as: The labels of the set X ⊂ X × L can be represented by L (X) = {L (x) : x ∈ X}, where L : X × L → L is defined by L ((x, l)) = l. The distinct label indicator is defined by ∆ (X) = δ |X| (|L (X)|).
The parameter of a labeled multi-Bernoulli(LMB) RFS can be described as a set {(r (ζ) , p (ζ) ) : ζ ∈ Ψ} with index set Ψ. We extend the problem on spaceX to spaceX × L, thus, the probability density distribution of labeled LMB-RFS is given by [34]: The following simplified alternative form of the LMB can be simplified as:

GLMB RFS
A generalized label multi-Bernoulli RFS under the state space X and the label space L has the following distribution [34]: where ξ= (θ 1:k ) ∈ Θ is a historical association maps. The non-negative ω (ξ) (L) and a probability density p (ξ) satisfy:
Assuming a single target state of position is x i = (x i , y i ), each set of sensors consists of two acoustic sensors.

Time Difference of Arrival
The time difference τ is expressed as: where both x i and s q,j are in Cartesian coordinates, · is the Euclidean-norm and v is the velocity of sound.

Angle of Arrival
For sensor 1, AOA can be expressed as: For ease of understanding, the TDOA and AOA are illustrated in the positioning system of

TDOA Measurement
The signals observed by a pair of sensors can be mathematically described as: where z 1 (t) and z 2 (t) are the signals received by the pair of sensor array, s (t) is the true signal, n 1 (t) and n 2 (t) are noise signals, τ 1:2 is the time difference between two sensors detecting the signal, α 1 and α 2 are signal amplitudes [8].
The time difference can be estimated by the generalized cross correlation (GCC) method [4]: Here, R GCC τ q is the GCC, where Z 1 (ω) is the Fourier transforms of z 1 (t) and Z * 2 (ω) is the Fourier transform conjugate of z 2 (t). ψ 1,2 (ω) is the weight function of GCC. To reduce environmental noise and reverberation interference, we choose the phase transform (PHAT) as our weight function, the formula is given by

AOA Measurement
The AOA is calculated by the positional relationship between the sensor and the target position. The difference in the AOAs of the sensor array is calculated by combining the measurements of the two sensors:

Motion Model
We take the CV model as an example of a linear model, also known as a non-maneuver model: where x is the location of the target,ẋ is the velocity of the target,ẍ is the acceleration of the target, w (t) is zero mean white noise. Let T denotes the sampling interval, then the discrete-time model is given by:

Target State Estimation
The purpose of multi-target tracking is to jointly estimate the target cardinality and target states based on the observations. Multi-target tracking can be transformed into the recursive Bayesian estimation problem by modeling the state and measurement of multi-target using RFS. The RFS approach can effectively deal with the uncertainty of data association between the target and the measurement and the state probability density function of the set of targets. We use π k ( ·| Z 1:k ) to indicate the RFS posterior probability density of multi-target state; f k|k−1 (· |· ) to represent the multi-target transition density; g k (· |· ) to represent the likelihood function. The posterior probability density of multiple targets is recursively calculated by the following prediction and update steps [11][12][13][14]: where the set integral on F (X × L) → R is defined as: All information about the multiple targets states are included in the multi-target posterior, for example, the number and state of the target at the current time.
We experience with q pairs of sensors, therefore the above update step and prediction step can be rewritten as [59][60][61]: π k X k Z The recursive process is not easy to deal with exactly due to the non-linear of the observation equation. The sequential Monte Carlo(SMC) methods are a viable approach to approximate the integrals of interest using random samples.

Particle Filter Implementation
Since the multi-target posterior probability density recursion requires the calculation of multi-set integral (24) and (25), its computational complexity is much larger than that of the single-target filtering process [16,61,62]. By the SMC method, the weighted particles can be estimated by recursive approximation to estimate the posterior probability density.
At the current time k, the particles are sampled by SMC, obtained from the spatial distribution of the target.X represents the set of importance weighted particles at time k − 1 and the multi-target posterior probability density can be expressed as: Algorithm 1 Particle Filter 1: for k = 1, · · · , T do 2: Normalise weights ω as the estimated posterior probability density;

The Multi-Sensor Likelihood
Given the multi-target state X, each (x, l) ∈ X is either detected with probability p D,m (x, l) and generates observation z with likelihood function g (z |x, l ). For S sensors, the multi-sensor and multi-target mapping [56] is defined by θ (m) : L → 0, 1 · · · , Z (m) , m = 1, ..., S. The set Θ represents the space of vector maps θ = (θ (1) , ..., θ (S) ). Assuming that the target and clutter generation are independent and the multi-sensor likelihood function is given by [34]: where p D,m (x, l) is the probability detection, K m is Poisson clutter, for sensor m.

GLMB Filter
The GLMB filter is a Bayesian recursion from the multi-Bernoulli distribution, which satisfies the following formula [34]: The forward propagation expression of GLMB Filter is as follows: The distribution of multi-target prior probability is given by the Equation (29), thus, the multi-target prediction is still the multi-Bernoulli distribution and the prediction step can be expressed as: where Here, the ω B (I + ∩ B) and ω s (ξ) (I + ∩ L) are weights of the birth labels (I + ∩ B) and surviving labels (I + ∩ L), respectively. p B (x, l) is density of a new-born target, p (ξ) s (x, l) is the probability density of the surviving target obtained from the prior probability p (ξ) (·, l). f (x|·, l) means density weighted by the probability of survival p s (·, l).
Given the predicted density as Equation (29), the update step can be expressed in the form of a truncated estimate:

Gibbs-GLMB Filter
Gibbs sampling is a special case of continuous Markov Chain Monte Carlo (MCMC), which can transform sampling from high-dimensional space to low-dimensional one [63]. Assuming that the target state is X k = x k,1 , · · · , x k,N(k) , which obeys the probability distribution π, the probability distribution of the target state is π x k,1 , · · · , x k,N(k) .

6: end for
In Algorithm 2, x k,1:n−1 is the samples x k,1 , · · · , x k,n1 that have generated at current time, x k−1,n+1:N(k) is associations x k,n+1 , · · · , x k,N(k) at previous time. The Gibbs sampling algorithm reduces the joint density estimation problem to conditional probability to reduce the sampling difficulty and finally updates all parameters by the iterative process of each parameter.
In the calculation process of the update step and the prediction step of the GLMB filter, the number of weights and data quantities of the update and the prediction step are exponentially increasing. By using optimal assignment implementation and the kth shortest path, the complexity of the measurement quantity is cubic [35]. The GLMB filter is truncated by Gibbs sampling, thereby joint prediction and update reduce the complexity of the measurements to linear.

Simulation Environment Settings
We use six sensors consisted of three arrays to observe three acoustic targets as shown in Figure 2. Three pairs of sensors track the target, each sensor's observation distance is 150 m, the simulated sound velocity is 344 m/s, the surviving probability is P S = 0.99 and the clutter intensity of Poisson distribution is λ c = 2. The scenario last 100 s, maximum number of targets is 3. The motion model is a linear state space equation(CV motion model) and the state of target is expressed as: where A is the target state transition matrix, B is the noise matrix and ω k is process noise and follows a standard Gaussian distribution. The sampling period is T = 1. In Figure 2, the sensors are displayed by the blue circle, the black circle represents the starting point and the triangle represents the end position. The target location is unknown and two scenarios were compared in this section. The position and velocity vector of target and sensor are represented as x k = p k,x ,ṗ k,x , p k,y ,ṗ k,y T and s i k = q i k,x ,q i k,x , q i k,y ,q i k,y T , respectively. The range dependent detection probability is defined as: where P D,max = 0.98, Σ D = diag(500, 500) 2 and C = 1 0 0 0 0 0 1 0 .
In the scenario 1, We build a parallel model and the survival period of three targets are 1 s-100 s, 10 s-90 s and 20 s-80 s, respectively. The initial states of the three targets are an LMB-RFS with parameters r B,k (l i ) , p B,k (l i ) In the scenario 2, the survival period of three targets are 1 s-90 s, 1 s-80 s and 30 s-100 s, respectively. Target 1 and target 2 are born in the same position at the same time. The initial parameters The experiment uses the three Matlab audio files sample1.wav, sample2.wav and sample3.wav as the acoustic signals of the three targets in the Figure 3.
where, z k is nonlinear. At time k,τ q is the time difference observed by a pair of sensors, δ q is the angle difference between a pair of sensors receiving signals, σ τ = 0.001 s and σ δ = (π/720) rad are the standard deviations of the Gaussian distributed measurement noise. In the scenario 1, three pairs of sensor arrays detected the measurements data of target 1 as show in the Figure 5.

Scenario 1
The simulation time is 100 s. The black line in the Figure 6 is the real trajectory of the target and the Red circle and the Color points are the estimated target location. It can be seen from the two pictures that the target tracks obtained from the Gibbs-GLMB filtering and the GLMB filtering are basically consistent with the true trajectory of the targets. From the simulation results in Figures 6-8, it can be seen that the tracking performance of both algorithms is better. In Figures 7 and 8, the cross points are all measurements in the simulation of 1 s-100 s and the points generated outside the target track are false alarms caused by clutter interference. When a target is born, the random clutter may cause false alarms at the position. Nevtheless, in the subsequent tracking and localization, most of these false alarms will be eliminated. Through 100 times Monte Carlo(MC) simulations, the number of targets is estimated as shown in the Figure 9. The red line is algorithm Gibbs-GLMB and the black dotted line is algorithm GLMB. We can see from the comparison of the two algorithms that the number of Gibbs-GLMB estimates is more accurate overall but the estimated number of targets has a large deviation when the actual number of targets changes. The cardinality estimates of targets based on GLMB Filter is always a slightly higher than the true comparison.   We use the Optimal Subpattern Assignment (OSPA) distance [64] (c = 100, p = 1) to analyze tracking performance. Figure 10 shows the simulation result over 100 MC runs. We can see that OSPA-Loc of two algorithms are very small in the whole process, indicating good estimation performance of the tracker. As shown in the results, the number of targets increases in 0 s, 10 s and 20 s, the OSPA of the GLMB fluctuated, however, the Gibbs-GLMB fluctuated more strongly than GLMB. When the number of targets decreases in 80 s and 90 s, the Gibbs-GLMB results have a large fluctuation, because the three parallel targets are relatively close, the target 3 is false detection; In 87 s-90 s, false cardinality estimates occurs due to the symmetrical geometric relationship between the sensors and the targets. while the GLMB stays a little high but remains relatively stable.

Scenario 2
The target simulation time is 100 s in the Figure 11. The legend in the pictures is the same as that in scenario 1. It also can be seen from the two pictures that the target tracks obtained from the Gibbs-GLMB filtering and the GLMB filtering are basically consistent with true trajectory of the targets. In the Figure 12, the blue box and the red circle are the estimated position of the particle point and the true position of the target point, respectively. We can see that the targets can be effectively detected through the particles after some steps even the two targets overlap in the same position at beginning.
From the simulation results in Figures 11-14, it can be seen that Gibbs-GLMB and GLMB are able to track the target by acoustic and quickly detect the new-born target. But target 3(x 3 = [0 m, 0.8 m/s, 95 m, −0.5 m/s] T ) with a special starting position which on a sensor position (0 m, 95 m), there will be a measurement error and the tracking will be missing detection. Through 100 times Monte Carlo(MC) simulations, the cardinality of targets is estimated as shown in Figure 15. Both algorithms can accurately estimate the cardinality of targets but the cardinality estimates has a large error at the beginning of Gibbs-GLMB. In this scenario, there is no fluctuation when the cardinality of targets changes.     Figure 16 is the OSPA distance over 100 MC runs. The statistical results can further illustrate that the proposed Gibbs-GLMB is more accurate than the GLMB throughout the tracking process, although the error is larger at the beginning of the experiment.

Performance
The filtering algorithms are run in the same PC. The configuration is as follows: CPU: Inter Core MLi5-4200H@2.80 GHz, RAM: 8 GB, using the software matlab2017b, the computation load and accuracy is as described in the following table: As shown in Table 1, the proposed Gibbs-GLMB filtering is significantly reduced in time-consuming and the cardinality estimates is more accurate compared to the GLMB filtering. The Gibbs-GLMB filtering improves the speed of operation greatly and reduces the complexity of the procedure by integrating updating and prediction of the GLMB filter into one step and by combining the Gibbs sampling algorithm to evaluate the updated target combination, eliminating the target combination with smaller possibility and retaining the target combination with larger weight so the tracking accuracy will be better, and the computational complexity is also significantly reduced due to the reduction of the target combination.

Conclusions
By introducing the multi-sensor acoustic array and signal detection model, we proposes to use TDOA and AOA measurements, combined with the Gibbs-GLMB filter to track multiple acoustic sources. In this paper, we use RFS theory which can solve the loss of correspondence between set elements with labels and use PHAT combined with the GCC algorithm which improves the result of TDOA calculation through real acoustic signals. The feasibility of the method is verified by tracking multiple nonlinear moving targets. The experimental simulation results show that the Gibbs-GLMB filter can effectively track multi-target but the sensor position will affect the results of the tracking. Compared with the GLMB filter, Gibbs-GLMB filter runs faster and the results are more accurate. The method proposed in this paper is only implemented under ideal simulation conditions. In the future, we will consider applying it to real audio experiments and design an effective sensor array distribution.