Analysis of Filtering Methods for Satellite Autonomous Orbit Determination Using Celestial and Geomagnetic Measurement

Satellite autonomous orbit determination OD is a complex process using filtering method to integrate observation and orbit dynamic equations effectively and estimate the position and velocity of a satellite. Therefore, the filtering method plays an important role in autonomous orbit determination accuracy and time consumption. Extended Kalman filter EKF , unscented Kalman filter UKF , and unscented particle filter UPF are three widely used filtering methods in satellite autonomous OD, owing to the nonlinearity of satellite orbit dynamic model. The performance of the system based on these three methods is analyzed under different conditions. Simulations show that, under the same condition, the UPF provides the highest OD accuracy but requires the highest computation burden. Conclusions drawn by this study are useful in the design and analysis of autonomous orbit determination system of satellites.


Introduction
Orbit determination OD of satellite plays a significant role in satellite missions, aiming at estimating the ephemeris of a satellite at a chosen epoch accurately.To date, the conventional OD system is dominated by measurements based on 1 ground tracking approaches 1 such as range, range rate, and angle, and 2 Global Position System measurement 2, 3 .The orbit determination technologies have shown fair performance on various space missions.However, its high cost, lack of robustness to loss of contact, space segment degradation, and other factors promote the application of autonomous OD system, which is less costly and less vulnerable in hostile environment 4 .
In general, orbit determination is the process of estimating the satellite's state variables position and velocity by comparing in statistical sense the difference between the measurement data and the estimated data.Orbit determination system, as shown in Figure 1, usually includes sensor subsystem, model subsystem, and filter subsystem.Sensor subsystem contains sensing instruments, such as star sensor, earth sensor, and magnetometer, in order to measure and process the original measurements which are functions of state variables.Model system generates estimated data including state model and measurement model.In the filter subsystem, the optimal algorithms filtering methods process both data from sensor subsystem and from model subsystem and then estimate state variables.
Owing to the nonlinear dynamic model of satellite orbit motion, the filtering method applied in OD system should be appropriate for nonlinear system 5, 6 .Extended Kalman filter EKF , unscented Kalman filter UKF , and unscented particle filter UPF are three main methods used in satellite OD system.The EKF is based on the analytical Taylor series expansion of the nonlinear systems and measurement equations.It works on the principle that the state distribution is approximated by a Gaussian random variable.However, the Taylor series approximations in EKF introduce large errors due to the neglected nonlinearities 7 .The UKF uses the true nonlinear model and a set of sigma sample points produced by the unscented transformation to capture the mean and covariance of state, but the UKF has the limitation that it does not apply to general non-Gaussian distribution 8, 9 .The particle filter PF is a computer-based method for implementing a recursive Bayesian filter by Monte Carlo simulations.The performance of the PF largely depends on the choices of importance sampling density and resampling scheme 10, 11 .Among many improved PF methods, UPF is a hybrid of the UKF and the particle filter which uses the UKF to get better importance sampling density 12, 13 .It combines the merits of unscented transformation and particle filtering and avoids their limitations.
A variety of autonomous orbit determination methods have been proposed and explored, including a magnetometer-based OD method 14, 15 , a celestial OD method 16, 17 , a landmark OD method 18, 19 , and an X-ray pulsar OD method 20, 21 .The first two methods can be used in low earth orbit LEO satellite autonomous orbit determination system.Thus, in this paper, these two OD methods are selected for analysis.This paper is divided into five sections.After this introduction, the basic descriptions of three filtering methods in autonomous OD system are given in Section 2.Then, the state model and measurement models in OD model subsystem are described in detail in Section 3. In Section 4, simulations are shown for analyzing and comparing three filtering methods.Finally, conclusions are drawn in Section 5.

Filtering Methods
The best known algorithms to solve the problem of autonomous satellite orbit determination are the EKF, UKF, and UPF.In this section, we shall present the theories of the three filter algorithms.These algorithms will be incorporated into the filtering framework based on the dynamic state-space model as follows: where x k−1 denotes the state of the system at time k − 1, z k denotes the observations at step k, w k denotes the process noise, and v k denotes the measurement noise.The mappings f and h represent the process and measurement models.
for all k, j, and Q k is the process noise covariance at step k, R k is the measurement noise covariance at step k.

Extended Kalman Filter
A Kalman filter that linearizes about the current mean and covariance is referred to as an extended Kalman filter or EKF.The EKF is the minimum mean-square-error estimator based on the Taylor series expansion of the nonlinear functions.For example, Using only the linear expansion terms, it is easy to derive the update equations for the mean and covariance of the Gaussian approximation to the distribution of the states 12 .
The equations for the extended Kalman filter fall into two groups: time update equations and measurement update equations.The specific equations for the time and measurement updates are presented below as shown in 2.3 ∼ 2.8 22 .
1 Time Update Predicted state estimate: Predicted estimate covariance:

2.4
The time update equations project the state, x k , and covariance, P k , estimates from the previous time step k − 1 to the current time step k, Φ k is the state transition matrix at step k, which is defined to be the following Jacobians:

2.6
Updated state estimate: Updated estimate covariance: where K k is known as the Kalman gain.The measurement update equations correct the state and covariance estimates with the measurement z k .H k is the observation matrix at step k, which is defined to be the following Jacobians:

2.9
The major drawback of EKF is that it only uses the first order terms in the Taylor series expansion.Sometimes it may introduce large estimation errors in a nonlinear system and lead to poor representations of the nonlinear functions and probability distributions of interest.As a result, this filter can diverge 23 .

Unscented Kalman Filter
The unscented Kalman filter UKF 8, 24 uses the unscented transformation to capture the mean and covariance estimates with a minimal set of sample points.The UKF process is identical to the standard EKF process with the prediction-estimation recursive loop.The exception is that the UKF uses the sigma points and the nonlinear equations to compute the predicted states and measurements and the associated covariance matrices.If the dimension of state is n × 1, the 2n 1 sigma point and their weight are computed by 9 where τ ∈ R, P k | k i is the ith column of the matrix square root.The UKF process can be described as follows.
1 Time 0, initialize the UKF with x 0 and P 0 as follows: 2 Time k, define 2n 1 sigma points from

2.12
The equations for the UKF fall into two groups the same as EKF: time update equations and measurement update equations.The specific equations for the time and measurement updates are presented below.
1 Time Update

Unscented Particle Filter
The unscented particle filter UPF is a hybrid of the UKF and the particle filter which uses the UKF to get better importance sampling density.A pseudo-code description of UPF is as follows 11-13 . 1 Initialization: Time 0. Generate N samples x i 0 , i 1, 2, . . ., N from the prior p x 0 , and set the importance weight w i 0 of each sample 1/N: 2 Time k.
I a Update the particles with the UKF: i calculate sigma points from { x i k−1 , P i k−1 } using 2.12 , ii propagate particle into future by 2.13 , iii incorporate new observation to update the measurement by 2.14 and obtain b Sample a new particle x i k and make k and normalize the importance weights w i k : where p z k | x i k is likelihood probability distribution, which is given by measurement model The basic idea of resampling is to eliminate particles with small weights and to concentrate on particles with large weights.Multiply/suppress particles { x i k , P i k } with high/low importance weights w i k , respectively, to obtain N random particles The overall state estimation and covariance are

2.17
Mathematical Problems in Engineering 7

State Model
The state model dynamical model of the celestial OD system for a near-Earth satellite based on the orbital dynamics in the Earth-Centered Inertial ECI frame J2000.0 is 3.1 Equation 3.1 can be written in a general state equation as where X x y z v x v y v z T is the state vector.x, y, z, v x , v y , v z are satellite positions and velocities of the three axes, respectively, μ is the gravitational constant of earth, J 2 is the second zonal coefficient and has the value 0.0010826269 25 , and R e is the earth's radius.ΔF x , ΔF y , ΔF z are the perturbations including high order nonspherical earth perturbations, third-body perturbations, atmospheric drag perturbations, solar radiation perturbations, and other perturbations, which are considered as process noises w t .

Celestial Orbit Determination and Its Measurement
The celestial OD method is based on the fact that the position of a celestial body in the inertial frame at a certain time is known and that its position measured in the spacecraft body frame is a function of the satellite's position.To earth satellite, stars are distributed all over the sky, and the positions of Earth are fixed at a certain time.The geometric relationship among stars, the Earth, and satellite enables us to determine the position of the satellite 26 .Satellite celestial OD methods can be broadly separated into two major approaches: directly sensing horizon method and indirectly sensing horizon method.In this paper, the directly sensing horizon method is used.
The angle between a star and the earth, α, as shown in Figure 2, is a kind of directly sensing horizon measurement of satellite celestial OD system, which is measured by star sensor and earth sensor.The measurement model using the star-earth angle is given by 27 where r is the position vector of the satellite, which is the same as that in 3.2 , s is the position vector of the star in the earth-centered inertial frame, ν α is the measurement noise.Assuming a measurement Z 1 α and measurement noise V 1 v α , 3.2 can be written as a general measurement equation 3.4

Geomagnetic Orbit Determination and Its Measurement
Geomagnetic OD system relies on measurements from a three-axis magnetometer to determine satellite position and orbit.It uses a model of Earth's magnetic field and a model of orbital dynamics to predict the time-varying magnitude of Earth's magnetic field vector at the space.OD system compares the time history of the predicted magnitude and the measured magnitude time history in filter sense to obtain the optimal estimated state position and velocity 14 .

Magnetic Model
Two main models used for describing Earth's magnetic vector in the geodetic reference frame are World Magnetic Model WMM and International Geomagnetic Reference Field IGRF 28 .The WMM 2005 is selected in this paper for geomagnetic orbit determination 29 .
According to the WMM model 2005, the vector field B can be written as the gradient of a potential function B r, λ, θ, t −∇V r, λ, θ, t , 3.5 where r, λ, θ represent the radius, the longitude, and the colatitude in a spherical, geocentric reference frame, respectively.This potential V can be expanded in terms of spherical harmonics: V r, λ, θ, t where N 12 is the degree of the expansion of the WMM, a is the standard Earth's magnetic reference radius, g m n t and h m n t are the time-dependent Gauss coefficients of degree n and order m, and P m n cos θ are the Schmidt normalized associated Legendre polynomials.

Magnetic Measurement Model
Based on the relationship between magnetic vector, which is obtained by the magnetometer, and the earth magnetic model, the measurement model can be written as where B s is the magnetic vector of local position in sensor coordinates, which can be obtained from vector magnetometer system consisting of three mutually orthogonal, singleaxis magnetometers.B t B n B e B v T is the magnetic vector of local position in geocentric coordinates, and it can be obtained from WMM according to local longitude, latitude, and height, as shown in Figure 3; A sb , A bi , and A it are the transformation matrices from satellite body coordinates to sensor coordinates, from earth inertial coordinates to satellite body coordinates, and from earth inertial coordinates to geocentric coordinates, respectively.v B is the measurement noise.Assuming a measurement Z 2 B s and measurement noise V 2 v B , 3.7 can be written as a general measurement equation as 3.8

Simulation Condition
The trajectory used in the following simulation is a LEO satellite whose orbital parameters are semimajor axis a 7136.and attitude data of the satellite are produced by the Satellite Tool Kit STK software 30 .
The accuracy of star sensor and earth sensor is selected 3 and 0.02 • , respectively.The stellar database used in simulation is the Tycho stellar catalog 31 .The magnetometer measurement and geomagnetic model accuracy is considered as 100 nT 32 .

Performances under Different Sampling Intervals
Figures 4 and 5 show the performances comparison among the EKF, UKF, and UPF methods of celestial OD system and geomagnetic OD system, respectively.Data is obtained with a 3 s sampling interval during the 600 min period 6 orbits .Tables 1 and 2 present the details of the simulation results of celestial OD system and geomagnetic OD system under different sampling intervals, respectively.The simulations in Figures 4 and 5 suggest that the EKF-based OD system performance is the worst.In contrast, UPF-based OD system provides the highest OD accuracy.As the details in Tables 1 and 2, regardless the celestial OD and geomagnetic OD system, the different sampling intervals can strongly affect the OD accuracy.OD performance is degraded remarkably with increasing sample interval.However, under the same sampling interval, the EKF method is the most sensitive to the sampling interval, for the nonlinear error increases rapidly with the longer sampling interval.In contrast, the UKF and UPF perform distinctly better.

Performance under Different Noise Distributions
This subsection reports how different noise distributions affect the OD performances using three filters.We selected three common noise distributions in navigation, and they are normal distribution, student's t distribution, and uniform distribution 33 .Figures 6 and 7 show the OD results of celestial OD system and geomagnetic OD system using three filters under student's t noise distributions, respectively.All performance curves were obtained with 15 s sampling interval during the 600 min period 6 orbits .Tables 3 and 4 present the details of the simulation results of celestial OD system and geomagnetic OD system under three different noise distributions, respectively.
As the results in Figures 6 and 7 showed, the UPF-based geomagnetic OD system provides the highest OD accuracy.As the details in Tables 3 and 4 demonstrated, OD performance under different noise distribution is similar.In general, the EKF performance

Computation Cost of Three Methods
Besides the accuracy, the computation cost is another essential requirement to evaluate the performance of filtering methods.Table 5 gives the computation cost of the three methods for the celestial orbit determination system and the geomagnetic orbit determination system, respectively.As in the theoretical value of computation cost, where Φ is the process Jacobian, n is the order of the Φ, and in the simulation n equals 6.The simulation results presented here were run on a 2.66 GHz Inter Core2 Duo CPU with 32-bit Windows 7 system.The simulation time of celestial OD system in Table 5 demonstrates that the UPF demands the highest computation time, which is almost twenty times sample number higher than UKF, and EKF requires almost a quarter of the computation time of UKF.However, the simulation time of geomagnetic OD system is not the same amount as celestial OD system, and the EKFbased geomagnetic OD system takes significantly longer time, since the time for computing measurement Jacobians takes a lot of computer resource.

Conclusion
The problem of choosing a suitable filtering method for the orbit determination application has been studied here.Three filtering methods for the autonomous orbit determination using either celestial or geomagnetic measurements have been studied and their performances have been compared for the estimation problem.
The algorithms are tested with STK satellite orbit data, and the simulation results demonstrate that UPF yields the best OD accuracy and the EKF yields the worst under the same condition.The main reason is that the state equations and measurement equations for autonomous orbit determination system are significantly nonlinear as well as the non-Gaussian errors.
In addition, the paper analyzed the computation cost of the three filtering methods, and UPF-based OD system can provide the highest OD accuracy, though it requires the largest computation time.However, the UPF can finally meet the real-time requirements, as with the development of computer technology.

Figure 1 :
Figure 1: The process of orbit determination.

Figure 2 :
Figure 2: The measurement of celestial OD system.

Figure 3 :
Figure 3: The measurement of geomagnetic OD system.

Figure 4 :Figure 5 :
Figure 4: Three filtering methods results of celestial OD system T 3 s .

Figure 6 :Figure 7 :
Figure 6: Celestial OD results of three filtering methods under Student's t noise distributions.

Table 1 :
Performance of celestial OD system under different sampling intervals.

Table 2 :
Performance of geomagnetic OD system under different sampling intervals.

Table 3 :
Performance of celestial OD system under different noise distributions.

Table 4 :
Performance of geomagnetic OD system under different noise distributions.

Table 5 :
Comparison of computation cost.