Sensor Fusion of Gaussian Mixtures for Ballistic Target Tracking in the Re-Entry Phase

A sensor fusion methodology for the Gaussian mixtures model is proposed for ballistic target tracking with unknown ballistic coefficients. To improve the estimation accuracy, a track-to-track fusion architecture is proposed to fuse tracks provided by the local interacting multiple model filters. During the fusion process, the duplicate information is removed by considering the first order redundant information between the local tracks. With extensive simulations, we show that the proposed algorithm improves the tracking accuracy in ballistic target tracking in the re-entry phase applications.


Introduction
The problem of ballistic target tracking in the re-entry phase has attracted much attention due to both of its theoretical and practical significance. Technically speaking we need to develop a stochastic nonlinear filter for state estimation with respect to the ballistic target dynamics. Practical applications are in the fields of surveillance for safety against the re-entry of space debris produced by old satellites or spacecraft at the end of their lifetime [1][2][3]. Tracking a ballistic target has been identified as a stressful filtering problem due to the strong nonlinearity exhibited by the forces acting on the target [2,4,5]. Besides, the nonlinearity also results from the bearing and range measurements which are given with nonlinear measurement functions in relation to the target [6]. In the literature several nonlinear filtering methods have been applied to this tracking problem under the assumption that the ballistic coefficient is already known [7], including extended Kalman filter [8], unscented Kalman filter [9], ensemble Kalman filter [10], cubature Kalman filter [11] and particle filter [3]. The tracking problem becomes more difficult when the ballistic coefficient is known crudely or totally unknown. In such a case, an approach which has been investigated to be effective is the multiple model approach [12][13][14][15]. The structure of multiple model filter consists of a bank of local filters, with each filter modelled with a different value of ballistic coefficient. The output of the multiple model filter is in terms of Gaussian mixtures obtained by weighted mixing of the estimate from each local filter.
The technology of track-to-track fusion (T2TF) makes it possible to achieve further improvement of tracking accuracy than in the single sensor case. A basic T2TF system [16] consists of two sensors and a fusion center. Each sensor generates a local track by observing the target of interest and estimating with filtering. The fusion center fuses local tracks to obtain a fused track with improved accuracy. There has been a great deal of work in developing T2TF algorithms [17][18][19][20], however, most of these algorithms have been designed for systems with Gaussian density based filters (e.g., Kalman filter, extended Kalman filter and unscented Kalman filter) and may not be applicable to systems with multiple model filters. Demand for tracking ballistic targets with unknown ballistic coefficients poses a new requirement on using multiple model filters, which makes it necessary to fuse local tracks in where x k−1 = x p,k−1 , x v,k−1 , y p,k−1 , y v,k−1 T is the target state vector at time step k − 1. x p,k−1 and y p,k−1 represent positions, x v,k−1 and y v,k−1 represent velocities in Cartesian coordinates (x, y). w k is a zero-mean white Gaussian process noise vector with covariance Q given by: where q is a parameter related to the process noise intensity. The nonlinear function Ψ (·) in Equation (1) is given by: where matrices F and G are given by: Sensors 2016, 16, 1289 3 of 17 and f (·) is the aerodynamic drag given by where β is the ballistic coefficient, g is the gravity acceleration, ρ (·) is the air density [25] which decays exponentially with the altitude y p,k−1 as

Sensor Configurations and Measurement Model
The sensor configuration involves two ground-radar sensors S 1 and S 2 , which are located at x S l R , y S l R for l = 1, 2 respectively. The measurement equation is given by: where z is the measurement noise for sensor S l . The sensor measures the range r S l k and the bearing θ k in relation to the target. Thus the measurement model can be written as where

Local Filter
Since the actual ballistic coefficient is unknown to the sensors, an interactive multiple model (IMM) filter [26] with r = 2 models is applied at the single sensor level. The models are based on the extended Kalman filter (EKF), but having different value of ballistic coefficient β S l i , i = 1, · · · , r. In general, the IMM filter merges state estimates computed under each possible model using local filters, with each filter using a different combination of the previous model-conditioned estimates [27]. The recursion equations for sensor S l consist of four major steps in each cycle as follows.
Step 1: Calculation of the mixing probabilities for i, where c Step 2: Mixing for j = 1, · · · , r:  (5) and (8) are shown in Appendixes A and B respectively.
Step 4: Mode probability update for j = 1, · · · , r A block diagram for one cycle of the IMM filter with two models (r = 2) for sensor S l is given in Figure 1.

Fusion Architecture
Among the various choices for multiple model filters, the IMM filter has been shown to be one of the most effective schemes for hybrid systems [28]. In this study we assume that the local sensors run IMM filters for handling unknown ballistic coefficient, therefore the output of each local sensor is a Gaussian mixtures model. In this regard, the track-to-track fusion (T2TF) problem we consider here involves the fusion of Gaussian mixtures densities. As shown in Figure 2, the fusion architecture presented includes two ground-radar sensors S 1 and S 2 with local measurements obtained periodically and local tracks updated synchronously. In addition, the fusion center may transmit the latest fused track back to one local sensor. When the local sensor receives the fused track, it will be used to replace the local one. In practice, this condition usually happens when the fusion center is collocated with one local sensor [29]. In this study, we assume that after each fusion, the local track of sensor S 1 is replaced with the fused track from the fusion center, while sensor S 2 operates independently without information feedback from the fusion center.

Basic Fusion Process and Redundant Information
The basic fusion process [19] is formulated as: where     c p x , which is to be subtracted out. While the removal of duplicate information is straightforward in the theoretical formulation [30], identification of duplicate information for distributed estimation system can be difficult in practical implementation. In this paper, the first order redundant information [16,31] is considered. In this way, the fusion center only needs to keep track of the previous data received from sensor S2 at the previous communication time step and remove it when fusing the current estimates from sensor S1 and S2. Specifically, when both    Sensors 2016, 16, 1289 5 of 18 Figure 1. Structure of the IMM algorithm for sensor S1.

Basic Fusion Process and Redundant Information
The basic fusion process [19] is formulated as: where     c p x , which is to be subtracted out. While the removal of duplicate information is straightforward in the theoretical formulation [30], identification of duplicate information for distributed estimation system can be difficult in practical implementation. In this paper, the first order redundant information [16,31] is considered. In this way, the fusion center only needs to keep track of the previous data received from sensor S2 at the previous communication time step and remove it when fusing the current estimates from sensor S1 and S2. Specifically, when both

Basic Fusion Process and Redundant Information
The basic fusion process [19] is formulated as: where p f (x) is the fused estimate, p S 1 (x) and p S 2 (x) are the estimates from local sensors S 1 and S 2 respectively, c is a normalization constant. The common information between the local estimates is given in the denominator by p c (x), which is to be subtracted out. While the removal of duplicate information is straightforward in the theoretical formulation [30], identification of duplicate information for distributed estimation system can be difficult in practical implementation. In this paper, the first order redundant information [16,31] is considered. In this way, the fusion center only needs to keep track of the previous data received from sensor S 2 at the previous communication time step and remove it when fusing the current estimates from sensor S 1 and S 2 . Specifically, when both p S 1 (x) and p S 2 (x) are Gaussian mixtures, namely: where r = 2 since the IMM filter is with two models for each sensor. Then the common information at time step k can be obtained as: where: and f S 2 J,i,k (·) is the Jacobian calculated at the estimated statex S 2 0i,k−1|k−1 .

Fusion of Gaussian Mixtures
In a distributed fusion problem, assuming two Gaussian mixtures, ,k|k are to be fused with a common prior distribution p c (x). With the standard Bayesian fusion formula, i.e., Equation (15), the fused probability density function (PDF) can be obtained as: To obtain an analytical fused result and avoid the potential complexity, one idea is to approximate the denominator p c (x) with a single Gaussian PDF using moment matching [27], namely: where: With this approximation, the fusion expressions to be applied in IMM filter which is the major contribution of the current work are derived as follows: where and with The derivation of Equations (23)-(25) is presented in Appendix C.

Gaussian Mixtures Reduction
As one can see from Equation (22), the fused Gaussian mixtures PDF has an exponentially growing number of components as more Gaussian mixtures are multiplied in the long run. Thus it is necessary to manage the components growth with a method of Gaussian mixtures reduction. Suppose that we are given a Gaussian mixtures model with n components, and we wish to approximate it with a mixture of r components (r < n). In general, the Gaussian mixtures reduction algorithm can be operated in the following manner.
While more than r components remain, choose two components that in a sense to be least dissimilar and replace them by their moment matching merge as: where (w i ,x i , P i ) and w j ,x j , P j are two weighted Gaussian components to be merged, and (w m ,x m , P m ) is their moment matching approximation. For the dissimilarity measure between two components of a Gaussian mixtures model, we adopt a metric that is proposed in Reference [32] as an upper bound on the discrimination of the Gaussian mixtures after the merge from the Gaussian mixtures before the merge. The dissimilarity measure is given as follows: Thus in each iteration, a Gaussian mixtures model with n components is constructed after the Gaussian mixtures fusion, then we operate the procedure of Gaussian mixtures reduction in an iterative manner until n = r, so as to meet the requirement that the IMM filter is set to be with r models. In our case, the fused Gaussian mixtures consist of n = 4 components after fusion in Equation (22), namely: where: Note that the subscript k|k is omitted for the sake of brevity. We set (w 11 ,x 11 , P 11 ) and (w 22 ,x 22 , P 22 ) as major components, since they are consistent with model 1 and 2 in a straightforward manner. For the rest two cross-components, namely (w 12 ,x 12 , P 12 ) and (w 21 ,x 21 , P 21 ), they are supposed to be merged with one of the major components according to the dissimilarity measure given in Equation (29). More specifically, for (w 12 ,x 12 , P 12 ), if D ((w 12 ,x 12 , P 12 ) , (w 11 ,x 11 , P 11 )) < D ((w 12 ,x 12 , P 12 ) , (w 22 ,x 22 , P 22 )), then we merge (w 12 ,x 12 , P 12 ) with (w 11 ,x 11 , P 11 ), otherwise we merge (w 12 ,x 12 , P 12 ) with (w 22 ,x 22 , P 22 ); for (w 21 ,x 21 , P 21 ), if D ((w 21 ,x 21 , P 21 ) , (w 11 ,x 11 , P 11 )) < D ((w 21 ,x 21 , P 21 ) , (w 22 ,x 22 , P 22 )), then we merge (w 21 ,x 21 , P 21 ) with (w 11 ,x 11 , P 11 ), otherwise we merge (w 21 ,x 21 , P 21 ) with (w 22 ,x 22 , P 22 ). In the end, a flow chart of the proposed fusion methodology is given in Figure 3 as follows: where: Note that the subscript | k k is omitted for the sake of brevity. We set   11 11 11 , , w x P and   22 22 22 , , w x P as major components, since they are consistent with model 1 and 2 in a straightforward manner. For the rest two cross-components, namely   12 12 12 , , w x P and   21 21 21 , , w x P , they are supposed to be merged with one of the major components according to the dissimilarity measure given in Equation (29 , , w x P . In the end, a flow chart of the proposed fusion methodology is given in Figure 3 as follows:

Simulation Results
The proposed algorithm is applied in a re-entry ballistic target tracking scenario to verify the performance through a series of simulation runs. To the best author's knowledge, the Gaussian

Simulation Results
The proposed algorithm is applied in a re-entry ballistic target tracking scenario to verify the performance through a series of simulation runs. To the best author's knowledge, the Gaussian mixtures fusion methodology has not been applied in a distributed ballistic target tracking scenario for track-to-track fusion before. The system dynamic model is represented as Equation (1), where we use T = 2 s, q = 1 m 2 ·s −3 and g = 9.81 m·s 2 . The initial state x 0 is: with initial covariance: For the air density ρ (·), we have c 1 = 1.227, c 2 = 1.093 × 10 −4 for y p < 9144 m and, and c 1 = 1.754, c 2 = 1.49 × 10 −4 for y p > 9144 m. The actual target ballistic coefficient is β = 40,000 kg·m −1 ·s −2 . Figure 4 shows the target trajectory in the X-Y plane. Figures 5 and 6 show the velocity of the ballistic target, and the aerodynamic drag acceleration against time, respectively. It can be observed that the velocity decreases with the increment of aerodynamic drag f (·). with initial covariance: For the air density     , we have c1 = 1.227, c2 = 1.093 × 10 −4 for yp < 9144 m and, and c1 = 1.754, c2 = 1.49 × 10 −4 for yp > 9144 m. The actual target ballistic coefficient is β = 40,000 kg·m −1 ·s −2 . Figure 4 shows the target trajectory in the X-Y plane. Figures 5 and 6 show the velocity of the ballistic target, and the aerodynamic drag acceleration against time, respectively. It can be observed that the velocity decreases with the increment of aerodynamic drag   f  .  For the local filters, it is assumed that the actual ballistic coefficient of the target is unknown to both sensors, so that = = 60000 kg · m · s and = = 10000 kg · m · s are used as the IMM models for both sensor 1 S and 2 S . Besides, the Markov chain transition matrix was taken as: For the measurement model, we assume that the two sensors are homogeneous, in the sense that R S 1 r = R S 2 r = 100 2 m 2 and R S 1 θ = R S 2 θ = 0.05 2 rad 2 . Besides, we consider that the sensors are located at x S 1 R , y S 1 R = (0, 0) and x S 2 R , y S 2 R = (50, 000, 0), respectively. For the local filters, it is assumed that the actual ballistic coefficient of the target is unknown to both sensors, so that β S 1 1 = β S 2 1 = 60, 000 kg·m −1 ·s −2 and β S 1 2 = β S 2 2 = 10, 000 kg·m −1 ·s −2 are used as the IMM models for both sensor S 1 and S 2 . Besides, the Markov chain transition matrix was taken as:  The proposed Gaussian mixtures fusion method is verified by comparing the results between sensors S1 and S2. Note that the estimate from sensor S1 is identical to the fused result, for the reason that after each fusion the local track from sensor S1 is replaced with the fused track from the fusion The proposed Gaussian mixtures fusion method is verified by comparing the results between sensors S 1 and S 2 . Note that the estimate from sensor S 1 is identical to the fused result, for the reason that after each fusion the local track from sensor S 1 is replaced with the fused track from the fusion center. Unlike sensor S 1 , the estimate from sensor S 2 is solely obtained from one singe IMM filter, since sensor S 2 operates by itself without information feedback from the fusion center. Figures 7 and 8 show the estimated model probabilities for sensors S 1 and S 2 , respectively. It can be seen that during 0-50 s the probabilities for models 1 and 2 are equal to each other, for the reason that the aerodynamic drag is zero which makes it not able to judge which model matches better. The ballistic coefficient becomes observable from t = 50 s, at that time the aerodynamic drag begins to emerge. It can be observed that the model probabilities for sensor S 1 is more stable than that for sensor S 2 , namely the proposed fusion method improves the stability of the estimated model probabilities. center. Unlike sensor S1, the estimate from sensor S2 is solely obtained from one singe IMM filter, since sensor S2 operates by itself without information feedback from the fusion center. Figures 7 and 8 show the estimated model probabilities for sensors S1 and S2, respectively. It can be seen that during 0-50 s the probabilities for models 1 and 2 are equal to each other, for the reason that the aerodynamic drag is zero which makes it not able to judge which model matches better. The ballistic coefficient becomes observable from t = 50 s, at that time the aerodynamic drag begins to emerge. It can be observed that the model probabilities for sensor S1 is more stable than that for sensor S2, namely the proposed fusion method improves the stability of the estimated model probabilities.    We further calculate the estimated ballistic coefficient with the total probability theorem [25] as: We further calculate the estimated ballistic coefficient with the total probability theorem [25] as: where µ S l i is the estimated model probability for β S l i . Figures 9 and 10 plot the estimated value of ballistic coefficient from sensors S 1 and S 2 , respectively. It can be seen that after the drag force emerges, the estimated value of the ballistic coefficient for sensor S 1F , namely the fused one, appears to be more accurate and stable than that for sensor S 2 , in the sense that the first one ranges around [3.6,4.4] × 10 4 , and the latter one ranges around [1.5,6.0] × 10 4 . ballistic coefficient from sensors S1 and S2, respectively. It can be seen that after the drag force emerges, the estimated value of the ballistic coefficient for sensor S1F, namely the fused one, appears to be more accurate and stable than that for sensor S2, in the sense that the first one ranges around [3.6,4.4] × 10 4 , and the latter one ranges around [1.5,6.0] × 10 4 . Figure 9. Estimated value of the ballistic coefficient for sensor S1. Figure 9. Estimated value of the ballistic coefficient for sensor S 1 . Figure 9. Estimated value of the ballistic coefficient for sensor S1. In the following, we compare the position estimation accuracy in terms of root mean square error (RMSE) between the proposed Gaussian mixtures fusion method, the single IMM filter and the covariance intersection (CI) method. The CI method [18] is the most well-known fusion technique which yields consistent estimates by optimizing a nonlinear cost function associated with the fused covariance. For the reason that CI cannot be applied to Gaussian mixtures fusion directly, it requires one to approximate the output of IMM filter with the single Gaussian density using moment matching before track fusion. Figures 11 and 12 and Table 1 show the RMSE of the algorithms in X and Y axes, respectively, by performing 1000 runs of Monte Carlo simulation. It can be seen that both of the fusion methods, namely the proposed Gaussian mixtures fusion, and CI provide performance In the following, we compare the position estimation accuracy in terms of root mean square error (RMSE) between the proposed Gaussian mixtures fusion method, the single IMM filter and the covariance intersection (CI) method. The CI method [18] is the most well-known fusion technique which yields consistent estimates by optimizing a nonlinear cost function associated with the fused covariance. For the reason that CI cannot be applied to Gaussian mixtures fusion directly, it requires one to approximate the output of IMM filter with the single Gaussian density using moment matching before track fusion. Figures 11 and 12 and Table 1 show the RMSE of the algorithms in X and Y axes, respectively, by performing 1000 runs of Monte Carlo simulation. It can be seen that both of the fusion methods, namely the proposed Gaussian mixtures fusion, and CI provide performance improvement over the single IMM filter. Furthermore, the performance of the proposed Gaussian mixtures fusion is better than CI. The first reason is that we remove the first order redundant information between the local tracks during Gaussian mixtures fusion process, however, the CI only provides a conservative estimate due to its ignorance of the cross correlations between the local tracks. improvement over the single IMM filter. Furthermore, the performance of the proposed Gaussian mixtures fusion is better than CI. The first reason is that we remove the first order redundant information between the local tracks during Gaussian mixtures fusion process, however, the CI only provides a conservative estimate due to its ignorance of the cross correlations between the local tracks.       The second reason is that instead of fusing Gaussian mixtures directly as for the proposed fusion method, for CI one needs to approximate the original local Gaussian mixtures estimate with mean and covariance before track fusion, thus the resulting approximation error could degrade the estimation performance.

Conclusions
A Gaussian mixtures fusion algorithm for the track-to-track fusion problem is proposed in this study. The common information to be reduced during the fusion process is approximated with the first order redundant information between the local tracks. The proposed fusion algorithm is applied to tracking a ballistic target with unknown ballistic coefficient using IMM filters. A series of Monte Carlo simulations are conducted to evaluate the sensor fusion performance. The results indicate that the proposed algorithm improves the estimation accuracy in terms of the root mean square error. There are two issues to be addressed in the future. First, a more detailed analysis of the effects of the first order approximation must be carried out. Second, experiments should be conducted to explore the effect of using higher order approximation to formulate the common information between the local estimates.

Appendix C
The multivariate Gaussian PDF N (x;x, P) can be written in canonical notation [33] as: where: d is the dimensionality of x, |Y| det(Y). Define: i,k|k + Y S 2 j,k|k − Y c k|k x = exp ζ S 1 i,k|k + ζ S 2 j,k|k − ζ c k|k − ζ ij,k|k + ζ ij,k|k + y T ij,k|k x − 1 2 x T Y ij,k|k x = exp ζ S 1 i,k|k + ζ S 2 j,k|k − ζ c k|k − ζ ij,k|k · exp ζ ij,k|k + y T ij,k|k x − 1 2 x T Y ij,k|k x (C6) and the scaling factor is: with: