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.

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

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 10βˆ’4≀𝛼≀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βˆ’π±βˆ’π‘˜+1ξ‚„T.(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ξ€·ππœπœπ‘˜+1ξ€Έβˆ’1𝐏,withπœˆπœˆπ‘˜+1=2π‘›π‘Žπ‘–=0π‘Š(c)𝑖y𝑖,π‘˜+1βˆ’π²βˆ’π‘˜+1y𝑖,π‘˜+1βˆ’π²βˆ’π‘˜+1ξ€»T,𝐏π‘₯π‘¦π‘˜+1=2π‘›π‘Žπ‘–=0π‘Š(c)π‘–ξ‚ƒπŒπ‘₯𝑖,π‘˜+1βˆ’π±βˆ’π‘˜+1ξ‚„ξ€Ίy𝑖,π‘˜+1βˆ’π²βˆ’π‘˜+1ξ€»T.(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.