Table of Contents Author Guidelines Submit a Manuscript
Mathematical Problems in Engineering
Volume 2009 (2009), Article ID 140963, 12 pages
http://dx.doi.org/10.1155/2009/140963
Research Article

A Discussion Related to Orbit Determination Using Nonlinear Sigma Point Kalman Filter

1INPE (National Institute for Space Research), DMC, Avenue dos Astronautas, 1.758. Jd. Granja, CEP 12227-010, São José dos Campos, SP, Brazil
2UNESP (Univ Estadual Paulista), Guaratinguetá, CEP 12516-410, SP, Brazil

Received 30 July 2009; Accepted 9 November 2009

Academic Editor: Tadashi Yokoyama

Copyright © 2009 Paula Cristiane Pinto Mesquita Pardal et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Herein, the purpose is to present a Kalman filter based on the sigma point unscented transformation development, aiming at real-time satellite orbit determination using GPS measurements. First, a brief review of the extended Kalman filter will be done. After, the sigma point Kalman filter will be introduced as well as the basic idea of the unscented transformation, in which this filter is based. Following, the unscented Kalman filter applied to orbit determination will be explained. Such explanation encloses formulations about the orbit determination through GPS; the dynamic model; the observation model; the unmodeled acceleration estimation; also an application of this new filter approaches on orbit determination using GPS measurements discussion.

1. Introduction

In orbit determination of artificial satellites, the dynamic system and the measurements equations are of nonlinear nature. It is a nonlinear problem in which the disturbing forces are not easily modeled. The problem consists of estimating variables that completely specify the body trajectory in the space, processing a set of information (pseudorange measurements) related to this body. A tracking network on Earth or through sensors, like the GPS receiver onboard Topex/Poseidon (T/P) satellite, can collect such observations.

The Global Positioning System (GPS) is a powerful and low cost means to allow computation of orbits for artificial Earth satellites by means of redundant measurements. The T/P is an example of using GPS for space positioning.

The Extended Kalman Filter (EKF) implementation in orbit estimation, under inaccurate initial conditions and scattered measurements, can lead to unstable or diverging solutions. For solving the problem of nonlinear nature, convenient extensions of the Kalman filter have been sought. In particular, the unscented transformation was developed as a method to propagate mean and covariance information through nonlinear transformations. The Sigma Point Kalman Filter (SPKF) appears as an emerging estimation algorithm applied to nonlinear system, without needing linearization steps.

2. Extended Kalman Filter

The real-time estimators are a class of estimators which fulfill the real time requirements. They are recursive algorithms and produce sequentially the state to be estimated. Among them, the Kalman filter and its variations are outstanding.

The Kalman filter is the recursive estimator most used nowadays because it is easy to implement and to use on digital computers. Its recursiveness leads to lesser memory storage, which makes it ideal for real-time applications. The EKF is a nonlinear version of the Kalman filter that generates reference trajectories which are updated at each measurement processing, at the corresponding instant. Details can be found in Brown and Hwang [1].

Due to the complexity of modeling the artificial satellites orbit dynamics accurately, the EKF is generally used in works of such nature. The EKF algorithm always brings up to date reference trajectory around the most current available estimate.

The KF filter consists of phases of time and measurement updates. In the first, state and covariance are propagated from one previous instant to a later one, meaning that they are propagated between discrete instants of the system dynamics model. In the second one, state and covariance are corrected for the later instant corresponding to the measurement time, through the observations model. This method has, therefore, recursive nature and it does not need to store the measurements previously in large matrices.

Exploiting the assumption that all transformations are quasilinear, the EKF simply linearizes all nonlinear transformations and substitutes the Jacobian matrices for the linear transformations in the Kalman filter equations. Although the EKF maintains the elegant and computationally efficient recursive update form of the Kalman Filter, it suffers a number of serious limitations.

The first limitation is that linearized transformations are only reliable if the error propagation can be matched approximated by a linear function. If this condition does not hold, the linearized approximation can be extremely poor. At best, this undermines the performance of the filter. The second is that linearization can be applied only if the Jacobian matrix exists. However, this is not always the case. Some systems contain discontinuities, singularities, and the states themselves are inherently discrete. And the last is that calculating Jacobian matrices can be a very difficult and error-prone process. The Jacobian equations frequently produce many pages of dense algebra that must be converted to code. This introduces numerous opportunities for human coding errors that may degrade the performance of the final system in a manner that cannot be easily identified and debugged. Regardless of whether the obscure code associated with a linearized transformation is or is not correct, it presents a serious problem for subsequent users who must validate it for use in any high integrity system.

Summarizing, the Kalman filter can be applied to nonlinear systems if a consistent set of predicted quantities can be calculated. These quantities are derived by projecting a prior estimate through a nonlinear transformation. Linearization, as applied in the EKF, is widely recognized to be inadequate, but the alternatives incur substantial costs in terms of derivation and computational complexity. Therefore, there is a strong need for a method that is probably more accurate than linearization but does not incur the implementation nor computational costs of other higher-order filtering schemes. The Unscented Transformation was developed to meet these needs.

3. Sigma Point Kalman Filters

If the dynamics system and the observation model are linear, the conventional Kalman filter can be used fearlessly. Although, not rarely, the dynamic systems and the measurement equations are nonlinear, convenient extensions of the Kalman Filter have been sought.

The SPKF is a new estimator that allows similar performance than the Kalman filter for linear systems and it elegantly extends to nonlinear systems, without the linearization steps. These filters are a new approach to generalize the Kalman filter for nonlinear process and observation models. A set of weighted samples, called sigma points, is used for normalizing mean and covariance of a probability distribution. This technique is claimed to lead to a filter more accurate and easier to implement than the EKF or a second-order Gaussian filter.

The SPKF approach is described by Van Merwe et al. [2] as follows.

(1)A set of weighted samples is deterministically calculated, based on mean and covariance decomposition of a random variable. One minimum need is that the first- and second-order momentums are known.(2)The sigma points are propagated through the real nonlinear function, using only functional estimation, that is, analytical derivatives are not used to generate a posteriori set of sigma points.(3)The later statistics are calculated using propagated sigma points functions and weights. In general, they assume the form of a simple weighted average of the mean and the covariance.
3.1. Basic Idea: The Unscented Transformation

A recent method to calculate the statistics of a random variable that passes through a nonlinear transformation is the unscented transform (UT). The UT builds on the principle that it is easier to approximate a probability distribution than it is to approximate an arbitrary nonlinear function, which is detailed in Julier and Uhlmann [3]. The approach is simple: a set of points (sigma points) are chosen so that their mean and covariance are x and 𝐏𝐱𝐱, according to Julier et al. [4, 5], and Julier and Uhlmann [6]. The nonlinear function is applied to each point, in turn, to yield a cloud of transformed points. The statistics of the transformed points, mean 𝐲, and covariance 𝐏𝐲𝐲 predicted, can then be calculated to form an estimate of the nonlinearly transformed mean and covariance.

There are several fundamental differences to the particle filters, although this method bears resemblance it. First, the sigma points are deterministically chosen so that they exhibit certain specific properties (like a given mean and covariance), and are not drawn at random. Second, sigma points can be weighted in ways that are inconsistent with the distribution interpretation of sample points in a particle filter.

The 𝑛-dimensional random variable 𝐱, where n is the state vector dimension, with 𝐱 mean and 𝐏𝐱𝐱 covariance, is approximated by 2𝑛+1 weighted points, given by

𝝌0=𝝌𝐱,𝑖=𝐱+(𝑛+𝜅)𝐏𝐱𝐱𝑖,𝝌𝑖+𝑛=𝐱(𝑛+𝜅)𝐏𝐱𝐱𝑖(3.1) in which 𝜅, ((𝑛+𝜅)𝐏𝐱𝐱)𝑖 is the 𝑖th row or column of the root square matrix of (𝑛+𝜅)𝐏𝐱𝐱, and 𝑊𝑖 is the weight associated to the 𝑖th point

𝑊0=𝜅,𝑊(𝑛+𝜅)𝑖=1𝑊2(𝑛+𝜅),𝑖=1,,𝑛,𝑖+𝑛=12(𝑛+𝜅),𝑖=1,,𝑛.(3.2)

The transformation in the prediction step of the EKF occurs as follows.

(1)Transform each point through the function to yield the set of transformed sigma points 𝑦𝑖𝝌=𝐟𝑖.(3.3)(2)The observations mean is given by the weighted average of the transformed points 𝐲=2𝑛𝑖=0𝑊𝑖𝑦𝑖.(3.4)(3)The covariance is the weighted outer product of the transformed points

𝐏𝑦𝑦=2𝑛𝑖=0𝑊𝑖𝑦𝑖𝑦𝐲𝑖𝐲T.(3.5)

As the algorithm works with a finite number of sigma points, it naturally lends itself to being used in a “black box” filtering library. Given a model with defined inputs and outputs, a standard routine can be used to calculate the predicted quantities as necessary for any given transformation. The computational cost of the algorithm is the same order of magnitude as the EKF. The most expensive operations are to calculate the matrix square root and the outer products which are required to compute the covariance of the projected sigma points.

Any set of sigma points that encodes the mean and covariance correctly calculates the projected mean and covariance correctly to the second order. Therefore, the estimate implicitly includes the second-order “bias correction” term of the truncated second-order filter, but without the need to calculate any derivatives.

The algorithm can be used with discontinuous transformations. Sigma points can pass over a discontinuity and, thus, can approximate the effect of a discontinuity on the transformed estimate.

3.2. The Unscented Kalman Filter

Using UT, the Kalman filter processes, as summarized in the following steps.

The steps shown in Figure 1 are detailed next.

140963.fig.001
Figure 1: Modified EKF, leading to UKF.
(1)To predict the new state system ̂𝐱(𝑘+1𝑘) and its associated covariance 𝐏(𝑘+1𝑘), taking into account the effects of the process white gaussian noise.(2)To predict the expected observation ̂𝐳(𝑘+1𝑘)and its residual covariance (innovation) 𝐏𝜈𝜈(𝑘+1𝑘), considering the effects of the observation noise.(3)To predict the cross correlation matrix 𝐏𝑥𝑧(𝑘+1𝑘).

These steps are put in order in the Kalman filter with the restructuring of dynamics, state vector, and observations models. First, the state vector is added of the noise vector 𝐰𝑘, with dimension 𝑞×1, in order to obtain a vector of dimension 𝑛𝑎=𝑛+𝑞,

𝐱𝑎𝐱(𝑘)=𝑘𝐰𝑘.(3.6)

The dynamics model is rewritten in function of 𝐱𝑎(𝑘) as

𝐱𝐱(𝑘+1)=𝐟𝑎(𝑘),(3.7) and the UT uses 2𝑛𝑎+1 sigma points, generated by

̂𝐱𝑎̂𝟎(𝑘𝑘)=𝐱(𝑘𝑘)𝑞×1,𝐏𝑎(𝑘𝑘)=𝐏(𝑘𝑘)𝐏𝑥𝜈(𝐏𝑘𝑘)𝑥𝜈(𝑘𝑘)𝐐𝑥𝜈(𝑘𝑘).(3.8)

The matrices in the principal diagonal of 𝐏𝑎(𝑘𝑘) are the variances, and the ones out of it are the correlations between the state dynamic errors and the Gaussian process noises.

There are several extensions and modifications that can be done in this basic method to consider specific details for one given application. In the next section, it will be presented a discussion of the orbit determination, in real time, using UKF.

3.3. Comparing EKF and UKF

The conventional nonlinear filters, such as the linearized or the extended Kalman Filter, many times have a poor performance when applied to nonlinear problems, due to two known difficulties.

The linearization (of the dynamic and the measurements models) can lead to a highly instable performance of the filter if the time discretization is not enough small.

The derivation of the Jacobian is not simple in most applications, and usually it makes the implementation difficult.

The UKF has more advantages, when compared to the EKF, as Lee and Alfriend [7, 8] wrote, in the following aspects.

(i)It allows more stable and accurate estimates of mean and covariance.(ii) It can estimate discontinuous functions.(iii)No explicit derivation of the Jacobian and/or Hessian matrix is necessary.(iv)It is suitable for parallel processing.

4. Using UKF on Orbit Determination

Before presenting a discussion of the orbit determination, in real time, using UKF, some points need to be outlined: the orbit determination; the dynamic model; the observations model.

4.1. The Orbit Determination

The orbit determination will be based on GPS technology, whose working principle is based on the geometric method. In such method, the observer knows the set of satellites position in the reference system, obtaining its own position in the same reference frame. Figure 2 presents the basic parameters used by GPS for user position determination.

140963.fig.002
Figure 2: The geometric method.

Here, 𝑅GPSi is the position of the 𝑖th GPS satellite in the reference system; 𝜌𝑖 is the pseudorange, the user satellite position in respect to the 𝑖th GPS satellite; 𝑟𝑢 is the user satellite position in the reference system.

4.2. The Dynamic Model

In the case of orbit determination via GPS, the ordinary differential equations that represent the dynamic model are as follows:

̇̇𝑟=𝑣,𝑣=𝜇𝑟𝑟3̇̇+𝑎+𝑤,𝑏=𝑑,𝑑=0+𝑤𝑑(4.1) with variables given in the inertial reference frame. In the equations above, 𝑟 is the vector containing the position components (𝑥,𝑦,𝑧); 𝑣 is the vector of velocity components; 𝑎 represents the modeled perturbations; 𝑤 is the white noise vector with covariance 𝑄; 𝑏 is the user clock bias; 𝑑 is the user clock drift; 𝑤𝑑 is the white noise on the drift rate with variance 𝑄𝑑.

4.3. The Observations Model

The nonlinear equation of the observations model is given by:

𝐳𝑘=𝐡𝑘𝐱𝑘,𝑡+𝝂𝑘,(4.2) where 𝐳𝑘 is the vector of m observations; 𝐡𝑘(𝐱𝑘) is the nonlinear function of state 𝐱𝑘, with dimension m; 𝝂𝑘 is the vector of observation errors with dimension m.

4.4. Estimation of the Unmodeled Accelerations

Some spacecraft missions require precise orbit knowledge to support payload experiments. Sometimes after launch, ground-based orbit determination solutions do not provide the level of accuracy expected. After verifying all known dynamic models, there may be a residual signature in the orbit as result of unmodeled accelerations. This leads to attempt to estimate anomalous accelerations during the orbit fit, if sufficient data exist. If successful, the acceleration estimates can improve the fit residuals, and also results in better orbital position estimates, as can be seen in Soyka and Davis [9].

Unmodeled accelerations may have many reasons: truncation of geopotential field; limitations of modeling solar pressure, Earth albedo, Earth infrared radiation, drag, and others. Some of these accelerations can be corrected through the use of higher fidelity dynamic and physical modeling, while others require postlaunch calibration.

The use of periodic accelerations, with a period near once per revolution of the satellite orbit, has been used within precision orbit determination programs to improve the accuracy of the derived ephemeris.

4.4.1. Anomalous Accelerations Modeling

When defining an anomalistic or periodic acceleration, one must consider three aspects: the subarc interval, the type of function, and the coordinate frame.

Subarc Interval
The subarc interval is the time of duration or number of revolutions for a given acceleration to be active. As its name implies, it is usually a subset of the total arc. A reason to break an arc into a subarc is to allow for better overall fits.

Type of Function
The underlying mathematical function of an acceleration function is usually a constant, a sine, or a cosine function.
The constant function is the most basic: a constant force in a specific direction. And the periodic functions (sine or cosine) have amplitude, frequency, and phase associated with them. The periodic functions are written as acc=𝐴sin𝜔𝑡+𝜙𝐴oracc=𝐵cos𝜔𝑡+𝜙𝐵,(4.3) where 𝐴 and 𝐵 are amplitudes, ω is the frequency, t is the time elapsed since the start of the periodic function reference point or subarc interval, 𝜙𝐴 and 𝜙𝐵 are the phase offsets. Either of these accelerations can be rewritten as: acc=𝐴sin(𝜔𝑡)+𝐵cos(𝜔𝑡),(4.4) where for a sine acceleration with phase, 𝐴=+𝐴cos𝜙𝐴;𝐵=+𝐵sin𝜙𝐴, and, for a cosine acceleration with phase, 𝐴=𝐵cos𝜙𝐵;𝐵=+𝐵sin𝜙𝐵. When estimated, the amplitudes 𝐴 and 𝐵 will adjust themselves to produce an effective phase offset.

Coordinate Frame
The selection of the start of the subarc can be important, especially for noncircular orbits. Conventionally, equator crossings, argument of perigee, mean anomaly, or orbit angle have been used as reference point.

4.5. Discussion of the Application

Now that the purpose of the discussion is known, it is possible to present a discussion about the subject.

In order to generate the UKF, it is necessary to rewrite the Kalman filter from UT. First, the state vector is increased, with the measurements noise vector 𝐰𝑘, yielding a vector with dimension 𝑛𝑎=𝑛+𝑞. The increased versions of the state and the covariance are

𝐱𝑎𝐱(𝑘)=𝑘𝐰𝑘,𝐏𝑎𝑘=𝐏𝑘𝟎𝟎𝐐𝑘.(4.5)

This increase can still contain 𝝊𝑘, the Gaussian process noise, of dimension 𝑙×1. The new covariance matrix would have 𝐑𝑘 in the principal diagonal, the covariance of such noise, and the new state vector dimension would be 𝑛𝑎=𝑛+𝑞+𝑙, according to Lee and Alfriend [10].

Next, the increased set of sigma points is built:

𝝌𝑎0,𝑘=𝐱𝑎𝑘,𝝌𝑎𝑖,𝑘=𝐱𝑎𝑘+(𝑛𝑎+𝜆)𝐏𝑎𝑘𝑖,𝑖=1,,𝑛𝑎,𝝌𝑎𝑖,𝑘=𝐱𝑎𝑘(𝑛𝑎+𝜆)𝐏𝑎𝑘𝑖,𝑖=𝑛𝑎+1,,2𝑛𝑎(4.6) with 𝜆=𝛼2(𝑛𝑎+𝜅)𝑛𝑎, where 𝛼 control the sigma points scattered about the mean 𝐱𝑎𝑘, and it is usually chosen small, in the interval 104𝛼1, as Jwo and Lai [11] wrote; 𝜅 provides an extra degree of freedom; ((𝑛𝑎+𝜆)𝐏𝑎𝑘)𝑖 is the 𝑖th row or column of the root square matrix of (𝑛𝑎+𝜆)𝐏𝑎𝑘.

In the propagation step, the state vector and covariance predicted are calculated based on the mean and the covariance of the propagated sigma points, transformed from the state vector and dynamical noises 𝐱𝑘+1=2𝑛𝑎𝑖=0𝑊𝑖𝝌𝑥𝑖,𝑘+1,𝐏𝑘+1=2𝑛𝑎𝑖=0𝑊𝑖𝝌𝑥𝑖,𝑘+1𝐱𝑘+1𝝌𝑥𝑖,𝑘+1𝐱𝑘+1T.(4.7)

The prediction of the observation vector and the innovation matrix, 𝐏𝜈𝜈𝑘+1, is done the same way. That means, the observation and the innovation are predicted from mean and covariance of the transformed sigma points:

𝐲𝑘+1=2𝑛𝑎𝑖=0𝑊𝑖(𝑚)y𝑖,𝑘+1,(4.8) where y𝑖,𝑘+1 represents the sigma vectors propagated through the nonlinear equation of the observation model, yielding the transformed sigma points from the state vector and the dynamics noise, shown before.

In the up to date (correction) step of measurement, the Kalman gain, 𝒦𝑘+1, is calculated based on the correlation matrix between the measurement and the observation, 𝐏𝑥𝑦𝑘+1, and the innovation matrix, both predicted 𝒦𝑘+1=𝐏𝑥𝑦𝑘+1𝐏𝜐𝜐𝑘+11𝐏,with𝜈𝜈𝑘+1=2𝑛𝑎𝑖=0𝑊(c)𝑖y𝑖,𝑘+1𝐲𝑘+1y𝑖,𝑘+1𝐲𝑘+1T,𝐏𝑥𝑦𝑘+1=2𝑛𝑎𝑖=0𝑊(c)𝑖𝝌𝑥𝑖,𝑘+1𝐱𝑘+1y𝑖,𝑘+1𝐲𝑘+1T.(4.9) Finally, the up to date state and covariance are

𝐱+𝑘+1=𝐱𝑘+1+𝒦𝑘+1𝐲𝑘+1𝐲𝑘+1,𝐏+𝑘+1=𝐏𝑘+1𝒦𝑘+1𝐏𝜈𝜈𝑘+1𝒦T𝑘+1,(4.10) where 𝐲 is the vector effectively measured in the instant 𝑘+1.

The process is repeated for the next instant, and the up to date mean (from state 𝐱+𝑘+1) and covariance will be used to deterministically generate the sigma points of the next instant.

5. Final Comments

The SPKF estimation technique has been investigated in very many different applications. After analyzing the investigations, some comments may be done.

(i)The main advantages of the SPKF are: easy to implement; computationally strong; high accuracy.(ii)The motivation for applying such technique herein: recent nonlinear state estimate techniques will be applied to the specific problem of orbit determination using real data from GPS, instead of simulated data. This can improve the results, when compared to the EKF, and get better the estimated state variables accuracy.(iii)alman filter may be used as the estimation algorithm. However, not rarely, the dynamic systems and the measurements equations are of nonlinear nature. For solving such a problem, convenient extensions of the Kalman filter have been sought. The EKF implementation in orbit estimation, under inaccurate initial conditions and scattered measurements, can lead to unstable solutions. Nevertheless, the unscented transformation was developed as a method to propagate mean and covariance information through nonlinear transformations. The SPKF appears as an emerging estimation algorithm applied to nonlinear system, without linearization steps.

Acknowledgment

The authors wish to express their appreciation for the support provided by FAPESP (The State of São Paulo Research Foundation), under Contract 07/53256-1.

References

  1. R. G. Brown and P. Y. C. Hwang, Introduction to Random Signals and Applied Kalman Filtering, John Wiley & Sons, New York, NY, USA, 3rd edition, 1985. View at Zentralblatt MATH
  2. R. D. Van Merwe, E. A. Wan, and S. I. Julier, “Sigma-point kalman filters for nonlinear estimation and sensor-fusion—applications to integrated navigation,” in Collection of Technical Papers—AIAA Guidance, Navigation, and Control Conference, vol. 3, pp. 1735–1764, Providence, RI, USA, August 2004.
  3. S. J. Julier and J. K. Uhlmann, “New extension of the Kalman filter to nonlinear systems,” in Signal Processing, Sensor Fusion, and Target Recognition VI, vol. 3068 of Proceedings of SPIE, pp. 182–193, Orlando, Fla, USA, April 1997.
  4. S. J. Julier, J. K. Uhlmann, and H. F. Durrant-Whyte, “New approach for filtering nonlinear systems,” in Proceedings of the American Control Conference, vol. 3, pp. 1628–1632, Seattle, Wash, USA, June 1995. View at Scopus
  5. S. Julier, J. Uhlmann, and H. F. Durrant-Whyte, “A new method for the nonlinear transformation of means and covariances in filters and estimators,” IEEE Transactions on Automatic Control, vol. 45, no. 3, pp. 477–482, 2000. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  6. S. J. Julier and J. K. Uhlmann, “Unscented filtering and nonlinear estimation,” Proceedings of the IEEE, vol. 92, no. 3, pp. 401–422, 2004. View at Publisher · View at Google Scholar
  7. D.-J. Lee and K. T. Alfriend, “Precise real-time orbit estimation using the unscented Kalman filter,” Advances in the Astronautical Sciences, vol. 114, supplement, pp. 1835–1854, 2003. View at Google Scholar
  8. D.-J. Lee and K. T. Alfriend, “Adaptive sigma point filtering for state and parameter estimation,” in Collection of Technical Papers—AIAA/AAS Astrodynamics Specialist Conference, vol. 2, pp. 897–916, Providence, RI, USA, August 2004.
  9. M. T. Soyka and M. A. Davis, “Estimation of periodic accelerations to improve orbit ephemeris accuracy,” Advances in the Astronautical Sciences, vol. 108, pp. 1123–1140, 2001. View at Google Scholar
  10. D.-J. Lee and K. T. Alfriend, “Sigma point filtering for sequential orbit estimation and prediction,” Journal of Spacecraft and Rockets, vol. 44, no. 2, pp. 388–398, 2007. View at Publisher · View at Google Scholar
  11. D.-J. Jwo and C.-N. Lai, “Unscented Kalman filter with nonlinear dynamic process modeling for GPS navigation,” GPS Solutions, vol. 12, no. 4, pp. 249–260, 2008. View at Publisher · View at Google Scholar