Improved Multitarget Tracking in Clutter Using Bearings-Only Measurements

Multitarget tracking in clutter using bearings-only measurements is a challenging problem. In this paper, a performance improved nonlinear filter is proposed on the basis of the Random Finite Set (RFS) theory and is named as Gaussian mixture measurements-based cardinality probability hypothesis density (GMMbCPHD) filter. The GMMbCPHD filter enables to address two main issues: measurement-origin-uncertainty and measurement nonlinearity, which constitutes the key problems in bearings-only multitarget tracking in clutter. For the measurement-origin-uncertainty issue, the proposed filter estimates the intensity of RFS of multiple targets as well as propagates the posterior cardinality distribution. For the measurement-origin-nonlinearity issue, the GMMbCPHD approximates the measurement likelihood function using a Gaussian mixture rather than a single Gaussian distribution as used in extended Kalman filter (EKF). The superiority of the proposed GMMbCPHD are validated by comparing with several state-of-the-art algorithms via intensive simulation studies.


Introduction
Multitarget tracking in clutter [1] is an interesting but difficult problem needed to be investigated, especially when only bearings-only measurements are available. Two main issues should be addressed: measurement-origin-uncertainty [2] and measurement nonlinearity [3]. The measurement-origin-uncertainty shows the common situation that one cannot tell a measurement originated from a target or clutter. As the coordinates of the tracker are usually different from those of the measurements, different levels of measurement nonlinearity arises, such as, the bearings-only measurement shows high-level nonlinearity while the range-bearing measurement gives relative low-level nonlinearity. Besides, targets may not be detected by the passive sensors thus miss-detection problem needs to be considered. Furthermore, the target tracking using bearings-only measurements is further complicated by the observability problem, in which the motion of the observer is suggested to outmaneuver the targets in order to satisfy the observability condition [4,5].
For multitarget tracking in clutter, the measurement-origin-uncertainty problem has been addressed by many classical and emerging methods. One popular traditional approach is the joint probability data association (JPDA) [1,6]. It enumerates and probabilistically evaluates every feasible measurement-to-target association event, and the target states are then estimated by using the marginal association probability. However, one main limitation of JPDA is that it can only track known and fixed number of targets. To override this limitation, joint integrated probability data association is proposed in this paper. Inspired from the idea of the Gaussian mixture measurement methodology, the GMMbCPHD algorithm approximates the likelihood function of the bearings-only measurements by a Gaussian mixture of numbers of refined Gaussian probability density functions, thereafter, rederives the update procedure of the cardinality and intensity density of the multitarget state set, eventually, the improved formulations on both the localization and position intensity update are obtained. For the sake of verifying the effectiveness and superiority of the proposed filter, numbers of classical experiments are simulated to compare the GMMbCPHD with several state-of-the-art algorithms, i.e., GMbPHDwEKF, GMbPHDwUKF, GMMbPHD, GM-based cardinality PHD with EKF (GMbCPHDwEKF) and GM-based cardinality PHD with UKF (GMbCPHDwUKF), the tracking performance is shown in terms of the optimal subpattern assignment (OSPA) [22] as well as OSPA localization and OSPA cardinality.
In Section 2, the target motion model, measurement model and clutter model are presented. The proposed GMMbCPHD filter is given in detail in Section 3. Simulation studies are presented in Section 4. Finally, a conclusion is proposed in Section 5.

Target Motion Model
In this paper, we assume that the motion of each target follows the continuous white noise acceleration motion model [1,12] given by where the state propagation matrix Φ is time-invariant, and ν k−1 is a sequence of zero mean, white Gaussian process noises with covariance Besides, T is the sampling time, I 2 is the 2 × 2 identity matrix, and q is the power spectral density (PSD) [1]. Please note that the target t with state x t k is presented as where (x t k , y t k ) and velocity (ẋ t k ,ẏ t k ) are position and velocity at time k, respectively. In a target state set X k = {x t k,1 , . . . , x t k,N k }, the number of targets N k can be random value and the sequence of target states in the set can also be random. In this paper, we assume that χ denote the single target state space and F (χ) denote the set of all finite subsets of χ, i.e., x t k,1 , . . . , x t k,N k ∈ χ with X k ∈ F (χ).

Bearings-Only Measurement Model
Let x s k denote the sensor state shown as x s k = [x s k , y s k ,ẋ s k ,ẏ s k ] . At each sampling time k, we assume that each target t can be detected with the probability P D,k . Using θ t k to denote the bearings-only measurement of target t if it is detected. The measurement generated by the passive sensor is represented as where and k is white Gaussian measurement noise with zero mean and covariance which is uncorrelated with ν k . Obviously, the measurement set Θ k (x t k ) generated by each target t can be either {θ t k } or ∅ where denote that the target t is detectable or not detected, respectively.

Clutter Measurement Model
In practical situations, some false measurements originated from clutter may also be detected by the sensor and form an RFS K k at each time k. These false measurements are always assumed to be uniformly distribute in the sensor detection space. Finally, the bearings-only measurements RFS Z k at time k can be expressed as Furthermore, we use Z k = {Z 1 , . . . , Z k } to denote the measurement collection from time 1 to time k.

The Proposed GMMbCPHD Filter
In the framework of the RFS-based methods, the Bayesian recursion of multitarget posterior density propagates in time as where f k|k−1 (X k |X) denotes the transition density of multitarget state and g k (Z k |X k ) represents the measurement likelihood. Furthermore, p k X k |Z k denotes the multitarget state posterior density, with µ s (dX) denoting an proper reference measure on F (χ) [18]. Obviously, the recursion in (9) and (10) is difficult to be calculated. Thus a first moment approximation was proposed in [18] and named as PHD filter. The PHD propagates the intensity density of the multitarget state density and the general form is given by (without target spawning): where v k|k−1 (x) and v k|k (x) denote the predicted intensity and the updated intensity, respectively, P S,k (η) is the target survival probability from time k − 1 to k. f k|k−1 (x|η) and g k (z|x) represent the single target transition density and measurement likelihood function, respectively. γ k (x) denotes the prior intensity of spontaneous target births and κ k (z) is the clutter intensity at time k.
To improve the estimation accuracy of the number of targets, except the intensity density, the CPHD [23] also propagates the cardinality density: where: numbers with E 0 (Z) = 1 by convention; p k−1|k−1 denotes the posterior cardinality distribution at time k − 1; p γ,k is the prior cardinality distribution of spontaneous births at time k; p κ,k is the cardinality distribution of clutter at time k; P S,k (x) is the probability of survival at time k; f k|k−1 (x k |x k−1 ) is the target transition density from time k − 1 to time k; γ k (x) is the prior intensity of spontaneous target births at time k; P D,k (x) is the probability of detection at time k; g k (z k |x k ) is the likelihood function; κ k (z) is the clutter intensity at time k.

Cardinality and Intensity Prediction of GMMbCPHD
The cardinality distribution of multitarget is propagated from time k − 1 to time k, given by Based on the Gaussian mixture assumption, the posterior intensity at time k − 1 can be represented as Assume that θ, r, c, and s denote the bearing, range, course and speed in polar coordinates, respectively. The function Ψ (x; γ k (θ, r, c, s)) is utilized here to interpret the transformation of the density γ k (θ, r, c, s) from polar coordinates to Cartesian coordinates. Then, the predicted intensity at time k is given by where In Equation (23), P S,k (x) denotes the survived probability and is usually independent of the target state, for simplicity, it is abbreviated as P S,k . In Equation (24), ω b k denotes the expected number of targets that are born at time k, U (θ; H θ ) denotes a uniform distribution with respect to θ over the region H θ [24],r and σ 2 r are the prior known mean range and it corresponded variance, respectively. Besides,s and σ 2 s denotes the prior known mean speed and it corresponded variance, respectively, with σ 2 c denoting the prior known course variance.x k|k−1 are the predicted target state estimates and its corresponded covariance, which are calculated based on the prediction step of the Kalman filter:x A binary variable β is utilized to augment the target state so as to distinguish the surviving components from the birth components.
where ω As mentioned in [25], in order to avoid the cardinality estimation bias, the surviving component and birth components are suggested to be separately considered.

GMM Model in GMMbCPHD
In the proposed GMMbCPHD filter, the measurement likelihood function g k (z|x) is modeled by Gaussian mixtures but not a single Gaussian distribution. As the measurement uncertainty of bearings-only measurement is non-Gaussian in Cartesian coordinates, thus the key idea of the GMM method is to model this non-Gaussian measurement uncertainty area by several refined components, where each component can be approximated by a Gaussian distribution We use the range interval [r k,min , r k,max ] and received measurement θ t k along with the standard deviation σ θ to illustrate the measurement uncertainty in Cartesian coordinates. Then the range interval is divided into A k subintervals, given by [21,26] r k,a+1 r k,a = τ k ; a = 1, . . . , A k , where τ k = r k,max r k,min Please note that each component a is determined by r k,a , r k,a+1 , θ t k − σ θ , θ t k + σ θ in polar coordinates. Based on the information ofr = (r k,a + r k,a+1 ) /2 and ∆r = (r k,a+1 − r k,a ) /2, the probability density function of measurement component a is approximated by a Gaussian distribution with meanẑ k,a and covariance R k,a defined in Cartesian coordinates, obtained bŷ where the transformation matrix As the area of each component is different, thus the component weight is chosen to be proportional to the area, given by [26] and Finally, the likelihood function (β = 0) is approximated using a Gaussian mixture where the observation matrix H = 1 0 0 0 0 1 0 0 , with the constant C k calculated by [27] From the measurement model defined in equation (3), we can know the measurement noise follows a Gaussian distribution, therefore the likelihood function defined in the polar coordinates (β = 1) can be obtained by

Cardinality and Intensity Update of GMMbCPHD
The new targets are always assumed to be detected when they are birth. For the sake of simplicity, the detection probability for surviving target is assumed to be independent from their state, which results in the definition as follows To make the mathematical derivation procedure friendlier, the notations defined below are necessary where The inner product with respect to v k|k−1 and ψ k is defined and obtained by where q k|k−1 denotes the likelihood of measurement componentẑ k,a with respect to component i, and the innovation covariance: and the updated target states and corresponding covariance are given by: with the Kalman gain: In the third step of Equation (47), as both likelihood function g k (z|x, β = 1) and the target birth model γ k (θ, r, c, s) are given in the same polar coordinates, we have: Then, where 1 H (θ) denotes the indicator function of the bearings-only measurement space region H θ and V H is the volume of H θ . Besides, the following approximation is made: The transformation function Ψ (x; ϕ k (θ, r, c, s)) is used to transform the function ϕ k (θ, r, c, s) from polar coordinate to Cartesian coordinate, which is obtained by: where the weight λ k,a is calculated by Equation (33), the position component ofx k θ t k andJ k θ t k are obtained by Equations (56) and (57), respectively, and the velocity component can be approximately calculated in the way in [28,29]:x where: P xx = (∆r) 2 sin 2 θ t k +r 2 σ 2 θ cos 2 θ t k , P yy = (∆r) 2 cos 2 θ t k +r 2 σ 2 θ sin 2 θ t k , Hence, the elementary symmetric functions are calculated over the set defined by The inner products involved in the CPHD recursions are obtained by Finally, the posterior PHD is approximated by and where In GMMbCPHD, the component management is same as that of GMMbPHD, while the state extraction is quite different as that of the GMMbPHD. Since the cardinality distribution has been estimated in the GMMbCPHD, the estimated target states can be extracted based on this cardinality distribution. The detailed extraction implementation procedure is: if the cardinality with N has the biggest probability in the distribution, N target states are extracted from the posterior intensity with the N biggest weights.

Simulation Experiments
For the sake of fair comparison as well as not losing generality, the simulation scenarios follow the exact same as those in [17]. The improvements of cardinality and localization estimation of GMMbCPHD are presented by comparing with GMbPHDwEKF, GMbPHDwUKF, GMMbPHD, GMbCPHDwEKF as well as GMbCPHDwUKF.

Experiment 1
The initial position of the bearings-only sensor installed in a maneuvering moving platform is [−4200 m, 3500 m]. To satisfy the observability condition, the sensor platform firstly moves at a speed of 5 knots and changes its course twice: the first one changes from 220 • to 60 • from the time 840th seconds to 1360th seconds, and the second time changes from 60 • to 220 • from the time 1860th seconds to 2040th seconds. In this paper, we define the positive direction of measurements is that the clockwise rotation from the positive Y-axis. In different periods, the number of targets changes along with time and Table 1 Figure 1 presents the moving condition of the sensor and the targets.  The simulation lasts 3000 s and the sampling interval is 10 s. Suppose that the measurement noise has the standard deviation σ θ = 1 • , and the survival probability of each target is P S,k = 0.98. To simulate the practical situations, we assume that the clutter distributes uniformly in the measurement space and the number follows the Poisson distribution with mean 15. As each target cannot be detected at each sampling time, the detection probability is assumed to be 0.95. For the new born targets, the prior target course and speed are set to be θ t k − π ands = 10 knots with standard deviation σ c = 50 • and σ s = 4 knots, respectively. Usually, the initial target range information is barely known, for simplicity, the initial target range information required in the UKF-based and EKF-based algorithm is assumed to follow a Gaussian distribution in this simulation, with known meanr = 12,000 m and standard deviation σ r = 4000 m and, the minimum and maximum target range information required in the GMMbPHD and GMMbCPHD are set to be [300 m, 18, 000 m], respectively. The number of measurement components at each scan is predefined to be A k = 8. In each filter, the birth intensity is assumed to be ω b k = 0.05. and the biggest number of Gaussian components is M k = 100.

Experiment 2
To show the advantages of proposed filter in more challenging situation, this experiment considers a higher clutter measurement density as well as a lower target detection probability compared to the experiment 1, i.e., the averaged number of clutter measurements increases to 30 while the target detection probability decreases to 0.85. All remained parameters are exactly same as those in experiment 1. Figure 2 presents the demonstration of the measurements for experiment 2.

Experiment 3
Some preselected parameters for each filter are changed in this simulation experiment. The survival probability P S,k and the birth intensity ω b k are changed to 0.95 and 0.01, respectively. Furthermore, the biggest number of Gaussian components M k increases to 150. The remained parameters in this simulation experiment are same as those in experiment 1.

Experiment 4
To illuminate the effectiveness and superiority of our proposed method in the scenario of maneuvering target tracking, this experiment extra adds another two maneuvering targets (Targets #6 and #7) to the scenario in experiment 1, as depicted in Figure 3. Targets #6 and #7 are birth in the point [1000 m, −8000 m] and [3000 m, 3000 m], respectively. The traveling courses of Target #6 is changed from 350 • to 270 • at 1500th seconds and from 180 • to 240 • at 1680th seconds for target #7. Targets #6 survives from 200th seconds to 2800th seconds, and target #7 survives from 400th seconds to 3000th seconds. The speeds of both maneuvering targets are 8 knots.

Simulation Results
The average OSPA with order 2 and cutoff 400 is used as a metric to evaluate performances of each filter over 500 Monte Carlo runs. The simulation results of experiment 1, experiment 2, experiment 3 and experiment 4 are shown in Figures 4-7, respectively. Basically, the OSPA error consists of two parts (OSPA localization error and OSPA cardinality error). To fully show the benefits of the filters, both the OSPA distance and the OSPA components (OSPA localization and cardinality) are also shown in each experiment result. The execution time to show the computational load of each filter in experiment 1 is given in Figure 8.  In all simulation experiments, as can be seen obviously from Figures 4a, 5a, 6a and 7a, the proposed GMMbCPHD filter delivers much lower OSPA distance error than other filters. As the nonlinear bearings-only measurement is approximated by Gaussian mixtures using GMM technique which is more accurate than that is approximated by a single Gaussian distribution, the localization estimation of targets in GMMbPHD and GMMbCPHD is much more accurate than EKF and UKF-based filters, such as GMbPHDwEKF, GMbPHDwUKF, GMbCPHDwEKF and GMbCPHDwUKF as shown in Figures 4b, 5b, 6b and 7b. From Figures 4c, 5c, 6c and 7c, we can see that the CPHD-based filters deliver much lower OSPA cardinality error than those of PHD-based filters such as GMbPHDwEKF, GMbPHDwUKF and GMMbPHD as the cardinality distribution has been modeled and propagated properly in GMbCPHDwEKF, GMbCPHDwUKF and GMMbCPHD.
Compared results of experiment 1 and experiment 2, we can see that the OSPA distance, OSPA localization and OSPA cardinality in experiment 1 are more accurate than those in experiment 2. As experiment 2 has much higher clutter density and lower target detection probability than those in experiment 1, the tracking situation in experiment 2 is more challenging than that in experiment 1. Fortunately, the proposed GMMbCPHD also performs much better than other filters. From results of experiment 3, we can see that GMMbCPHD has stable performance than other filters as the changed preselected parameters do not have much influence on the the proposed filter. However, these preselected parameters have negative effects on filters such as GMbCPHDwEKF and GMbCPHDwUKF. When two extra maneuvering targets are considered, the OSPA errors of all filters are increased as shown in Figure 7. However, the proposed GMMbCPHD also shows the best performance compared with other filters.
In all simulation experiments, all compared filters present a large OSPA distance error and OSPA localization error before 1000 s as the target states are not observable. After the maneuvering motion of sensor, the observability condition is satisfied and the OSPA distance and localization errors decrease sharply. Based on the discussion above, the improvements of proposed GMMbCPHD reflect in both target localization estimation and target cardinality estimation. As GMM technique can improve the localization tracking accuracy, the GMMbPHD and GMMbCPHD show advantages on the OSPA localization. Besides, estimating cardinality distribution can improve the cardinality estimation of multitarget states and the GMMbCPHD gives the benefits on OSPA cardinality. However, the GMM technique and cardinality distribution estimation will increase the computational load of the filter. Figure 8 shows the execution times of all filters and we can see that GMMbPHD and GMMbCPHD filters have a little bit heavier load than other filters. As the measurement likelihood function is modeled by single Gaussian distribution in GMbCPHDwEKF and the likelihood function of each measurement is modeled using eight Gaussian measurement components in proposed GMMbCPHD, the computational time of GMMbCPHD is almost eight times heavier than that of GMbCPHDwEKF as shown in Figure 8. Please note that each measurement component is used to update predicted intensity component using Kalman filter. Similarly, the GMbCPHDwUKF approximates the measurement likelihood function using nine sigma points, but these sigma points are utilized in only one Kalman filter cycle. Thus, the computational time of GMbCPHDwUKF just shows a little bit heavier than that of GMbCPHDwEKF. For simulation scenarios with more number of targets or clutter measurements, the computational loads of all methods obviously increase as more intensity components will be generated and more Kalman filter cycles are performed in algorithms. Though the GMM technique takes much computational sources in the GMM-based filters, fortunately, the execution time of proposed GMMbCPHD shows much smaller than the real time. Please note that all simulation cases are executed on the platform with a 2.4 GHz Intel Core i5 CPU, Windows 7, and Matlab.

Conclusions
In this paper, we proposed a random finite set (RFS)-based filter, Gaussian mixture measurements-based cardinality probability hypothesis density (GMMbCPHD), which delivers significant benefits from two aspects: Gaussian mixture measurements (GMM) and cardinality distribution estimation. The GMM models the nonlinear bearings-only measurement likelihood using Gaussian mixtures but not single Gaussian distribution. The GMM is used in the proposed GMMbCPHD and each measurement component is used to update the predicted intensity components. Thus, the number of updated intensity components is also larger than that in GMbCPHDwEKF and GMbCPHDwUKF. In traditional GMbCPHD, the likelihood function is modeled by single Gaussian distribution; however, that is approximated using Gaussian mixtures in GMMbCPHD. Each measurement is associated with each predicted intensity component with state updated using Kalyan filter. Furthermore, the predicted cardinality distribution is also updated using measurement components of all received measurement. As the GMM is more precise to model the likelihood function, the estimation of updated intensity and cardinality show significant improvement. The GMMbCPHD also propagates the cardinality distribution estimation which can improve the target number estimation in multitarget tracking. In different simulation scenarios with different parameters, simulation results show that the proposed filter performs much better than other filters on both optimal subpattern assignment (OSPA) localization and OSPA cardinality.