# Onboard and Real-Time Artificial Satellite Orbit Determination Using GPS

Authors: Ana Paula Marins Chiaradia; Hélio Koiti Kuga; Antonio Fernando Bertachini de Almeida Prado
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

An algorithm for real-time and onboard orbit determination applying the Extended Kalman Filter (EKF) method is developed. Aiming at a very simple and still fairly accurate orbit determination, an analysis is performed to ascertain an adequacy of modeling complexity versus accuracy. The minimum set of to-be-estimated states to reach the level of accuracy of tens of meters is found to have at least the position, velocity, and user clock offset components. The dynamical model is assessed through several tests, covering force model, numerical integration scheme and step size, and simplified variational equations. The measurement model includes only relevant effects to the order of meters. The EKF method is chosen to be the simplest real-time estimation algorithm with adequate tuning of its parameters. In the developed procedure, the obtained position and velocity errors along a day vary from 15 to 20 m and from 0.014 to 0.018 m/s, respectively, with standard deviation from 6 to 10 m and from 0.006 to 0.008 m/s, respectively, with the SA either on or off. The results, as well as analysis of the final adopted models used, are presented in this work.


---

## Body

## 1. Introduction

The Global Positioning System (GPS) is a satellite navigation system that has been used to determine the position, velocity, and time with high accuracy of an artificial satellite, such as TOPEX/POSEIDON (T/P) [1], Jason-1, Jason-2, CHAMP, and GRACE satellites. Current samples of satellites having an onboard GPS receiver are Jason [2–4], GRACE [5, 6], CHAMP [7, 8], and KOMPSAT [9, 10].The TOPEX/POSEIDON (T/P) mission is jointly conducted by the United States National Aeronautics and Space Administration (NASA) and the French Space Agency, Centre National d’Etudes Spatiales (CNES). The main goal of this mission is to improve the knowledge of the global ocean circulation. Other applications include the ocean tides, geodesy and geodynamics, ocean wave height, and wind speed [1]. The T/P spacecraft orbits the Earth at an altitude of 1336 km, inclination of 66°, and near-zero eccentricity. The period of the orbit is 1.87 hr. The satellite orbit must be determined with an RMS radial accuracy of 13 cm. This is an extremely stringent accuracy requirement for a satellite of this shape and altitude [11]. The T/P GPS receiver can track up to 6 GPS satellites at once in both frequencies if antispoofing is inactive.Its successor, Jason-1, continues the time series of centimeter-level ocean topography observations as it follows on to the highly successful TOPEX/Poseidon (T/P) radar altimeter satellite. Jason has the same orbital parameters of T/P satellite. Preliminary tests of orbits computed using Jason-1 tracking data in a reduced dynamic strategy suggest that the RMS radial accuracies are already better than 2 cm, and that the 1 cm goal is within reach [2].The GRACE (Gravity Recovery and Climate Experiment) satellites, launched in March 2002 with about 470 km in altitude, are each equipped with a BlackJack GPS onboard receiver for precise orbit determination and gravity field recovery [5]. With primary objective of determining the Earth’s gravity field and its temporal variations with unprecedented accuracy, accurate orbits for GRACE are required. The orbit determination has been obtained with accuracy better than 5 cm in each direction for GRACE orbits. The relative accuracy of the two GRACE satellites is about 1 cm in position and 10 micrometers per second in velocity.The CHAllenging Minisatellite Payload (CHAMP) is a German small satellite used for geophysical research and application. This mission is scheduled to last five years in order to provide a sufficiently long observation time to resolve long-term temporal variations primarily in the magnetic field, in the gravity field, and within the atmosphere. CHAMP was launched on 15 July 2000 and reentered the Earth’s atmosphere on September 20, 2010 [12]. The position accuracy is about 0.8 m for pseudorange in single frequency [13].Two types of orbit determination are performed at the KOMPSAT using GPS data. One is the operational orbit determination using GPS navigation solutions providing the position accuracy around 4.45 m RMS and the velocity accuracy around 0.005 m/s RMS. The other is the precise orbit determination using GPS raw measurement data by DGPS technique. In this case, the difference between NASA GSFC Precise Orbit Ephemeris and precise orbit determination results using TEC are radial 6.9 cm RMS, along track 19.4 cm RMS, and cross track 6.9 cm RMS [8].Gomes et al. [14] analyzed a real-time orbit estimator using the raw navigation solution provided by GPS receivers. The estimation algorithm considers a Kalman filter with a rather simple orbit dynamic model and random walk modeling of the receiver clock bias and drift. Using the Topex/Poseidon satellite as test bed, characteristics of model truncation, sampling rates, and degradation of the GPS receiver (Selective Availability) were analyzed. The position precision obtained was around 20–30 m in 3D and around 0.5 m/s in velocity.Pardal et al. [15] determined the orbit of an artificial satellite and analyzed its implications, using least squares algorithms through sequential given rotations as the method of estimation, and data of the GPS receivers. This approach has the goal of improving the performance of the orbit estimation process and, at the same time, at the same time, of minimizing the computational procedure cost. Perturbations due to high order geopotential and direct solar radiation pressure and the position of the GPS antenna on the satellite body were taken into account.In this paper, an algorithm to determine onboard the satellite orbit in real-time using the GPS system and Kalman filtering is developed. It used a simplified and compact model with low computational cost. The extended Kalman filter (EKF) estimates the state vector, composed of the position and velocity components, bias, drift, and drift rate of the GPS receiver clock. A simple fixed step-size-fourth order Runge-Kutta numerical integrator is found to be suitable to integrate the differential equations of the orbital motion. The force model in the equations of motion may consider the perturbations due to the geopotential up to any order and degree of the spherical harmonics. The state error covariance matrix is computed through the transition matrix, which is calculated analytically in an optimized way. The raw single frequency pseudorange GPS measurements are used as observations by the Kalman filter. They are modeled taking into account most of the GPS satellite and receiver clock offsets. The main aim is to develop a GPS-based real-time orbit determination system for any artificial Earth satellite. The target is orbit accuracy in the order of tens of meters with simple modeling of the dynamics or measurement systems, and still keeping low requirements of speed and computational burden.


## 2. Orbit Determination

The orbit determination process consists of obtaining values of the parameters that completely specify the motion of an orbiting body, like artificial satellites, based on a set of observations of the body. It involves nonlinear dynamical and nonlinear measurement systems, which depends on the tracking system and estimation technique (e.g., Kalman filtering or least squares [16–21]). The dynamical system model consists of orbital motion dynamics, measurements models, Earth’s rotation effects, and perturbation models. As Montenbruck and Gill [22], besides the state variables that define the initial conditions, these models depend on a variety of parameters that affects either the dynamical motion or the measurement process. Due to the complexity of the applied models, it is hardly possible to solve directly for any of these parameters from a given set of observations.The observation may be obtained from the ground station networks using laser, radar, Doppler, and so forth, or by space navigation systems, as the Global Positioning System (GPS). The choice of the tracking system depends on a compromise between the goals of the mission and the tools available.In the case of the GPS, the advantages are global coverage, high precision, low cost, and autonomous navigation resources. The GPS may provide orbit determination with accuracy at least as good as methods using ground-tracking networks. The latter provides standard precision around hundred meters and the former can provide precision as tight as some centimeters.The GPS provides, at a given instant, a set of many redundant measurements, which makes the orbit position observable geometrically.After some advances of technology, the single frequency GPS receivers provide a good basis to achieve fair precision at relatively low cost, still attaining the accuracy requirements of the mission.The GPS allows the receiver to determine its position and time geometrically anywhere at any instant with data from at least four satellites. The principle of navigation by satellites is based in sending signals and data from the GPS satellites to a receiver that is inside the satellite that needs to have its orbit determined. This receiver measures the travel time of the signal and then calculates the distance between the receiver and the GPS satellite. If the clocks are synchronized, measurements from three GPS satellites are enough to obtain its position. If the clocks are not synchronized, four measurements are required. Those measurements of distances are called pseudorange, as shown in Figure1.Figure 1
Orbit determination using GPS.


## 3. Description of the Algorithm

The algorithm is developed to determine artificial satellite orbits according the dynamical model (the force model).  The method is used to integrate the differential equations of a satellite motion, the computation of the state transition matrix, the measurement model, and the estimation technique.For that sake, the following aspects are verified:(a)
concerning the dynamical model: what is better order and degree truncation of the spherical harmonic coefficients of the JGM-2 model in order to reach the objectives in terms of low computational cost, relative accuracy (tens of meters), and the lowest possible order and degree of the harmonical coefficients;(b)
concerning the integration method: with step size of the fixed step RK4 integrator to use, considering low computational cost, accuracy, and the selection of appropriate reference systems and transformations;(c)
considering the measurement model: where the errors source is included in the model; and(d)
considering the filtering: the EKF is used and the use of the simplified transition matrix considering either the pure Keplerian motion or including at least the J2 effects.


### 3.1. Force Model and Numerical Integration

The main forces acting on an artificial satellite may be classified as gravitational (e.g., geopotential, third-body perturbations, and tides) and nongravitational forces (e.g., atmospheric drag, solar radiation pressure). These forces are rather difficult to model, so a mathematically inaccurate model of the motion with respect to true motion is sometimes adopted in Kuga [23, 24].To choose a simplified force model to be adopted in this work, the following factors are considered: generality, orbit accuracy, and computational cost for orbit determination in a real-time and onboard environment. Therefore, only forces due to the Earth gravitational field should be implemented.The considered harmonic coefficients of the Earth geopotential field are set up to order and degree 10 of JGM-2 model, according to studies developed in Chiaradia et al. [25], without overloading the processing time.The acceleration and the related partial derivatives matrices are computed through the recurrence relations, according to Pines [26], in Earth-fixed (EF) coordinates. The coordinate transformation from True of Date (ToD) to Pseudo-Earth fixed equator and prime meridian (PEF) takes into account the Earth sidereal rotation, but the polar motion is neglected in this work.The integration of the satellite motion equation is carried out using the simple 4th order of the Runge-Kutta (RK4) algorithm. It is implemented without any mechanism of step adjustment or error control. The RK4 is considered an adequate numerical integrator, due to its simplicity, fair accuracy, low truncation error, and low computational burden. An initialization procedure is not necessary and the step size is quite easy to change. The 30-second integration step interval is used.


### 3.2. Measurement Model

With the advances of technology, the single frequency GPS receivers provide a good basis to achieve fair precision at a relatively low cost, still attaining the accuracy requirements of general missions. In this way, only the code pseudorange in L1 frequency is used to determine the satellite orbit, providing the accuracy required by the mission, similar to the one shown in this work.The errors considered in the C/A code pseudorange measurements are the GPS satellite and receiver clock offsets. Chiaradia et al. [25, 27] studied the need of considering the ionospheric correction model. It has been concluded that the ionospheric effect (around 3–6 meters) does not affect significantly the final accuracy of the orbit. In this way, it will be neglected once the goal is not to achieve high precision (around centimeters). In high precision and onboard applications, some ionospheric correction for single frequency should be applied, but this may add some complexity if we need to implement a better ionospheric model.The equation of the C/A code pseudorange in L1 frequency is given by(1)ρc=ρ+c[ΔtGPS(t)-Δtu(t)],
where ρc is the C/A code pseudorange in L1, c is the speed of light, ΔtGPS is the GPS satellite clock offset, Δtu is the receiver clock offset, t is the observation instant in GPS time, and ρ is the geometric range given by
(2)ρ=(xGPS-x)2+(yGPS-y)2+(zGPS-z)2,
where x, y, and z are the positional states of the user satellite at reception time (in GPS time) and xGPS, yGPS, and zGPS are the positional states of the GPS satellite at transmission time (in GPS time), corrected for light time delay.The second term on the right side of (1) is the clock bias that represents the combined clock offsets of the satellite and of the receiver with respect to GPS time. Each GPS satellite contributes with one clock bias. The information for GPS satellite clocks is known and transmitted via a navigation message in the form of three polynomial coefficients with a reference time toc. The clock correction of the GPS satellite for epoch t is
(3)t=tGPS-ΔtGPS,(4)ΔtGPS=af0+af1(t-toc)+af2(t-toc)2+δRGPS,(5)δRGPS=2c2aμesinE,
where tGPS is the spacecraft code phase time at message transmission time in seconds; af0, af1, and af2 are bias, drift, and drift rate of GPS satellite clock, respectively; toc is the reference time, in seconds, measured from the GPS time weekly epoch; δR is a small relativistic correction of the GPS satellite clock caused by orbital eccentricity e; a is the semimajor axis of the orbit; μ is the Earth gravitational constant; and E is the eccentric anomaly.The relativistic effect is included in the clock polynomial broadcast via navigation message. If a more accurate equation is required, the relativistic effects must be subtracted from the clock polynomial coefficients. The polynomial coefficientaf0, af1, and af2 are transmitted in units of sec, nondimensional, and 1/sec, respectively. The clock data reference time toc is also broadcasted. The value of tGPS must account for the beginning or the end-of-week crossovers. The user may approximate t by tGPS in (3).The user clock offset is part of the estimated state vector given by [28]
(6)cΔtu=b0+b1ΔT+b2ΔT2+δRUSER,
where b0, b1, and b2 are the bias, drift, and drift rate of the user clock and ΔT is the elapsed time since the instant of the first measurement. The relativistic effect in the user clock is calculated by using the best available estimated state vector in the epoch.


### 3.3. GPS Travel and Reception Times

The computation of geometric rangeρ is required for both pseudorange and carrier phase processing. The GPS position coordinates are used at the instant of emission. However, the reception time is used to compute it. Then, this time is corrected by subtracting the signal travel time to obtain the emission time. Therefore, the travel time is computed through an iterative process that starts assuming an average value for travel time τ. Next, the GPS position for the epoch (t-τ) is interpolated and then the geometric range is computed, which can be used to reconstruct the travel time by [29]
(7)τ=ρc.
The light-time iteration is usually performed in an inertial system with position vector r=(x,y,z) and the GPS satellite position vector rGPS=(xGPS,yGPS,zGPS). Then, the positions from inertial to the Earth-fixed WGS84 system or vice versa are transformed. So, the signal path is given as [22]
(8)ρ=RWGS(t)(rGPS(t-τ)-r(t)),
where RWGS is the rotation matrix from inertial to Earth-fixed WGS-84 system. Making use of the approximation
(9)RWGS(t)≈Rz(ωeτ)RWGS(t-τ),
where ωe is the Earth rotation rate, the inertial position of the GPS satellite may be substituted by the corresponding Earth-fixed position:
(10)rGPSWGS(t-τ)=RWGS(t-τ)rGPS(t-τ).
This yields in the Earth-fixed reference frame
(11)ρ=cτ=|Rz(ωeτ)rGPS(t-τ)-rWGS(t)|.
If the discrepancy between the first and second approximation of τ is greater than a specified criterion, the iteration is repeated; that is, a new satellite position is interpolated and a new distance is computed and so on. Generally, a couple of iterations are sufficient.This computation is not affected by the receiver clock errorΔtu and the satellite clock error Δtsv. However, the ionosphere troposphere, hardware, and multipath corrupt the computed nominal emission time. All of these effects are neglected in the present context. The same process corrects the reception time as well.


### 3.4. Estimation Technique

The standard Kalman filter algorithm is used to estimate the spacecraft orbit onboard. This filter is robust, with easy implementation for applications in real-time and no requirement of iterating the previously collected data, and is able to provide the current orbit in real-time. The Extended Kalman Filter (EKF) is a version which is applicable to nonlinear problems, and it is composed of time-update and measurement-update cycles. The time-update phase updates the state and the covariance matrix along the time using dynamical equations. In this work, the estimated state vector is given by(12)x=(r,v,b)T,
where r=(x,y,z)T and v=(x˙,y˙,z˙)T are the spacecraft’s position and velocity vectors, b=(b0,b1,b2)T, where b0, b1, and b2 are bias, drift, and drift rate of the receiver clock, respectively, and t is the instant of integration. All coordinates refer to ToD system.The algorithm used can be found in references [18, 19]. The details of how it is implemented here can be found in reference [30].The state error covariance matrix is propagated through the transition matrix. The conventional filter is implemented by calculating only the upper triangular part of the state error covariance matrix. The other elements are obtained imposing symmetry to the matrix. The lack of symmetry is one of the highest sources of truncation error per Kalman filter cycle [31].One method to avoid the problem of high computational cost and extended analytical expressions of the transition matrix consists of propagating the state vector using the complete force model and then, computing the transition matrix using a simplified force model. The analytical calculation of the transition matrix of the Keplerian motion is a reasonable approximation when only short time intervals of the observations and reference instant are involved [31]. On the other hand, the inclusion of J2 (Earth flattening) effect in the transition matrix can approximately be done by adopting the Markley’s method [32].Chiaradia et al. [33] studied these two methods. The first method is an approximation of the Keplerian motion, providing an analytical solution, which is then calculated numerically by solving the Kepler’s equation. The second one is a local numerical approximation that includes the effect of J2. They were evaluated according to accuracy, processing time, and handling complexity of the equations for two kinds of orbits. It has been concluded that the Keplerian motion model, which was developed by Kuga [31], is still a better suited force model to be adopted in the evaluation of a general-purpose state transition matrix. Some other papers studying similar or related topics can be found in references [34–36].


## 3.1. Force Model and Numerical Integration

The main forces acting on an artificial satellite may be classified as gravitational (e.g., geopotential, third-body perturbations, and tides) and nongravitational forces (e.g., atmospheric drag, solar radiation pressure). These forces are rather difficult to model, so a mathematically inaccurate model of the motion with respect to true motion is sometimes adopted in Kuga [23, 24].To choose a simplified force model to be adopted in this work, the following factors are considered: generality, orbit accuracy, and computational cost for orbit determination in a real-time and onboard environment. Therefore, only forces due to the Earth gravitational field should be implemented.The considered harmonic coefficients of the Earth geopotential field are set up to order and degree 10 of JGM-2 model, according to studies developed in Chiaradia et al. [25], without overloading the processing time.The acceleration and the related partial derivatives matrices are computed through the recurrence relations, according to Pines [26], in Earth-fixed (EF) coordinates. The coordinate transformation from True of Date (ToD) to Pseudo-Earth fixed equator and prime meridian (PEF) takes into account the Earth sidereal rotation, but the polar motion is neglected in this work.The integration of the satellite motion equation is carried out using the simple 4th order of the Runge-Kutta (RK4) algorithm. It is implemented without any mechanism of step adjustment or error control. The RK4 is considered an adequate numerical integrator, due to its simplicity, fair accuracy, low truncation error, and low computational burden. An initialization procedure is not necessary and the step size is quite easy to change. The 30-second integration step interval is used.


## 3.2. Measurement Model

With the advances of technology, the single frequency GPS receivers provide a good basis to achieve fair precision at a relatively low cost, still attaining the accuracy requirements of general missions. In this way, only the code pseudorange in L1 frequency is used to determine the satellite orbit, providing the accuracy required by the mission, similar to the one shown in this work.The errors considered in the C/A code pseudorange measurements are the GPS satellite and receiver clock offsets. Chiaradia et al. [25, 27] studied the need of considering the ionospheric correction model. It has been concluded that the ionospheric effect (around 3–6 meters) does not affect significantly the final accuracy of the orbit. In this way, it will be neglected once the goal is not to achieve high precision (around centimeters). In high precision and onboard applications, some ionospheric correction for single frequency should be applied, but this may add some complexity if we need to implement a better ionospheric model.The equation of the C/A code pseudorange in L1 frequency is given by(1)ρc=ρ+c[ΔtGPS(t)-Δtu(t)],
where ρc is the C/A code pseudorange in L1, c is the speed of light, ΔtGPS is the GPS satellite clock offset, Δtu is the receiver clock offset, t is the observation instant in GPS time, and ρ is the geometric range given by
(2)ρ=(xGPS-x)2+(yGPS-y)2+(zGPS-z)2,
where x, y, and z are the positional states of the user satellite at reception time (in GPS time) and xGPS, yGPS, and zGPS are the positional states of the GPS satellite at transmission time (in GPS time), corrected for light time delay.The second term on the right side of (1) is the clock bias that represents the combined clock offsets of the satellite and of the receiver with respect to GPS time. Each GPS satellite contributes with one clock bias. The information for GPS satellite clocks is known and transmitted via a navigation message in the form of three polynomial coefficients with a reference time toc. The clock correction of the GPS satellite for epoch t is
(3)t=tGPS-ΔtGPS,(4)ΔtGPS=af0+af1(t-toc)+af2(t-toc)2+δRGPS,(5)δRGPS=2c2aμesinE,
where tGPS is the spacecraft code phase time at message transmission time in seconds; af0, af1, and af2 are bias, drift, and drift rate of GPS satellite clock, respectively; toc is the reference time, in seconds, measured from the GPS time weekly epoch; δR is a small relativistic correction of the GPS satellite clock caused by orbital eccentricity e; a is the semimajor axis of the orbit; μ is the Earth gravitational constant; and E is the eccentric anomaly.The relativistic effect is included in the clock polynomial broadcast via navigation message. If a more accurate equation is required, the relativistic effects must be subtracted from the clock polynomial coefficients. The polynomial coefficientaf0, af1, and af2 are transmitted in units of sec, nondimensional, and 1/sec, respectively. The clock data reference time toc is also broadcasted. The value of tGPS must account for the beginning or the end-of-week crossovers. The user may approximate t by tGPS in (3).The user clock offset is part of the estimated state vector given by [28]
(6)cΔtu=b0+b1ΔT+b2ΔT2+δRUSER,
where b0, b1, and b2 are the bias, drift, and drift rate of the user clock and ΔT is the elapsed time since the instant of the first measurement. The relativistic effect in the user clock is calculated by using the best available estimated state vector in the epoch.


## 3.3. GPS Travel and Reception Times

The computation of geometric rangeρ is required for both pseudorange and carrier phase processing. The GPS position coordinates are used at the instant of emission. However, the reception time is used to compute it. Then, this time is corrected by subtracting the signal travel time to obtain the emission time. Therefore, the travel time is computed through an iterative process that starts assuming an average value for travel time τ. Next, the GPS position for the epoch (t-τ) is interpolated and then the geometric range is computed, which can be used to reconstruct the travel time by [29]
(7)τ=ρc.
The light-time iteration is usually performed in an inertial system with position vector r=(x,y,z) and the GPS satellite position vector rGPS=(xGPS,yGPS,zGPS). Then, the positions from inertial to the Earth-fixed WGS84 system or vice versa are transformed. So, the signal path is given as [22]
(8)ρ=RWGS(t)(rGPS(t-τ)-r(t)),
where RWGS is the rotation matrix from inertial to Earth-fixed WGS-84 system. Making use of the approximation
(9)RWGS(t)≈Rz(ωeτ)RWGS(t-τ),
where ωe is the Earth rotation rate, the inertial position of the GPS satellite may be substituted by the corresponding Earth-fixed position:
(10)rGPSWGS(t-τ)=RWGS(t-τ)rGPS(t-τ).
This yields in the Earth-fixed reference frame
(11)ρ=cτ=|Rz(ωeτ)rGPS(t-τ)-rWGS(t)|.
If the discrepancy between the first and second approximation of τ is greater than a specified criterion, the iteration is repeated; that is, a new satellite position is interpolated and a new distance is computed and so on. Generally, a couple of iterations are sufficient.This computation is not affected by the receiver clock errorΔtu and the satellite clock error Δtsv. However, the ionosphere troposphere, hardware, and multipath corrupt the computed nominal emission time. All of these effects are neglected in the present context. The same process corrects the reception time as well.


## 3.4. Estimation Technique

The standard Kalman filter algorithm is used to estimate the spacecraft orbit onboard. This filter is robust, with easy implementation for applications in real-time and no requirement of iterating the previously collected data, and is able to provide the current orbit in real-time. The Extended Kalman Filter (EKF) is a version which is applicable to nonlinear problems, and it is composed of time-update and measurement-update cycles. The time-update phase updates the state and the covariance matrix along the time using dynamical equations. In this work, the estimated state vector is given by(12)x=(r,v,b)T,
where r=(x,y,z)T and v=(x˙,y˙,z˙)T are the spacecraft’s position and velocity vectors, b=(b0,b1,b2)T, where b0, b1, and b2 are bias, drift, and drift rate of the receiver clock, respectively, and t is the instant of integration. All coordinates refer to ToD system.The algorithm used can be found in references [18, 19]. The details of how it is implemented here can be found in reference [30].The state error covariance matrix is propagated through the transition matrix. The conventional filter is implemented by calculating only the upper triangular part of the state error covariance matrix. The other elements are obtained imposing symmetry to the matrix. The lack of symmetry is one of the highest sources of truncation error per Kalman filter cycle [31].One method to avoid the problem of high computational cost and extended analytical expressions of the transition matrix consists of propagating the state vector using the complete force model and then, computing the transition matrix using a simplified force model. The analytical calculation of the transition matrix of the Keplerian motion is a reasonable approximation when only short time intervals of the observations and reference instant are involved [31]. On the other hand, the inclusion of J2 (Earth flattening) effect in the transition matrix can approximately be done by adopting the Markley’s method [32].Chiaradia et al. [33] studied these two methods. The first method is an approximation of the Keplerian motion, providing an analytical solution, which is then calculated numerically by solving the Kepler’s equation. The second one is a local numerical approximation that includes the effect of J2. They were evaluated according to accuracy, processing time, and handling complexity of the equations for two kinds of orbits. It has been concluded that the Keplerian motion model, which was developed by Kuga [31], is still a better suited force model to be adopted in the evaluation of a general-purpose state transition matrix. Some other papers studying similar or related topics can be found in references [34–36].


## 4. Simulations

To validate and to analyze the algorithm, the Topex/Poseidon (T/P) satellite data is used. This satellite carries a dual frequency receiver GPS onboard to experimentally test the ability of the GPS to provide Precise Orbit Determination (POD). The following files are used:(i)
the T/P observation files that broadcast the code and carrier pseudorange measurements in two frequencies in 10-second GPS time steps and are provided by the GPS Data Processing Facility of the Jet Propulsion Laboratory (JPL) in RINEX format;(ii)
the T/P Precise Orbit Ephemeris (POE) files that are generated by JPL; and(iii)
the broadcast GPS navigation message file in RINEX format provided by Crystal Dynamics Data Information System (CDDIS) of the Goddard Space Flight Center (GSFC).The estimated position and velocity are compared against the T/P-POE. It is claimed to estimate the T/P position with accuracy, which is better than 15 cm. The states in the POE are provided in one minute UTC time steps in ToD system. The T/P GPS measurements are provided in 10 seconds of GPS time. According to IERS, the difference between the UTC and GPS times is approximately an integer number of seconds, increasing timely with the introduction of leap seconds. For example, the difference is 9 seconds for 1993.However, in the process of orbit determination, the state is estimated in 30-second intervals in UTC time with transmission and reception time corrections. Then, these data instants are not coincident. Therefore, the states are interpolated through an interpolation (Polint) subroutine [37]. With this approach, the mean error of the interpolated states is 0.068 m and 2.5 × 10-4m/s for position and velocity, respectively, which does not add any significant bias to the accuracy evaluation of the results.Some parameters to be used in evaluating the algorithm in the tests are defined. The real position error is given by(13)Δr≡[∑i=13(xi-x^i)2]1/2,
where xi and x^i, i=1,2,3, are the reference components (or real) and estimated position vectors, respectively. The estimated position error is given by
(14)Δr^=[∑i=13Pii]1/2,
where Pii, i=1,2,3, represents the values of the diagonal elements of the state covariance matrix of the estimated position. The real velocity error is given by
(15)Δv≡[∑i=46(xi-x^i)2]1/2,
where xi and x^i,  i=4,5,6, are the components of the reference (or real) and estimated velocity vectors, respectively. The estimated velocity error is given by
(16)Δv^=[∑i=46Pii]1/2,
where Pii, i=4,5,6, represents the values of the diagonal elements of the state covariance matrix of the estimated velocity. And, the residual is given by:
(17)Δρ=z-ρc,
where z and ρc are the observed and calculated pseudorange measurements, respectively.


## 5. Results

The algorithm is analyzed and tested using 8 complete days of data. Figures2 and 3 show the real and filtered position and velocity errors in m and m/s, with SA (Select Availability) off and on, respectively. These figures show the results for one day of filtering.Figure 2
Estimated and real error with SA off.Figure 3
Estimated and real error with SA on.Figure4 shows the typical behavior of the pseudorange residuals for all tests.Figure 4
Typical behavior of the residuals.Figures5 and 6 show the histogram and normal distribution of the residuals.Figure 5
Histogram of the residuals.Figure 6
Normal distribution of the residuals.Table1 shows the statistical error for position and velocity vectors and the residual for all tests.Table 1
Statistical errors with SA on or off.
Date
Δ
r (m)
Δ
v (m)
Residuals (m)
11/18/1993
15.5
±
6.8
0.014
±
0.006
0.027
±
13.2
11/19/1993
17.4
±
6.7
0.016
±
0.006
-
0.130
±
13.4
11/20/1993
17.6
±
8.4
0.017
±
0.007
-
0.160
±
22.6
01/03/1994
16.5
±
8.5
0.015
±
0.008
0.004
±
22.1
01/04/1994
19.5
±
9.6
0.018
±
0.009
-
0.160
±
22.8
01/05/1994
19.6
±
7.8
0.018
±
0.008
0.040
±
22.9
01/21/1994
19.0
±
7.7
0.018
±
0.008
-
0.050
±
23.9
01/22/1994
16.5
±
8.6
0.015
±
0.008
-
0.007
±
23.3It can be noted that, for all days tested, the real position and velocity errors are fewer than with the estimated position and velocity errors by the filter. It shows the conservative behavior of the filter. In all cases, the filtering takes around one hour to converge. Before achieving the convergence, the onboard computer can use the GPS navigation solution provided by GPS receiver with 100-meter position error. The position accuracy with the SA off (11/18/1993 and 11/19/1993) or on (other days) is from 15 to 20 m, with standard deviation from 6 to 10 m. The velocity accuracy varies from 0.014 to 0.018 m/s, with standard deviation going from 0.006 to 0.009 m/s, as shown in Table1.Through Figures5 and 6 and Table 1, it can be noted that the residuals present normal distribution with mean zero and standard deviation of around 23 m. On November 18 and 19, 1993, the SA is off and the standard deviation of the residuals is around 13.3 m. This fact shows how the SA affects the statistics of the residuals, although the estimated state had not been significantly degraded.


## 6. Conclusions

The main goal of this work is to achieve accuracy around tens of meters along with minimum computational cost when determining, in real-time, artificial satellite orbits onboard considering a simplified model. For its development, single frequency GPS measurements are used, considering the effects of the clock offsets of the GPS and user satellites and user relativistic effects.This algorithm is tested with real data from a satellite with a GPS receiver onboard. The satellite orbit is estimated using the developed algorithm with a good accuracy and minimum computational cost. The algorithm uses a large 30-second step size of propagation (10 seconds step size can be used as well), the geopotential model up to order and degree 10, and the computation method of the transition matrix considering the Keplerian motion. The position accuracy obtained is better than 20 m with the SA on or off. The results are according to expected.



---

*Source: 530516-2013-02-27.xml*