Table of Contents Author Guidelines Submit a Manuscript
Mathematical Problems in Engineering
VolumeΒ 2012, Article IDΒ 207381, 15 pages
Research Article

An Economic Hybrid 𝐽𝟐 Analytical Orbit Propagator Program Based on SARIMA Models

1Departamento de MatemΓ‘ticas y ComputaciΓ³n, Universidad de La Rioja, 26004 LogroΓ±o, Spain
2Departamento 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.


We present a new economic hybrid analytical orbit propagator program based on SARIMA models, which approximates to a 4Γ—4 tesseral analytical theory for a Quasi-Spot satellite. The J2 perturbation is described by a first-order closed-form analytical theory, whereas the effects produced by the higher orders of J2 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: π±π’œπ‘‘ξ€·=𝐹𝑑,𝐱𝑑0ξ€Έ,(1.1) where 𝐱𝑑0 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: πœ€π‘‘π‘–=π±π‘‘π‘–βˆ’π±π’œπ‘‘π‘–,(1.2) 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 ̂𝐱𝑑𝑖=π±π’œπ‘‘π‘–+Μ‚πœ€π‘‘π‘–.(1.3)

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 1β„‹=2πœ‡(𝐗⋅𝐗)βˆ’π‘Ÿξ‚Έ1βˆ’π½2ξ‚€π›Όπ‘Ÿξ‚2𝑃2ξ‚€π‘§π‘Ÿξ‚ξ‚Ή,(2.1) where π‘Ÿ=‖𝐱‖=√π‘₯2+𝑦2+𝑧2, πœ‡ is the gravitational constant, 𝛼 the equatorial radius of the Earth, 𝐽2 the oblateness coefficient, and 𝑃2 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 𝑙=𝑀, 𝑔=πœ”, β„Ž=Ξ©, √𝐿=πœ‡π‘Ž, 𝐺=βˆšπœ‡π‘Ž(1βˆ’π‘’2), 𝐻=βˆšπœ‡π‘Ž(1βˆ’π‘’2)cos𝑖, 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 πœ‡β„‹=βˆ’22𝐿2βˆ’πœ–2πœ‡π‘Ÿξ‚€π›Όπ‘Ÿξ‚2ξ€·1βˆ’3𝑠2sin2ξ€Έ,(𝑓+𝑔)(2.2) where πœ–=𝐽2 is a small parameter, 𝑠=sin𝑖, 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 𝒦0=β„‹0,(2.3)𝒦1=β„‹1βˆ’πœ‡2πΏξ…ž3πœ•π’²πœ•π‘™ξ…ž.(2.4)

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 𝑙′: 𝒦1=3𝛼2πœ‡4𝑠′24𝐿′6πœ‚β€²3βˆ’π›Ό2πœ‡42πΏξ…ž6πœ‚ξ…ž3,(2.5) and then 𝒲1 is computed as 𝒲1=πΏξ…ž3πœ‡2ξ€œξ€·β„‹1βˆ’π’¦1ξ€Έ=πœ‡π‘‘π‘™2𝛼2ξ€·3π‘ ξ…ž2ξ€Έπœ™βˆ’2ξ…ž4πΏξ…ž3πœ‚ξ…ž3+πœ‡2𝛼2π‘’ξ…ž(3𝑠′2βˆ’2)4πΏξ…ž3πœ‚ξ…ž3βˆ’sin𝑓′3πœ‡2𝛼2π‘’ξ…žπ‘ ξ…ž28πΏξ…ž3πœ‚ξ…ž3𝑓sinξ…ž+2π‘”ξ…žξ€Έβˆ’3πœ‡2𝛼2π‘ ξ…ž28πΏξ…ž3πœ‚ξ…ž3ξ€·sin2π‘“ξ…ž+2π‘”ξ…žξ€Έβˆ’πœ‡2𝛼2π‘’ξ…žπ‘ ξ…ž28πΏξ…ž3πœ‚ξ…ž3ξ€·sin3π‘“ξ…ž+2π‘”ξ…žξ€Έ,(2.6) where βˆšπœ‚β€²=1βˆ’π‘’β€²2 and πœ™ξ…ž=π‘“β€²βˆ’π‘™β€².

Hence, up to the first order, the transformed Hamiltonian is given by πœ‡π’¦=βˆ’22πΏξ…ž2+πœ–3𝛼2πœ‡4π‘ ξ…ž24πΏξ…ž6πœ‚ξ…ž3βˆ’π›Ό2πœ‡42πΏξ…ž6πœ‚ξ…ž3ξƒͺ.(2.7)

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 d𝑙′d𝑑=πœ•π’¦=πœ‡πœ•πΏβ€²2𝐿′3ξ‚΅+πœ–3𝛼2πœ‡42𝐿′7πœ‚β€²3βˆ’9𝛼2πœ‡4𝑠′24𝐿′7πœ‚β€²3ξ‚Ά,d𝑔′d𝑑=πœ•π’¦ξ‚΅πœ•πΊβ€²=πœ–3𝛼2πœ‡4𝐿′7πœ‚β€²4βˆ’15𝛼2πœ‡4𝑠′24𝐿′7πœ‚β€²4ξ‚Ά,dβ„Žβ€²d𝑑=πœ•π’¦πœ•πΊβ€²=βˆ’πœ–3𝛼2πœ‡4π‘ξ…ž2πΏξ…ž7πœ‚ξ…ž4,d𝐿′d𝑑=d𝐺′d𝑑=d𝐻′d𝑑=0.(2.8)

By integrating (2.8) we can directly obtain that the values of the momenta 𝐿′, πΊξ…ž, and 𝐻′ are constants, whereas the variables 𝑙′, π‘”ξ…ž, and β„Žβ€² yield π‘™ξ…ž=ξƒ¬πœ‡2πΏξ…ž3+πœ–3𝛼2πœ‡42πΏξ…ž7πœ‚ξ…ž3βˆ’9𝛼2πœ‡4π‘ ξ…ž24πΏξ…ž7πœ‚ξ…ž3ξ€·ξƒͺξƒ­π‘‘βˆ’π‘‘0ξ€Έ+π‘™ξ…ž0,π‘”ξ…ž=ξƒ¬πœ–ξƒ©3𝛼2πœ‡4πΏξ…ž7πœ‚ξ…ž4βˆ’15𝛼2πœ‡4π‘ ξ…ž24πΏξ…ž7πœ‚ξ…ž4ξ€·ξƒͺξƒ­π‘‘βˆ’π‘‘0ξ€Έ+π‘”ξ…ž0,β„Žξ…ž=ξƒ¬βˆ’πœ–3𝛼2πœ‡4π‘ξ…ž2πΏξ…ž7πœ‚ξ…ž4ξƒ­ξ€·π‘‘βˆ’π‘‘0ξ€Έ+β„Žξ…ž0,(2.9) where π‘™ξ…ž0, π‘”ξ…ž0, β„Žξ…ž0, πΏξ…ž0, πΊξ…ž0, π»ξ…ž0 are the transformed initial conditions 𝑙0, 𝑔0, β„Ž0, 𝐿0, 𝐺0, 𝐻0 at the epoch 𝑑0.

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 93 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 𝑑0. Next, it transforms the initial conditions into the Delaunay variables (𝑙0,𝑔0,β„Ž0,𝐿0,𝐺0,𝐻0) and transports them across the inverse transformation of the Delaunay normalization (π‘™ξ…ž0,π‘”ξ…ž0,β„Žξ…ž0,πΏξ…ž0,πΊξ…ž0,π»ξ…ž0). 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 (π‘Ž=7148, 𝑒=0.001, 𝑖=98∘).

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 𝐽2 analytical theory when compared with a more complex perturbation model is about 360km.

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 12.4∘ and 9.5∘, respectively.

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: πœ™(𝐡)(1βˆ’π΅)π‘‘Μƒπœ€π‘‘=πœƒ(𝐡)πœˆπ‘‘,(3.1) where Μƒπœ€π‘‘=πœ€π‘‘βˆ’πœ‡, πœ‡ is the mean of the original time series, and πœˆπ‘‘ is a white noise residual. 𝐡 is the backward shift, such that π΅Μƒπœ€π‘‘=Μƒπœ€π‘‘βˆ’1, 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: πœ™(𝐡)=1βˆ’πœ™1π΅βˆ’πœ™2𝐡2βˆ’β‹―βˆ’πœ™π‘π΅π‘,(3.2) where πœ™π‘–, 𝑖=1,…,𝑝, 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: πœƒ(𝐡)=1βˆ’πœƒ1π΅βˆ’πœƒ2𝐡2βˆ’β‹―βˆ’πœƒπ‘žπ΅π‘ž,(3.3) where πœƒπ‘–, 𝑖=1,…,𝑝, 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: Ξ¦(𝐡𝑠)πœ™(𝐡)(1βˆ’π΅π‘ )𝐷(1βˆ’π΅)π‘‘Μƒπœ€π‘‘=Θ(𝐡𝑠)πœƒ(𝐡)πœˆπ‘‘,(3.4) where Ξ¦(𝐡𝑠)=1βˆ’Ξ¦1π΅π‘ βˆ’Ξ¦2𝐡2π‘ βˆ’β‹―βˆ’Ξ¦π‘ƒπ΅π‘ƒπ‘ ,Θ(𝐡𝑠)=1βˆ’Ξ˜1π΅π‘ βˆ’Ξ˜2𝐡2π‘ βˆ’β‹―βˆ’Ξ˜π‘žπ΅π‘„π‘ ,(3.5) 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 10 satellite cycles. This number of cycles was experimentally calculated; however the models obtained from fewer than 10 cycles were less accurate, but from above 10 cycles the increase in accuracy was not significant either. The force model used to generate the simulated data is a 4Γ—4 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: πœ€π‘₯𝑑=π‘₯π‘‘βˆ’π‘₯π’œπ‘‘,(4.1) 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 𝐽3 and 𝐽4), as well as the higher orders of the analytical solution, that is, the error of the analytical theory π’ͺ(𝐽22), during 10 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.

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 βˆ’1 (βˆ’0.9607), as well as between their respective conjugate momenta time series errors, πœ€πΏπ‘‘ and πœ€πΊπ‘‘, where their correlation coefficient is near 1 (0.9982).

Table 1: Correlation matrix of the πœ€π‘₯𝑑 series.

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 (βˆšπ‘’=1βˆ’(𝐺/𝐿)2) 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 βˆ’1 value in the correlation coefficient, whilst πœ€πΏπ‘‘ and πœ€πΊπ‘‘ are almost the same, and therefore the sign in the correlation coefficient is positive with a near 1 value.

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 10 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 0.63km 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.

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 0.6921>0.05). Moreover, its periodogram (see Figure 5) shows high peaks at low frequencies (𝑓=0.01 and 0.02). 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(𝑝,𝑑,π‘ž)(𝑃,𝐷,𝑄)10 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(6,1,7)(3,1,3)10 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 0.182 and 0.993, respectively.

Table 2: SARIMA(6,1,7)(3,1,3)10 models.

We must note that the model used to approximate the πœ€π‘™π‘‘ time series is also a SARIMA(6,1,7)(3,1,3)10, 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.


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 30 days, involves several mathematical expressions of more than 10000 terms, whereas the whole Z2DN1 analytical theory only needs to evaluate 93 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 10 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 9.4∘ and 8.2∘, respectively. These errors have been reduced to 3∘ in the case of the mean anomaly and to 1.3∘ in the argument of the perigee, with respect to the Z2DN1 AOPP.

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 352.076km while for Z2DN1-SARIMA it is only 23.7489km. 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.

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, 10 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.


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.


  1. D. Brouwer, β€œSolution of the problem of artificial satellite theory without drag,” The Astronomical Journal, vol. 64, pp. 378–397, 1959. View at Publisher Β· View at Google Scholar
  2. K. Aksnes, β€œA second-order artificial satellite theory based on an intermediate orbit,” Astronomical Journal, vol. 75, no. 9, pp. 1066–1076, 1970. View at Publisher Β· View at Google Scholar
  3. H. Kinoshita, β€œThird-order solution of an artificial satellite theory, Smithsonian Astrophysical Observatory,” Special Report No. 379, Cambridge, Mass, USA, 1977. View at Google Scholar
  4. 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.
  5. G. E. P. Box and G. M. Jenkins, Time Series Analysis: Forecasting and Control, Holden-Day, San Francisco, Calif, USA, 1976.
  6. A. Deprit, β€œDelaunay normalisations,” Celestial Mechanics and Dynamical Astronomy, vol. 26, no. 1, pp. 9–21, 1982. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH
  7. A. A. Kamel, β€œPeturbation method in the theory of nonlinear oscillations,” Celestial Mechanics, vol. 3, no. 1, pp. 90–106, 1970. View at Publisher Β· View at Google Scholar
  8. 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. View at Publisher Β· View at Google Scholar Β· View at Scopus
  9. 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. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH
  10. 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. View at Google Scholar Β· View at Zentralblatt MATH
  11. 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. View at Publisher Β· View at Google Scholar
  12. 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. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH Β· View at Scopus
  13. R Development Core Team, R: A Language and Environment For Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria, 2011,
  14. K. Chan, β€œTSA: Time Series Analysis, R package version 0. 98,” 2010,
  15. R. J. Hyndman, β€œForecast: Forecasting functions for time series, R package age version 2.09,” 2010,
  16. A. Trapletti and K. Hornik, β€œTseries: Time Series Analysis and Computational Finance, R package version 0.10-27,” 2011,
  17. B. Garfinkel, β€œTesseral harmonic perturbations of an artificial satellite,” The Astronomical Journal, vol. 70, pp. 784–786, 1965. View at Publisher Β· View at Google Scholar
  18. 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. View at Google Scholar
  19. 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.