Abstract
A novel fifth-degree cubature Kalman filter is proposed to improve the accuracy of real-time orbit determination by ground-based radar. A fully symmetric cubature rule, approaching the Gaussian weighted integral of a nonlinear function in general form, is introduced, and the specific points and weights are calculated by matching the monomials of degree not greater than five with the exact values. On the basis of the above rule, a novel fifth-degree cubature Kalman filter, which can achieve a higher accuracy than UKF and CKF, is derived under the Bayesian filtering framework. Then, to describe the nonlinear system more accurately, the orbital dynamics equation with J2 perturbation is used as the state equation, and the nonlinear relationship between the radar measurement elements and orbital states is built as the measurement equation. The simulation results show that, compared with the traditional third-degree algorithm, the proposed fifth-degree algorithm has a higher accuracy of orbit determination.
1. Introduction
With the increasing number of satellites launched into orbit every year, the monitoring and cataloguing of satellites play an important role in improving the rate of utilisation of space resources and alleviating the pressure on orbit resources. As a type of sensor in space surveillance systems, ground-based radar is equipped without considering the influence of the weather and other special circumstances, and the use of its measurement data for real-time orbit determination is a key technology in space target tracking [1, 2]. Due to the nonlinearity of the satellite orbital dynamics model with the influence of orbital perturbation, the essence of orbit determination in real-time is to achieve the optimal estimation of the orbital state by means of nonlinear filtering technology under the Bayesian framework using the measured ranging, velocity, and angle data with measurement noise, which has important research value.
The core issue in nonlinear Kalman filtering is to calculate the intractable multidimensional vector integral such as the “nonlinear function Gaussian probability density function (pdf),” for which it is difficult to achieve the analytical solution [3, 4]. At present, two methods, including the approximation of the nonlinear function and the approximation of the Gaussian pdf, are mainly taken. In the former method, the nonlinear function is approximated by the polynomial and results in an extended Kalman filter (EKF) [5, 6], divided difference Kalman filter (DDKF) [7], and polynomial Kalman filter (PKF) [8, 9], where the first-order Taylor expansion, the multidimensional Stirling interpolation, and polynomials including Chebyshev and Fourier-Hermit are adopted to approximate the nonlinear function, respectively, in EKF, DDKF, and PKF. However, the aforementioned methods tend to be restricted when the system has strong nonlinearity with high dimensionality. For the latter, the Gaussian pdf is approximated using the deterministic sampling approach, which mainly includes the unscented transform (UT) and spherical-radial rule (SRR). Then, the unscented Kalman filter (UKF) [10, 11] and cubature Kalman filter (CKF) [12–14] are obtained by embedding UT and SRR into the Bayesian filtering framework, respectively, these have a wide range of applications in engineering [15–20], but these two types of algorithm have only third-degree filtering accuracy, which is required to be further improved. In this paper, a novel fifth-degree cubature Kalman filter is proposed without differential operation to improve the filtering accuracy from third degree to fifth degree. The integral points with corresponding weights in the general cubature rule are calculated by matching the monomials of degree no more than five with their exact values in the fully symmetric region. Then, the proposed filtering method, which can achieve a higher accuracy compared to that with UKF and CKF, is deduced by embedding the novel fifth-degree cubature rule into the Bayesian filtering framework. The proposed filtering algorithm is applied in real-time orbit determination, and a more accurate orbit estimate is achieved.
The remainder of the paper is organised as follows: the traditional cubature Kalman filter is described in Section 2, the novel fifth-degree cubature rule and related Kalman filter are derived in Section 3, the mathematical models used for orbit determination are described in Section 4, the numerical experiment and results are presented in Section 5, and the conclusions are drawn in Section 6.
2. The Traditional Cubature Kalman Filter
The core problem in any nonlinear Kalman filter is to calculate the integral , for which, in general, it is difficult to find the analytical solution, where denotes the arbitrary function, denotes the Gaussian pdf. Specifically, an integral of the form in the Cartesian coordinate system is considered. Let with , where denotes the direction vector and denotes the radius, so that and then the integral can be rewritten in a spherical-radial coordinate system as follows:where is the surface of the sphere defined by and is the area element on . Thus, the integral is decomposed into spherical integral and radial integral , respectively, and approximately represented using numerical integration as follows:where denote the integral points and weights of the spherical integral and denotes the number of integral points. Similarly, denote the integral points and weights of the radial integral, and denotes the number of points. From the third-degree spherical-radial cubature rule used in [12], we obtain thatwhere is the surface area of the unit sphere, is the Gamma function, and denotes the unit vector with the th element being 1. is substituted into , to get
Due to identity equation
it may be seen that
The calculation process used in the traditional CKF is listed as follows.
Time Update. Evaluate the cubature points .
Evaluate the propagated cubature points .
Estimate the predicted state .
Estimate the predicted error covariance .
Measurement Update. Evaluate the cubature points .
Evaluate the propagated cubature points .
Estimate the predicted measurement .
Estimate the measurement covariance matrix .
Estimate the cross-covariance matrix .
Estimate the Kalman gain .
Estimate the updated state .
Estimate the corresponding error covariance .
From the algorithm we see that 2n points are adopted when approximating the Gaussian pdf. To improve the filtering accuracy, more points, with corresponding weights, are needed.
3. Fifth-Degree Cubature Rule and Cubature Filtering Algorithm
3.1. Fifth-Degree Cubature Rule
The integral , for which it is difficult to find the analytical solution, can be approximated using the cubature rule by selecting the appropriate cubature points and corresponding weights, where denotes the cubature points, denotes the weights that do not depend on the integrand, and L denotes the number of cubature points. We will write to denote an arbitrary point in real n-dimensional space. By a monomial of degree d, we mean a function of the form , where the indices are nonnegative integers such that . The following definitions and lemma are introduced.
Definition 1 (see [21]). is a region in -dimensional space; given , the fully symmetric set of , , is the set of all points , where is any permutation of .
Definition 2 (see [21]). A region is said to be fully symmetric if and only if implies .
Definition 3 (see [21]). An integration rule R is said to be fully symmetric if and only if, whenever is an abscissa of the rule R, every element of is an abscissa of R and the same weight corresponds to all abscissas belonging to a given fully symmetric set.
Lemma 4 (see [21]). A fully symmetric rule applied to a fully symmetric -dimensional region is of degree if and only if it is exact for all monomials of degree of the form
The following cubature rule is considered:
The above rule is fully symmetric; therefore, it will be of degree five if it is exact for the monomials 1, , , and ; thus the following equations are obtained:
For the case , we have , where the Gamma function satisfies and ; thus we obtain
Formula (21) is combined with formula (22) to solve the following parameters as
Thus the specific form of rule is achieved by substituting formula (23) into formula (20), and the total number of cubature points required for the rule is . Furthermore, the integral is written using the rule in the following form:
The following vectors are defined:where denotes the th column of the matrix. From formula (5), the fifth-degree cubature rule approximating the Gaussian weighted integral of the nonlinear function is obtained by combining formula (24) with the vectors defined in formulae (25) as follows:
3.2. Fifth-Degree Cubature Kalman Filter
The following discrete nonlinear dynamic system is considered:where denotes the state vector, denotes the measurement vector, and are known nonlinear functions, and the process noise and measurement noise are uncorrelated zero mean Gaussian white noise with covariance matrixes and , respectively.
The cubature points and corresponding weights are obtained from the fifth-degree cubature rule shown in formula (26) as follows:
The proposed fifth-degree cubature Kalman filter is deduced by using the points and weights, shown, respectively, in formulae (28), under the Bayesian filtering framework with reference to the third-degree algorithm in [12], and the specific calculation steps are listed as follows.
Step 1 (filter initialisation). Cycle , and complete the following steps.
Step 2 (time update). Calculate the points using the a posteriori state estimation and a posteriori error covariance matrix .Calculate the nonlinear propagation of using .Calculate the a priori state estimation by weighted merging .Calculate the a priori error covariance matrix .
Step 3 (measurement update). Calculate the points using the a priori state estimation and a priori error covariance matrix .Calculate the nonlinear propagation of using .Calculate the predicted measurement value by weighted merging .Calculate the predicted measurement covariance matrix .Calculate the cross-covariance matrix .Calculate the Kalman filtering gain .Calculate the a posteriori state estimation .Calculate the a posteriori error covariance matrix .The pseudocode representing the proposed algorithm is given in Algorithm 1.
| ||||||||||||||||||||||||||||||||||||||||||||||
Remark 5. The proposed method is differential free; that is, there is no need to calculate the Jacobian matrix.
Remark 6. Compared with CKF of third degree, the filtering accuracy of the proposed method is improved to fifth degree.
Remark 7. The general method of computation of the cubature rule is given in the proposed filtering method, without dividing the intractable integral into spherical integral and radial integral.
4. Mathematical Model for Orbit Determination
4.1. Orbital Dynamics Model
Satellites in orbit are subjected to various perturbations, mainly including nonspherical gravitational perturbation, third body gravitational perturbation, atmospheric drag perturbation, and solar radiation pressure perturbation, among which the J2 nonspherical gravitational perturbation is the most influential perturbation. To describe the in-orbit motion of the satellite more accurately, the orbital dynamics model with J2 perturbation in the earth central fixed (ECF) coordinate system is used as follows to describe the orbit of the satellite [22]:where denotes the position and velocity of satellite in ECF, denotes the earth gravity constant, denotes the harmonic coefficient, denotes the radius of the earth, denotes the angular velocity of the earth, and is the sum of other perturbations, which can be approximated as zero mean Gaussian white noise in the study. Formula (42) can be written in discrete state equation form as follows:where denotes the orbital state at time and denotes the process noise.
4.2. Radar Measurement Model
The radar measurement model is described in local horizontal (LH) coordinate system, and the transformation matrix from ECF to LH is as follows:where and denote the geocentric longitude and the geocentric latitude of radar, respectively, which can also be represented by geocentric coordinates . We define as the satellite orbital states in LH, and the following equations are obtained:
The geometric relationship between orbital states and radar measurement values is obtained as follows:where denotes the ranging value, denotes the velocity value, denotes the azimuth angle, and denotes the elevation angle. Formula (46) is written in the following measurement equation form:where denotes the measurement values at time and denotes the measurement noise.
5. The Numerical Experiment
The simulation platform is built using the Satellite Tool Kit (STK) and MATLAB, the satellite runs in low-earth sun-synchronous orbit, and the reference orbit data is generated by the High-Precision Orbit Prediction (HPOP) algorithm. In HPOP, the 21-order earth gravity model is taken into account, the atmospheric drag coefficient , the Jacchia-Roberts model is adopted as the atmospheric density model, the solar radiation pressure coefficient , the area-mass ratio is 0.02 m2/kg, the third body gravitational perturbations of sun and moon are considered, and the tidal perturbation is considered. The accuracy of ranging, velocity, and angle values is assumed to be 20 m, 0.1 m/s, and 0.015°, respectively. The six orbital elements, including semimajor axis (), eccentricity (), inclination (), Right Ascension of Ascending Node (RAAN), argument of perigee (), true anomaly (), and the latitude and longitude of the radar station, are listed in Table 1.
The initial filter state = .
The initial covariance matrix = .
The access time from radar station to satellite is from 1 July, 2015, 16:14:00 to 1 July, 2015, 16:21:00, and the root mean square error (RMSE) is adopted to evaluate the real-time orbit determination results. The filtering cycle is 1 s, and we ran 200 Monte-Carlo simulations. The UKF and CKF are compared in this experiment to validate the performance of the proposed algorithm. The RMSEs of the three algorithms are shown in Figure 1, and the statistical average RMSEs are summarised in Table 2. From the results it may be seen that the orbit determination accuracy obtained by UKF is almost consistent with that of CKF due to the two algorithms being made to adopt the third-degree deterministic sampling method. By contrast with CKF, the proposed fifth-degree cubature Kalman filter is capable of achieving a higher accuracy, and the position accuracy is increased by 3.204 m with velocity accuracy increased by 0.04 m/s. For a low-earth-orbit satellite, the atmospheric drag perturbation has an influence on the orbit, meaning that the state equation (formula (42)) cannot describe the orbit exactly, and if the high-precision orbit perturbation model is used, the excessive computation demand will impair the filtering algorithm’s real-time performance; however, the errors caused by the orbital model are generally acceptable due to the access time of the LEO satellite being short.
(a) Position RMSRs
(b) Velocity RMSEs
6. Conclusion
In this paper, a novel fifth-degree cubature Kalman filter is proposed to improve the accuracy of real-time orbit determination by ground-based radar. The integral points and weights in the general cubature rule are solved by matching the monomials with degree not greater than five with their exact values, and then the fifth-degree cubature rule is deduced. The proposed novel fifth-degree cubature Kalman filter, which can achieve a higher filtering accuracy than UKF and CKF, is derived by using the aforementioned rule based on the Bayesian filtering framework. The simulation results show that the position accuracies achieved by CKF and the proposed algorithm are 27.148 m and 23.944 m, respectively, with the velocity accuracies being 0.347 m/s and 0.307 m/s, respectively. Compared with the results of CKF, the position accuracy and velocity accuracy are improved by 3.204 m and 0.04 m/s, respectively, thus verifying the validity of the proposed algorithm.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work is supported partly by the National High Technology Research and Development Program of China (2015AA7026085).