# An Economic HybridJ2 Analytical Orbit Propagator Program Based on SARIMA Models

Authors: Juan Félix San-Juan; Montserrat San-Martín; Iván Pérez
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

We present a new economic hybrid analytical orbit propagator program based on SARIMA models, which approximates to a4×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.


---

## Body

## 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 callHybrid 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 velocityxt for which an approximate analytical solution is known:
(1.1)xtA=F(t,xt0),
where xt0 is the satellite’s initial time position and velocity. Moreover, at any moment ti, a precise observation xti can be obtained. This observation is related to xti𝒜 by the following linear relation:
(1.2)εti=xti-xtiA,
where εti 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 εti series, we apply statistical techniques in time series analysis.The firstn values of εti are used to estimate a model by means of these techniques. From this model a forecast of the ε^ti 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)x^ti=xtiA+ε^ti.In this paper, the orbit propagator Z2DN1 derived from a first-order closed-form analytical integration of themain 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 themain 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 (x, X) as
(2.1)H=12(X⋅X)-μr[1-J2(αr)2P2(zr)],
where r=∥x∥=x2+y2+z2, μ is the gravitational constant, α the equatorial radius of the Earth, J2 the oblateness coefficient, and P2 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-e2), H=μa(1-e2)cosi, 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
(2.2)H=-μ22L2-ϵ2μr(αr)2(1-3s2sin⁡2(f+g)),
where ϵ=J2 is a small parameter, s=sini, 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
(2.3)K0=H0,(2.4)K1=H1-μ2L′3∂W∂l′.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′:
(2.5)K1=3α2μ4s′24L′6η′3-α2μ42L′6η′3,
and then 𝒲1 is computed as
(2.6)W1=L′3μ2∫‍(H1-K1)dl=μ2α2(3s′2-2)ϕ′4L′3η′3+μ2α2e′(3s′2-2)4L′3η′3sin⁡f′-3μ2α2e′s′28L′3η′3sin⁡(f′+2g′)-3μ2α2s′28L′3η′3sin⁡(2f′+2g′)-μ2α2e′s′28L′3η′3sin⁡(3f′+2g′),
where η′=1-e′2 and ϕ′=f′-l′.Hence, up to the first order, the transformed Hamiltonian is given by(2.7)K=-μ22L′2+ϵ(3α2μ4s′24L′6η′3-α2μ42L′6η′3).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
(2.8)dl′dt=∂K∂L′=μ2L′3+ϵ(3α2μ42L′7η′3-9α2μ4s′24L′7η′3),dg′dt=∂K∂G′=ϵ(3α2μ4L′7η′4-15α2μ4s′24L′7η′4),dh′dt=∂K∂G′=-ϵ3α2μ4c′2L′7η′4,dL′dt=dG′dt=dH′dt=0.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
(2.9)l′=[μ2L′3+ϵ(3α2μ42L′7η′3-9α2μ4s′24L′7η′3)](t-t0)+l0′,g′=[ϵ(3α2μ4L′7η′4-15α2μ4s′24L′7η′4)](t-t0)+g0′,h′=[-ϵ3α2μ4c′2L′7η′4](t-t0)+h0′,
where l0′, g0′, h0′, L0′, G0′, H0′ are the transformed initial conditions l0, g0, h0, L0, G0, H0 at the epoch t0.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 evaluate93 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 epocht0. Next, it transforms the initial conditions into the Delaunay variables (l0,g0,h0,L0,G0,H0) and transports them across the inverse transformation of the Delaunay normalization (l0′,g0′,h0′,L0′,G0′,H0′). Then, the program provides Delaunay’s variables at epoch tf from integrated Hamilton equations (l′,g′,h′,L′,G′,H′). Finally, the direct transformation of the Delaunay normalization is applied, and therefore the osculating Keplerian elements (a,e,g,h,i,l) and the state vector (x,y,z,x˙,y˙,z˙) 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 (a=7148, e=0.001, i=98∘).Figure2 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 J2 analytical theory when compared with a more complex perturbation model is about 360km.Distance, along-track, cross-track, and radial errors for a Quasi-Spot satellite.
(a)
(b)
(c)
(d)Figure3 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.Relative errors of the orbital elements for a Quasi-Spot satellite.
(a)
(b)
(c)
(d)
(e)
(f)


## 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 ε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:
(3.1)ϕ(B)(1-B)dε~t=θ(B)νt,
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:
(3.2)ϕ(B)=1-ϕ1B-ϕ2B2-⋯-ϕpBp,
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:
(3.3)θ(B)=1-θ1B-θ2B2-⋯-θqBq,
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:
(3.4)Φ(Bs)ϕ(B)(1-Bs)D(1-B)dε~t=Θ(Bs)θ(B)νt,
where
(3.5)Φ(Bs)=1-Φ1Bs-Φ2B2s-⋯-ΦPBPs,Θ(Bs)=1-Θ1Bs-Θ2B2s-⋯-ΘqBQs,
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 Figure4. 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.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 during10 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)εtx=xt-xtA,
where x represents each of the Delaunay variables (l,g,h,L,G,H), xt is the simulated data at epoch t, and xt𝒜 is the data from the analytical theory at the same epoch. Therefore, these six time series (εtl,εtg,εth,εtL,εtG,εtH) 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 J3 and J4), as well as the higher orders of the analytical solution, that is, the error of the analytical theory 𝒪(J22), 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εtx reveals that all variables show cyclical patterns or periodicities and, moreover, there is very similar behaviour between the time series εtl and εtg and εtL and εtG. For example, this study for εtl and εtg can be seen in Figure 5. However εth and εtH do not show any similar behaviour between them or with the rest of the time series.Periodograms and autocorrelation functions ofεtl (left) and εtg (right).
(a)
(b)
(c)
(d)
(e)
(f)On the other hand, the correlation matrix of theεtx series is shown in Table 1. This matrix presents a strong relationship between εtl and εtg, as their correlation coefficient is near -1 (-0.9607), as well as between their respective conjugate momenta time series errors, εtL and εtG, where their correlation coefficient is near 1 (0.9982).Table 1
Correlation matrix of theεtx series.
ε
t
l
ε
t
g
ε
t
h
ε
t
L
ε
t
G
ε
t
H
ε
t
l
1.0000
ε
t
g
−0.9607
1.0000
ε
t
h
−0.0575
−0.1713
1.0000
ε
t
L
−0.1037
0.0934
0.0451
1.0000
ε
t
G
−0.1128
0.1168
−0.0070
0.9982
1.0000
ε
t
H
0.4033
−0.4114
0.0420
0.0066
0.0005
1.0000It 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 (e=1-(G/L)2) are almost constant. Figure 6 shows εtl, εtg and εtL, εtG time series. As can be observed, εtl and εtg are almost symmetric with respect to the x-axis, which explains the negative sign and the near -1 value in the correlation coefficient, whilst εtL and εtG are almost the same, and therefore the sign in the correlation coefficient is positive with a near 1 value.Blue representsεtl, and εtL errors. Dashed red represents εtg and εtG errors.
(a)
(b)Finally, this preliminary study can be completed by analyzing the results obtained when theεtx time series are combined with the data obtained from the analytical theory (lt𝒜,gt𝒜,ht𝒜,Lt𝒜,Gt𝒜,Ht𝒜) during the first 10 satellite cycles. This test allows us to consider several possibilities. The first consists of considering each series separately, for instance, (lt𝒜+εtl,gt𝒜,ht𝒜,Lt𝒜,Gt𝒜,Ht𝒜). 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)
(lt𝒜+εtl,gt𝒜+εtg,ht𝒜,Lt𝒜,Gt𝒜,Ht𝒜),(ii)
(lt𝒜,gt𝒜,ht𝒜,Lt𝒜+εtL,Gt𝒜+εtG,Ht𝒜),(iii)
(l+εtl,gt𝒜+εtg,ht𝒜,Lt𝒜+εtL,Gt𝒜+εtG,Ht𝒜),(iv)
(lt𝒜,gt𝒜,ht𝒜+εth,Lt𝒜,Gt𝒜,Ht𝒜+εtH).Figure7 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 εtl and εtg can be seen in the first plot; the distance error is reduced to 0.63km after ten satellite cycles, whereas εtL and εtG only remove part of the short-period variations, as can be seen in the second plot. The third case collects the corrections due to εtl, εtg, εtL, and εtG, which produce a distance error similar to the first case. Finally, the corrections due to εth and εtH do not have any effect on the distance error, as can be observed in the last plot.Blue represents the distance error betweenxti and xti𝒜, whereas dashed red is the distance error between xti and the corrected xti𝒜 from the exact εtx.
(a)
(b)
(c)
(d)Next we focus our attention on carrying out a hybrid-AOPP from(lt𝒜+εtl,gt𝒜+εtg,ht𝒜,Lt𝒜,Gt𝒜,Ht𝒜). 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εtl and εtg

To estimate the model of theεtg 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 (f=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εtg. 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 εtg 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.
Coefficients
ε
t
l
ε
t
g
ϕ
1
−0.8167142581776827
−0.8312754102958296
ϕ
2
0.5740668630692614
0.5649950239467217
ϕ
3
1.1200190576112987
1.1322145773449677
ϕ
4
−0.3087269639691513
−0.2946598587646277
ϕ
5
−0.8631472597887673
−0.8663595825987899
ϕ
6
−0.4970128310547080
−0.5012076262623468
θ
1
0.8789635328932761
0.8901438163231377
θ
2
−0.9965719584754135
−0.9834775736328421
θ
3
−1.6921042908035564
−1.7044561386678760
θ
4
0.5206249360614155
0.4906661687944900
θ
5
1.5190020299259381
1.5136974115875366
θ
6
0.3606298858846529
0.3820758183271771
θ
7
−0.4721289483358399
−0.4568107050477830
Φ
1
1.4676902604680244
1.4557054355067451
Φ
2
−1.3618276394078148
−1.3525881149216927
Φ
3
0.7076454044934742
0.7023619159769664
Θ
1
0.8412748367694940
0.8370722835862818
Θ
2
−0.8240640881053386
−0.8282527898256336
Θ
3
0.9939014392344844
0.9966285211365780We must note that the model used to approximate theεtl time series is also a SARIMA(6,1,7)(3,1,3)10, which confirms the similar behaviour previously detected to εtg, 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.


## 4.1. Previous Statistical Analysis

The time series analysis begins calculating the linear relations:(4.1)εtx=xt-xtA,
where x represents each of the Delaunay variables (l,g,h,L,G,H), xt is the simulated data at epoch t, and xt𝒜 is the data from the analytical theory at the same epoch. Therefore, these six time series (εtl,εtg,εth,εtL,εtG,εtH) 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 J3 and J4), as well as the higher orders of the analytical solution, that is, the error of the analytical theory 𝒪(J22), 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εtx reveals that all variables show cyclical patterns or periodicities and, moreover, there is very similar behaviour between the time series εtl and εtg and εtL and εtG. For example, this study for εtl and εtg can be seen in Figure 5. However εth and εtH do not show any similar behaviour between them or with the rest of the time series.Periodograms and autocorrelation functions ofεtl (left) and εtg (right).
(a)
(b)
(c)
(d)
(e)
(f)On the other hand, the correlation matrix of theεtx series is shown in Table 1. This matrix presents a strong relationship between εtl and εtg, as their correlation coefficient is near -1 (-0.9607), as well as between their respective conjugate momenta time series errors, εtL and εtG, where their correlation coefficient is near 1 (0.9982).Table 1
Correlation matrix of theεtx series.
ε
t
l
ε
t
g
ε
t
h
ε
t
L
ε
t
G
ε
t
H
ε
t
l
1.0000
ε
t
g
−0.9607
1.0000
ε
t
h
−0.0575
−0.1713
1.0000
ε
t
L
−0.1037
0.0934
0.0451
1.0000
ε
t
G
−0.1128
0.1168
−0.0070
0.9982
1.0000
ε
t
H
0.4033
−0.4114
0.0420
0.0066
0.0005
1.0000It 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 (e=1-(G/L)2) are almost constant. Figure 6 shows εtl, εtg and εtL, εtG time series. As can be observed, εtl and εtg are almost symmetric with respect to the x-axis, which explains the negative sign and the near -1 value in the correlation coefficient, whilst εtL and εtG are almost the same, and therefore the sign in the correlation coefficient is positive with a near 1 value.Blue representsεtl, and εtL errors. Dashed red represents εtg and εtG errors.
(a)
(b)Finally, this preliminary study can be completed by analyzing the results obtained when theεtx time series are combined with the data obtained from the analytical theory (lt𝒜,gt𝒜,ht𝒜,Lt𝒜,Gt𝒜,Ht𝒜) during the first 10 satellite cycles. This test allows us to consider several possibilities. The first consists of considering each series separately, for instance, (lt𝒜+εtl,gt𝒜,ht𝒜,Lt𝒜,Gt𝒜,Ht𝒜). 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)
(lt𝒜+εtl,gt𝒜+εtg,ht𝒜,Lt𝒜,Gt𝒜,Ht𝒜),(ii)
(lt𝒜,gt𝒜,ht𝒜,Lt𝒜+εtL,Gt𝒜+εtG,Ht𝒜),(iii)
(l+εtl,gt𝒜+εtg,ht𝒜,Lt𝒜+εtL,Gt𝒜+εtG,Ht𝒜),(iv)
(lt𝒜,gt𝒜,ht𝒜+εth,Lt𝒜,Gt𝒜,Ht𝒜+εtH).Figure7 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 εtl and εtg can be seen in the first plot; the distance error is reduced to 0.63km after ten satellite cycles, whereas εtL and εtG only remove part of the short-period variations, as can be seen in the second plot. The third case collects the corrections due to εtl, εtg, εtL, and εtG, which produce a distance error similar to the first case. Finally, the corrections due to εth and εtH do not have any effect on the distance error, as can be observed in the last plot.Blue represents the distance error betweenxti and xti𝒜, whereas dashed red is the distance error between xti and the corrected xti𝒜 from the exact εtx.
(a)
(b)
(c)
(d)Next we focus our attention on carrying out a hybrid-AOPP from(lt𝒜+εtl,gt𝒜+εtg,ht𝒜,Lt𝒜,Gt𝒜,Ht𝒜). 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εtl and εtg

To estimate the model of theεtg 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 (f=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εtg. 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 εtg 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.
Coefficients
ε
t
l
ε
t
g
ϕ
1
−0.8167142581776827
−0.8312754102958296
ϕ
2
0.5740668630692614
0.5649950239467217
ϕ
3
1.1200190576112987
1.1322145773449677
ϕ
4
−0.3087269639691513
−0.2946598587646277
ϕ
5
−0.8631472597887673
−0.8663595825987899
ϕ
6
−0.4970128310547080
−0.5012076262623468
θ
1
0.8789635328932761
0.8901438163231377
θ
2
−0.9965719584754135
−0.9834775736328421
θ
3
−1.6921042908035564
−1.7044561386678760
θ
4
0.5206249360614155
0.4906661687944900
θ
5
1.5190020299259381
1.5136974115875366
θ
6
0.3606298858846529
0.3820758183271771
θ
7
−0.4721289483358399
−0.4568107050477830
Φ
1
1.4676902604680244
1.4557054355067451
Φ
2
−1.3618276394078148
−1.3525881149216927
Φ
3
0.7076454044934742
0.7023619159769664
Θ
1
0.8412748367694940
0.8370722835862818
Θ
2
−0.8240640881053386
−0.8282527898256336
Θ
3
0.9939014392344844
0.9966285211365780We must note that the model used to approximate theεtl time series is also a SARIMA(6,1,7)(3,1,3)10, which confirms the similar behaviour previously detected to εtg, 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

Figure8 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 (ε^tl,ε^tg) and the osculating Keplerian elements (a,e,g,h,i,l) and state vector (x,y,z,x˙,y˙,z˙) 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ω/n, where n 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 first10 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.Figure9 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.Z2DN1-SARIMA hybrid-AOPP. Relative errors of the mean anomaly and argument of the perigee for a Quasi-Spot satellite.
(a)
(b)Figure10 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.Z2DN1-SARIMA hybrid-AOPP. Distance, along-track, cross-track, and radial errors for a Quasi-Spot satellite.
(a)
(b)
(c)
(d)


## 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 first10 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.Figure9 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.Z2DN1-SARIMA hybrid-AOPP. Relative errors of the mean anomaly and argument of the perigee for a Quasi-Spot satellite.
(a)
(b)Figure10 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.Z2DN1-SARIMA hybrid-AOPP. Distance, along-track, cross-track, and radial errors for a Quasi-Spot satellite.
(a)
(b)
(c)
(d)


## 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 εtx 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.



---

*Source: 207381-2012-08-26.xml*