An Economic Hybrid J 2 Analytical Orbit Propagator Program Based on SARIMA Models

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


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 x t for which an approximate analytical solution is known: where x t 0 is the satellite's initial time position and velocity. Moreover, at any moment t i , a precise observation x t i can be obtained. This observation is related to x A t i by the following linear relation: where ε t i 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 ε t i series, we apply statistical techniques in time series analysis. The first n values of ε t i are used to estimate a model by means of these techniques. From this model a forecast of the ε t i 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.

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.

Mathematical Problems in Engineering 3
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 x, X as where r x x 2 y 2 z 2 , μ is the gravitational constant, α the equatorial radius of the Earth, J 2 the oblateness coefficient, and P 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 l, g, h, L, G, H . This set of canonical action-angle variables can be defined in terms of the orbital elements such as l M, g ω, h Ω, L √ μa, G μa 1 − e 2 , H μa 1 − e 2 cos i, where M, ω, Ω, a, e, i 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 J 2 is a small parameter, s sin i, and f is the true anomaly. Next, we normalize the Hamiltonian 2.2 by applying the Lie transform ϕ : l, g, h, L, G, H → l , g , h , L , G , H , 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 l : and then W 1 is computed as where η 1 − e 2 and φ f − l .

Mathematical Problems in Engineering
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 L , G , and H , and so therefore the equations of motion are obtained as dl dt dL dt dG dt dH dt 0.

2.8
By integrating 2.8 we can directly obtain that the values of the momenta L , G , and H are constants, whereas the variables l , g , and h yield 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.
Z2DN1 begins by initializing the physical parameters and the initial conditions at epoch t 0 . Next, it transforms the initial conditions into the Delaunay variables l 0 , g 0 , h 0 , L 0 , G 0 , H 0 and transports them across the inverse transformation of the Delaunay normalization l 0 , g 0 , h 0 , L 0 , G 0 , H 0 . Then, the program provides Delaunay's variables at Mathematical Problems in Engineering  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 J 2 analytical theory when compared with a more complex perturbation model is about 360 km. 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

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 ε t be a discrete time series, in an ARIMA p, d, q 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 ε t ε t − μ, μ is the mean of the original time series, and ν t is a white noise residual. B is the backward shift, such that B ε t ε t−1 , whilst d is the number of times that ε t 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 ε t are time invariants. φ B represents the autoregressive AR part, where each ε t is made up of a linear combination from prior observations p, which can be expressed as a polynomial in B of degree p in the following form:  where φ i , i 1, . . . , p, are the AR parameters. θ B represents the moving average MA part, which describes the relation of ε t with past residuals and can also be expressed as a polynomial in B of degree q in the following form: where θ i , i 1, . . . , p, are the MA parameters.
In the case that ε t series shows seasonal behaviour, it can be included in 3.1 . The extended model is known as a Seasonal ARIMA model or SARIMA p, d, q P, D, Q s and takes the following form: represent the seasonal part with periodicity s.
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 d, D, p, q, P , and Q. 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.

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.

Previous Statistical Analysis
The time series analysis begins calculating the linear relations: where x represents each of the Delaunay variables l, g, h, L, G, H , x t is the simulated data at epoch t, and x A t is the data from the analytical theory at the same epoch. Therefore, these six time series ε l t , ε g t , ε h t , ε L t , ε G t , ε H t 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 J 3 and J 4 , as well as the higher orders of the analytical solution, that is, the error of the analytical theory O J 2 2 , 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 ε x t reveals that all variables show cyclical patterns or periodicities and, moreover, there is very similar behaviour between the time series ε l t and ε g t and ε L t and ε G t . For example, this study for ε l t and ε g t can be seen in Figure 5. However ε h t and ε H t do not show any similar behaviour between them or with the rest of the time series.
On the other hand, the correlation matrix of the ε x t series is shown in Table 1. This matrix presents a strong relationship between ε l t and ε g t , as their correlation coefficient is near −1 −0.9607 , as well as between their respective conjugate momenta time series errors, ε L t and ε G t , where their correlation coefficient is near 1 0.9982 . 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   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 ε g t . Consequently, the tentative models should incorporate both seasonal and nonseasonal parameters.
We analyzed different SARIMA p, d, q P, D, Q 10 models in order to approximate the ε g t 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.
We must note that the model used to approximate the ε l t time series is also a SARIMA 6, 1, 7 3, 1, 3 10 , which confirms the similar behaviour previously detected to ε g t , 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 ε l t , ε g t and the osculating Keplerian elements a, e, g, h, i, l and state vector x, y, z,ẋ,ẏ,ż are calculated.

Z2DN1-SARIMA Hybrid-AOPP
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 ω/n, where n is the mean motion of the satellite, and the usual Garfinkel assumptions 17 with a precision of about one kilometer after

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 10 shows the distance, along-track, cross-track, and radial errors. The maximum distance error obtained from Z2DN1 AOPP is 352.076 km while for Z2DN1-SARIMA it is only 23.7489 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.

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 ε x t 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.