Simplex Cubature Kalman-Consensus Filter for Distributed Space Target Tracking

A simplex cubatureKalman-consensus filter, which is suitable for distributed space target tracking usingmultiple radars, is proposed to improve the target tracking accuracy. The detailed orbital dynamics model and radar measurement model are given as the system filtering models.The intractable nonlinear Gaussian weighted integral in the filter is decomposed into the spherical integral and radial integral, which are calculated using the spherical simplex rule and the second-order Gauss-Laguerre quadrature rule, respectively. In this way, a new simplex cubature rule is derived. By means of the statistical linear regression method, the posterior mean, covariance, and the transmitted messages in the extended Kalman-consensus filter are approximated using the deduced simplex cubature rule, which results in the proposed simplex cubature Kalman-consensus filter. No data fusion center exists in the filter, and each radar only needs to exchange the information with its neighbors to reach a consensus estimate. The simulation results show that the proposed filter can achieve more accurate results compared with the cubature Kalman-consensus filter.


Introduction
In recent years, the orbital resources are becoming increasingly strained on account of the increase in the number of space targets, which should be monitored to improve the utilization efficiency of the space resources.Ground-based radar is a significant sensor in the space surveillance system, it provides round the clock working capability, and it is a key technology to use its measurement data for space target tracking [1].
The space target tracking can be considered as a nonlinear system state estimation problem, and two methods, including the particle filter and the nonlinear Kalman filter, are mainly taken.Particle filter is a type of Monte Carlo based filter, which is theoretically applicable to arbitrary nonlinear non-Gaussian system without any assumptions on the posterior probability density function (PDF) [2].However, in practical applications, there appear some problems containing particle degradation and depletion, large amount of calculation, and the choice of the importance function that will affect the filtering accuracy as well as the computational efficiency [3].Therefore it is not suitable for application in the occasion that requires high real-time performance [4].
In the nonlinear Kalman filter, the posterior PDF is assumed to be Gaussian distribution, and the suboptimal estimate of the nonlinear system state is obtained.The extended Kalman filter (EKF) [5] is the most widely used nonlinear Kalman filter in the past several decades, it uses the multidimensional Taylor series expansion to linearize the nonlinear function, and then the standard Kalman filter is applied.For arbitrary nonlinear transformation, it can be seen using the Taylor series that the linearized mean can only match to the first order of the true one; hence the EKF is considered to achieve the first-order filtering accuracy.Moreover, the accuracy and the numerical stability of EKF are reduced for the strong nonlinear system, and the calculation of Jacobian matrix imposes further restrictions on the system models.
Based on the assumption that the approximation to the posterior PDF is easier than that of arbitrary nonlinear function, Julier [6] proposed the unscented Kalman filter (UKF).UKF is a deterministic sampling filter; that is, the posterior mean and covariance are calculated using the sampling points, which are generated using the nonlinear propagation of sigma points selected as certain criterion, and the thirdorder accuracy can be achieved [7].UKF is a derivationfree filter, which can commendably overcome the defects of EKF.However, there exist some tunable parameters in sigma points, the selection of which lacks rigorous mathematical basis, and the negative weight on the center point may reduce the numerical stability for high-dimensional system [8,9].
By means of the coordinate transformation, Arasaratnam [10] proposed the cubature Kalman filter (CKF), in which the key intractable integral is decomposed into the spherical integral and radial integral.These two integrals are approximated using different numerical methods, respectively, and result in the spherical-radial cubature rule, which is used to calculate the posterior mean and covariance in the nonlinear Kalman filter framework.CKF can be regarded as a special case of UKF with the parameter  = 0 [11], whereas CKF gives the rigorous mathematical reason why  should be chosen zero for the first time.Furthermore, CKF has higher numerical stability than UKF and has a wide range of applications in engineering [12,13].
In order to further improve the estimation accuracy of CKF, Wang [14] introduced the transformation group of the regular simplex to compute the spherical integral and proposed the spherical simplex-radial cubature Kalman filter (SSRCKF).The simulation results show that SSRCKF can achieve more accurate results compared with CKF for high-dimensional system.Bhaumik [15,16] adopts -order generalized Gauss-Laguerre quadrature rule to calculate the radial integral and puts forward the cubature quadrature Kalman filter (CQKF).It is pointed out that CQKF is the generalized form of CKF; that is, CKF is the simplified form of CQKF with the first-order Gauss-Laguerre quadrature rule.Jia [17] proposes the higher-order cubature Kalman filter; it improves the filtering accuracy with much more points needed.
However, the space target tracking accuracy obtained by a single radar cannot satisfy the demands in some practical applications; thus the information fusion by multiple radars should be considered.In tradition, the information is centralized fused; namely, all radars send their measurement data to a data fusion center, on which the centralized filter is carried out.In this mode, the whole system will collapse once the center failure occurs.With the development of distributed sensor network technology, distributed filter has become a subject undergoing intense study [18,19].Olfati [20,21] gives a computational framework of the Kalmanconsensus filter (KCF) for linear system, in which multiple sensors synchronously observe the same target, and the estimation errors among local nodes are eliminated through information collaborative interaction, such that the state estimates of various sensors reach a consensus.Furthermore, a consensus-based distributed Kalman filter for state estimation in a sensor network considering the random failures of the connections is proposed in [22], and a Kalman filter type consensus+innovations distributed linear estimator of the linear time-invariant dynamical systems is developed in [23].These two distributed filters have shown better performance compared with KCF in some cases; however, they are designed for linear systems.For nonlinear system, Pellett [24] proposes the extended Kalman-consensus filter; nevertheless, the filter holds the inherent defects of EKF, that is, to calculate the Jacobian matrix and achieve the first-order accuracy.
In this paper, a simplex cubature Kalman-consensus filter (SCKCF) is put forward to improve the distributed space target tracking accuracy by multiple radars.First, for the intractable nonlinear Gaussian weighted integral, the spherical simplex rule and the second-order Gauss-Laguerre quadrature rule are utilized to calculate the spherical integral and radial integral, respectively, and a new simplex cubature rule is derived.Then, combined with the statistical linear regression method, the SCKCF is proposed in the extended Kalman-consensus filter framework.Different from the centralized filter, there is no data fusion center in the SCKCF, and each radar only exchanges information with its neighbors, which can effectively improve the system's fault tolerance and scalability.The simulation results show that the proposed SCKCF can achieve higher orbit determination accuracy than cubature Kalman-consensus filter.
The rest of this paper is organized as follows: the mathematical models for space target tracking are provided in Section 2. The new simplex cubature rule is derived in Section 3. The simplex cubature Kalman-consensus filter is proposed in Section 4. The simulation results and analysis are presented in Section 5.The conclusion is given in Section 6.

Mathematical Models for Space Target Tracking
In this section, the space target orbital dynamics model and the radar measurement model, which are considered as the state equation and measurement equation in the filter, respectively, are given below.

The Orbital Dynamics Model.
The space target orbit is described in the J2000 earth inertial coordinate system ( − , shown in blue in Figure 1), and the orbital dynamics model with  2 nonspherical gravitational perturbation is given as follows: where   = (,,) Τ and  V = (V  , V  , V  ) Τ denote the position and velocity of the space target, respectively.The parameters  2 , , and   represent the harmonic coefficient, the earth gravitational constant, and the radius of the earth, respectively.The perturbation (  ,   ,   ) Τ is the sum of the high-order nonspherical perturbation, threebody gravitational perturbation, and solar radiation pressure perturbation, which can be considered as the zero mean white Gaussian noise in this study.
Equation ( 1) can be written in the following discrete form by using the fourth-order Runge-Kutta method: where ) Τ denotes the orbit state at time ,  −1 represents the process noise at time  − 1, and f(⋅) is the nonlinear function in (1).

The Measurement Model.
The measurement model is established in the radar horizon coordinate system ( −  ℎ  ℎ  ℎ , the yellow one in Figure 1), in which the orbital state is expressed as And the geometrical relationships between the range   , velocity Ṙ  , azimuth angle   , elevation angle   of radar, and orbital state are given as follows: However, the orbital model and the measurement model are established in  −  and  −  ℎ  ℎ  ℎ , respectively, and it is required to perform a coordinate transformation, the detailed process of which is listed in Appendix.Consequently, (3) can be written in the following measurement equation form.
where Z  = (  , Ṙ  ,   ,   ) Τ denotes the measurement element at time , k  represents the measurement noise at time , and h(⋅) is the nonlinear function in (3).

Simplex Cubature Rule
For nonlinear systems, as shown in ( 2) and ( 4), the posterior mean value cannot be propagated directly using the nonlinear function, and it results in the crucial problem in nonlinear Kalman filter being the calculation of the integral of the form "nonlinear function × Gaussian PDF".Specifically, that is to calculate the integral   = ∫ R  g(x)(x; x, P  )x, where x ∈ R  , g(x) denotes the arbitrary nonlinear function, and (x; x, P  ) represents the Gaussian distribution with mean x and covariance P  .Owing to it being difficult to achieve the analytical solution, finding the high accuracy integral rule for numerical approximation becomes momentous, and the cubature rule is taken into account in this section.

Spherical Simplex Rule for Spherical
Integral.The integral of the form  g = ∫ R  g(x) exp(−x Τ x)x is considered first, and it needs to be transformed into the spherical-radial coordinate system to compute the integral  g numerically.For this purpose, let x = s, where s = ( 1 ,  2 , ⋅ ⋅ ⋅ ,   ) Τ represents the direction vector such that s Τ s = 1 and  = √ x Τ x ∈ [0, ∞) denotes the radius.Then the integral  g can be decomposed into the spherical integral () and the radial integral R as follows [10], which are easier to approximate numerically compared with  g . where denotes the spherical surface, and (⋅) represents the area element on   .
The spherical integral () is approximated using the following third-order spherical simplex rule consisting of 2+ 2 points [14].
where   and ω represent the quadrature points and the corresponding weights, respectively. denotes the number of the points.Points   can be obtained by solving the following -order Chebyshev-Laguerre polynomials.
Weights ω are calculated using the following: 3.3.Simplex Cubature Rule for Approximation.Through substituting (7) and S() = (√) into ( 9), the radial integral  is written as follows: It can be proved that   has the following equivalent form: Combined with (13), we achieve the following generalized simplex cubature rule to approximate the nonlinear Gaussian weighted integral.
It can be seen from ( 14) that 2( + 1) points with corresponding weights are needed.For  ≥ 3, the difficulty in finding a general expression between  and   , ω will bring some trouble; meanwhile, with the third-order spherical simplex rule, the improvement in accuracy caused by the higher-order Gauss-Laguerre quadrature rule is limited in practical applications.Hence, we choose  = 2, and ( 10) and (11) turn into the following forms, respectively.
The solutions of  1 ,  2 , ω1 , ω2 can be solved from ( 15) and ( 16) as follows: The above rule consists of 4(+1) points, and the cubature points x() and weights   in (18) are given as follows: where matrix a = [a 1 , a 2 , ⋅ ⋅ ⋅ , a +1 ] and [⋅]  indicates the -th column of the matrix.Take  = 3 for example; a is the matrix given as follows. [

Distributed Simplex Cubature Kalman-Consensus Filter
The centralized method and the distributed one in the space target tracking using multiple radars are generally illustrated in Figures 2(a) and 2(b), respectively.It is obvious that the main difference of these two methods is the existence of the fusion center, which decides the information fusion mode.Suppose the radars in Figure 2(b) are in a wireless network, and if a message is transmitted between two radars, we use  to denote the sending radar and  to denote the receiving and calculating radar.The radars that have the ability to communicate with the radar  are defined as its neighbor radars, which are denoted as   , and the set   =   ∪  represents   with the radar  itself included.Here, the following two assumptions are given as in [20].
(1) The radars have the ability to communicate with each other, and the links are stable and reliable.
(2) The information can be transmitted and processed within one filter cycle.
The state equation and the measurement equation, which are given in ( 2) and ( 4), respectively, are rewritten in the general nonlinear forms, and the following nonlinear system is taken into account.
where  denotes the models running on the radar .x , ∈ R  and z , ∈ R  are the state and measurement vectors, respectively, while w , and k , denote the white Gaussian noise with the covariance being Q , and R , , respectively.For linear system, the measurement model in ( 22) reduces to z , = H  x , + k , , where H  denotes the measurement matrix.The information vector and matrix that exchanged among the neighbor radars are defined using the measurement matrix H  in the Kalman-consensus filter [21].Hence, for nonlinear system, the nonlinear function z , = h  (x , ) + k , needs to be linearized to find the equivalent expression of H  , and the proposed simplex cubature rule should be introduced to improve the filtering accuracy without derivation.
The statistical linear regression method is adopted, and the resulting linearized function is more accurate in a statistical sense than simply using a first-order truncated Taylor series [25].For clarity, the subscripts are omitted in the absence of ambiguity, and the objective is to find the linear estimator of z = h(x) + k, such that ẑ = Hx + b + k.The matrix  H and vector b are determined by minimizing the criterion function {H, b} = arg min Ε{‖h(x) − Hx − b‖ 2  2 }, through the partial derivative with respect to  to zero and the gradient with respect to  to zero, respectively, we obtain H = P Τ  P −1 , where P  denotes the cross covariance, and P represents the covariance of x.
By means of the proposed simplex cubature rule and the equation H , = P Τ ,, (P − , ) −1 (with the subscript), the proposed SCKCF is derived based on the KCF framework, and the specific calculation steps are given below.
Step 2 (time update).The posterior state x+ ,−1 and covariance P + ,−1 at time  − 1 are used instead of x and P  in (19), respectively, to calculate the cubature points x() ,−1 , which are propagated by nonlinear function f(⋅) as follows.
The prior state estimate x− , and the covariance P − , are calculated using the propagated points X ()  , and weights   , where   are given in (20).
Step 3 (measurement update).The prior state x− , and covariance P − , at time  are used instead of x and P  in (19), respectively, to calculate the cubature points x() , , which are propagated by nonlinear function h  (⋅) as follows.
The predicted measurement ẑ, and cross covariance P ,, are approximated using the propagated points Z ()  , and weights   , where   are given in (20).
Compute the information vector and matrix of the sending radar: where  , = z , − ẑ, denotes the residual.
Step 4 (broadcast and receive message).The message m , = (u , , U , , x− , ) is broadcasted to the neighbor radars, meanwhile receiving the same defined messages m , = (u , , U , , x− , ) from its neighbors.Step 5 (information fusion).The received messages are fused to obtain matrix S , and vector g , . Step where  = /(1 + ‖P − , ‖  ),  is a small constant, and ‖ ⋅ ‖  is the Frobenius norm of the matrix.
By means of the information exchange among the neighbor radars, the estimate of each radar, that is, x+ , , can reach a asymptotical consensus through the above distributed filter [21].

Simulation Results and Analysis
A space target tracking simulation is given in this section to show the performance of the proposed SCKCF.The space target simulator, which runs the high-precision orbit propagation algorithm and provides the orbit simulation data, is shown in Figure 3.The perturbations, mainly including the high-order nonspherical gravitational perturbation, the atmospheric drag, the solar radiation pressure, the threebody gravity, and the tide perturbation, are taken into account, where the atmospheric drag coefficient   = 2.2, and the area-to-mass ratio of satellite is 0.02 m 2 /kg.Furthermore, the Jacchia-Roberts model is used as the atmospheric density model and the solar radiation pressure coefficient   = 1.The space target runs in sun-synchronous low earth orbit, the orbit epoch is 1 July, 2016, 12:00:00 (UTC), and the six orbit elements are  = 6978.14km, = 0,  = 97.79∘ , Ω = 279.76∘ ,  = 0 ∘ , and  = 0 ∘ , respectively.
The latitudes and longitudes as well as the communication topology of the six radars are shown in Figure 4.The measurement errors of the range, velocity, and the angles are 60m, 0.1m/s, and 0.02 ∘ , respectively.The minimum elevation angle is 10 ∘ , and the access time between the space target and all the radars is from 04:41:50 to 04:48:00 on 2 July, 2016, for a total of 370s.As shown in Figure 4, each radar only exchanges information with its two neighbor radars, so there is no information fusion center in the system and the failures of any one cannot result in the collapse of the space target tracking system.
In this simulation, the proposed SCKCF is compared with the CKF performed by single radar, the cubature Kalman-consensus filter (CKCF), and the centralized filter.The parameter  = 0.01, the initial filter state x+ ,0 , and the initial covariance matrix P + ,0 are given as follows: The filtering interval is 1s, that means the data should be collected and processed in 1s.The metrics used to compare the space target tracking accuracy of various filters is the root mean square error (RMSE).The Monte Carlo simulation is conducted 200 times, and the RMSEs of the position and velocity are shown in Figures 5-10.The logarithmic scale is used on the y-axis in order to improve the readability.From the figures, we can see that the RMSEs of the CKCF and the proposed SCKCF are significant smaller than that of the CKF performed by single radar, indicating that the multiple radars information fusion could obtain higher space target tracking accuracy.Moreover, the RMSEs of the centralized filter are the smallest, indicating that the distributed filters could acquire the suboptimal estimates compared with the centralized one.
The mean RMSEs from 200s to 370s are calculated and summarized in Tables 1 and 2. As shown in the figures, the CKF and the centralized filter achieve the worst and best accuracy, respectively.As for the two distributed filters, the proposed SCKCF achieves more accurate results compared with the CKCF, because SCKCF adopts the proposed simplex cubature rule while the CKCF adopts the conventional cubature rule.In the new simplex cubature rule, the spherical integral and radial integral are calculated using the spherical simplex rule and second-order Gauss-Laguerre quadrature rule, respectively; that is, two points are used in the radial integral while only one point is used in that of the conventional cubature rule, and the spherical simplex rule has       better performance than spherical rule in this application; thus the proposed simplex cubature rule outperforms in approximating the nonlinear orbital states.

Conclusions
To improve the distributed space target tracking accuracy by multiple radars, a simplex cubature Kalman-consensus filter is proposed in this paper.The intractable integral is approximated using a new simplex cubature rule, which is more accurate compared with the traditional cubature rule.By means of the statistical linear regression method, the posterior mean, covariance, and the transmitted messages in the Kalman-consensus filter are approximated using the deduced simplex cubature rule, which results in the proposed SCKCF.There is no data fusion center in the filter, and each radar needs to exchange the information only with its neighbors to reach a consensus estimate.The simulation results show that the proposed SCKCF can achieve more

Figure 5 :
Figure 5: RMSEs of the space target tracking of radar 1.

Figure 6 :
Figure 6: RMSEs of the space target tracking of radar 2.

Figure 7 :
Figure 7: RMSEs of the space target tracking of radar 3.