# Simplex Cubature Kalman-Consensus Filter for Distributed Space Target Tracking

Authors: Zhaoming Li; Wenge Yang; Dan Ding; Yurong Liao
Journal: Wireless Communications and Mobile Computing
 (2018)

---

## Abstract

A simplex cubature Kalman-consensus filter, which is suitable for distributed space target tracking using multiple 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.


---

## Body

## 1. 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 third-order accuracy can be achieved [7]. UKF is a derivation-free 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 n-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 Kalman-consensus 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 Section2. 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.


## 2. 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.


### 2.1. The Orbital Dynamics Model

The space target orbit is described in the J2000 earth inertial coordinate system (O-XYZ, shown in blue in Figure 1), and the orbital dynamics model with J2 nonspherical gravitational perturbation is given as follows:(1)x˙=vxy˙=vyz˙=vzv˙x=μxr3J2Rer27.5z2r2-1.5-1+fxv˙y=μyr3J2Rer27.5z2r2-1.5-1+fyv˙z=μzr3J2Rer27.5z2r2-4.5-1+fzr=x2+y2+z2where Xp=x,y,zΤ and Xv=vx,vy,vzΤ denote the position and velocity of the space target, respectively. The parameters J2, μ, and Re represent the harmonic coefficient, the earth gravitational constant, and the radius of the earth, respectively. The perturbation fx,fy,fzΤ is the sum of the high-order nonspherical perturbation, three-body gravitational perturbation, and solar radiation pressure perturbation, which can be considered as the zero mean white Gaussian noise in this study.Figure 1
The coordinate systems.Equation (1) can be written in the following discrete form by using the fourth-order Runge- Kutta method:(2)Xk=fXk-1+wk-1where Xk=Xp,kΤ,Xv,kΤ=xk,yk,zk,vx,k,vy,k,vz,kΤ denotes the orbit state at time k, wk-1 represents the process noise at time k-1, and f· is the nonlinear function in (1).


### 2.2. The Measurement Model

The measurement model is established in the radar horizon coordinate system (O-XhYhZh, the yellow one in Figure 1), in which the orbital state is expressed as X⌢k=X⌢p,kΤ,X⌢v,kΤΤ, where X⌢p,k=x⌢k,y⌢k,z⌢kΤ and X⌢v,k=v⌢x,k,v⌢y,k,v⌢z,kΤ. And the geometrical relationships between the range Rk, velocity R˙k, azimuth angle Ak, elevation angle Ek of radar, and orbital state are given as follows:(3)Rk=x⌢k2+y⌢k2+z⌢k2R˙k=x⌢kv⌢x,k+y⌢kv⌢y,k+z⌢kv⌢z,kx⌢k2+y⌢k2+z⌢k2Ak=arc⁡tan⁡y⌢kx⌢kEk=arc⁡tan⁡z⌢kx⌢k2+y⌢k2However, the orbital model and the measurement model are established inO-XYZ and O-XhYhZh, 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.(4)Zk=hXk+vkwhere Zk=Rk,R˙k,Ak,EkΤ denotes the measurement element at time k, vk represents the measurement noise at time k, and h· is the nonlinear function in (3).


## 2.1. The Orbital Dynamics Model

The space target orbit is described in the J2000 earth inertial coordinate system (O-XYZ, shown in blue in Figure 1), and the orbital dynamics model with J2 nonspherical gravitational perturbation is given as follows:(1)x˙=vxy˙=vyz˙=vzv˙x=μxr3J2Rer27.5z2r2-1.5-1+fxv˙y=μyr3J2Rer27.5z2r2-1.5-1+fyv˙z=μzr3J2Rer27.5z2r2-4.5-1+fzr=x2+y2+z2where Xp=x,y,zΤ and Xv=vx,vy,vzΤ denote the position and velocity of the space target, respectively. The parameters J2, μ, and Re represent the harmonic coefficient, the earth gravitational constant, and the radius of the earth, respectively. The perturbation fx,fy,fzΤ is the sum of the high-order nonspherical perturbation, three-body gravitational perturbation, and solar radiation pressure perturbation, which can be considered as the zero mean white Gaussian noise in this study.Figure 1
The coordinate systems.Equation (1) can be written in the following discrete form by using the fourth-order Runge- Kutta method:(2)Xk=fXk-1+wk-1where Xk=Xp,kΤ,Xv,kΤ=xk,yk,zk,vx,k,vy,k,vz,kΤ denotes the orbit state at time k, wk-1 represents the process noise at time k-1, and f· is the nonlinear function in (1).


## 2.2. The Measurement Model

The measurement model is established in the radar horizon coordinate system (O-XhYhZh, the yellow one in Figure 1), in which the orbital state is expressed as X⌢k=X⌢p,kΤ,X⌢v,kΤΤ, where X⌢p,k=x⌢k,y⌢k,z⌢kΤ and X⌢v,k=v⌢x,k,v⌢y,k,v⌢z,kΤ. And the geometrical relationships between the range Rk, velocity R˙k, azimuth angle Ak, elevation angle Ek of radar, and orbital state are given as follows:(3)Rk=x⌢k2+y⌢k2+z⌢k2R˙k=x⌢kv⌢x,k+y⌢kv⌢y,k+z⌢kv⌢z,kx⌢k2+y⌢k2+z⌢k2Ak=arc⁡tan⁡y⌢kx⌢kEk=arc⁡tan⁡z⌢kx⌢k2+y⌢k2However, the orbital model and the measurement model are established inO-XYZ and O-XhYhZh, 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.(4)Zk=hXk+vkwhere Zk=Rk,R˙k,Ak,EkΤ denotes the measurement element at time k, vk represents the measurement noise at time k, and h· is the nonlinear function in (3).


## 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 IN=∫RngxNx;x^,Pxdx, where x∈Rn, gx denotes the arbitrary nonlinear function, and Nx;x^,Px represents the Gaussian distribution with mean x^ and covariance Px. 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.


### 3.1. Spherical Simplex Rule for Spherical Integral

The integral of the formIg=∫Rngxexp⁡-xΤxdx is considered first, and it needs to be transformed into the spherical-radial coordinate system to compute the integral Ig numerically. For this purpose, let x=rs, where s=s1,s2,⋯,snΤ represents the direction vector such that sΤs=1 and r=xΤx∈0,∞ denotes the radius. Then the integral Ig can be decomposed into the spherical integral Sr and the radial integralR as follows [10], which are easier to approximate numerically compared with Ig.(5)Sr=∫Ungrsdσs(6)R=∫0∞Srrn-1exp⁡-r2drwhere Un=s∈Rn:s12+s22+⋯+sn2=1 denotes the spherical surface, and σ· represents the area element on Un.The spherical integralSr is approximated using the following third-order spherical simplex rule consisting of 2n+2 points [14].(7)Sr≈Tn2n+1∑i=1n+1grai+g-raiwhere Tn=2πn/Γn/2 denotes the surface area of the unit sphere and Γx=∫0∞exp⁡-ttx-1dt is the Gamma function. The vectors ai=ai,1,ai,2,⋯,ai,nΤ,i=1,2,⋯,n+1 are the vertices of the n-dimensional simplex, and the elements are defined as follows.(8)ai,j=-n+1nn-j+2n-j+1,j<in+1n-i+1nn-i+2,i=j0,j>i


### 3.2. Generalized Gauss-Laguerre Quadrature Rule for Radial Integral

As for the radial integralR=∫0∞Srrn-1exp⁡-r2dr, we need to transform it into the standard form. For this, let r=t and R=1/2∫0∞Sttn/2-1exp⁡-tdt is deduced. Further, let S~t=St and β=n/2-1, then R is written in the form of R=1/2∫0∞S~ttβexp⁡-tdt, which is approximated numerically using the following generalized Gauss-Laguerre quadrature rule [15].(9)∫0∞S~ttβexp⁡-tdt≈∑i=1pω~iS~tiwhere ti and ω~i represent the quadrature points and the corresponding weights, respectively. p denotes the number of the points. Points ti can be obtained by solving the following m-order Chebyshev-Laguerre polynomials.(10)Lpβt=-1pt-βexp⁡tdpdtptβ+pexp⁡-t=0Weightsω~i are calculated using the following:(11)ω~i=p!Γβ+p+1tiL˙pβti2


### 3.3. Simplex Cubature Rule for Approximation

Through substituting (7) and S~t=St into (9), the radial integral R is written as follows:(12)R=12∫0∞S~ttβexp⁡-tdt≈12∑j=1pω~jStj=Tn4n+1∑j=1pω~j∑i=1n+1gtjai+g-tjai=πn2n+1Γn/2∑j=1pω~j∑i=1n+1gtjai+g-tjaiIt can be proved thatIN has the following equivalent form:(13)IN=1πn∫Rng2Pxx+x^exp-xΤxdxCombined with (13), we achieve the following generalized simplex cubature rule to approximate the nonlinear Gaussian weighted integral.(14)IN≈12n+1Γn/2∑j=1pω~j∑i=1n+1g2tjPxai+x^+g-2tjPxai+x^It can be seen from (14) that 2n+1p points with corresponding weights are needed. For p≥3, the difficulty in finding a general expression between n and ti, ω~i 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 p=2, and (10) and (11) turn into the following forms, respectively.(15)L2βt=t-βexp⁡td2dt2tβ+2e-t=0(16)ω~i=2Γβ+3tiL˙2βti2The solutions oft1, t2, ω~1, ω~2 can be solved from (15) and (16) as follows:(17)t1=n2+1+n2+1t2=n2+1-n2+1,ω~1=nΓn/22n+4+22n+4ω~2=nΓn/22n+4-22n+4The new simplex cubature rule is derived by substituting (17) into (14) as follows:(18)IN≈n4n+1n+2+2n+4∑i=1n+1gx^+n+2+2n+4Pxai+gx^-n+2+2n+4Pxai+n4n+1n+2-2n+4×∑i=1n+1gx^+n+2-2n+4Pxai+gx^-n+2-2n+4Pxai=∑i=14n+4ωigx^iThe above rule consists of4n+1 points, and the cubature points x^i and weights ωi in (18) are given as follows:(19)x^i=x^+n+2+2n+4Pxa,-ai,i=1,2,⋯,2n+2x^+n+2-2n+4Pxa,-ai-2n-2,i=2n+3,⋯,4n+4(20)ωi=n4n+1n+2+2n+4,i=1,2,⋯,2n+2n4n+1n+2-2n+4,i=2n+3,⋯,4n+4where matrix a=a1,a2,⋯,an+1 and ·i indicates the i-th column of the matrix. Take n=3 for example; a is the matrix given as follows.(21)1-13-13-13083-23-230063-63


## 3.1. Spherical Simplex Rule for Spherical Integral

The integral of the formIg=∫Rngxexp⁡-xΤxdx is considered first, and it needs to be transformed into the spherical-radial coordinate system to compute the integral Ig numerically. For this purpose, let x=rs, where s=s1,s2,⋯,snΤ represents the direction vector such that sΤs=1 and r=xΤx∈0,∞ denotes the radius. Then the integral Ig can be decomposed into the spherical integral Sr and the radial integralR as follows [10], which are easier to approximate numerically compared with Ig.(5)Sr=∫Ungrsdσs(6)R=∫0∞Srrn-1exp⁡-r2drwhere Un=s∈Rn:s12+s22+⋯+sn2=1 denotes the spherical surface, and σ· represents the area element on Un.The spherical integralSr is approximated using the following third-order spherical simplex rule consisting of 2n+2 points [14].(7)Sr≈Tn2n+1∑i=1n+1grai+g-raiwhere Tn=2πn/Γn/2 denotes the surface area of the unit sphere and Γx=∫0∞exp⁡-ttx-1dt is the Gamma function. The vectors ai=ai,1,ai,2,⋯,ai,nΤ,i=1,2,⋯,n+1 are the vertices of the n-dimensional simplex, and the elements are defined as follows.(8)ai,j=-n+1nn-j+2n-j+1,j<in+1n-i+1nn-i+2,i=j0,j>i


## 3.2. Generalized Gauss-Laguerre Quadrature Rule for Radial Integral

As for the radial integralR=∫0∞Srrn-1exp⁡-r2dr, we need to transform it into the standard form. For this, let r=t and R=1/2∫0∞Sttn/2-1exp⁡-tdt is deduced. Further, let S~t=St and β=n/2-1, then R is written in the form of R=1/2∫0∞S~ttβexp⁡-tdt, which is approximated numerically using the following generalized Gauss-Laguerre quadrature rule [15].(9)∫0∞S~ttβexp⁡-tdt≈∑i=1pω~iS~tiwhere ti and ω~i represent the quadrature points and the corresponding weights, respectively. p denotes the number of the points. Points ti can be obtained by solving the following m-order Chebyshev-Laguerre polynomials.(10)Lpβt=-1pt-βexp⁡tdpdtptβ+pexp⁡-t=0Weightsω~i are calculated using the following:(11)ω~i=p!Γβ+p+1tiL˙pβti2


## 3.3. Simplex Cubature Rule for Approximation

Through substituting (7) and S~t=St into (9), the radial integral R is written as follows:(12)R=12∫0∞S~ttβexp⁡-tdt≈12∑j=1pω~jStj=Tn4n+1∑j=1pω~j∑i=1n+1gtjai+g-tjai=πn2n+1Γn/2∑j=1pω~j∑i=1n+1gtjai+g-tjaiIt can be proved thatIN has the following equivalent form:(13)IN=1πn∫Rng2Pxx+x^exp-xΤxdxCombined with (13), we achieve the following generalized simplex cubature rule to approximate the nonlinear Gaussian weighted integral.(14)IN≈12n+1Γn/2∑j=1pω~j∑i=1n+1g2tjPxai+x^+g-2tjPxai+x^It can be seen from (14) that 2n+1p points with corresponding weights are needed. For p≥3, the difficulty in finding a general expression between n and ti, ω~i 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 p=2, and (10) and (11) turn into the following forms, respectively.(15)L2βt=t-βexp⁡td2dt2tβ+2e-t=0(16)ω~i=2Γβ+3tiL˙2βti2The solutions oft1, t2, ω~1, ω~2 can be solved from (15) and (16) as follows:(17)t1=n2+1+n2+1t2=n2+1-n2+1,ω~1=nΓn/22n+4+22n+4ω~2=nΓn/22n+4-22n+4The new simplex cubature rule is derived by substituting (17) into (14) as follows:(18)IN≈n4n+1n+2+2n+4∑i=1n+1gx^+n+2+2n+4Pxai+gx^-n+2+2n+4Pxai+n4n+1n+2-2n+4×∑i=1n+1gx^+n+2-2n+4Pxai+gx^-n+2-2n+4Pxai=∑i=14n+4ωigx^iThe above rule consists of4n+1 points, and the cubature points x^i and weights ωi in (18) are given as follows:(19)x^i=x^+n+2+2n+4Pxa,-ai,i=1,2,⋯,2n+2x^+n+2-2n+4Pxa,-ai-2n-2,i=2n+3,⋯,4n+4(20)ωi=n4n+1n+2+2n+4,i=1,2,⋯,2n+2n4n+1n+2-2n+4,i=2n+3,⋯,4n+4where matrix a=a1,a2,⋯,an+1 and ·i indicates the i-th column of the matrix. Take n=3 for example; a is the matrix given as follows.(21)1-13-13-13083-23-230063-63


## 4. 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 Figures2(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 s to denote the sending radar and c to denote the receiving and calculating radar. The radars that have the ability to communicate with the radar c are defined as its neighbor radars, which are denoted as Nc, and the set Jc=Nc∪c represents Nc with the radar c itself included. Here, the following two assumptions are given as in [20].Figure 2
Multiple radar space target tracking.
(a)
Multiple radars centralized space target tracking (b)
Multiple radars distributed space target tracking(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.(22)xc,k=fxc,k-1+wc,kzc,k=hcxc,k+vc,kwhere c denotes the models running on the radar c. xc,k∈Rn and zc,k∈Rq are the state and measurement vectors, respectively, while wc,k and vc,k denote the white Gaussian noise with the covariance being Qc,k and Rc,k, respectively.For linear system, the measurement model in (22) reduces to zc,k=Hcxc,k+vc,k, where Hc denotes the measurement matrix. The information vector and matrix that exchanged among the neighbor radars are defined using the measurement matrix Hc in the Kalman-consensus filter [21]. Hence, for nonlinear system, the nonlinear function zc,k=hcxc,k+vc,k needs to be linearized to find the equivalent expression of Hc, 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=hx+v, such that z^=Hx+b+v. The matrix H and vector b are determined by minimizing the criterion function H,b=arg⁡min⁡Εhx-Hx-b22, through the partial derivative with respect to b to zero and the gradient with respect to A to zero, respectively, we obtain H=PxzΤP-1, where Pxz denotes the cross covariance, and P represents the covariance of x.By means of the proposed simplex cubature rule and the equationHc,k=Pc,xz,kΤPc,k--1 (with the subscript), the proposed SCKCF is derived based on the KCF framework, and the specific calculation steps are given below.Step 1 (filter initialization).
Give the initial filter valuesx^c,0+ and Pc,0+.
Cyclek=1,2,⋯, and finish the following steps.Step 2 (time update).
The posterior statex^c,k-1+ and covariance Pc,k-1+ at time k-1 are used instead of x^ and Px in (19), respectively, to calculate the cubature points x^c,k-1i, which are propagated by nonlinear function f· as follows.(23)Xc,ki=fx^c,k-1i
The prior state estimatex^c,k- and the covariance Pc,k- are calculated using the propagated points Xc,ki and weights ωi, where ωi are given in (20).(24)x^c,k-=∑i=14n+4ωiXc,ki(25)Pc,k-=∑i=14n+4ωiXc,ki-x^c,k-Xc,ki-x^c,k-Τ+Qc,kStep 3 (measurement update).
The prior statex^c,k- and covariance Pc,k- at time k are used instead of x^ and Px in (19), respectively, to calculate the cubature points x^c,ki, which are propagated by nonlinear function hc· as follows.(26)Zc,ki=hcx^c,ki
The predicted measurementz^c,k and cross covariance Pc,xz,k are approximated using the propagated points Zc,ki and weights ωi, where ωi are given in (20).(27)z^c,k=∑i=14n+4ωiZc,ki(28)Pc,xz,k=∑i=12n+2ωix^c,ki-x^c,k-Zc,ki-z^c,kΤ
Compute the information vector and matrix of the sending radar:(29)uc,k=Hc,kΤRc,k-1εc,k=Pc,k--1Pc,xz,kRc,k-1εc,k(30)Uc,k=Hc,kΤRc,k-1Hc,k=Pc,k--1Pc,xz,kRc,k-1Pc,xz,kΤPc,k--1where εc,k=zc,k-z^c,k denotes the residual.Step 4 (broadcast and receive message).
The messagemc,k=uc,k,Uc,k,x^c,k- is broadcasted to the neighbor radars, meanwhile receiving the same defined messages ms,k=us,k,Us,k,x^s,k- from its neighbors.Step 5 (information fusion).
The received messages are fused to obtain matrixSc,k and vector gc,k.(31)Sc,k=∑s∈JcUs,k(32)gc,k=∑s∈Jcus,k-Us,kx^c,k--x^s,k-Step 6 (state update).
The posterior covariancePc,k+ and state estimate x^c,k+ are calculated.(33)Pc,k+=Pc,k--1+Sc,k-1(34)x^c,k+=x^c,k-+Pc,k+gc,k+γPc,k-∑s∈Ncx^s,k--x^c,k-where γ=α/1+Pc,k-F, α is a small constant, and ·F 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^c,k+, can reach a asymptotical consensus through the above distributed filter [21].


## 5. 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 Figure3. The perturbations, mainly including the high-order nonspherical gravitational perturbation, the atmospheric drag, the solar radiation pressure, the three-body gravity, and the tide perturbation, are taken into account, where the atmospheric drag coefficient cd=2.2, and the area-to-mass ratio of satellite is 0.02 m2/kg. Furthermore, the Jacchia-Roberts model is used as the atmospheric density model and the solar radiation pressure coefficient cr=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 a=6978.14km, e=0, i=97.79°, Ω=279.76°, ω=0°, and f=0°, respectively.Figure 3
Space target simulator.The latitudes and longitudes as well as the communication topology of the six radars are shown in Figure4. 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.Figure 4
Locations and communication topology.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^c,0+, and the initial covariance matrix Pc,0+ are given as follows:(35)x^c,0+=-1485607,4323895,5272734,-327,5803,-4830ΤPc,0+=diag⁡10002,10002,10002,102,102,102The covariancesQc,k and Rc,k are given as follows:(36)Qc,k=diag⁡0.0012,0.0012,0.0012,0.000012,0.000012,0.000012Rc,k=diag⁡602,0.12,0.02π1802,0.02π1802The 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 Figures5–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.Figure 5
RMSEs of the space target tracking of radar 1.
(a)
Position RMSEs (b)
Velocity RMSEsFigure 6
RMSEs of the space target tracking of radar 2.
(a)
Position RMSEs (b)
Velocity RMSEsFigure 7
RMSEs of the space target tracking of radar 3.
(a)
Position RMSEs (b)
Velocity RMSEsFigure 8
RMSEs of the space target tracking of radar 4.
(a)
Position RMSEs (b)
Velocity RMSEsFigure 9
RMSEs of the space target tracking of radar 5.
(a)
Position RMSEs (b)
Velocity RMSEsFigure 10
RMSEs of the space target tracking of radar 6.
(a)
Position RMSEs (b)
Velocity RMSEsThe mean RMSEs from 200s to 370s are calculated and summarized in Tables1 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.Table 1
Mean position RMSEs.
Filters Radar 1/m Radar 2/m Radar 3/m Radar 4/m Radar 5/m Radar 6/m CKF 23.7473 35.4105 37.2135 23.8687 21.1308 17.2393 CKCF 6.1371 6.6633 6.9016 6.8462 6.6887 6.7597 Proposed SCKCF 5.4843 5.3668 5.1280 4.6526 4.3772 4.4636 Centralized filter 3.8374Table 2
Mean velocity RMSEs.
Filters Radar 1/(m⋅s−1) Radar 2/(m⋅s−1) Radar 3/(m⋅s−1) Radar 4/(m⋅s−1) Radar 5/(m⋅s−1) Radar 6/(m⋅s−1) CKF 0.1750 0.2317 0.2317 0.1801 0.2077 0.1999 CKCF 0.0713 0.0433 0.0411 0.0449 0.0397 0.0439 Proposed SCKCF 0.0496 0.0343 0.0317 0.0334 0.0342 0.0339 Centralized filter 0.0256


## 6. 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 accurate results compared with the CKCF, thus verifying the validity of the proposed filter. The neighboring radars are required to communicate mutually in the proposed SCKCF. However, if the link between the radars can be simplified to one-way mode without reducing accuracy, that is, the communication topology is described using the directed graph instead of the undirected graph, it can further reduce the communication burden to a certain extent. Therefore, more advanced consensus algorithms should be investigated and combined with the CKF in the next study.



---

*Source: 1476426-2018-06-24.xml*