Forecasting satellite trajectories by interpolating hybrid orbit propagators

TL;DR

This paper introduces a method for forecasting satellite trajectories by interpolating parameters of a hybrid orbit propagator, enhancing prediction accuracy for nearby initial conditions using Holt-Winters interpolation.

Contribution

It demonstrates that interpolating hybrid propagator parameters can produce accurate trajectory forecasts for nearby initial conditions, reducing the need for developing new propagators.

Findings

Interpolation yields similar accuracy to non-interpolated methods for J2 effect modeling.

Parameter interpolation effectively extends hybrid propagator applicability.

The approach simplifies trajectory prediction for similar initial conditions.

Abstract

A hybrid orbit propagator based on the analytical integration of the Kepler problem is designed to determine the future position and velocity of any orbiter, usually an artificial satellite or space debris fragment, in two steps: an initial approximation generated by means of an integration method, followed by a forecast of its error, determined by a prediction technique that models and reproduces the missing dynamics. In this study we analyze the effect of slightly changing the initial conditions for which a hybrid propagator was developed. We explore the possibility of generating a new hybrid propagator from others previously developed for nearby initial conditions. We find that the interpolation of the parameters of the prediction technique, which in this case is an additive Holt-Winters method, yields similarly accurate results to a non-interpolated hybrid propagator when modeling the J2 effect in the main problem propagation.

Tables

Table 1. Table 1: Distance error (km) after propagating the studied satellite.

Propagation span	Analytic method	Hybrid method
	(Kepler)	(Kepler + $J_2$)
$1$ day	$1197.10$	$0.45$
$2$ days	$2379.94$	$0.83$
$7$ days	$7900.47$	$3.63$
30 days	$14504.69$	$13.73$

Table 2. Table 2: Position error (km) after the interpolated hybrid propagation of the intermediate initial conditions I e ​ 2 subscript 𝐼 𝑒 2 I_{e2} through different interpolation methods.

Propagation	Weighted	Linear	Lagrange	Spline
span	average	regression	polynomial
$1$ day	$9.070$	$1.772$	$2.021$	$0.469$
$2$ days	$18.440$	$3.563$	$4.070$	$0.836$
$7$ days	$66.222$	$12.460$	$14.237$	$3.498$
30 days	$272.777$	$48.653$	$55.715$	$13.598$

Table 3. Table 3: Position error (km) after propagating the intermediate initial conditions I e ​ 2 subscript 𝐼 𝑒 2 I_{e2} .

Propagation	Analytic method	Hybrid method	Spline-interpolated hybrid method
span	(Kepler)	(Kepler + $J_2$)	(Kepler + $J_2$)
$1$ day	$1244.061$	$0.600$	$0.469$
$2$ days	$2472.668$	$0.840$	$0.836$
$7$ days	$8165.104$	$3.711$	$3.498$
30 days	$14504.581$	$14.982$	$13.598$

Table 4. Table 4: Position error (km) after propagating the intermediate initial conditions I 4 ​ i subscript 𝐼 4 𝑖 I_{4i} .

Propagation	Analytic method	Hybrid method	Spline-interpolated hybrid method
span	(Kepler)	(Kepler + $J_2$)	(Kepler + $J_2$)
$1$ day	$1232.193$	$0.533$	$0.451$
$2$ days	$2449.259$	$0.819$	$0.821$
$7$ days	$8098.784$	$3.828$	$3.601$
30 days	$14510.092$	$15.090$	$13.691$

Table 5. Table 5: Position error (km) after propagating the intermediate initial conditions I e ​ i subscript 𝐼 𝑒 𝑖 I_{ei} .

Propagation	Analytic method	Hybrid method	Spline-interpolated hybrid method
span	(Kepler)	(Kepler + $J_2$)	(Kepler + $J_2$)
$1$ day	$1240.220$	$0.561$	$0.411$
$2$ days	$2465.050$	$0.823$	$0.711$
$7$ days	$8143.626$	$3.724$	$2.936$
30 days	$14506.802$	$14.950$	$11.619$

Equations

$${\mathbf{x}}_f^{\mathcal{I}}=\mathcal{I}(t_f,{\mathbf{x}}_1).$$

$${\mathbf{\epsilon }}_i={\mathbf{x}}_i-{\mathbf{x}}_i^{\mathcal{I}},$$

$${\hat{\mathbf{x}}}_f={\mathbf{x}}_f^{\mathcal{I}}+{\hat{\mathbf{\epsilon }}}_f.$$

$${\epsilon }_t={\mu }_t+S_t+{\nu }_t.$$

$${\mu }_t=A+Bt.$$

$${\hat{\epsilon }}_t=A_{t-1}+B_{t-1}+S_{t-s},$$

$$A_t$$

$$B_t$$

$$S_t$$

Full text

1 Introduction

The propagation of perturbed orbits is a well-known problem which implies having to tackle a set of three second-order or six first-order differential equations, so as to determine the position and velocity of an orbiter at a given final time $t_{f}$ from its situation at an initial instant $t_{1}$.

As these equations are not directly integrable, there are three well-established techniques aimed at providing a solution to the problem. Each of these methods can be characterized in terms of the formulation of the equation of motion, the integration method used to obtain the solution to this equation, which can be numerical or analytical, and, finally, the perturbation model taken into account.

General perturbation theories apply perturbation methods to the determination of an analytical solution. Such solution, which is an explicit function of time and some physical constants, allows for a fast determination of the coordinates at $t_{f}$. In addition, being an analytical expression, it embeds the dynamics of the problem. Nevertheless, in order to avoid extreme complexity, the analytical solution is usually a low-order approximation in which only the most relevant forces are considered.

Special perturbation theories, in contrast, perform a numerical integration of the problem. They have the advantage of allowing for the consideration of any effect into the model, even the complex ones, thus leading to highly accurate solutions. Nonetheless, the disadvantage lies in the necessity to take small integration steps, which implies long computational time.

Semianalytical techniques take advantage of both theories. They allow for the consideration of complex perturbing effects into the model, which is simplified by means of analytical methods so as to remove the short-period component. Consequently, the new equations of motion can be numerically integrated through longer steps, resulting in reduced computational time.

More recently, the hybrid propagation methodology has been presented. It is based on the combination of any of the aforementioned integration methods and a forecasting technique. The former generates an initial solution, which is approximate because of the assumed simplifications and inaccuracies in the perturbation models. The latter makes use of forecasting techniques, based on either statistical time series models [7, 8] or machine learning methods [5], in order to provide, once adjusted with a set of real observations that include the dynamics neglected in the initial approximation, a prediction of its error. The sum of this error prediction and the initial solution generates the final result.

The forecasting component of a hybrid propagator needs a set of control data, deduced from precise observations or accurately computed coordinates, so that the statistical or machine learning technique can model dynamics not present in the first stage of the method.

Nevertheless, a grid of hybrid propagators for a set of relatively close initial conditions can be constructed, in such a way that hybrid propagators for intermediate cases can be directly deduced from the grid, with no need for control data. By doing so, the study of initial conditions in the surroundings of an orbiter can be easily handled with no need to recompute the parameters of the forecasting component of the hybrid method.

In this paper we will consider the so-called main problem of the artificial satellite theory, that is, the Kepler problem only perturbed by the flattening of the Earth. We will create a hybrid propagator, composed of a general perturbation theory derived from the Kepler problem plus an additive Holt-Winters method modeling the $J_{2}$ effect, for a certain orbiter. In order to handle both eccentricity and inclination small variations, we will construct a grid of hybrid propagators around the studied satellite. After that, we will prove that the forecasting component of the hybrid propagator, when eccentricity and/or inclination slightly vary, can be directly derived from the grid by simply interpolating the parameters of the Holt-Winters method.

The outline of the paper is divided into seven sections. Section 2 presents the principles of the hybrid propagation methodology, whereas Section 3 focuses on the use of an exponential smoothing technique, the Holt-Winters method, as the forecasting stage of hybrid propagators. The described concepts are applied to the creation of a hybrid propagator for a certain satellite in Section 4. With the aim of studying its surroundings, a grid of initial conditions, together with its corresponding hybrid propagators, is created around the studied satellite in Section 5. Section 6 illustrates how to develop new hybrid propagators for nearby initial conditions through the interpolation of parameters from other propagators in the grid. Finally, Section 7 summarizes the conclusions of the study and future lines of research.

2 Hybrid propagation methodology

The hybrid propagation methodology is aimed at estimating the position and velocity of an orbiter, which can be an artificial satellite or a fragment of space debris, at a given final time $t_{f}$, $\hat{\bm{x}}_{f}$, starting from the known position and velocity at an initial instant $t_{1}$, $\bm{x}_{1}$. It is worth noting that any set of canonical or non-canonical variables can be used for this purpose.

In a first stage, an integration method $\mathcal{I}$ is used to calculate an initial approximation of $\hat{\bm{x}}_{f}$:

[TABLE]

The integration method is applied to a mathematical model that not always describes the physical phenomena exactly. Moreover, when the general perturbation theory or semianalytical techniques are used, only the most important forces and low-order approximations are usually considered; otherwise cumbersome expressions would be obtained. Due to all these facts, $\bm{x}_{f}^{\mathcal{I}}$ is an initial approximation that needs to be complemented in a second stage in order to obtain $\hat{\bm{x}}_{f}$.

The information that this second stage needs to model and reproduce, that is, the missing dynamics, has to be deduced from a control interval $[t_{1},t_{T}]$, with $t_{T}<t_{f}$. Throughout this interval both the initial approximation $\bm{x}_{i}^{\mathcal{I}}$ and the exact position and velocity $\bm{x}_{i}$ are assumed to be known, for example by means of precise observations or intensive and accurate numerical propagations. Therefore, the error due to the missing dynamics for any instant in the control interval can be expressed as

[TABLE]

and the time series of the errors of each of the six variables during the control interval, which we will call control data, can be constructed as $\mathcal{\bm{\varepsilon}}_{1},\ldots,\mathcal{\bm{\varepsilon}}_{T}$.

The processing of this time series, by means of either statistical techniques or machine learning methods, allows for the modeling of its behavior and, more importantly, its prediction at any time outside the control interval. Therefore, an estimation of the error at the final instant $t_{f}$, $\hat{\mathcal{\bm{\varepsilon}}}_{f}$, can be determined, and thus the desired value of $\hat{\bm{x}}_{f}$ can be calculated as

[TABLE]

3 Exponential smoothing method for time series forecasting

Exponential smoothing methods consider a time series $\varepsilon_{t}$ as the combination of three components: the trend $\mu_{t}$, or secular variation, the seasonal component $S_{t}$, or periodic oscillation, and the irregular or non-predictable component $\nu_{t}$. In the case of an additive composition, $\varepsilon_{t}$ can be expressed as

[TABLE]

In particular, the Holt-Winters method [10] considers a linear trend with level $A$ and slope $B$:

[TABLE]

According to this method, and taking into account that $\nu_{t}$ cannot be predicted, the next value of a time series can be estimated from past values as

[TABLE]

where $s$ is the period of the seasonal component, and $A$, $B$, and $S$ can be determined from previous values according to the following recurrences

[TABLE]

in which $\alpha$, $\beta$, and $\gamma$ are three smoothing parameters with values in the interval $[0,1]$.

Algorithm 1 shows how to apply the Holt-Winters method to the prediction of future time series values. The inputs to the algorithm are the amount of data per revolution, $s$, the number of revolutions in the control interval, $c$, the number of time steps after the control interval for which the time series value has to be predicted, $h$, and the control data, $\{\varepsilon_{t}\}_{t=1}^{T}$, with $T=s\times c$. The output is $\hat{\varepsilon}_{T+h|T}$, that is, the forecast of the time series at the final instant $t_{f}=t_{T+h}$, based on the last control data, $\varepsilon_{T}$.

The algorithm starts by estimating the initial parameters $A_{0}$, $B_{0}$, $S_{-s+1},\ldots$, $S_{-1}$, and $S_{0}$, which is accomplished through a classical additive decomposition into trend and seasonal variation over the three first revolutions. A linear regression over the trend provides the initial level $A_{0}$ and slope $B_{0}$, whereas the seasonal component yields the values of $S_{-s+1},\ldots$, $S_{-1}$, and $S_{0}$.

Then, an iterative process takes place by applying Eqs. (6) and (7) to the control interval (lines 2–7). As a result, the expressions of the parameters $A_{t}$, $B_{t}$, $S_{t}$, and the single-step time series prediction $\hat{\varepsilon}_{t}$ are obtained as functions of the smoothing parameters $\alpha$, $\beta$, and $\gamma$.

After that, an error measure is selected among mean square error, MSE, mean absolute error, MAE, and mean absolute percentage error, MAPE.

The selected error measure applied to the control interval yields an expression which is a function of the smoothing parameters. Then, an optimization method is necessary to determine the values of the smoothing parameters that minimize this error measure. The limited memory algorithm L-BFGS-B [4], which is a variation of the BFGS method [9], allows to impose restrictions on the smoothing parameters, and hence is the algorithm that has been used.

Once the optimal smoothing parameters have been found, the time series parameters $A_{T}$, $B_{T}$, $S_{T-s+1},\ldots,S_{T-1},S_{T}$ are determined for the last period of the control interval, from which the forecasted time series value at the final instant, that is, $h$ epochs ahead, $\hat{\varepsilon}_{f}=\hat{\varepsilon}_{T+h|T}$, can be calculated (line 11).

4 Application of the hybrid propagation methodology

In this section, the described hybrid methodology is applied to the propagation of an orbit with the following initial conditions: semi-major axis $a=7228$ km, eccentricity $e=0.06$, and inclination $i=49^{\circ}$. The first stage of the method is an analytical theory derived from the Kepler problem, that is, considering no perturbations at all, whereas the second part is an additive Holt-Winters method, designed to model the perturbation caused by the flattening of the Earth, which corresponds to the $J_{2}$ term in its gravitational potential. Therefore, the complete hybrid propagator is adapted to the main problem of the artificial satellite theory; consequently, its results will be compared with those obtained from a highly accurate numerical integration of the main problem by means of a high-order Runge-Kutta method.

The solution to the Kepler problem provided by the analytical expression in the first stage of the hybrid propagator is characterized by constant values in all the classical orbital elements, except in the mean anomaly, whose values evolve following the orbiter angular position. In contrast, when the $J_{2}$ effect is considered, no orbital element remains constant, so that, in general, secular, short-period, and long-period effects can be found in the evolution of orbital elements. The goal of the Holt-Winters method in the second stage of the hybrid propagator is the modeling and reproduction of such dynamics. The difference between the initial Kepler solution and the desired main problem solution translates into a position error of about 14500 km after 20 days of propagation, which represents approximately the distance between the apogee and perigee of the orbit.

The hybrid methodology can be applied to any set of variables, although Delaunay variables $(l,g,h,L,G,H)$ will be used in this case. The first step consists in preparing the control data, which is composed of two time series: the initial approximations generated by the analytical expression derived from the Kepler problem, and the accurate solutions calculated by means of a high-order Runge-Kutta method. The last time series could be substituted for a set of precise observations in case they were available. The subtraction of both data sets yields the time series of the error, which contains the dynamics missing from the initial approximation. It is worth noting that the control data set should be large enough so as to include any pattern to be modeled by the second stage of the method. In this case, a control interval of ten revolutions has been chosen, which represents a time span of nearly 17 hours, taking into account that the aforementioned orbital elements correspond to an orbital period of 101.926 minutes. The sampling rate for the time series has been taken equal to 12 samples per orbiter revolution, which corresponds to a sampling period of $101.926/12=8.494$ minutes.

Before processing data, angular variables are homogenized to the interval $(-\pi,\pi]$ by adding or subtracting complete spins to values outside this interval. An univariate Holt-Winters model is considered for the time series of the error of each Delaunay variable, except for $\varepsilon_{t}^{H}$, which is 0 in this case, which means that the analytical approximation is perfect for this variable, and hence there is no need to complement it in the second stage of the hybrid propagator.

Then, a preliminary analysis of the five remaining time series is performed through the study of their sequence graphics, periodograms, and autocorrelation functions (ACF). This analysis reveals the existence of three main seasonal components, with periods a third, a half, and one Keplerian period, that is, 33.976, 50.964, and 101.926 minutes, although the last one is the most remarkable and includes the others.

Next, Algorithm 1 is applied, selecting MSE as the error measure needed to determine the optimal values for the smoothing parameters $\alpha$, $\beta$, and $\gamma$.

Once the five Holt-Winters models corresponding to Delaunay variables $l$, $g$, $h$, $L$, and $G$ have been created, they are integrated into the hybrid propagator so as to evaluate its accuracy through the comparison with a precise numerical propagation by means of a high-order Runge-Kutta method.

Table 1 compares the position error, after different propagation spans, between the analytical approximation, which only considers the Kepler problem, and the hybrid propagation, which models the main problem. As can be seen, the latter presents reduced errors, even after 30 days of propagation, which implies that the forecasting part of the hybrid method has been able to model most of the $J_{2}$ effect.

5 Creation of a grid from control data

After developing a hybrid propagator for the studied satellite, the effect of a slight change in the initial conditions will be analyzed. For that purpose, small variations in eccentricity and inclination will be considered. We construct a grid of initial conditions around the studied satellite, modifying its eccentricity in $0.5\times 10^{-2}$ steps and its inclination in $1^{\circ}$ steps, as shown in Figure 1.

Next, we develop a new hybrid propagator for each initial condition $S_{ei}$ in the grid, following the steps described in the previous section for the studied satellite $S_{33}$. It is worth noting that control data are necessary for that process. Our final objective, in the next section, will be to verify the possibility to develop new hybrid propagators for initial conditions within the margins of the grid without having to follow the complete process, and hence with no need for control data.

We finish the creation of the grid hybrid propagators by analyzing their position errors with respect to the accurate numerical integration of their initial conditions. Figure 2 shows their distribution after different propagation spans by means of boxplot graphics. It can be seen that their average values agree with those shown in Table 1 for the studied satellite $S_{33}$.

In general, the distributions of the position errors are symmetrical, showing little dispersion and only a few outliers in the case of a 30-day propagation horizon. All the initial conditions have similar dynamic behavior, which leads to the homogeneity of the obtained position errors. Such results constitute an appropriate scenario for the adaptation of the developed hybrid propagators to nearby initial conditions.

6 Propagation of new orbits

At this point, hybrid propagators for the studied satellite $S_{33}$ and its surrounding grid $S_{ei}$ have been developed. Now, we want to propagate nearby initial conditions $I$ which occupy intermediate positions within the limits of the grid (Figure 1).

The analytical theory in the first stage of the hybrid propagators is the same for all the cases. However, each set of initial conditions requires an individual Holt-Winters model in the forecasting stage of its hybrid propagator, aimed at modeling and predicting the effect of the $J_{2}$ perturbation under its particular conditions. In order to take advantage of the nearby hybrid propagators developed in advance, and also to avoid the need for control data, a new strategy is proposed: the interpolation of the parameters $A_{T}$, $B_{T}$, $S_{T-s+1},\ldots,S_{T-1},S_{T}$ of the intermediate $I$ Holt-Winters models from those corresponding to the studied satellite $S_{33}$ and its surrounding grid $S_{ei}$.

Several interpolation methods have been compared. Some of them only allow for one-dimensional interpolation, while others permit multi-dimensional interpolation. We perform comparisons on $I_{e2}$, which only needs one-dimensional interpolation because only one of its elements, the eccentricity, differs from the values in the grid.

In the first place, a weighted average technique is used. We take the inverse of the difference in eccentricity as weight, and interpolate $I_{e2}$ Holt-Winters parameters from those of $S_{12}$, $S_{22}$, $S_{32}$, $S_{42}$, and $S_{52}$, which share the same inclination with $I_{e2}$. In the second place, the linear regression method is applied, deducing $I_{e2}$ parameters from the nearest straight lines to $S_{12}$, $S_{22}$, $S_{32}$, $S_{42}$, and $S_{52}$ parameters. The third interpolation approach is performed through Lagrange polynomials, by deducing $I_{e2}$ parameters from the fourth-order polynomials passing through $S_{12}$, $S_{22}$, $S_{32}$, $S_{42}$, and $S_{52}$ parameters. As it is known, the order of the Lagrange polynomials would increase if more initial conditions were available on the grid. Finally, spline interpolation is used. This is the only considered method that permits multi-dimensional interpolation. The two-dimensional spline interpolation implemented in the Akima package [3] of the R programming language [6], which will be the method to be applied to the case of $I_{ei}$ because both its eccentricity and inclination differ from all the initial conditions present on the grid, is based on References [1] and [2].

Table 2 presents the results obtained for each of the four aforementioned interpolation methods by means of the position error of the interpolated hybrid propagators developed for $I_{e2}$. As can be seen, spline interpolation leads to the best results for all the propagation spans, followed by linear regression, Lagrange polynomial, and, finally, the weighted average technique, which yields the worst results.

The analysis of these interpolation methods applied to the other set of intermediate initial conditions that requires one-dimensional interpolation, $I_{4i}$, yields the same conclusions; therefore spline is selected as the interpolation method to be used. Then, an interpolated hybrid propagator is also developed for $I_{ei}$, making use of two-dimensional spline interpolation, as mentioned previously.

Figure 3 represents the position errors obtained for the spline-interpolated hybrid propagation of the three intermediate initial conditions, and compares them with the distributions of the corresponding position errors for the hybrid propagation of the grid initial conditions (Figure 2). As could be expected, due to the homogeneous behavior of all the initial conditions in the grid, the position errors of the intermediate cases are very similar to the grid average. The case of the two-dimensional interpolation in eccentricity and inclination, $I_{ei}$, is remarkable because of its especially low errors, to the extent that it constitutes a low-error outlier for a propagation horizon of 30 days.

Tables 3, 4, and 5 compare the results of propagating the three sets of intermediate initial conditions $I_{e2}$, $I_{4i}$, and $I_{ei}$ through the mere analytic, the hybrid, and the spline-interpolated hybrid methods. In general, it can be verified that the latter propagators outperform the non-interpolated hybrid ones, especially in the case of the two-dimensionally spline-interpolated hybrid propagator for $I_{ei}$.

7 Conclusion and future work

In this work, we have presented an advance in the hybrid propagation methodology. Hybrid propagators are composed of an integration theory plus a forecasting technique. The latter is developed from control data so as to complement the approximation generated by the former by modeling and reproducing the missing dynamics. We have explored the possibility of deducing the forecasting stage directly from other hybrid propagators developed for surrounding initial conditions. This approach avoids the need for control data, and makes it possible to have a grid of hybrid propagators prepared in advance for a region of initial conditions of interest. We have verified that the spline interpolation of the parameters of an additive Holt-Winters forecasting method from nearby hybrid propagators yields similar accuracy, or even better, to a non-interpolated hybrid propagator. The study has been conducted using the main problem of the artificial satellite theory as the propagation model, with the forecasting stage modeling the complete $J_{2}$ effect.

At present, we are testing the hybrid propagation methodology considering neural networks instead of the Holt-Winters algorithm as time series forecasters. During the second semester of 2017, a competition organized by the European Space Agency will be launched through the Advanced Concepts Team competition website, Kelvins,111https://kelvins.esa.int/ in order to encourage the machine learning community to get involved and participate in the problem.

Acknowledgments

This work has been funded by the Spanish State Research Agency and the European Regional Development Fund under Project ESP2016-76585-R (AEI/ERDF, EU). Support from the European Space Agency through Project Ariadna Hybrid Propagation (ESA Contract No. 4000118548/16/NL/LF/as) is also acknowledged.