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

Abstract

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.

1. 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.

Figure 1 

The process of orbit determination.

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.

2. 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:𝐱_𝑘=𝑓𝐱_𝑘−1,𝑘−1+𝐰_𝑘,𝐳_𝑘=ℎ𝐱_𝑘,𝑘+𝐯_𝑘, (2.1) where 𝐱_𝑘−1 denotes the state of the system at time 𝑘−1, 𝐳_𝑘 denotes the observations at step 𝑘, 𝐰_𝑘 denotes the process noise, and 𝐯_𝑘 denotes the measurement noise. The mappings 𝑓 and ℎ represent the process and measurement models. 𝐄(𝐰_𝑘𝐰^𝑇_𝑗)=𝐐_𝑘, 𝐄(𝐯_𝑘𝐯^𝑇_𝑗)=𝐑_𝑘, for all 𝑘, 𝑗, and 𝐐_𝑘 is the process noise covariance at step 𝑘, 𝐑_𝑘 is the measurement noise covariance at step 𝑘.

2.1. 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,𝑓𝐱_𝑘=𝑓𝐱_{𝑘∣𝑘−1}+𝜕𝑓𝐱_𝑘𝜕𝐱_𝑘|_𝐱𝑘=𝐱_{𝑘∣𝑘−1}𝐱_𝑘−𝐱_{𝑘∣𝑘−1}+⋯. (2.2) 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:𝐱_{𝑘∣𝑘−1}=𝑓𝐱_𝑘−1,𝑘−1. (2.3)

Predicted estimate covariance:𝐏⁻_𝑘=𝚽_𝑘𝐏_𝑘−1𝚽^𝑇_𝑘+𝐐_𝑘−1. (2.4)

The time update equations project the state, 𝐱_𝑘, and covariance, 𝐏_𝑘, estimates from the previous time step 𝑘−1 to the current time step 𝑘, Φ_𝑘 is the state transition matrix at step 𝑘, which is defined to be the following Jacobians:𝚽_𝑘=𝜕𝑓𝜕𝐱_𝑘|_𝐱𝑘=𝐱_{𝑘∣𝑘−1}. (2.5)

(2)Measurement Update

Near-Optimal Kalman gain:𝐊_𝑘=𝐏_𝑘⁻¹𝐇^𝑇_𝑘𝐇_𝑘𝐏_𝑘⁻¹𝐇^𝑇_𝑘+𝐑_𝑘⁻¹⁻¹. (2.6)

Updated state estimate:𝐱_𝑘=𝐱_{𝑘∣𝑘−1}+𝐊_𝑘𝐳_𝑘−ℎ𝐱_{𝑘∣𝑘−1},𝑘. (2.7)

Updated estimate covariance:𝐏_𝑘=𝐈−𝐊_𝑘𝐇_𝑘𝐏⁻_𝑘, (2.8) where 𝐊_𝑘 is known as the Kalman gain. The measurement update equations correct the state and covariance estimates with the measurement 𝐳_𝑘. 𝐇_𝑘 is the observation matrix at step 𝑘, which is defined to be the following Jacobians:𝐇_𝑘=𝜕ℎ𝜕𝐱_𝑘|_𝐱𝑘=𝐱_{𝑘∣𝑘−1}. (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].

2.2. 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 𝑛×1, the 2𝑛+1 sigma point and their weight are computed by [9]𝐗_0,𝑘=𝐱_𝑘,𝑊₀=𝜏(),𝐗𝑛+𝜏_𝑖,𝑘=𝐱_𝑘+√𝑛+𝜏𝐏(𝑘∣𝑘)_𝑖,𝑊_𝑖=1[2(,𝐗𝑛+𝜏)]_{𝑖+𝑛,𝑘}=𝐱_𝑘−√𝑛+𝜏𝐏(𝑘∣𝑘)_𝑖,𝑊_𝑖+𝑛=1[2(,𝑛+𝜏)]𝑖=1,2,…,𝑛, (2.10) where 𝜏∈𝐑, (𝐏(𝑘∣𝑘))_𝑖is the 𝑖th column of the matrix square root. The UKF process can be described as follows.(1)Time = 0, initialize the UKF with 𝐱₀ and 𝐏₀ as follows:𝐱₀=𝐸𝐱₀,𝐏₀=𝐸𝐱₀−𝐱₀𝐱₀−𝐱₀^𝑇. (2.11)(2)Time = 𝑘, define 2𝑛+1 sigma points from𝐗_𝑘−1=𝐗_0,𝑘𝐗_𝑖,𝑘𝐗_{𝑖+𝑛,𝑘},𝑖=1,2,…,𝑛. (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𝐗_{𝑘∣𝑘−1}=𝑓𝐗_𝑘−1,𝐱,𝑘−1_−𝑘=^2𝑛_𝑖=0𝑊_𝑖𝐗_{𝑖,𝑘∣𝑘−1},𝐏⁻_𝑘=^2𝑛_𝑖=0𝑊_𝑖𝐗_{𝑖,𝑘∣𝑘−1}−𝐱⁻_𝑘⋅𝐗_{𝑖,𝑘∣𝑘−1}−𝐱⁻_𝑘^𝑇+𝐐_𝑘,𝐙_{𝑘∣𝑘−1}=ℎ𝐗_{𝑘∣𝑘−1},𝐳̂,𝑘_−𝑘=^2𝑛_𝑖=0𝑊_𝑖𝐙_{𝑖,𝑘∣𝑘−1}. (2.13)(2)Measurement Update𝐏𝑧_𝑘𝑧_𝑘=^2𝑛_𝑖=0𝑊_𝑖𝐙_{𝑖,𝑘∣𝑘−1}−𝐳̂⁻_𝑘𝐙_{𝑖,𝑘∣𝑘−1}−𝐳̂⁻_𝑘^𝑇+𝐑_𝑘,𝐏𝑥_𝑘𝑧_𝑘=^2𝑛_𝑖=0𝑊_𝑖𝐗_{𝑖,𝑘∣𝑘−1}−𝐱⁻_𝑘𝐙_{𝑖,𝑘∣𝑘−1}−𝐳̂⁻_𝑘^𝑇,𝐊_𝑘=𝐏𝑥_𝑘𝑧_𝑘𝐏𝑧⁻¹_𝑘𝑧_𝑘,𝐱_𝑘=𝐱⁻_𝑘+𝐊_𝑘𝐳_𝑘−𝐳̂⁻_𝑘,𝐏_𝑘=𝐏⁻_𝑘−𝐊_𝑘𝐏𝑧_𝑘𝑧_𝑘𝐊^𝑇_𝑘. (2.14)

2.3. 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 𝑁 samples 𝐱^𝑖₀, (𝑖=1,2,…,𝑁) from the prior 𝑝(𝐱₀), and set the importance weight 𝑤^𝑖₀ of each sample 1/𝑁:𝐱^𝑖₀=𝐸𝐱^𝑖₀,𝐏^𝑖₀=𝐸𝐱^𝑖₀−𝐱^𝑖₀𝐱^𝑖₀−𝐱^𝑖₀^𝑇,𝑤^𝑖₀=1𝑁. (2.15)

(2) Time = 𝑘. (I)(a) Update the particles with the UKF:(i)calculate sigma points from {𝐱^𝑖_𝑘−1,𝐏^𝑖_𝑘−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 𝐱^𝑖_𝑘 and make𝐱^𝑖_𝑘∼𝑞(𝐱^𝑖_𝑘∣𝐱^𝑖_𝑘−1,𝐳_𝑘)=𝑁(𝐱^𝑖_𝑘,𝐏^𝑖_𝑘).(II) Compute the importance weight 𝑤^𝑖_𝑘 and normalize the importance weights 𝑤^𝑖_𝑘:𝑤^𝑖_𝑘=𝑤^𝑖_𝑘−1⋅𝑝𝐳_𝑘∣𝐱^𝑖_𝑘𝑝𝐱^𝑖_𝑘∣𝐱^𝑖_𝑘−1𝑞𝐱^𝑖_𝑘∣𝐱^𝑖_𝑘−1,𝐳_𝑘−1,𝑤^𝑖_𝑘=𝑤^𝑖_𝑘∑^𝑁_𝑖=1𝑤^𝑖_𝑘, (2.16) where 𝑝(𝐳_𝑘∣𝐱^𝑖_𝑘) is likelihood probability distribution, which is given by measurement model 𝐳_𝑘=ℎ(𝐱_𝑘,𝑘)+𝐯_𝑘, 𝑝(𝐱^𝑖_𝑘∣𝐱^𝑖_𝑘−1) is the forward transition probability distribution, which is given by process model 𝐱_𝑘=𝑓(𝐱_𝑘−1,𝑘−1)+𝐰_𝑘, 𝑞(𝐱^𝑖_𝑘∣𝐱^𝑖_𝑘−1,𝐳_𝑘−1)is the proposal distribution [12].(III) Resampling step:

The basic idea of resampling is to eliminate particles with small weights and to concentrate on particles with large weights. Multiply/suppress particles {𝐱^𝑖_𝑘,𝐏^𝑖_𝑘} with high/low importance weights 𝑤^𝑖_𝑘, respectively, to obtain 𝑁 random particles {𝐱^𝑖_𝑘,𝐏^𝑖_𝑘}.(IV) Output step:

The overall state estimation and covariance are𝐱_𝑘=^𝑁_𝑖=1𝑤^𝑖_𝑘𝐱^𝑖_𝑘,𝐏_𝑘=^𝑁_𝑖=1𝑤^𝑖_𝑘𝐏^𝑖_𝑘=^𝑁_𝑖=1𝑤^𝑖_𝑘𝐱^𝑖_𝑘−𝐱_𝑘𝐱^𝑖_𝑘−𝐱_𝑘^𝑇. (2.17)

3. System Models

3.1. 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𝑑𝑥𝑑𝑡=𝑣_𝑥,𝑑𝑦𝑑𝑡=𝑣_𝑦,𝑑𝑧𝑑𝑡=𝑣_𝑧,𝑑𝑣_𝑥𝑑𝑡=−𝜇𝑥𝑟³⋅1−𝐽₂𝑅_𝑒𝑟²7.5𝑧²𝑟²−1.5+Δ𝐹_𝑥,𝑑𝑣_𝑦𝑑𝑡=−𝜇𝑦𝑟³⋅1−𝐽₂𝑅_𝑒𝑟²7.5𝑧²𝑟²−1.5+Δ𝐹_𝑦,𝑑𝑣_𝑧𝑑𝑡=−𝜇𝑧𝑟³⋅1−𝐽₂𝑅_𝑒𝑟²7.5𝑧²𝑟²−4.5+Δ𝐹_𝑧,𝑟=𝑥²+𝑦²+𝑧². (3.1)

Equation (3.1) can be written in a general state equation as𝐗̇(𝑡)=𝑓(𝐗(𝑡),𝑡)+𝐰(𝑡), (3.2) where 𝐗=[𝑥𝑦𝑧𝑣_𝑥𝑣_𝑦𝑣_𝑧]^𝑇 is the state vector. 𝑥, 𝑦, 𝑧, 𝑣_𝑥, 𝑣_𝑦, 𝑣_𝑧 are satellite positions and velocities of the three axes, respectively, 𝜇 is the gravitational constant of earth, 𝐽₂ is the second zonal coefficient and has the value 0.0010826269 [25], and 𝑅_𝑒 is the earth’s radius. Δ𝐹_𝑥,Δ𝐹_𝑦,Δ𝐹_𝑧 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 𝐰(𝑡).

3.2. 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]𝛼=arccos−𝐬⋅𝐫𝑟+𝜈_𝛼, (3.3) where 𝐫 is the position vector of the satellite, which is the same as that in (3.2), 𝐬 is the position vector of the star in the earth-centered inertial frame, 𝜈_𝛼 is the measurement noise.

Figure 2 

The measurement of celestial OD system.

Assuming a measurement 𝑍₁=[𝛼] and measurement noise 𝑉₁=[𝑣_𝛼], (3.2) can be written as a general measurement equation𝐙₁(𝑡)=ℎ₁[𝐗(𝑡),𝑡]+𝐕₁(𝑡). (3.4)

3.3. 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].

3.3.1. 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 𝐁 can be written as the gradient of a potential function𝐁(𝑟,𝜆,𝜃,𝑡)=−∇𝑉(𝑟,𝜆,𝜃,𝑡), (3.5) where (𝑟,𝜆,𝜃) represent the radius, the longitude, and the colatitude in a spherical, geocentric reference frame, respectively.

This potential 𝑉 can be expanded in terms of spherical harmonics:𝑉(𝑟,𝜆,𝜃,𝑡)=^𝑁𝑛_𝑛=1𝑚=0𝑉^𝑚_𝑛=𝑎^𝑁_𝑛=1𝑎𝑟^𝑛𝑛+1_𝑚=0𝑔^𝑚_𝑛(𝑡)cos(𝑚𝜆)+ℎ^𝑚_𝑛(𝑡)sin(𝑚𝜆)𝑃̆^𝑚_𝑛(cos𝜃), (3.6) where 𝑁=12 is the degree of the expansion of the WMM, 𝑎 is the standard Earth’s magnetic reference radius, 𝑔^𝑚_𝑛(𝑡) and ℎ^𝑚_𝑛(𝑡) are the time-dependent Gauss coefficients of degree 𝑛 and order 𝑚, and 𝑃̆^𝑚_𝑛(cos𝜃) are the Schmidt normalized associated Legendre polynomials.

3.3.2. 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𝐁_𝑠=𝐀_sb𝐀_bi𝐀_it𝐁_𝑡+𝐯_𝐵, (3.7) where 𝐁_𝑠 is the magnetic vector of local position in sensor coordinates, which can be obtained from vector magnetometer system consisting of three mutually orthogonal, single-axis magnetometers. 𝐁_𝑡=[𝐵_𝑛𝐵_𝑒𝐵_𝑣]^𝑇 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; 𝐀_sb, 𝐀_bi, and 𝐀_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. 𝐯_𝐵 is the measurement noise.

Figure 3 

The measurement of geomagnetic OD system.

Assuming a measurement 𝐙₂=𝐁_𝑠 and measurement noise 𝐕₂=𝐯_𝐵, (3.7) can be written as a general measurement equation as𝐙₂(𝑡)=ℎ₂[𝐗(𝑡),𝑡]+𝐕₂(𝑡). (3.8)

4. Analysis and Comparison

4.1. Simulation Condition

The trajectory used in the following simulation is a LEO satellite whose orbital parameters are semimajor axis 𝑎=7136.635444km, eccentricity 𝑒=1.809×10⁻³, inclination 𝑖=65^∘, right ascension of the ascending node Ω=30^∘, and the argument of perigee 𝜔=30^∘. The orbit 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].

4.2. 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.

Table 1 

Performance of celestial OD system under different sampling intervals.

Table 2 

Performance of geomagnetic OD system under different sampling intervals.

(a)

(b)

Figure 4 

Three filtering methods results of celestial OD system (𝑇=3s).

(a)

(b)

Figure 5 

Three filtering methods results of geomagnetic OD system (𝑇=3s).

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.

4.3. 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 𝑡 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 𝑡 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.

Table 3 

Performance of celestial OD system under different noise distributions.

Table 4 

Performance of geomagnetic OD system under different noise distributions.

(a)

(b)

Figure 6 

Celestial OD results of three filtering methods under Student’s 𝑡 noise distributions.

(a)

(b)

Figure 7 

Geomagnetic OD results of three filtering methods under Student’s 𝑡 noise distributions.

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 is the worst and the UPF performance is the best, no matter what measurement errors are chosen.

4.4. 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, 𝑛 is the order of the Φ, and in the simulation 𝑛 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 EKF-based geomagnetic OD system takes significantly longer time, since the time for computing measurement Jacobians takes a lot of computer resource.

Table 5 

Comparison of computation cost.

5. 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.

Acknowledgments

The work described in this paper was supported by the National Natural Science Foundation of China (60874095) and Hi-Tech Research and Development Program of China. The authors would like to thank all members of Science and Technology on Inertial Laboratory and Fundamental Science on Novel Inertial Instrument & Navigation System Technology Laboratory, for their useful comments regarding this work effort.