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

Authors: Xiaolin Ning; Xin Ma; Cong Peng; Wei Quan; Jiancheng Fang
Journal: Mathematical Problems in Engineering
 (2012)

---

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


---

## Body

## 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 Figure1, 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 Section2. 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:(2.1)xk=f(xk-1,k-1)+wk,zk=h(xk,k)+vk,
where xk-1 denotes the state of the system at time k-1, zk denotes the observations at step k, wk denotes the process noise, and vk denotes the measurement noise. The mappings f and h represent the process and measurement models. E(wkwjT)=Qk, E(vkvjT)=Rk, for all k, j, and Qk is the process noise covariance at step k, Rk is the measurement noise covariance at step k.


### 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,(2.2)f(xk)=f(x̂k∣k-1)+∂f(xk)∂xk|xk=x̂k∣k-1(xk-x̂k∣k-1)+⋯.
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 UpdatePredicted state estimate:(2.3)x̂k∣k-1=f(x̂k-1,k-1).Predicted estimate covariance:(2.4)Pk-=ΦkPk-1ΦkT+Qk-1.The time update equations project the state,xk, and covariance, Pk, 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.5)Φk=∂f∂xk|xk=x̂k∣k-1.(2)
Measurement UpdateNear-Optimal Kalman gain:(2.6)Kk=Pk-1HkT(HkPk-1HkT+Rk-1)-1.Updated state estimate:(2.7)x̂k=x̂k∣k-1+Kk(zk-h(x̂k∣k-1,k)).Updated estimate covariance:(2.8)Pk=(I-KkHk)Pk-,
where Kk is known as the Kalman gain. The measurement update equations correct the state and covariance estimates with the measurement zk. Hk is the observation matrix at step k, which is defined to be the following Jacobians:(2.9)Hk=∂h∂xk|xk=x̂k∣k-1.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 n×1, the 2n+1 sigma point and their weight are computed by [9](2.10)X0,k=x̂k,W0=τ(n+τ),Xi,k=x̂k+n+τ(P(k∣k))i,Wi=1[2(n+τ)],Xi+n,k=x̂k-n+τ(P(k∣k))i,Wi+n=1[2(n+τ)],i=1,2,…,n,
where τ∈R, (P(k∣k))iis the ith column of the matrix square root. The UKF process can be described as follows.(1)
Time = 0, initialize the UKF withx̂0 and P0 as follows:(2.11)x̂0=E[x0],P0=E[(x0-x̂0)(x0-x̂0)T].(2)
Time =k, define 2n+1 sigma points from(2.12)Xk-1=[X0,kXi,kXi+n,k],i=1,2,…,n.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(2.13)Xk∣k-1=f(Xk-1,k-1),x̂k-=∑i=02nWiXi,k∣k-1,Pk-=∑i=02nWi[Xi,k∣k-1-x̂k-]⋅[Xi,k∣k-1-x̂k-]T+Qk,Zk∣k-1=h(Xk∣k-1,k),ẑk-=∑i=02nWiZi,k∣k-1.(2)
Measurement Update(2.14)Pẑkẑk=∑i=02nWi[Zi,k∣k-1-ẑk-][Zi,k∣k-1-ẑk-]T+Rk,Px̂kẑk=∑i=02nWi[Xi,k∣k-1-x̂k-][Zi,k∣k-1-ẑk-]T,Kk=Px̂kẑkPẑkẑk-1,x̂k=x̂k-+Kk(zk-ẑk-),Pk=Pk--KkPẑkẑkKkT.


### 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.GenerateN samples x0i, (i=1,2,…,N) from the prior p(x0), and set the importance weight w0i of each sample 1/N:(2.15)x̂0i=E[x0i],P0i=E[(x0i-x̂0i)(x0i-x̂0i)T],w0i=1N.(2) Time =k.(I)
(a) Update the particles with the UKF:
(i)
calculate sigma points from{x̂k-1i,Pk-1i} using (2.12),(ii)
propagate particle into future by (2.13),(iii)
incorporate new observation to update the measurement by (2.14) and obtain {x̂ki,Pki}.(b) Sample a new particlex̂ki and makex̂ki~q(xki∣xk-1i,zk)=N(x̂ki,Pki).(II)
Compute the importance weightw̅ki and normalize the importance weights wki:(2.16)w̅ki=wk-1i⋅p(zk∣x̂ki)p(x̂ki∣xk-1i)q(x̂ki∣xk-1i,zk-1),wki=w̅ki∑i=1Nw̅ki,
where p(zk∣x̂ki) is likelihood probability distribution, which is given by measurement model zk=h(xk,k)+vk, p(x̂ki∣xk-1i) is the forward transition probability distribution, which is given by process model xk=f(xk-1,k-1)+wk, q(x̂ki∣xk-1i,zk-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{x̂ki,Pki} with high/low importance weights wki, respectively, to obtain N random particles {x̃ki,P̃ki}.(IV)
Output step:The overall state estimation and covariance are(2.17)x̂k=∑i=1Nwkix̃ki,Pk=∑i=1NwkiP̃ki=∑i=1Nwki(x̃ki-x̂k)(x̃ki-x̂k)T.


## 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,(2.2)f(xk)=f(x̂k∣k-1)+∂f(xk)∂xk|xk=x̂k∣k-1(xk-x̂k∣k-1)+⋯.
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 UpdatePredicted state estimate:(2.3)x̂k∣k-1=f(x̂k-1,k-1).Predicted estimate covariance:(2.4)Pk-=ΦkPk-1ΦkT+Qk-1.The time update equations project the state,xk, and covariance, Pk, 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.5)Φk=∂f∂xk|xk=x̂k∣k-1.(2)
Measurement UpdateNear-Optimal Kalman gain:(2.6)Kk=Pk-1HkT(HkPk-1HkT+Rk-1)-1.Updated state estimate:(2.7)x̂k=x̂k∣k-1+Kk(zk-h(x̂k∣k-1,k)).Updated estimate covariance:(2.8)Pk=(I-KkHk)Pk-,
where Kk is known as the Kalman gain. The measurement update equations correct the state and covariance estimates with the measurement zk. Hk is the observation matrix at step k, which is defined to be the following Jacobians:(2.9)Hk=∂h∂xk|xk=x̂k∣k-1.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 n×1, the 2n+1 sigma point and their weight are computed by [9](2.10)X0,k=x̂k,W0=τ(n+τ),Xi,k=x̂k+n+τ(P(k∣k))i,Wi=1[2(n+τ)],Xi+n,k=x̂k-n+τ(P(k∣k))i,Wi+n=1[2(n+τ)],i=1,2,…,n,
where τ∈R, (P(k∣k))iis the ith column of the matrix square root. The UKF process can be described as follows.(1)
Time = 0, initialize the UKF withx̂0 and P0 as follows:(2.11)x̂0=E[x0],P0=E[(x0-x̂0)(x0-x̂0)T].(2)
Time =k, define 2n+1 sigma points from(2.12)Xk-1=[X0,kXi,kXi+n,k],i=1,2,…,n.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(2.13)Xk∣k-1=f(Xk-1,k-1),x̂k-=∑i=02nWiXi,k∣k-1,Pk-=∑i=02nWi[Xi,k∣k-1-x̂k-]⋅[Xi,k∣k-1-x̂k-]T+Qk,Zk∣k-1=h(Xk∣k-1,k),ẑk-=∑i=02nWiZi,k∣k-1.(2)
Measurement Update(2.14)Pẑkẑk=∑i=02nWi[Zi,k∣k-1-ẑk-][Zi,k∣k-1-ẑk-]T+Rk,Px̂kẑk=∑i=02nWi[Xi,k∣k-1-x̂k-][Zi,k∣k-1-ẑk-]T,Kk=Px̂kẑkPẑkẑk-1,x̂k=x̂k-+Kk(zk-ẑk-),Pk=Pk--KkPẑkẑkKkT.


## 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.GenerateN samples x0i, (i=1,2,…,N) from the prior p(x0), and set the importance weight w0i of each sample 1/N:(2.15)x̂0i=E[x0i],P0i=E[(x0i-x̂0i)(x0i-x̂0i)T],w0i=1N.(2) Time =k.(I)
(a) Update the particles with the UKF:
(i)
calculate sigma points from{x̂k-1i,Pk-1i} using (2.12),(ii)
propagate particle into future by (2.13),(iii)
incorporate new observation to update the measurement by (2.14) and obtain {x̂ki,Pki}.(b) Sample a new particlex̂ki and makex̂ki~q(xki∣xk-1i,zk)=N(x̂ki,Pki).(II)
Compute the importance weightw̅ki and normalize the importance weights wki:(2.16)w̅ki=wk-1i⋅p(zk∣x̂ki)p(x̂ki∣xk-1i)q(x̂ki∣xk-1i,zk-1),wki=w̅ki∑i=1Nw̅ki,
where p(zk∣x̂ki) is likelihood probability distribution, which is given by measurement model zk=h(xk,k)+vk, p(x̂ki∣xk-1i) is the forward transition probability distribution, which is given by process model xk=f(xk-1,k-1)+wk, q(x̂ki∣xk-1i,zk-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{x̂ki,Pki} with high/low importance weights wki, respectively, to obtain N random particles {x̃ki,P̃ki}.(IV)
Output step:The overall state estimation and covariance are(2.17)x̂k=∑i=1Nwkix̃ki,Pk=∑i=1NwkiP̃ki=∑i=1Nwki(x̃ki-x̂k)(x̃ki-x̂k)T.


## 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(3.1)dxdt=vx,dydt=vy,dzdt=vz,dvxdt=-μxr3⋅[1-J2(Rer)2(7.5z2r2-1.5)]+ΔFx,dvydt=-μyr3⋅[1-J2(Rer)2(7.5z2r2-1.5)]+ΔFy,dvzdt=-μzr3⋅[1-J2(Rer)2(7.5z2r2-4.5)]+ΔFz,r=x2+y2+z2.Equation (3.1) can be written in a general state equation as(3.2)Ẋ(t)=f(X(t),t)+w(t),
where X=[xyzvxvyvz]T is the state vector. x, y,z, vx, vy, vz are satellite positions and velocities of the three axes, respectively, μ is the gravitational constant of earth, J2 is the second zonal coefficient and has the value 0.0010826269 [25], and Re is the earth’s radius. ΔFx,ΔFy,ΔFz 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).


### 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](3.3)α=arccos(-s⋅rr)+να,
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.Figure 2
The measurement of celestial OD system.Assuming a measurementZ1=[α] and measurement noise V1=[vα], (3.2) can be written as a general measurement equation(3.4)Z1(t)=h1[X(t),t]+V1(t).


### 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 fieldB can be written as the gradient of a potential function(3.5)B(r,λ,θ,t)=-∇V(r,λ,θ,t),
where (r,λ,θ) represent the radius, the longitude, and the colatitude in a spherical, geocentric reference frame, respectively.This potentialV can be expanded in terms of spherical harmonics:(3.6)V(r,λ,θ,t)=∑n=1N∑m=0nVnm=a∑n=1N(ar)n+1∑m=0n[gnm(t)cos⁡(mλ)+hnm(t)sin⁡(mλ)]P̆nm(cos⁡θ),
where N=12 is the degree of the expansion of the WMM, a is the standard Earth’s magnetic reference radius, gnm(t) and hnm(t) are the time-dependent Gauss coefficients of degree n and order m, and P̆nm(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(3.7)Bs=AsbAbiAitBt+vB,
where Bs 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. Bt=[BnBeBv]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; Asb, Abi, and Ait 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. vB is the measurement noise.Figure 3
The measurement of geomagnetic OD system.Assuming a measurementZ2=Bs and measurement noise V2=vB, (3.7) can be written as a general measurement equation as(3.8)Z2(t)=h2[X(t),t]+V2(t).


## 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(3.1)dxdt=vx,dydt=vy,dzdt=vz,dvxdt=-μxr3⋅[1-J2(Rer)2(7.5z2r2-1.5)]+ΔFx,dvydt=-μyr3⋅[1-J2(Rer)2(7.5z2r2-1.5)]+ΔFy,dvzdt=-μzr3⋅[1-J2(Rer)2(7.5z2r2-4.5)]+ΔFz,r=x2+y2+z2.Equation (3.1) can be written in a general state equation as(3.2)Ẋ(t)=f(X(t),t)+w(t),
where X=[xyzvxvyvz]T is the state vector. x, y,z, vx, vy, vz are satellite positions and velocities of the three axes, respectively, μ is the gravitational constant of earth, J2 is the second zonal coefficient and has the value 0.0010826269 [25], and Re is the earth’s radius. ΔFx,ΔFy,ΔFz 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).


## 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](3.3)α=arccos(-s⋅rr)+να,
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.Figure 2
The measurement of celestial OD system.Assuming a measurementZ1=[α] and measurement noise V1=[vα], (3.2) can be written as a general measurement equation(3.4)Z1(t)=h1[X(t),t]+V1(t).


## 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 fieldB can be written as the gradient of a potential function(3.5)B(r,λ,θ,t)=-∇V(r,λ,θ,t),
where (r,λ,θ) represent the radius, the longitude, and the colatitude in a spherical, geocentric reference frame, respectively.This potentialV can be expanded in terms of spherical harmonics:(3.6)V(r,λ,θ,t)=∑n=1N∑m=0nVnm=a∑n=1N(ar)n+1∑m=0n[gnm(t)cos⁡(mλ)+hnm(t)sin⁡(mλ)]P̆nm(cos⁡θ),
where N=12 is the degree of the expansion of the WMM, a is the standard Earth’s magnetic reference radius, gnm(t) and hnm(t) are the time-dependent Gauss coefficients of degree n and order m, and P̆nm(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(3.7)Bs=AsbAbiAitBt+vB,
where Bs 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. Bt=[BnBeBv]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; Asb, Abi, and Ait 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. vB is the measurement noise.Figure 3
The measurement of geomagnetic OD system.Assuming a measurementZ2=Bs and measurement noise V2=vB, (3.7) can be written as a general measurement equation as(3.8)Z2(t)=h2[X(t),t]+V2(t).


## 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 fieldB can be written as the gradient of a potential function(3.5)B(r,λ,θ,t)=-∇V(r,λ,θ,t),
where (r,λ,θ) represent the radius, the longitude, and the colatitude in a spherical, geocentric reference frame, respectively.This potentialV can be expanded in terms of spherical harmonics:(3.6)V(r,λ,θ,t)=∑n=1N∑m=0nVnm=a∑n=1N(ar)n+1∑m=0n[gnm(t)cos⁡(mλ)+hnm(t)sin⁡(mλ)]P̆nm(cos⁡θ),
where N=12 is the degree of the expansion of the WMM, a is the standard Earth’s magnetic reference radius, gnm(t) and hnm(t) are the time-dependent Gauss coefficients of degree n and order m, and P̆nm(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(3.7)Bs=AsbAbiAitBt+vB,
where Bs 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. Bt=[BnBeBv]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; Asb, Abi, and Ait 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. vB is the measurement noise.Figure 3
The measurement of geomagnetic OD system.Assuming a measurementZ2=Bs and measurement noise V2=vB, (3.7) can be written as a general measurement equation as(3.8)Z2(t)=h2[X(t),t]+V2(t).


## 4. Analysis and Comparison

### 4.1. Simulation Condition

The trajectory used in the following simulation is a LEO satellite whose orbital parameters are semimajor axisa=7136.635444km, eccentricity e=1.809×10-3, inclination i=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

Figures4 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.
Sampling intervalRMS (after convergence)Maximum (after convergence)Position error/mVelocity error/m/sPosition error/mVelocity error/m/sT=3sEKF203.3182110.196622532.9872720.486705UKF161.3127230.162900357.5456000.407503UPF159.7560790.160780354.5553590.402378T=15sEKF271.6409530.287641703.0005810.613708UKF245.9393020.219864593.8290980.563749UPF245.2296830.219566593.8020480.563396T=60sEKF934.2389390.9766412457.8951702.348656UKF736.8762880.6999422322.2918962.302674UPF735.1669320.6988082314.7612492.294924Table 2
Performance of geomagnetic OD system under different sampling intervals.
Sampling intervalRMS (after convergence)Maximum (after convergence)Position error/mVelocity error/m/sPosition error/mVelocity error/m/sT=3sEKF591.2536280.6019461129.3655101.061735UKF376.8943720.371566877.9094030.799659UPF376.5168630.366538861.5139750.783449T=15sEKF1481.6737521.3588672851.6601772.607626UKF705.7654500.6482281325.5012671.276198UPF705.1612630.6478761323.9761071.274585T=60sEKF4343.7831624.30892114643.27574111.712547UKF3904.8905443.63379810403.52854110.328902UPF3904.7478923.63346010401.27913910.328554Three filtering methods results of celestial OD system(T=3s).
(a)(b)Three filtering methods results of geomagnetic OD system(T=3s).
(a)(b)The simulations in Figures4 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’st distribution, and uniform distribution [33].Figures6 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.Table 3
Performance of celestial OD system under different noise distributions.
Noise distributionRMS (after convergence)Maximum (after convergence)Position error/mVelocity error/m/sPosition error/mVelocity error/m/sNormaldistributionEKF271.6409530.287641703.0005810.613708UKF245.9393020.219864593.8290980.563749UPF245.2296830.219566593.8020480.563396Student’st distributionEKF387.6180440.411669827.3583000.972190UKF257.1043600.261611604.3252650.664447UPF246.2941970.250511587.1542290.639735UniformdistributionEKF385.8036090.398842994.7142640.902353UKF350.0795160.339869888.7906770.788940UPF349.5373250.338473880.4172060.782641Table 4
Performance of geomagnetic OD system under different noise distributions.
Noise distributionRMS (after convergence)Maximum (after convergence)Position  error/mVelocity error/m/sPosition error/mVelocity error/m/sNormal distributionEKF1481.6737521.3588672851.6601772.607626UKF705.7654500.6482281325.5012671.276198UPF705.1612630.6478761323.9761071.274585Student’st distributionEKF1059.8748490.8288952900.8001783.184135UKF670.3216610.6048351594.0606921.436411UPF669.9590040.6045131593.0858091.436468Uniform distributionEKF1246.5291301.2077284029.1877743.253167UKF776.6224860.7675621985.8705771.781823UPF776.5261440.7674521985.7046991.781594Celestial OD results of three filtering methods under Student’st noise distributions.
(a)(b)Geomagnetic OD results of three filtering methods under Student’st noise distributions.
(a)(b)As the results in Figures6 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. Table5 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 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.
FiltermethodTheoretical value of computation costSimulation value of computation costΦ is full matrixΦ is diagonal matrixComputation time per orbit of celestial OD system(s)Computation time per orbit of geomagnetic OD system (s)EKFCEKF=3n3+3n2+4nCEKF′=6n2+4n0.39.5UKFCUKF=2n3+12n2+14n+5CUKF′=2n3+12n2+14n+51.611.7UPF(sample num = 20)CUPF=samplenum·CUKFCUPF′=samplenum·CUKF′39.5313.9


## 4.1. Simulation Condition

The trajectory used in the following simulation is a LEO satellite whose orbital parameters are semimajor axisa=7136.635444km, eccentricity e=1.809×10-3, inclination i=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

Figures4 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.
Sampling intervalRMS (after convergence)Maximum (after convergence)Position error/mVelocity error/m/sPosition error/mVelocity error/m/sT=3sEKF203.3182110.196622532.9872720.486705UKF161.3127230.162900357.5456000.407503UPF159.7560790.160780354.5553590.402378T=15sEKF271.6409530.287641703.0005810.613708UKF245.9393020.219864593.8290980.563749UPF245.2296830.219566593.8020480.563396T=60sEKF934.2389390.9766412457.8951702.348656UKF736.8762880.6999422322.2918962.302674UPF735.1669320.6988082314.7612492.294924Table 2
Performance of geomagnetic OD system under different sampling intervals.
Sampling intervalRMS (after convergence)Maximum (after convergence)Position error/mVelocity error/m/sPosition error/mVelocity error/m/sT=3sEKF591.2536280.6019461129.3655101.061735UKF376.8943720.371566877.9094030.799659UPF376.5168630.366538861.5139750.783449T=15sEKF1481.6737521.3588672851.6601772.607626UKF705.7654500.6482281325.5012671.276198UPF705.1612630.6478761323.9761071.274585T=60sEKF4343.7831624.30892114643.27574111.712547UKF3904.8905443.63379810403.52854110.328902UPF3904.7478923.63346010401.27913910.328554Three filtering methods results of celestial OD system(T=3s).
(a)(b)Three filtering methods results of geomagnetic OD system(T=3s).
(a)(b)The simulations in Figures4 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’st distribution, and uniform distribution [33].Figures6 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.Table 3
Performance of celestial OD system under different noise distributions.
Noise distributionRMS (after convergence)Maximum (after convergence)Position error/mVelocity error/m/sPosition error/mVelocity error/m/sNormaldistributionEKF271.6409530.287641703.0005810.613708UKF245.9393020.219864593.8290980.563749UPF245.2296830.219566593.8020480.563396Student’st distributionEKF387.6180440.411669827.3583000.972190UKF257.1043600.261611604.3252650.664447UPF246.2941970.250511587.1542290.639735UniformdistributionEKF385.8036090.398842994.7142640.902353UKF350.0795160.339869888.7906770.788940UPF349.5373250.338473880.4172060.782641Table 4
Performance of geomagnetic OD system under different noise distributions.
Noise distributionRMS (after convergence)Maximum (after convergence)Position  error/mVelocity error/m/sPosition error/mVelocity error/m/sNormal distributionEKF1481.6737521.3588672851.6601772.607626UKF705.7654500.6482281325.5012671.276198UPF705.1612630.6478761323.9761071.274585Student’st distributionEKF1059.8748490.8288952900.8001783.184135UKF670.3216610.6048351594.0606921.436411UPF669.9590040.6045131593.0858091.436468Uniform distributionEKF1246.5291301.2077284029.1877743.253167UKF776.6224860.7675621985.8705771.781823UPF776.5261440.7674521985.7046991.781594Celestial OD results of three filtering methods under Student’st noise distributions.
(a)(b)Geomagnetic OD results of three filtering methods under Student’st noise distributions.
(a)(b)As the results in Figures6 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. Table5 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 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.
FiltermethodTheoretical value of computation costSimulation value of computation costΦ is full matrixΦ is diagonal matrixComputation time per orbit of celestial OD system(s)Computation time per orbit of geomagnetic OD system (s)EKFCEKF=3n3+3n2+4nCEKF′=6n2+4n0.39.5UKFCUKF=2n3+12n2+14n+5CUKF′=2n3+12n2+14n+51.611.7UPF(sample num = 20)CUPF=samplenum·CUKFCUPF′=samplenum·CUKF′39.5313.9


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



---

*Source: 267875-2012-02-06.xml*