Two methods to evaluate the state transition matrix are implemented and analyzed to verify the computational cost and the accuracy of both methods. This evaluation represents one of the highest computational costs on the artificial satellite orbit determination task. The first method is an approximation of the Keplerian motion, providing an analytical solution which is then calculated numerically by solving Kepler's equation. The second one is a local numerical approximation that includes the effect of . The analysis is performed comparing these two methods with a reference generated by a numerical integrator. For small intervals of time (1 to 10 s) and when one needs more accuracy, it is recommended to use the second method, since the CPU time does not excessively overload the computer during the orbit determination procedure. For larger intervals of time and when one expects more stability on the calculation, it is recommended to use the first method.

1. Introduction

The orbit determination process consists of obtaining values of the parameters which 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., the Kalman filtering or least squares [16]). The dynamical system model consists of the description for the dynamics of the orbital motion of a satellite, measurement models, Earth’s rotation effects, and perturbation models. According to Montenbruck and Gill [7], these models depend on, besides the state variables that define the initial conditions, a variety of parameters that either affect the dynamical motion or the measurement process. Due to the complexity of the applied models, it is hardly possible to have a direct solution for any of these parameters from a given set of observations. To linearize the relation between the observables and the independent parameters, one should obtain simplified expressions that can be handled more easily. A linearization of the trajectory and measurement model requires a large number of partial derivatives, among them, the state transition matrix [79].

The function of the transition matrix is to relate the rectangular coordinate variations for the times and . The evaluation of the state transition matrix presents one of the highest computational costs on the artificial satellite orbit determination procedure, because it requires the evaluation of the Jacobian matrix and the integration of the current variational equations. This matrix can pose cumbersome analytical expressions when using a complex force model [10]. Binning [11] suggests one method to avoid the problem of the high computational cost and extended analytical expressions of the transition matrix. This method consists of propagating the state vector using complete force model and, then, to compute the transition matrix using a simplified force model. For sure, this simplification provides some impact on the accuracy of the orbit determination. 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. On the other hand, the inclusion of the Earth-flattening () effect in the transition matrix can be performed adopting Markley’s method [12].

In Chiaradia [13], Chiaradia et al. [14], and Gomes et al. [15], simplified and compact algorithms with low computational cost are developed for artificial satellite orbit determination in real time and onboard, using the global positioning system (GPS). The state vector, composed of the position, velocity, bias, drift, and drift rate of the GPS receiver clock, has been estimated by the extended Kalman filter in all three works. The fourth-order Runge-Kutta numerical integrator has been used to integrate the state vector. The equations of motion have considered only the perturbations due to the geopotential. The state error covariance matrix has been propagated through the transition matrix, which has been calculated considering the pure Keplerian motion. To improve the accuracy of those algorithms, two methods have been compared for calculating the transition matrix considering circular (1000 km of altitude) and elliptical (Molniya) orbits. Those methods considered the pure Keplerian motion and perturbed only by the flattening effect of the Earth. Markley's method is used to include the flattening of the Earth in the transition matrix. It is a method that allows the inclusion of more perturbations in a simple way.

In the present work, the spherical harmonic coefficients of degree and order up to 50 and the drag effect are included in the reference orbit provided by the RK78 numerical integrator (Runge-Kutta with Fehlberg coefficients of order 7-8). For this study, the atmospheric density is considered constant. For the simulations made here, circular and elliptical low (up to 300 km) orbits were used. To analyze the results, the reference orbit is compared with the two methods implemented in this work.

2. State Transition Matrix

The differential equation for the Keplerian motion is expressed by with , where is the magnitude of the satellite position vector and is the Earth gravitational constant. The equation that relates the state deviations in different times is given by where and are the position and velocity vectors at the time , respectively, and are the initial position and velocity vectors at the time , respectively, and is the transition matrix given by The submatrices , , , and are calculated in accordance with the two methods that will be shown in the next topics.

3. First Method: The Analytical Transition Matrix Solution Considering the Pure Keplerian Motion

Goodyear [16] published a method for the analytical calculation of a transition matrix for the two-body problem. This method is valid for any kind of orbit. Kuga [10] implemented this method using the same elegant and adequate formulation optimized for the Keplerian elliptical orbit problem. He performed some simplifications in this method to increase its numerical efficiency (processing time, memory, and accuracy). However, this method is used to propagate the position and velocity of the satellite and can be used in any kind of two-body orbits. Such equations are simple and easy to be developed.

According to Goodyear [16] and Shepperd [17], the four submatrices are obtained developing (2.4). Then, they are written as where is the identity matrix , is matrix , and the , , are the components of a matrix , which will be shown later, and , , , are calculated in the next topic [8, 9].

3.1. Calculating the Functions , , ,

Given , , and the propagation interval , where is the semimajor axis. The eccentric anomaly and the eccentricity for the initial orbit are calculated by with reduced to the interval 0 to 2π. The mean anomalies for the initial and propagated orbits, and , are given by with and reduced to the interval 0 to 2π. The eccentric anomaly for the propagated orbit is calculated by Kepler’s equation. Kepler’s equation is solved from an initial guess based on an approximation series, then iterated by Newton-Raphson’s method, until convergence to the level of 10−12 is achieved. The variation of the eccentric anomaly is calculated and reduced to the interval 0 to 2π: The transcendental functions , , for the elliptical orbit, according to Goodyear [16], are calculated by Therefore, the functions , , are valid only for elliptical orbits and are calculated by The propagated vectors e are given by There is no singularity problem and Kepler’s equation is solved through Newton Raphson’s method in double precision. The classical parameters a, e, i, Ω, ω, are constant in the Keplerian motion, and, therefore, there is no different subscript for them.

3.2. The Evaluation of the Matrix

First of all, it is necessary to calculate the secular component including the effect of multi-revolutions in the case where the orbit propagation time, , is larger than one orbital period: where provides the truncated integer number of the real argument x, and the variables and are related with the transcendental functions s4 and s5 (Goodyear [16]). Therefore, the components of the matrix are related to the partial derivatives of (3.8)–(3.11) with respect to the and , given by

3.3. The Property of the Transition Matrix

Sometimes the inverse matrix is required, such as in backward filters. It is also easily accomplished as follows. The inverse matrix , which propagates deviations backward from to , is given by which results from the canonic nature of the original equations [18, 19]. This also applies to the second method (Markley’s) if the perturbations are derived from a potential (e.g., effect).

4. Second Method: Markley’s Method

Markley's method uses two states, one at the time and the other at the time, and calculates the transition matrix between them by using , , , the radius of the Earth, and the two states. In this case, the effect of the Earth’s flattening is the most influent factor in the process. Markley’s method consists of making one approximation for the transition matrix of the state vector based on the Taylor series expansion for short intervals of propagation, . This method can be used by any kind of orbit and the equations are simple and easily implemented, as shown next. The state transition's differential equation is defined by where is the initial condition, and are the Cartesian state at the instant t, and are the Cartesian state at the instant t0, the matrix of zeros, the identity matrix ,   and the accelerations of the satellite.

Performing successive derivatives of (4.1), followed by substitutions, gives the derivative of the transition matrix: where The dot represents the derivative with respect to the time. Developing in Taylor's series at , using the matrices for and the initial condition , the transition matrix of the position and velocity obtained after some simplifications is given by [12] where The calculations of these matrices pose no problems, since the gradient matrix in the end of the propagating interval is a function of the final Cartesian state which should be calculated during the data processing and is available at no extra cost. The and, therefore, are symmetric if the perturbation derives from a potential. The gradient matrix including only the central force and the is given by The accelerations due the central force and the Earth’s flattening are given by The partial derivatives are [20]

5. Analyses of the Methods

First of all, the reference orbit is integrated using the numerical integrator Runge-Kutta of eighth order with automatic step size control (RK78). This integrator is implemented to integrate simultaneously the stated vectors considering the spherical harmonic coefficients up to 50th order and degree and the transition matrix considering the Earth's flattening and the atmospheric drag effects.

The transition matrix generated by the numerical integrator is used as the reference to compare the transition matrix generated by the two methods. To compare the accuracy of them, let us define the global relative error as [10] where is the component , of the transition matrix calculated by the RK78 considered as reference and is the component , of the transition matrix calculated by one of the methods. It gives a rough idea of the number of the common significant figures retained after the computations.

The first method considers the pure Keplerian motion and the second one considers the Keplerian motion and the effect. A whole day of integration is divided in intervals of 1 to 60 seconds, which produces from 86,400 to 1,440 steps of integration, respectively, as shown in Table 1. The total processing time of each method is also shown in Table 1, although it depends on the code, the computer, and the programmer skills (the evaluation used a regular PC Pentium II processor, 512 MB, running FORTRAN codes). However, it illustrates roughly the comparative expected performance. One can note that the CPU time difference at any step of integration and for any method is very small. The first method is slightly faster than the second one for small steps of time, like 1 and 10 seconds.

To analyze the accuracy of the methods, six comparisons are done, using two kinds of orbits: one circular and one elliptical. The transition matrix generated by each of the methods is compared with the reference generated by the RK78. The circular orbit is from the Topex/Poseidon satellite [21] and the elliptical one is from the Molniya satellite [12]. All those tests are performed for a period of 24 hours. The data of these orbits are shown in Table 2.

Besides, the RK4 (Runge-Kutta of fixed 4th order) numerical method for computing the transition matrix is also implemented and the transition matrix generated by the RK4 is also included in the comparison. The RK4 used the same dynamical model of the reference RK78. The results of those three comparisons for the circular orbit are shown in Table 3, and the results for the elliptical orbit are shown in Table 4. The errors are shown in terms of per step mean and standard deviations considering the whole day sample. The comparison is done for only these two kinds of orbits, where the Keplerian approximation for the transition matrix provides already reasonable accuracy. For large time interval, the second method is always disadvantageous because it is a local numerical approximation, as depicted in Tables 3 and 4.

6. Conclusions

The goal of this research was to compare and choose the most suited method to calculate the transition matrix used to propagate the covariance matrix of the position and velocity of the state estimator, (e.g., Kalman filtering or least squares), within the procedure of the artificial satellite orbit determination. The methods were evaluated according to accuracy, processing time, and handling complexity of the equations for two kinds of orbits: circular and elliptical. The processing time of the second method is around 30% larger than the first one; however, this difference is not considered enough to harm the computer load in the orbit determination tasks. The second method is more accurate for short intervals, as and 10 seconds, for the orbits considered. For other intervals of propagation, the first method shows to be more stable in the sense of keeping the same accuracy regardless of the step size. The equations of the second method are easier to be handled, which means that it is possible to include more perturbations easily. However, the first method is a closed analytical solution for the two-body problem only. The second method has no singularity and no restriction about the type of orbit; the first one is optimized for elliptical orbits (not optimized for parabolic and hyperbolic orbits). The second method performs an approximation in Taylor’s series, whereas in the first there is an analytical solution for the Keplerian motion. Therefore, for small intervals of time (1 to 10 seconds) and when one expects more accuracy, it is recommended to use the second method, since the CPU time does not overload excessively the computer in the orbit determination procedure. For larger intervals of time and when one expects more stability on the calculation, it is recommended to use the first method.


The authors wish to express their appreciation for the support provided by UNESP (Universidade Estadual Paulista “Júlio de Mesquita Filho”) of Brazil and INPE (Brazilian Institute for Space Research).