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

Abstract

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 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: 𝐱^𝒜𝑡=𝐹𝑡,𝐱_𝑡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 ℋ=12(𝐗⋅𝐗)−𝜇𝑟1−𝐽²𝛼𝑟²𝑃₂𝑧𝑟,(2.1) 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 𝑙=𝑀, 𝑔=𝜔, ℎ=Ω, 𝐿=√𝜇𝑎, 𝐺=√𝜇𝑎(1−𝑒²), 𝐻=√𝜇𝑎(1−𝑒²)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 ℋ=−𝜇²2𝐿²−𝜖2𝜇𝑟𝛼𝑟²1−3𝑠²sin²(𝑓+𝑔),(2.2) where 𝜖=𝐽₂ 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 𝒦₀=ℋ₀,(2.3)𝒦₁=ℋ₁−𝜇²𝐿^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 𝑙′: 𝒦₁=3𝛼²𝜇⁴𝑠′²4𝐿′⁶𝜂′³−𝛼²𝜇⁴2𝐿^6𝜂^3,(2.5) and then 𝒲₁ is computed as 𝒲₁=𝐿^3𝜇²ℋ¹=𝜇−𝒦₁𝑑𝑙₂𝛼²3𝑠^2−2𝜙^4𝐿^3𝜂^3+𝜇²𝛼²𝑒^{(3𝑠′2−2)}4𝐿^3𝜂^3−3𝜇₂sin𝑓′𝛼²𝑒^𝑠^28𝐿^3𝜂^3sin𝑓^+2𝑔^−3𝜇²𝛼²𝑠^28𝐿^3𝜂^3−𝜇₂sin2𝑓^+2𝑔^𝛼²𝑒^𝑠^28𝐿^3𝜂^3sin3𝑓^+2𝑔^,(2.6) where 𝜂′=√1−𝑒′² and 𝜙^=𝑓′−𝑙′.

Hence, up to the first order, the transformed Hamiltonian is given by 𝒦=−𝜇²2𝐿^2+𝜖3𝛼²𝜇⁴𝑠^24𝐿^6𝜂^3−𝛼²𝜇⁴2𝐿^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𝑡=𝜕𝜕𝒦𝐿′=𝐿𝜇′²³+𝜖3𝛼²𝜇⁴2𝐿′⁷𝜂′³−9𝛼²𝜇⁴𝑠′²4𝐿′⁷𝜂′³,d𝑔′d𝑡=𝜕𝜕𝐺𝒦′=𝜖3𝛼²𝜇⁴𝐿′⁷𝜂′⁴−15𝛼²𝜇⁴𝑠′²4𝐿′⁷𝜂′⁴,dℎ′d𝑡=𝜕𝜕𝐺𝒦′=−𝜖3𝛼²𝜇⁴𝑐^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 𝑙^=𝜇²𝐿^3+𝜖3𝛼²𝜇⁴2𝐿^7𝜂^3−9𝛼²𝜇⁴𝑠^24𝐿^7𝜂^3𝑡−𝑡⁰𝑔_+𝑙^0,=𝜖3𝛼²𝜇⁴𝐿^7𝜂^4−15𝛼²𝜇⁴𝑠^24𝐿^7𝜂^4𝑡−𝑡⁰ℎ_+𝑔^0,=−𝜖3𝛼²𝜇⁴𝑐^2𝐿^7𝜂^4𝑡−𝑡⁰+ℎ^0,(2.9) where 𝑙^0, 𝑔^0, ℎ^0, 𝐿^0, 𝐺^0, 𝐻^0 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 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 𝑡₀. Next, it transforms the initial conditions into the Delaunay variables (𝑙₀,𝑔₀,ℎ₀,𝐿₀,𝐺₀,𝐻₀) 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 𝐽₂ analytical theory when compared with a more complex perturbation model is about 360km.

(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 12.4^∘ and 9.5^∘, 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: 𝜙(𝐵)(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−𝜙₁𝐵−𝜙₂𝐵²−⋯−𝜙_𝑝𝐵^𝑝,(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−𝜃₁𝐵−𝜃₂𝐵²−⋯−𝜃_𝑞𝐵^𝑞,(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−Φ₁𝐵^𝑠−Φ₂𝐵^2𝑠−⋯−Φ_𝑃𝐵^𝑃𝑠,)=1−Θ₁𝐵^𝑠−Θ₂𝐵^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 𝐽₃ and 𝐽₄), as well as the higher orders of the analytical solution, that is, the error of the analytical theory 𝒪(𝐽²²), 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.

(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 −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−(𝐺/𝐿)²) 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.

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

(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 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(𝑝,𝑑,𝑞)(𝑃,𝐷,𝑄)₁₀ 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)₁₀ 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)₁₀ models.

We must note that the model used to approximate the 𝜀^𝑙𝑡 time series is also a SARIMA(6,1,7)(3,1,3)₁₀, 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 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.

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

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

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.