Improved Bearings-Only Multi-Target Tracking with GM-PHD Filtering.

In this paper, an improved nonlinear Gaussian mixture probability hypothesis density (GM-PHD) filter is proposed to address bearings-only measurements in multi-target tracking. The proposed method, called the Gaussian mixture measurements-probability hypothesis density (GMM-PHD) filter, not only approximates the posterior intensity using a Gaussian mixture, but also models the likelihood function with a Gaussian mixture instead of a single Gaussian distribution. Besides, the target birth model of the GMM-PHD filter is assumed to be partially uniform instead of a Gaussian mixture. Simulation results show that the proposed filter outperforms the GM-PHD filter embedded with the extended Kalman filter (EKF) and the unscented Kalman filter (UKF).

1. Introduction

Bearings-only multi-target tracking (MTT) [1,2,3] with clutter and missed detections is a challenging nonlinear problem. The filter for such a problem should not only deal with the nonlinearity of measurements, but also solve the measurement origin uncertainty of multiple targets. Bearings-only measurements [4] are often obtained by passive sensors, such as sonars, and contain a low level of information about detected targets, leading to the observability problem [5]. For bearings-only single-target tracking, the observable level can be evaluated by the Cramer–Rao lower bound (CRLB) [6,7]. However, it is difficult to evaluate the level of observability in a theoretical context for the bearings-only multi-target tracking with clutter and missed detections [8]. In order to overcome the non-observable problem, the sensor needs to outmaneuver the targets [5,9]. As the origin of a measurement is not known for target tracking in clutter, the tracking filter needs to have a proper track-to-measurement-association method [1] or something equivalent, such as symmetric measurement equations [10] and random finite sets [11].

For tracking targets using nonlinear measurements, the most popular method is the extended Kalman filter (EKF) [12]. It linearizes the nonlinear measurement function around the predicted target state under the assumption that the true target state is close enough to the predicted state. Then, the measurement update of the filter can be performed using the Kalman filter (KF) [12], which has a good performance for linear systems. The unscented Kalman filter (UKF) [13] applies the unscented transform instead of linearization to the nonlinear measurement function. The probability density functions (pdfs) are represented and propagated using several sigma points. The particle filter (PF) [14] represents the nonlinear pdfs by considering the amount of particles and can obtain better performances than the EKF and the UKF. However, particle filters have a heavy computational load. Apart from these nonlinear filters, many other filters exist for nonlinear systems, such as the cubature Kalman filter (CKF) [15], the shifted Rayleigh filter (SRF) [16] and the ensemble Kalman filter (EnKF) [17].

To handle the measurement origin uncertainty problem in MTT with clutters, many techniques have been developed. Many of these methods belong to the following categories: joint probability data association (JPDA) [18], multiple hypothesis tracking (MHT) [19,20] and random finite set (RFS) [11] based methods. The JPDA filters come from combing the probability data association (PDA) [1] filter with joint events and can only track fixed and known numbers of targets. To accommodate varied and unknown numbers of targets, the joint integrated PDA (JIPDA) [21] filter, which is able to estimate the probability of target existence, is proposed. To improve the tracking accuracy, a multi-scan multi-target tracking algorithm, the joint integrated track splitting (JITS), was developed in [22]. As the JIPDA and the JITS suffer from heavy computational load when tracking targets in mutual proximity, the linear multi-target IPDA (LMIPDA) and the linear multi-target ITS (LMITS) were proposed in [22,23], respectively. Besides, the iterative JIPDA (iJIPDA) [24] tracker is also a computationally-efficient algorithm for the MTT. The MHT filters attempt to maintain and evaluate a set of measurement hypotheses with high track scores. There are many versions of the MHT filter, and most of them can be grouped into two classes: the track-oriented class [25] and the measurement-oriented class [19]. The JPDA and MHT approaches are formulated via data association, which takes a large proportion of computational resources, while the RFS approach is an emerging paradigm and is established without data association. The RFS-based filters [26,27,28] treat multi-target states and measurements as the state finite set and the measurement finite set, respectively. The RFS-based filters try to estimate the target set based on the measurement set. In this way, the RFS-based filters are performed in a computational efficient way for the multi-target tracking.

Among RFS-based filters, the probability hypothesis density (PHD) [26] filter is a first moment approximation to the multi-target predicted and posterior densities. It propagates the target intensities without considering data associations between targets and measurements. As there is no close form to the PHD filter, a sequential Monte Carlo (SMC) or Gaussian mixture (GM) technique was implemented, which resulted in the SMC-PHD [29,30] filter or the GM-PHD [31] filter. For bearings-only MTT, it is difficult to apply the SMC-PHD filter because of the problem of extracting estimated target states. For this reason, the GM-PHD filter is considered in this paper. The GM-PHD filter associated with the EKF and the UKF are termed as GM-PHD-EKF and GM-PHD-UKF, respectively.

In the GM-PHD-EKF and the GM-PHD-UKF, the predicted intensity is modeled by a Gaussian mixture, and the likelihood function [12] is approximated using a single Gaussian distribution. The updated intensity of each filter is obtained after the predicted intensity is updated by the measurements. As bearings-only measurements suffer from severe nonlinearity, the likelihood function approximated by a single Gaussian is not accurate enough. Furthermore, the accuracy of the updated intensity cannot get enough improvements after the predicted intensity is updated by the inaccurate likelihood function of the measurements. In this paper, the Gaussian mixture measurements-PHD (GMM-PHD) filter is proposed to address bearings-only MTT with clutter and missed detections. To improve the tracking accuracy, the target intensities are approximated by Gaussian mixtures, and the likelihood function is also modeled by a Gaussian mixture in the GMM-PHD filter. In this way, the updated intensity of the GMM-PHD is much more accurate than those of the GM-PHD-EKF and the GM-PHD-UKF after the predicted intensity is updated by a more accurate likelihood function, which is modeled by Gaussian mixtures. The proposed filter is a nonlinear MTT algorithm that can address not only bearings-only measurements, but also other nonlinear measurements. For bearings-only MTT, the targets may appear at any place in the measurement space. In order to reduce the number of selected parameters and avoid the undesirable effect of poor parameter selection [8], the derivation of the GMM-PHD filter is processed with a partially-uniform target birth model [8,32], but not a Gaussian mixture birth model. In the simulation experiment, the performance of the GMM-PHD is compared with the GM-PHD-EKF and the GM-PHD-UKF in terms of the optimal subpattern assignment (OSPA) [33] distance, the OSPA localization, the OSPA cardinality and the average CPU time.

The rest of the paper is arranged as follows. The bearings-only MTT problem is presented in Section 2. The proposed GMM-PHD filter is derived in Section 3. A simulation study is given in Section 4. Finally, a conclusion is proposed in Section 5.

2. The Bearings-Only MTT Problem

In this paper, bearings-only MTT with one passive sensor on a maneuvering platform is studied. The targets obey the continuous white noise acceleration (CWNA) motion model [1,12]. The target motion model and the sensor measurement model in the 2D Cartesian coordinates case are considered in this section. For the target t, the target state with position and velocity at time k is expressed as:

An RFS of target states can contain any number of targets at time k and is given by:

Let denote the single target state space and denote the set of all finite subsets of , i.e., and .

As the target t is assumed to follow the CWNA motion model, its dynamic model can be expressed as: where the state propagation matrix Φ is time-invariant,

is a sequence of zero mean, white Gaussian process noises with covariance:

T is the sampling time, is the identity matrix and q is the power spectral density (PSD) [1]. The sensor state is given as:

Each target t can be detected with the probability at time k. The sensor can generate a measurement when the target is detected. We use to denote the target measurement at time k for target t. The measurement from the radar is: where: is zero mean, white Gaussian measurement noise with covariance: that it is uncorrelated with . Therefore, each target t can generate a measurement RFS , which can be either (the target t is detectable) or ∅ (the target t is not detected). In addition, the sensor also produces false measurements, which form an RFS at each time k. Then, the MTT measurement RFS at time k can be expressed as:

3. The Proposed GMM-PHD Filter

In RFS-based methods, the Bayesian recursion of multi-target posterior density is propagated in time as: where and denote the multi-target transition density and likelihood function, respectively. Here, denotes the multi-target posterior, density and is an appropriate reference measure on [31].

To reduce the computational intractability, a first moment approximation (the PHD filter) of the recursion is proposed in [26]. The general form of the PHD filter recursion (without target spawning) is given by: where and denote the predicted intensity from time to k and the updated intensity, respectively, is the survival probability, meaning the target still survives at time k, represents the single target transition density from time to k, denotes the prior intensity of spontaneous target births at time k, is the single target measurement likelihood function and is the clutter intensity.

3.1. Intensity Prediction

According to the Gaussian mixture assumption, the posterior intensity at time can be expressed as:

Let denote the bearing, range, course and speed in polar coordinates, respectively. The function represents the transformation of the density from polar coordinates to Cartesian coordinates. Then, the predicted intensity at time k is given by: where:

In Equation (17), the survived probability is assumed to be independent of the target state and is given as . In Equation (18), is the expected number of targets appearing at time k, is the uniform distribution in θ over the region , is the prior mean of the range and is its prior variance. Similarly, is the prior mean of the speed associated with its prior variance , and is the prior course variance. The predicted estimates and are obtained via the Kalman prediction:

The target state is augmented by a binary variable β, so we may distinguish between the surviving components and the birth components: where . As noted in [34], the surviving and birth components are separated to avoid biasing the cardinality estimates.

3.2. GMM Likelihood Approximation

In GM-PHD filters, the likelihood function is always modeled as a single Gaussian distribution. However, it is approximated by a Gaussian mixture in the proposed GMM-PHD filter.

Though the measurement noise is assumed to be Gaussian, the measurement uncertainty is non-Gaussian (non-ellipse) in Cartesian coordinates. The GMM measurement presentation first divides the non-elliptical measurement uncertainty into several segments. Then, each segment is modeled by one Gaussian distribution, and the whole measurement uncertainty can be approximated by a Gaussian mixture.

Suppose the measurement uncertainty in Cartesian coordinates is determined by the range interval and the measurement with standard deviation . The range interval is divided into subintervals, given by [35,36]: where:

Then, each segment a can be determined by in polar coordinates. Let and , then the segment a is approximated by a Gaussian distribution with mean and covariance in Cartesian coordinates: where:

Obviously, the area of each segment is different, and the weight of the segment a is proportional to its area, given by [35]: and

Then, the measurement likelihood function in Cartesian coordinate is approximated as a Gaussian mixture: where, is the observation matrix and the constant is calculated as [37]:

Since the measurement noise is modeled as Gaussian, the likelihood function in polar coordinates is expressed as:

Figure 1 gives an example of the GMM presentation. In the figure, the measurement uncertainty of the GMM-PHD is approximated by six measurement components (shown as the solid ellipses), and that of the GM-PHD-EKF is modeled by a single Gaussian distribution (shown as the dashed ellipse). Actually, the measurement uncertainty of the GM-PHD-UKF is also approximated by one Gaussian distribution. The true target is displayed by a cross. It is obvious that the measurement likelihood function approximated using Gaussian mixtures in the GMM-PHD is much more accurate than the approximation with one Gaussian in the GM-PHD-EKF. Thus, the GMM-PHD can perform with a much higher tracking accuracy than the GM-PHD-EKF and the GM-PHD-UKF.

Figure 1

An example of GMM presentation.

3.3. Intensity Update

New targets are assumed to be always detected at their time of birth. Furthermore, we also assume that the detection probability for a surviving target is independent of their state and:

According to Equation (14), the posterior intensity at time k is given by: and:

In the first step of Equation (34), the likelihood function is expressed by a Gaussian distribution Equation (31), but not a Gaussian mixture Equation (29), as the target birth model is given in polar coordinates (uniform across the bearing space and Gaussian in the range and velocity). As both of them are expressed in the same polar coordinates, we have:

In the above equation, where is the indicator function of the bearing space region and is the volume of , and we use the approximation: which is reasonable in practice by assuming that is small compared to the region . The transformation function transforms the function from polar coordinates to Cartesian coordinates, given by: where the weight is given by Equation (27), the position part of the mean and the covariance are given by Equations (24) and (25), respectively, and the velocity part is calculated by the approximation [38,39]: where:

In the denominator of Equations (33) and (34), we have the following: where is the likelihood of measurement component against the component i with the innovation covariance: and the updated target states and corresponding covariance are given by: with the Kalman gain:

The following approximation for the posterior intensity can be obtained: and: where:

3.4. Component Management and State Extraction

Just like the GM-PHD filter, the number of Gaussian components of the updated intensity exponentially increases in time. To reduce the computational load, same techniques (component merging and pruning) in [27,31] are also implemented in the GMM-PHD filter. If the weight of a Gaussian component is lower than the preset threshold, it will be discarded. When some of components are close enough, they will be merged into one Gaussian component. If the total number of Gaussian components is bigger than the maximum value , only components with high probability will be retained. More details about component management can be found in [31]. The state extraction is also the same as the method in [31]. If the weights are bigger than a threshold, the corresponding states are extracted as the outputs of the filter.

4. Simulation Experiments

To illustrate the improvements of the proposed GMM-PHD filter, we compare it with the GM-PHD-EKF and the GM-PHD-UKF, which also use a partially-uniform birth model. The GM-PHD-EKF and GM-PHD-UKF were introduced in [8] for a challenging bearings-only multi-target scenario.

4.1. Simulation Scenarios

4.1.1. Case 1

The sensor is located on a maneuvering platform with the initial position . The sensor platform moves at a speed of five knots and changes course two times to ensure the observability of the bearings-only target tracking. Note that the sensor should change course again to track the targets that appear after the second sensor maneuvering. The initial course of the sensor is and changes to from 14 min to 17 min. The second course change occurs from 31 min to 34 min, with the course changed from to . Note that the clockwise rotation from the positive Y-axis is considered to be positive. A time-varying number of targets is contained in the scenario, and the motion profile of targets is presented in Table 1. The starting positions of Targets are and , respectively. Similarly, the initial positions of Targets are and , respectively. Finally, the travel of Target starts at . An illustration of the sensor and the targets without process noise and clutter is presented in Figure 2.

Table 1

Motion profile of targets.

Target	Survival Time (s)	Course (degree)	Speed (knots)
#1	0,2400	95	8
#2	300,3000	20	7
#3	500,3000	280	8
#4	0,3000	275	7
#5	0,2700	215	10

Figure 2

The geometry of sensor and targets.

In the simulation, the sampling time of the sensor is 10 s, and the total simulation time is 3000 s. The sensor generates measurements with a standard deviation , and the survival probability of all targets is also the same, . The number of clutter measurements is assumed to follow a Poisson distribution with a mean of 15. The detection probability in Case 1 is . The clutter measurements are distributed uniformly across the whole bearing space . The PSD is .

In the target birth model, the prior selected parameters are the same for all filters. The prior target course is with standard deviation . The prior target speed is knots with the standard deviation knots. The prior target range for the GM-PHD-EKF and the GM-PHD-UKF is m with standard deviation m. The range interval for the GMM-PHD is . The number of measurement components in each scan is . The birth intensity is set to for all algorithms. The thresholds for the component merging and pruning are four and , respectively. Further, all filters have the maximum number of the Gaussian components, .

4.1.2. Case 2

In this case, a more challenging scenario with higher clutter intensity and lower detection probability is considered. Compared with Case 1, the mean number of clutter measurements is increased to 30, and the detection probability is decreased to 0.85. The other parameters of Case 2 are the same as those of Case 1. An example of the measurements of Case 2 is shown in Figure 3.

Figure 3

An example of measurements.

4.1.3. Case 3

In this case, some selected parameters of filters are changed to show whether the proposed GMM-PHD also has better performance than the GM-PHD-EKF and the GM-PHD-UKF. The survival probability in Case 3 is changed to 0.95. The birth intensity is set as 0.01 in all filters. Besides, the maximum number of the Gaussian components is increased to 150 in each filter. The other parameters of Case 3 are also the same as those of Case 1.

4.1.4. Case 4

In this case, two maneuvering Targets are added to Case 1. The starting positions of Targets are and , respectively. The course of Target is changed from to at 25 min. The initial course of Target is and changes to at 28 min. The survival times of Targets are and , respectively. The speed of two targets is eight knots. The other parameters of Case 4 are also the same as those of Case 1. An illustration of the sensor and the targets without process noise and clutter in Case 4 is presented in Figure 4.

Figure 4

The geometry of the sensor and targets.

4.2. Simulation Results and Analyses

The performances of all filters are compared with the average OSPA metric (with order two and cutoff 4000 m) over 500 Monte Carlo runs. The simulation results of Case 1, Case 2, Case 3 and Case 4 are presented in Figure 5, Figure 6, Figure 7 and Figure 8, respectively. The OSPA distance consists of two components: localization and cardinality. Figure 5a, Figure 6a, Figure 7a and Figure 8a show OSPA distances of three filters. The OSPA localizations are presented in Figure 5b, Figure 6b, Figure 7b and Figure 8b. The OSPA cardinalities are shown in Figure 5c, Figure 6c, Figure 7c and Figure 8c. The cardinality statistics are also considered and are shown in Figure 5d, Figure 6d, Figure 7d and Figure 8d. The execution time of different filters for Case 2 is also compared and shown in Figure 9.

Figure 5

Results of Case 1: (a) Optimal subpattern assignment (OSPA) distance; (b) OSPA localization; (c) OSPA cardinality; (d) Cardinality statistics.

Figure 6

Results of Case 2: (a) OSPA distance; (b) OSPA localization; (c) OSPA cardinality; (d) Cardinality statistics.

Figure 7

Results of Case 3: (a) OSPA distance; (b) OSPA localization; (c) OSPA cardinality; (d) Cardinality statistics.

Figure 8

Results of Case 4: (a) OSPA distance; (b) OSPA localization; (c) OSPA cardinality; (d) Cardinality statistics.

Figure 9

Execution time.

In Figure 5a, Figure 6a, Figure 7a and Figure 8a, the OSPA distance of the GM-PHD-UKF is slightly lower than that of the GM-PHD-EKF. The GMM-PHD performs significantly better than the other two filters. Before 1000 s, the target states are unobservable, and all of the filters have a large OSPA distance error. After the sensor maneuver, the target states become observable, and the OSPA distance decreases sharply. From 1500 s to 2500 s, the OSPA distances of all filters show a little bit of an increase as the measurements overlap, which makes it hard for the filters to estimate the target states.

In Figure 5b, Figure 6b, Figure 7b and Figure 8b, the OSPA localization of the GMM-PHD is much lower than those of the GM-PHD-EKF and the GM-PHD-UKF after the target states become observable. Especially in Case 4, which contains two maneuver targets, the OSPA localization errors of the GM-PHD-EKF and the GM-PHD-UKF have significant increases after target maneuvering, as shown in Figure 8b. However, the proposed GMM-PHD can keep a high tracking accuracy as the likelihood function is modeled by Gaussian mixtures, which can quickly responds to the target course change. The performance of the GM-PHD-UKF is also better than the GM-PHD-EKF on the OSPA localization metric in all cases.

In Case 1, Case 3 and Case 4, the three filters have almost the same OSPA cardinality error as shown in Figure 5c, Figure 7c and Figure 8c. Additionally, the three filters have accurate cardinality estimates in Figure 5d, Figure 7d and Figure 8d. For Case 2, the GMM-PHD has some improvement of the OSPA cardinality criteria compared to the GM-PHD-EKF and the GM-PHD-UKF in Figure 6c. Furthermore, the cardinality statistic of the GMM-PHD is more accurate than those of other two filters in Figure 6d. From Figure 5d, Figure 6d, Figure 7d and Figure 8d, three filters suffer from a delay in detecting a new target as the targets are not observable before the sensor maneuver.

According to the discussion above, the improvement of the GMM-PHD is mainly reflected in increasing the tracking accuracy of the target states, but not the cardinality estimation of the targets, for Case 1, Case 3 and Case 4. When the simulation situation becomes worse (the clutter intensity becomes bigger, and the probability of target detection decreases) in Case 2, the GMM-PHD shows improvement on both the track accuracy and the cardinality estimate. As the measurement likelihood is approximated by a Gaussian mixture, the number of track components in the GMM-PHD filter is bigger than for the GM-PHD-EKF and the GM-PHD-UKF. The execution time of the GMM-PHD filter is slightly larger than the other two filters, as shown in Figure 9. However, the execution times of these filters are much smaller than the real time. Note that the computational load of each filter becomes heavier if the number of targets increases, as the maximum number should be set bigger to ensure good performances.

5. Conclusions

For bearings-only multi-target tracking with clutter and missed detections, the PHD filter is popular and known to be effective. As there is no closed form solution to the original PHD recursion, the Gaussian mixture assumption is applied, and the GM-PHD filter was designed. However, the likelihood function of the GM-PHD filter is always modeled by a Gaussian distribution. In this paper, the GMM-PHD filter was proposed in a way that the likelihood function, as well as the posterior intensity are modeled by Gaussian mixtures. In this way, the likelihood function of the GMM-PHD can be approximated more accurately than those of the GM-PHD-EKF and the GM-PHD-UKF. The simulation results show that the proposed algorithm has significant improvement of target state estimation for challenging bearings-only multi-target tracking scenarios at the expense of computational resources. Besides, the targets appear at several fixed points in the standard form of the GM-PHD filter. The GMM-PHD was derived by using a partially uniform target birth model that can reduce the number of parameters that must be chosen by the end user.