Real time orbit determination using GPS navigation solution

TECHNICAL PAPERS

Real time orbit determination using GPS navigation solution

Vivian M. Gomes^I; Helio K. Kuga^II; Ana Paula M. Chiaradia^III

^I [e-mail]

^IIMember, ABCM [e-mail] Nacional de Pesquisas Espaciais - INPE Av. dos Astronautas 1758, Jd. Da Granja 12227-010 São José dos Campos, SP. Brazil

^[e-mail] São Paulo State University – UNESP Engineering Faculty of Guaratingueta Av. Dr. Ariberto Pereira da Cunha, 333 12516-410 Guaratingueta, SP. Brazil

ABSTRACT

This work analyses 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 analysed.

Keywords: GPS, orbit determination, Kalman filter, navigation solution

Introduction

The Global Positioning System (GPS) is a satellite navigation system for determining position, velocity and time with high accuracy. The GPS system allows a GPS receiver to determine its position and time at any place using data from at least four GPS satellites (Parkinson and Spilker, 1996). Using such a system, this work proposes to compute, in real time, a state vector composed of position, velocity, GPS receiver clock bias and drift of a satellite equipped with an on-board GPS receiver, by filtering the raw navigation solutions provided by the receiver. The Kalman filter is used to estimate the state vector based on the incoming observations from the receiver. The filter dynamic model includes perturbations due to geopotential and the bias is modeled as a random walk process. The observations include the raw navigation solution, composed by the position and time bias. Several tests are done using three days of observations of a GPS receiver on board TOPEX/POSEIDON (T/P) satellite, which are processed by the proposed algorithm. A comparison is done between the estimated state vector and the precise orbit ephemeris (POE) produced by JPL/NASA. Other characteristics are also analysed, including effects of a truncated dynamic model, sampling rate and Selective Availability (SA) effect, to depict the procedure performance in terms of accuracy and computational burden. Current samples of satellites having an on-board GPS receiver are Kompsat (Lee et al., 2004) and BIRD (Briess, 2000). Kompsat does not use the navigation solution for real time orbit determination, but for post processing orbit determination on ground. Unlikely BIRD, it has a receiver providing the navigation solution which is processed in real time on-board (Gill and Montenbruck, 2002; Montenbruck, 2000), however by an epoch state filter in conjunction with NORAD two-lines model (Hoots and Roehrich, 1980). With a growing trend of on-board GPS receivers flying on upcoming satellites this work investigates characteristics of a simple real time algorithm which could influence the performance of the system. Also, the next generation of Brazilian satellites with no exceptions are designed to have platforms carrying GPS receivers, such as the MMP (Multi-Mission Platform), EQUARS (Equatorial Scientific Satellite), and CBERS 2B, 3 and 4 series (China Brazil Earth Resources Satellite).

Methodology

The main purpose of this work is to study a rather simple but still fairly accurate algorithm to determine the artificial satellite orbits, in real time and with low computational burden, by using the raw navigation solutions provided by GPS receivers.

The orbit determination algorithm proposed in this work is implemented by Kalman filtering of the navigation solution from a GPS receiver. The states to be estimated are the position and velocity of the satellite, and bias and drift of the on board GPS receiver clock. The numerical integrator used is the fixed step size Runge-Kutta of fourth order, due to its numerical stability and simplicity of use. The Kalman filter is used due to its suitability to real time applications and recursive and sequential nature, without unnecessary storage of observations, which can be processed during the collecting phase. (Brown and Hwang, 1997).

Dynamical Model

The dynamic model implemented in the Kalman filter is a set of first order stochastic differential equations given by:

wherer is the position vector (x, y, z), v is the velocity vector, a is the acceleration, b is the bias of the GPS receiver clock, d is the drift. w_a is the dynamic noise on the acceleration coordinates and w_d is the dynamic noise on the drift, characterized by the following white noise characteristics:

where Q_ais the spectral power density of the white noise w_a, d(t-t) is the Dirac delta which is equal to 1 when (t=t), and zero otherwise, and Q_d is the spectral power density of the white noise w_d. Thus, the full state vector of the filter is composed by x = (r, v, b, d).

The acceleration vector a includes the geopotential effect which can be truncated to any order and degree. In real time environment it is noticed that the inclusion of high order terms of geopotential, as well as other sources of perturbation may lead to unnecessary computational burden without the corresponding improvement in accuracy. Also, it is known that, for short term predictions (30 minutes), the other sources of perturbation, including Sun and Moon third body effects, solar radiation pressure, and atmospheric drag (for LEO satellites) can be safely neglected (Gill and Montenbruck, 2002).

Measurement Model

The measurement model of the Kalman filter algorithm is given by:

where y = (x, y, z, b) is the observation vector obtained from the navigation solution provided by the GPS receiver at every sampling rate, b is the clock bias and the sensitivity matrix H is therefore given by:

The white noise vector v, which represents the random errors of the measurements, is modeled by:

where d_kjis the Kronecker delta which is equals to 1 when k = j and zero otherwise, and the observation error covariance matrix R is given by:

with s º 30m when the DOP (Dilution of Precision) information is not available. When the DOP information is available, then it is possible to form the matrix R, in a suited way, given by:

where s_r is the standard deviation of the pseudorange measurement, PDOP is the dilution of precision in position and TDOP is the dilution of precision in time. See Parkinson and Spilker (1996) e.g. to DOP definitions.

Trasition Matrix

Another potential source of computer burden is the evaluation of the transition matrix arising from the variational equations of the orbital motion. It was noticed that more complex computations, e.g. including J2 or other effects, does not benefit real time algorithms (Gill and Montenbruck, 2002; Chiaradia et al., 1999). Based on this, the calculation of the transition matrix considered only the Keplerian motion. The Keplerian transition matrix is quickly computable (Kuga, 1986) with the advantage of having analytical solution (Goodyear, 1965; Shepperd, 1985). The full transition matrix is given by:

The submatrices j₁₁, j₁₂, j₂₁ and j₂₂ correspond to the Keplerian motion, j13, j₁₄, j₂₃, j₂₄, j₃₁, j₃₂, j₄₁ and j₄₂ are null and the remaining elements are related to bias and drift of the GPS receiver clock. This partition is given by:

where Dt is the time elapsed between observations.

Simulations Data

For the test cases, it was used data from the TOPEX/POSEIDON satellite, from November 17, 18 and 19, 1993. The results were compared to POE (Precision Orbit Ephemeris) files, which are released by JPL/NASA.

For starting the procedure, it is used the two first navigation solutions provided by the GPS receiver to generate the initial states of the filter (position, velocity, bias and null drift). The standard Kalman filter algorithm was implemented (Brown and Hwang, 1997), where the time update cycle used the transition matrix, as proposed by Eqs. 7-8, to compute the predicted covariance. The measurement update cycle used the conventional Kalman form to process the navigation solution observations in a sequential manner.

Table 1 shows the filter initial state error covariance matrix, where s represents the standard deviation. These values were used throughout the tests. Table 2 shows the values of dynamic noise covariance matrices Q_a and Q_d used in the filter. The measurement errors were set as in Eqs. 5 or 6.

Results and Analyses

Initially, the influence of the dynamic model and sampling rate in the orbit determination accuracy were analyzed. These tests were done using data from T/P satellite, on November 19, 1993.

Next, it was analyzed the influence of SA (Selective Availability) on the T/P orbit determination. To do that, it was used data from November 17, 1993, when the SA was active until 12:00 UTC and inactive after this time.

Finally, it is presented the orbit determination of T/P by using an in-house navigation solution developed at INPE (Instituto Nacional de Pesquisas Espaciais), so as to compare the performance difference between different navigation solution algorithms.

All the statistics and details of the test cases reported may be found in Gomes (2004).

Dynamic Model and Sampling Rate Influence

Tests were done varying the degree and order of the harmonic coefficients, that represent the effect of the terrestrial gravitational potential. It covered since the ideal Keplerian model up to geopotential perturbation of degree and order 50, including the particular cases of J2, 10 x 10 and 50 x 50. Also the different sampling rates of 10 seconds, 1 minute, 5 minutes, 10 minutes were tested to verify the sensitivity of the algorithm concerning the accuracy of orbit determination. In this case, the integration step was varied assuming that the navigation solution was available every 10 seconds, 1 minute, 5 minutes and 10 minutes.

It was noticed that, when using a sampling rate of 10 seconds, there is no significant difference when either a simple model or a more complex one is used. For all models tested (Keplerian, J2, 10x10, 50x50), the filter residuals show mean values nicely fit around zero and similar RMS (e.g. 10.5m for x-position coordinate). The position and velocity errors (compared to POE) show differences smaller than 1 cm and 1 mm/s among the different models. So we can conclude that, with this 10 seconds sampling rate, a simple Keplerian model is enough to deliver results comparable to more complex models.

When using a sampling rate of 1 minute, there is an increase in the order of 30 m on the residuals RMS, however, the residuals still fit around zero for all tested models. Even with increased residuals, there was no meaningful difference on the errors, in relation to sampling of 10 seconds, with variations smaller than 0.5 m and 1 mm/s on the mean position and velocity errors.

For sampling rate of 5 minutes, it is noticed that the position residual RMS increases from meters to hundreds of meters when using the Keplerian model. However, it is possible to notice the improvement when the model considered J2. Using this model, the filter presented residuals mean value around zero and RMS (e.g. 62m for x-position coordinate) much smaller than when the Keplerian model is used (e.g. 470m for x-position coordinate). The further cases 10x10 and 50x50 did not present significant differences. In spite of the residuals increase when using the Keplerian model, there was no impact on the position and velocity error, showing that the filter performs very well, even with such sampling rate and a simple model.

When using an exaggerated sampling rate of 10 minutes with the Keplerian model, the residuals mean and standard deviations increase wildly, even when a more complex model is used. In this case, a big impact on position and velocity errors was observed, with a small improvement when the model with J2 is used.

From this range of tests one concludes that it is advisable to model at least the J2 effect in the satellite orbital motion model, when sampling rates from 10 seconds to 5 minutes are used. For the next tests it was fixed the harmonic coefficients up to 10x10 and sampling rate of 10 seconds, which is the nominal capacity of the GPS receiver on-board T/P.

Selective Availability (SA) Influence

For the first half of the day 17 November 1993, when SA was on, as well as for the second one, when SA was off, the filter residuals mean are around zero. However, the standard deviation is smaller when SA was off. With SA off the magnitude of standard deviation of residuals is around 10 m, and with SA on it is from 2 to 4 times higher.

In Table 3 we can see the filter residual statistics for this case. The 3D error of the estimates are also magnified with SA on. Table 4 shows how the error standard deviation behaves in those two situations. Position RMS errors increased from 30 to 60 m levels and velocity from 0.5 to 1 m/s levels.

INPE´s Navigation Solution

The navigation solutions used in this test were generated by a conventional least squares method developed at INPE. In this case, it is available the full DOP information to build the R matrix stepwise as in Eq (6).

In this case the standard deviations are slightly better than when the navigation solution from the GPS receiver is used. Figure 1 shows a typical residual history for the x- component, and Figure 2 shows the GDOP values along the day. It is noticed that the navigation solutions were generated using 4 to 6 GPS satellites, most of the GDOP values were around 2.5, and there were peaks of GDOP reaching 80 which were not shown in the figure. Figure 3 shows a typical history of the x-component error. It corresponds to the error estimated (lighter lines) by the filter.

We can notice, from Tables 5 and 6, that we achieve better results when using the INPE’s navigation solution. We speculate whether the more adequate weights of the measurement errors (through Eq. 6) is the responsible for this improvement. On the other hand, it also shows that the quality of results on orbit determination using navigation solutions can also depend on the algorithm used to calculate it. In this case, apparently, the algorithm used by INPE is better than the one used in the receiver on board Topex/Poseidon satellite.

Conclusions

The main purpose of this work was to implement and validate a procedure for artificial satellites orbit determination in real time, which is simple and yields a fair accuracy, by using navigation solutions as inputs to a Kalman filtering algorithm.

The cases tested showed that it is not necessary to use complex models to obtain good results when measurements are available at 10 seconds sampling rate. With a higher rate, greater than 5 minutes, it is still possible to get acceptable results using a model considering at least the J2 geopotential harmonic coefficient.

In the tests, it was obtained precision around 20-30 m in 3-D position and around 0.5 m/s in velocity. Such results are consistent with the ones expected, when solving the real time orbit determination problem using GPS.

An analyses of all the design parameters leads us to some guidelines when estimating in real time the orbit of satellites by using the navigation solutions provided by on-board GPS receivers. A sampling rate of 10 s (0.1 Hz) with 10x10 geopotential truncated model and Keplerian model for the transition matrix computation, plus a Runge-Kutta of 4-th order for numerical integration and usage of DOP information to weigh the measurement errors, provides a good balance of cost benefit ratio in terms of speed and computational burden for accuracies in the order of tens of meters.

Paper accepted January, 2007.

Technical Editor: Paulo E. Miyagi.

• Briess, K., 2000, "Satellitenbescheibung: BIRD Spacecraft Description", Vol. 1, DLR, Germany, TN-BIRD-1000-WP/084.
• Brown, R.G. & Hwang, P.Y.C., 1997, "Introduction to random signals and applied Kalman filtering". New York, John Wiley.
• Chiaradia, A.P.; Kuga, H.K. and Prado, A.F.B.A., 1999, "Investigation of simplified models for orbit determination using single frequency GPS measurements", Journal of the Brazilian Society of Mechanical Sciences, Vol. XXI, pp. 165-172.
• Gill, E. & Montenbruck, O., 2002, "On-board navigation system for the BIRD small satellite". DLR, Germany, Forschungsbericht 2002-06.
• Goodyear, W.H., 1965, "Completely general closed form solution for coordinates and partial derivatives of the two-body problem", The Astronomical Journal, Vol. 70, No. 3, pp. 189-192, Apr.
• Gomes, V.M., 2004, "Determinação de órbita de satélites artificiais em tempo real através da solução de navegação do GPS" (In Portuguese), Dissertação de Mestrado Engenharia e Tecnologias Espaciais, INPE, São Jose dos Campos.
• Hoots, F.R. & Roehrich, R.L., 1980, "Models for propagation of NORAD element sets. Aerospace Defense Command", USAF, Project Spacecraft Report 3, Dec.
• Kuga, H.K., 1986, "Transition matrix of the elliptical Keplerian motion" (In Portuguese), Internal report INPE-3779-NTE/250, INPE, São Jose dos Campos.
• Lee, B-S.; Lee, J-S.; Kim, J-H.; Lee, S-P.; Yoon, J-C.; Roh, K-M.; Park, E-S. and Choi, K-H., 2004, "Reconstruction of KOMPSAT-1 GPS navigation solutions using GPS data generation and preprocessing program", Acta Astronautica, Vol. 54, pp. 571-576.
• Montenbruck, O., 2000, "An epoch state filter for use with analytical orbit models of low Earth satellites", Aerospace Science and Technology, Vol. 4, pp. 277-287.
• Parkinson, B.W. & Spilker Jr., J.J., 1996, "Global Positioning System: theory and applications". AIAA, Vol. 1, (Progress in Astronautics and Aeronautics 163).
• Shepperd, S.W, 1985, "Universal Keplerian state transition matrix", Celestial Mechanics, Vol. 35, pp. 129-144.

Publication Dates

• Publication in this collection
  22 Oct 2007
• Date of issue
  Sept 2007

History

• Accepted
  Jan 2007
• Received
  Jan 2007