Mathematical Problems in Engineering

Volume 2012 (2012), Article ID 207381, 15 pages

doi:10.1155/2012/207381

Research Article

An Economic Hybrid Analytical Orbit Propagator Program Based on SARIMA Models

Juan Félix San-Juan,¹ Montserrat San-Martín,¹ and Iván Pérez²

¹Departamento de Matemáticas y Computación, Universidad de La Rioja, 26004 Logroño, Spain

²Departamento de Ingeniería Eléctrica, Universidad de La Rioja, 26004 Logroño, Spain

Received 28 April 2012; Accepted 8 July 2012

Academic Editor: Carlo Cattani

Copyright © 2012 Juan Félix San-Juan 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.

We present a new economic hybrid analytical orbit propagator program based on SARIMA models, which approximates to a tesseral analytical theory for a Quasi-Spot satellite. The J ₂ perturbation is described by a first-order closed-form analytical theory, whereas the effects produced by the higher orders of J ₂ and the perturbation of the rest of zonal and tesseral harmonic coefficients are modelled by SARIMA models. Time series analysis is a useful statistical prediction tool, which allows building a model for making future predictions based on the study of past observations. The combination of the analytical techniques and time series analysis allows an increase in accuracy without significant loss in efficiency of the new propagators, as a consequence of modelling higher-order terms and other perturbations are not taken into account in the analytical theory.

1. Introduction

An analytical orbit propagator program (AOPP) is an application which collects and arranges all mathematical expressions involved in an approximate analytical solution of the satellite equations of motion. The analytical solutions are known as General Perturbation Theories. It is noteworthy that the perturbation force model used and the order of the analytical approximation are closely related to the accuracy and computational efficiency of an AOPP.

In many situations, in order to improve the accuracy of the solution, it may be necessary to consider a more precise perturbation force model. The solution provided by the General Perturbation Theories may not be the best approach, because the calculating process generates unmanageably large mathematical expressions and, therefore, reduces the computational efficiency of its corresponding AOPP. Other alternatives, although computationally more expensive than an economic analytical approximation, are the Special Perturbation Theories, which directly integrate the equations of motion using numerical techniques, or by means of semi-analytical theories, which are a combination of General and Special Perturbation Theories.

In this paper we present a new methodology, which we will call Hybrid Perturbation Theories, to carry out new families of hybrid orbit propagator programs which combine a simplified analytical orbit propagator [1–4] with statistical time series models [5]. This combination allows an increase in accuracy for predicting the position of a satellite without significant loss in computational efficiency in the new hybrid propagators, as well as modelling higher-order terms and other perturbations not considered in the analytical theory.

Mathematically, the problem consists of estimating the satellite’s position and velocity for which an approximate analytical solution is known:

where is the satellite’s initial time position and velocity. Moreover, at any moment , a precise observation can be obtained. This observation is related to by the following linear relation:

where represents the errors produced by the perturbation forces not considered in the analytical theory and by the selfsame approximate analytical solution. In order to predict the future values of the series, we apply statistical techniques in time series analysis.

The first values of are used to estimate a model by means of these techniques. From this model a forecast of the error can be calculated. Finally these estimations are used to obtain the forecast of the satellite’s position and velocity by the relation

In this paper, the orbit propagator Z2DN1 derived from a first-order closed-form analytical integration of the main problem of the artificial satellite theory and the SARIMA time series models are described. Secondly, using the univariate Box-Jenkins time series analysis, a specific Z2DN1-SARIMA model is developed for a Quasi-Spot satellite so as to model the effects of some zonal and tesseral harmonics by means of the statistical part, where these influences have not been taken into consideration in the analytical part. The simulated data are obtained from the numerical integration for an Earth orbiter, which has only taken into account the perturbation due to the nonsymmetrical Earth gravity field up to the fourth degree and order. Finally, we compare the simulation with both the analytical propagator alone and the analytical-statistical hybrid propagator.

2. Z2DN1 Analytical Orbit Propagator Program

This AOPP has been derived from a first-order closed-form analytical theory of the main problem of the artificial satellite theory.

The main problem is defined as a Kepler problem perturbed by Earth’s oblateness. The Hamiltonian of this dynamical system can be written in a cartesian coordinate system (, ) as

where , is the gravitational constant, the equatorial radius of the Earth, the oblateness coefficient, and the second degree Legendre polynomial.

The first step to carry out the analytical theory consists of expressing the Hamiltonian (2.1) in terms of the Delaunay variables . This set of canonical action-angle variables can be defined in terms of the orbital elements such as , , , , , , where , , , , , are the mean anomaly, argument of the perigee, longitude of the ascending node, semimajor axis, eccentricity, and inclination, respectively. Then the transformed Hamiltonian is given as

where is a small parameter, , and is the true anomaly.

Next, we normalize the Hamiltonian (2.2) by applying the Lie transform , the so-called Delaunay Normalization [6], which up to first order reads

The Lie method solves (2.4) by choosing the form of the transformed Hamiltonian; the Delaunay Normalization takes the Hamiltonian as the average over the fastest angle :

and then is computed as

where and .

Hence, up to the first order, the transformed Hamiltonian is given by

We must remark that the Hamiltonian (2.7) is integrable. This Hamiltonian only depends on the momenta , , and , and so therefore the equations of motion are obtained as

By integrating (2.8) we can directly obtain that the values of the momenta , , and are constants, whereas the variables , , and yield

where , , , , , are the transformed initial conditions , , , , , at the epoch .

Finally, from (2.6) the first-order explicit equations of the direct and inverse transformations [7] are calculated.

From the above analytical theory an AOPP was derived, which has to evaluate terms. This AOPP has been called Z2DN1. The algebraic manipulations required to carry out this analytical theory and its corresponding AOPP were built using a set of Mathematica packages called MathATESAT [8]. Figure 1 shows the flowchart of the Z2DN1 analytical orbit propagator program.

Figure 1: Flowchart of the Z2DN1 orbit propagator program.

Z2DN1 begins by initializing the physical parameters and the initial conditions at epoch . Next, it transforms the initial conditions into the Delaunay variables and transports them across the inverse transformation of the Delaunay normalization . Then, the program provides Delaunay’s variables at epoch from integrated Hamilton equations . Finally, the direct transformation of the Delaunay normalization is applied, and therefore the osculating Keplerian elements and the state vector can be calculated.

This model has been compared with the numerical integration (8th-order Runge-Kutta method) of the equation of motion of a model, which includes the Earth’s zonal and tesseral harmonic coefficients of fourth degree and order in the case of a Quasi-Spot satellite (, , ).

Figure 2 shows the distance, along-track, cross-track, and radial errors in a time span interval of 30 days, which is about 430 satellite cycles. As can be observed, the distance error of the first-order analytical theory when compared with a more complex perturbation model is about km.

(a)

(b)

(c)

(d)

Figure 2: Distance, along-track, cross-track, and radial errors for a Quasi-Spot satellite.

Figure 3 shows the relative errors of the orbital elements for a Quasi-Spot satellite. The mean anomaly and argument of the perigee are the variables which present the worst performance. The maximum absolute errors in a time span interval of 30 days are about and , respectively.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3: Relative errors of the orbital elements for a Quasi-Spot satellite.

3. Statistical Time Series Analysis: SARIMA Model

Introduced by Box and Jenkins [5], the autoregressive integrated moving average (ARIMA) model has been one of the most popular approaches for time series forecasting. Let be a discrete time series, in an ARIMA model, in which the future value of a series is assumed to be a linear combination of its own past values and past residuals, expressed as follows:

where , is the mean of the original time series, and is a white noise residual. is the backward shift, such that , whilst is the number of times that needs to be differentiated to ensure its conversion to a stationary time series, that is, a time series in which the mean, variance, and autocorrelation functions of are time invariants. represents the autoregressive (AR) part, where each is made up of a linear combination from prior observations , which can be expressed as a polynomial in of degree in the following form:

where , , are the AR parameters. represents the moving average (MA) part, which describes the relation of with past residuals and can also be expressed as a polynomial in of degree in the following form:

where , , are the MA parameters.

In the case that series shows seasonal behaviour, it can be included in (3.1). The extended model is known as a Seasonal ARIMA model or SARIMA and takes the following form:

where

represent the seasonal part with periodicity .

To determine a suitable SARIMA model for a given series, we use the 3-stage Box-Jenkins methodology. This procedure is illustrated in Figure 4. At the identification stage, a preliminary SARIMA model is proposed from the analysis of the estimated autocorrelation function (ACF) and partial autocorrelation function (PACF), allowing us to determine the parameters , , , , , and . Then, the seasonal and nonseasonal AR and MA parameters are estimated at the second stage. The last stage, diagnostic checking, determines whether the proposed model is adequate or not. If the model is considered adequate, it can be used for forecasting future values; otherwise the process is repeated until a satisfactory model is found.

Figure 4: 3-stage Box-Jenkins methodology.

4. Time Series Analysis

In order to carry out the statistical part of the hybrid propagator, we consider a simulated data set, that is, position and velocity, taken from the numerical integration of the Quasi-Spot equations of motion during satellite cycles. This number of cycles was experimentally calculated; however the models obtained from fewer than cycles were less accurate, but from above cycles the increase in accuracy was not significant either. The force model used to generate the simulated data is a EGM-96 gravity field, whereas for the numerical integration a high-order Runge-Kutta method [9] is used.

It is noteworthy to mention that different sets of canonical and noncanonical variables can be used to develop the statistical part, such as cartesian variables, orbital elements, Delaunay variables, and polar-nodal variables. In this work, we will only take into account the Delaunay variables, which allow a direct visualization of the geometry of the orbit.

4.1. Previous Statistical Analysis

The time series analysis begins calculating the linear relations:

where represents each of the Delaunay variables , is the simulated data at epoch , and is the data from the analytical theory at the same epoch. Therefore, these six time series allocate all the information related to the perturbation forces not considered in the analytical theory (tesseral terms of fourth degree and order and the zonal coefficients and ), as well as the higher orders of the analytical solution, that is, the error of the analytical theory , during cycles (and 10 data points per cycle).

Then, the periodogram, a mathematical tool for examining cyclical behaviour in time series, and the autocorrelation functions, a measure of how a time series is correlated with itself at different time delays, are used to identify the time series models (see [5], for further details).

The study of the periodogram and autocorrelation functions of each reveals that all variables show cyclical patterns or periodicities and, moreover, there is very similar behaviour between the time series and and and . For example, this study for and can be seen in Figure 5. However and do not show any similar behaviour between them or with the rest of the time series.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 5: Periodograms and autocorrelation functions of (left) and (right).

On the other hand, the correlation matrix of the series is shown in Table 1. This matrix presents a strong relationship between and , as their correlation coefficient is near (), as well as between their respective conjugate momenta time series errors, and , where their correlation coefficient is near ().

Table 1: Correlation matrix of the series.

1.0000

−0.9607
1.0000

−0.0575
−0.1713
1.0000

−0.1037
0.0934
0.0451
1.0000

−0.1128
0.1168
−0.0070
0.9982
1.0000

0.4033
−0.4114
0.0420
0.0066
0.0005
1.0000

It is noteworthy to mention that although the intrinsic nature of the mean anomaly and the argument of the perigee is different, as mean anomaly is related to short-periodic terms and the argument of the perigee is related to long-periodic terms, the similar behaviours detected in the above statistical studies can be explained, because the Quasi-Spot satellite is near a repeat ground track orbit, in which the argument of the perigee and eccentricity () are almost constant. Figure 6 shows , and , time series. As can be observed, and are almost symmetric with respect to the -axis, which explains the negative sign and the near value in the correlation coefficient, whilst and are almost the same, and therefore the sign in the correlation coefficient is positive with a near value.

(a)

(b)

Figure 6: Blue represents , and errors. Dashed red represents and errors.

Finally, this preliminary study can be completed by analyzing the results obtained when the time series are combined with the data obtained from the analytical theory during the first satellite cycles. This test allows us to consider several possibilities. The first consists of considering each series separately, for instance, . In all these cases the accuracy is not as good as the approach given by the Z2DN1 AOPP. After considering other possibilities we show the relations obtained in previous statistical analyses:

(i), (ii), (iii), (iv).

Figure 7 shows the distance errors between the simulated data and analytical theory for the first ten cycles, and the simulated data and the above corrected analytical theories with the exact error added. The strong influence of and can be seen in the first plot; the distance error is reduced to km after ten satellite cycles, whereas and only remove part of the short-period variations, as can be seen in the second plot. The third case collects the corrections due to , , , and , which produce a distance error similar to the first case. Finally, the corrections due to and do not have any effect on the distance error, as can be observed in the last plot.

(a)

(b)

(c)

(d)

Figure 7: Blue represents the distance error between and , whereas dashed red is the distance error between and the corrected from the exact .

Next we focus our attention on carrying out a hybrid-AOPP from . The following step in the process of looking for the most suitable SARIMA models using the Box-Jenkins methodology is described below.

4.2. Time Series Estimation of and

To estimate the model of the time series, we use the Box-Jenkins methodology. In the first step, the stationary behaviour of the time series is analyzed. Figure 6 suggests that the variance is time-invariant, whereas for the mean value the plot is not conclusive. On the other hand, Figure 5 shows that the autocorrelation function (ACF) decreases slowly and the augmented Dickey-Fuller test [10] allows accepting the null hypothesis that the time series has a unit root (P value ). Moreover, its periodogram (see Figure 5) shows high peaks at low frequencies ( and ). Therefore, the time series does not seem stationary; thus differentiating the time series data may be necessary.

The second step analyzes the periodicity. The ACF shows a pronounced cyclical fluctuation with a strong correlation at lag 10. Besides, its periodogram shows a peak at the frequency of 0.2, which corresponds to a periodicity of 10. These patterns agree with the satellite cycle. This suggests that a seasonal model might be adequate to estimate . Consequently, the tentative models should incorporate both seasonal and nonseasonal parameters.

We analyzed different SARIMA models in order to approximate the time series, where the maximum likelihood method was used to estimate model parameters, as can be seen in Table 2. Finally, the diagnostic stage showed a good fit for the SARIMA model, in which the Jarque-Bera and Ljung-Box tests [11, 12] do not reject the null hypothesis of normality nor the no autocorrelation of residuals, with P values and , respectively.

Table 2: SARIMA models.

Coefficients

−0.8167142581776827
−0.8312754102958296

0.5740668630692614
0.5649950239467217

1.1200190576112987
1.1322145773449677

−0.3087269639691513
−0.2946598587646277

−0.8631472597887673
−0.8663595825987899

−0.4970128310547080
−0.5012076262623468

0.8789635328932761
0.8901438163231377

−0.9965719584754135
−0.9834775736328421

−1.6921042908035564
−1.7044561386678760

0.5206249360614155
0.4906661687944900

1.5190020299259381
1.5136974115875366

0.3606298858846529
0.3820758183271771

−0.4721289483358399
−0.4568107050477830

1.4676902604680244
1.4557054355067451

−1.3618276394078148
−1.3525881149216927

0.7076454044934742
0.7023619159769664

0.8412748367694940
0.8370722835862818

−0.8240640881053386
−0.8282527898256336

0.9939014392344844
0.9966285211365780

We must note that the model used to approximate the time series is also a SARIMA, which confirms the similar behaviour previously detected to , although the model parameters are slightly different, as can be seen in Table 2.

The GNU software R (version 2.14) [13] was the statistical tool used to perform all statistical analyses. In particular, the R packages TSA [14], forecast [15], and tseries [16] were used for all time series analyses.

5. Z2DN1-SARIMA Hybrid-AOPP

Figure 8 shows the flowchart of the Z2DN1-SARIMA hybrid-AOPP. We find the difference to the pure Z2DN1 AOPP (Figure 1) after applying the direct transformation of the Delaunay normalization. At this point, the Delaunay variables are combined with the new forecast and the osculating Keplerian elements and state vector are calculated.

Figure 8: Flowchart of the Z2DN1-SARIMA hybrid orbit propagator program.

At this point, it is noteworthy that a pure analytical theory which takes into account the perturbation of the fourth degree and order harmonic coefficients of the gravity field, considering the dimensionless parameter , where is the mean motion of the satellite, and the usual Garfinkel assumptions [17] with a precision of about one kilometer after days, involves several mathematical expressions of more than terms, whereas the whole Z2DN1 analytical theory only needs to evaluate terms. The technical details of the tesseral analytical theory have been developed in [18, 19].

5.1. Numerical Validations

Finally we analyze the behaviour of Z2DN1-SARIMA hybrid-AOPP designed for a Quasi-Spot satellite versus the numerical integration for an Earth orbiter, which has only taken into account the perturbation due to the nonsymmetrical Earth gravity field up to the fourth degree and order. The first cycles are considered for the estimation stage, whilst from the 10th and up to approximately the 430th cycle, which is about 30 days, are used in the forecasting stage.

Figure 9 shows the relative errors of the mean anomaly and argument of the perigee for a Quasi-Spot satellite. The maximum absolute errors of mean anomaly and argument of the perigee in a time span interval of 30 days are about and , respectively. These errors have been reduced to in the case of the mean anomaly and to in the argument of the perigee, with respect to the Z2DN1 AOPP.

(a)

(b)

Figure 9: Z2DN1-SARIMA hybrid-AOPP. Relative errors of the mean anomaly and argument of the perigee for a Quasi-Spot satellite.

Figure 10 shows the distance, along-track, cross-track, and radial errors. The maximum distance error obtained from Z2DN1 AOPP is km while for Z2DN1-SARIMA it is only km. We must remark that the accuracy obtained by the described hybrid-AOPP is only comparable to a higher-order analytical theory, which includes a more precise perturbation model.

(a)

(b)

(c)

(d)

Figure 10: Z2DN1-SARIMA hybrid-AOPP. Distance, along-track, cross-track, and radial errors for a Quasi-Spot satellite.

6. Conclusions and Future Works

A new methodology to carry out hybrid-AOPP families, based on the combination of an analytical orbit propagator program and statistical time series models, is presented. To illustrate this methodology, a hybrid-AOPP, named Z2DN1-SARIMA, has been developed, which combines an economic first-order closed-form analytical orbit propagator and two SARIMA time series models fitted to the case of the Quasi-Spot satellite. Although the increment in the computational time cost is not significant with respect to the pure analytical theory, the error of our theory is reduced in comparison to the pure Z2DN1 AOPP. The accuracy reached by our new hybrid model is similar to that obtained by a more complex zonal and tesseral analytical theory, but without the inconvenience of losing computational efficiency.

To calculate the SARIMA models, satellite cycles are considered and the univariate Box-Jenkins time series analysis is used to model the time series, using statistical software packages for R. Two of the six components were modelled, whilst, at present, we are working on the study of the bivariate SARIMA models in order to collect the similar behaviour found between the mean anomaly and the argument of the perigee. In the study of the argument of the node and the third component of angular momentum behaviour, we are performing an economic analytical theory, which includes tesseral coefficients.

The behaviour of the Z2DN1-SARIMA hybrid-AOPP with respect to other initial conditions near the Quasi-Spot conditions, as well as the adapted hybrid-AOPP, when other perturbations, like atmospheric drag, third body, and so on, are taken into account, is future works to be investigated.

Acknowledgments

This work has been supported by the Gobierno de La Rioja (Project Fomenta 2010/16). Special thanks are given to the reviewer for constructive advice and thoughtful reviews and comments pointed out and to K. Hatcher for the revision of this paper.

References

D. Brouwer, “Solution of the problem of artificial satellite theory without drag,” The Astronomical Journal, vol. 64, pp. 378–397, 1959.

K. Aksnes, “A second-order artificial satellite theory based on an intermediate orbit,” Astronomical Journal, vol. 75, no. 9, pp. 1066–1076, 1970.

H. Kinoshita, “Third-order solution of an artificial satellite theory, Smithsonian Astrophysical Observatory,” Special Report No. 379, Cambridge, Mass, USA, 1977.

F. R. Hoots and R. L. Roehrich, Spacetrack Report No. 3: Models for Propagation of the NORAD Element Sets, U.S. Air Force Aerospace Defense Command, Colorado Springs, Colo, USA, 1980.

G. E. P. Box and G. M. Jenkins, Time Series Analysis: Forecasting and Control, Holden-Day, San Francisco, Calif, USA, 1976.

A. Deprit, “Delaunay normalisations,” Celestial Mechanics and Dynamical Astronomy, vol. 26, no. 1, pp. 9–21, 1982.

A. A. Kamel, “Peturbation method in the theory of nonlinear oscillations,” Celestial Mechanics, vol. 3, no. 1, pp. 90–106, 1970.

J. F. San-Juan, L. M. López, and R. López, “MathATESAT: a symbolic-numeric environment in astrodynamics and celestial mechanics,” Lecture Notes in Computer Science, vol. 6783, no. 2, pp. 436–449, 2011.

J. R. Dormand and P. J. Prince, “Practical runge-kutta processes,” SIAM Journal on Scientific and Statistical Computing, vol. 10, no. 5, pp. 977–989, 1989.

D. A. Dickey and W. A. Fuller, “Distribution of the estimators for autoregressive time series with a unit root,” Journal of the American Statistical Association, vol. 74, no. 366, pp. 427–431, 1979.

C. M. Jarque and A. K. Bera, “Efficient tests for normality, homoscedasticity and serial independence of regression residuals,” Economics Letters, vol. 6, no. 3, pp. 255–259, 1980.

G. M. Ljung and G. E. P. Box, “On a measure of lack of fit in time series models,” Biometrika, vol. 65, no. 2, pp. 297–303, 1978.

R Development Core Team, R: A Language and Environment For Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria, 2011, http://www.Rproject.org/.

K. Chan, “TSA: Time Series Analysis, R package version 0. 98,” 2010, http://CRAN.R-project.org/package=TSA/

R. J. Hyndman, “Forecast: Forecasting functions for time series, R package age version 2.09,” 2010, http://CRAN.R-project.org/package=forecast/

A. Trapletti and K. Hornik, “Tseries: Time Series Analysis and Computational Finance, R package version 0.10-27,” 2011, http://CRAN.Rproject.org/package=tseries/

B. Garfinkel, “Tesseral harmonic perturbations of an artificial satellite,” The Astronomical Journal, vol. 70, pp. 784–786, 1965.

J. F. San-Juan and S. Serrano, “ATESAT: revisited and improvements. Analytical theory and numerical validation for a SPOT-like satellite,” Tech. Rep. No. DTS/MPI/MS/MN/2000-013, Centre National d'Etudes Spatiales, Toulouse, France, 2000.

S. Serrano, Teoras Analtícas del Movimiento de un Satélite Artificial alrededor de un Planeta. Ordenación Asintótica del Potencial en el Espacio [Fśico Ph.D.], Universidad de Zaragoza, 2003.