Refining satellite trajectories with celestial body features using neural networks

Contents lists available at ScienceDirect Expert Systems With Applications journal homepage: www.elsevier.com/locate/eswa Refining satellite trajectories with celestial body features using neural networks José Calderón a,c,∗, Daniel Ayala a,c, Rafael Ayala b, Luis Valencia-Cabrera c, Inma Hernández a,c, David Ruiz a,c aDEAL, ETSII, University of Seville, Av. Reina Mercedes, s/n, 41012, Seville, Spain bMolecular Cryo-Electron Microscopy Unit, OIST, 904-0411, Okinawa, Japan cSCORE Lab, I3US, University of Seville, Av. Reina Mercedes, s/n, 41012, Seville, Spain A R T I C L E I N F O Dataset link: SCD-ML Keywords: Satellites Orbit propagation Neural networks Features engineering A B S T R A C T Satellite orbit propagation involves predicting a satellite’s future position and velocity based on initial conditions. Traditional physical models, such as SGDP4, simplify the forces that act on the satellite to achieve high computational efficiency at the cost of reduced prediction accuracy, especially over longer time intervals where error accumulates. More sophisticated models like HPOP offer improved accuracy at the cost of high prediction times, rendering them unusable for realtime long-term predictions. Recent advancements have introduced machine learning techniques to refine these predictions and reduce errors. However, they often lack an analysis of model design choices, such as input feature selection and architectural configurations. Existing models do not incorporate features related to the state of celestial bodies, such as the positions of the Moon or Sun, which can influence the satellite’s trajectory. This paper proposes a novel model that integrates such features at both the initial time and throughout the prediction interval, leveraging their potential impact on the orbit of the satellite. The model is based on a neural network architecture employing GRU layers for encoding sequential data about the celestial conditions. Our results demonstrate that the inclusion of these sequential features significantly reduces prediction errors. Additionally, we have evaluated a variety of design choices such as independent sub-models for specific spatial coordinates and time intervals, further enhancing performance. These innovations lead to substantial improvements in both shortand long-term orbit predictions, providing a more robust and accurate alternative for satellite orbit propagation. 1. Introduction On October 4th, 1957, the Soviet Union launched the first artificial satellite, Sputnik-1. Since then, almost 20000 objects have been sent into space, of which more than 13000 are still in orbit, according to the UN (United Nations Office for Outer Space Affairs, 2024). This massive volume of objects makes it crucial to track them to ensure the efficiency, safety, and success of space missions and the proper functioning of satellites (Ayala, Ayala, Vidal, & Ruiz, 2023; Levit & Marshall, 2011; Refaat, Badawy, Ashry, & Adel, 2018). Orbit propagation consists in the prediction of the position and velocity of a satellite at some future date and time (known as epoch in the context of astronomy), given some initial conditions (Flores, Burhani, & Fantino, 2021). The trajectory of a satellite is influenced by various forces, mainly gravitational and centrifugal, as well as other perturbations such as the influence of the Sun, the Moon, and other ∗Corresponding author at: SCORE Lab, I3US, University of Seville, Av. Reina Mercedes, s/n, 41012, Seville, Spain. E-mail addresses: [email protected] (J. Calderón), [email protected] (D. Ayala), [email protected] (R. Ayala), [email protected] (L. Valencia-Cabrera), [email protected] (I. Hernández), [email protected] (D. Ruiz). celestial bodies (Shou, 2014). To make these predictions, two main types of models are used: analytical models and high-precision models. Analytical models, also known as simplified perturbation models, include models such as SGP4 (Simplified General Perturbations) and SDP4 (Simplified Deep-space Perturbations) (Hoods & Roehrich, 1988), which we collectively refer to as SGDP4. These models offer an approximate way to model the satellite orbit, allowing for very fast but less precise long-term predictions. They are useful for applications in which speed is crucial and accuracy can be sacrificed in favour of computational efficiency. High-precision models, such as HPOP (High Precision Orbit Propagator), use differential equations to model all the forces acting on the satellite. These models are extremely precise, but require high computing power, making them expensive and slow to apply, especially for long-term predictions, as can be seen in Table 1. https://doi.org/10.1016/j.eswa.2025.127453 Received 23 October 2024; Received in revised form 31 January 2025; Accepted 25 March 2025 Expert Systems With Applications 281 (2025) 127453 Available online 7 April 2025 0957-4174/© 2025 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/bync-nd/4.0/ ). J. Calderón et al. Table 1 Computational time for different prediction times with HPOP and SGDP4 propagators. Prediction Computational Computational time time time (hours) HPOP SGDP4 0,5 8,69 s 602.54 μs 1 13,57 s 617.03 μs 2 23,55 s 605.36 μs 5 34,93 s 637.44 μs 10 57,45 s 615.19 μs 24 4,09 min 608.23 μs 48 7,99 min 627.64 μs These models consider in detail gravitational perturbations, solar radiation pressure, atmospheric drag, among other factors, providing a more accurate prediction of the satellite orbit (Bate, Mueller, & White, 1971; Montenbruck, Gill, & Lutze, 2002; Vallado, 2001). In turn, they require the input of detailed information about the satellite’s properties, such as its mass, area, and others. In this context, Machine Learning emerges as a powerful tool to address this challenge: assist in providing a solution dramatically faster than HPOP, and with a significant improvement in accuracy with respect to SGDP4. Advances in Machine Learning algorithms and the availability of large volumes of data have the potential for the development of models that significantly improve some traditional forecasting methods (Lam et al., 2023; Xie, Yao, Li, Wang, Zheng, & Chen, 2024). Additionally, Machine Learning models can continuously improve as more data becomes available, ensuring that orbit propagation methods stay current and effective. To improve the accuracy and efficiency of orbit propagators, some authors have developed hybrid models, which combine a dynamic model with a machine learning method. For sequential satellite position data, approaches such as linear regression and Long Short-Term Memory (LSTM) neural networks (Ren et al., 2019), or combinations of autoencoders and random forests (Liu, Tarlow, Akbar, Donnellan, & Senkow, 2021), have been proposed. In cases where data lacks a sequential nature, alternative methods such as dense neural networks (San-Juana, Pérezb, Vergarac, San Martınd, Lópeze, Wittigf, & Izzog, 2018), support vector machines (Peng & Bai, 2017) or a combination of PCA and XGBoost (Zhai, Huyan, Hu, Jiang, & Li, 2022). Despite these advancements, they only test specific configurations, missing the chance to study the influence of certain architectures or data over the outcome of the models. Additionally, these works tend to focus solely on satellite-specific data, excluding the integration of external contextual information, such as environmental factors or external forces, which could further refine predictions. Among the different Machine Learning techniques, neural networks stand out for their ability to model complex non-linear relationships in data, making them particularly suitable for prediction and classification tasks. Neural networks have been shown to outperform other Machine Learning algorithms in terms of accuracy and efficiency in certain contexts (Ayala, Borrego, Hernández, & Ruiz, 2020; Ayala, Hernández, Ruiz, & Toro, 2019; LeCun, Bengio, & Hinton, 2015; Schmidhuber, 2015). For this reason, this work focuses on exploring how neural networks can be applied in conjunction with simplified perturbation models, developing a hybrid model with SGDP4 to optimize orbital prediction as illustrated in Fig. 1, reducing errors in predicted trajectories while keeping a low prediction time. Throughout this work, we present various models designed to address this specific task, evaluating the performance of each one compared to the reference SGDP4 model and analysing its impact on accuracy. Our most significant contribution is the study of the impact of including explicit features related to the state of the celestial bodies that may influence the satellite trajectory. Our working hypothesis is that integrating these novel features is crucial for the development of more effective prediction models. Fig. 1. Use of a model to refine orbit propagation. The output of the model corresponds to the green arrow. To achieve this goal, we created a dataset with data from Kosmos 2514, a satellite from the GLONASS constellation. Then, we developed and tested some neural network models with different input data, finding great error reduction in some of them. Our experiments show that the inclusion of features related to celestial bodies results in significant improvements, and is key when it comes to obtaining acceptable results in short-term predictions, in which the baseline prediction provided by the SGDP4 model is more accurate and therefore harder to improve or at least not to deteriorate. 2. Related work Several authors have already evaluated the outcomes of applying machine learning techniques to refine trajectory prediction. These studies highlight the potential of machine learning algorithms to enhance precision, showing improvements over only using traditional methods. Ren et al. (2019) proposed a combination of an LSTM neural network and linear regression. They use LSTM layers to predict orbital inclination, orbital eccentricity and average displacement, training one model to predict each variable. Besides, linear regression is used to estimate the ascending node, perigee angular distance and near point angle. After predicting the six elements of the TLE orbit, the position of the spacecraft can be obtained by a simple transformation that allows calculating the error. The inputs for the six models are the same: the six orbital elements, the timestamp to which these elements correspond and the time interval for the prediction. In these LSTM networks, each of the six orbital parameters is provided as a five-element sequence. To train these models, data from the IRIDIUM 118 satellite1 was used. The data was divided into training and testing sets, with 80% for training and 20% for testing. The authors consider the predicted position error to be within an acceptable range and its variation to stay relatively stable against the progress of time. Peng and Bai (2017) explored the viability of employing a Support Vector Machine (SVM) regression model to forecast the error in an assumed dynamical model of the ENVISAT satellite.2 Initially, they model three stations to generate discrete measurements, when the target satellite is visible to them, according to a ‘‘truth’’ dynamical model. Then they apply least squares estimation to get the state of the satellite. After obtaining estimations for all the tracks, the prediction process is straightforward, the SVM model is tasked with predicting the error in this forecast to refine the model estimates. The inputs for the SVM model are the duration of the prediction; position and velocity 1https://www.n2yo.com/satellite/?s=42807 (accessed 15 October 2024) 2https://www.n2yo.com/satellite/?s=27386 (accessed 15 October 2024) Expert Systems With Applications 281 (2025) 127453 2 J. Calderón et al. at the current epoch; estimated drag coefficient at the current epoch; maximal measured elevation in the current track and the corresponding range and azimuth; and predicted position and velocity. The output of the SVM is a six-component vector detailing the error in position and velocity across each axis. The conclusion drawn from their study was that the SVM model, once trained, had limited applicability for predictions extending too far into the future. Therefore, they advised that orbit predictions should be made within a relatively short timeframe to ensure accuracy. Peng and Bai (2018) later introduced a neural network approach, employing dense neural networks that used identical inputs from the same satellite as in the previously described SVM model. The outputs of the neural network mirrored those of the SVM, with a notable variation: the approach involved training a different model for each component, resulting in a total of six unique models. They studied also the ability of these models to generalize to future epochs and to different but nearby Resident Space Objects (RSOs) in subsequent epochs. In their analysis of this second task, the authors analysed the range of learning variables, excluding any variables whose ranges exhibited significant disparities between the training and testing datasets. Their findings suggested that the neural network models demonstrate strong generalization capabilities to future epochs. Furthermore, it was concluded that these neural networks could be generalized to a relatively broad spectrum of nearby RSOs not included in the training dataset, showcasing their versatility and potential for predictive accuracy in dynamic space environments. Liu et al. (2021) developed a hybrid model, combining SGP4 with autoencoders and random forests to reduce the propagation error. First, they used the encoder to obtain a representation of the distance error from SGP4. Then, a random forest model was applied to predict the embedding vector for the next time step and finally, a decoder obtained the distance error. This was done for each positional coordinate (𝑥, 𝑦, 𝑧) resulting in three different models. The input of the model was a 30-days time series of the SGDP4 distance error. The data they used come from three objects in low earth orbit (LEO): a research CubeSat (QuakeSat by Stanford University3), a satellite payload (COSMOS 20984), and a satellite debris (PEGASUS DEB debris5). These data were collected from the Space-Track API (SAIC, 2024). With this approach, they obtained an average 20%–30% improvement on 30-days orbit prediction. Zhai et al. (2022) proposed a combination of PCA and XGBoost model to improve the orbit prediction accuracy. The inputs of the model are prediction duration, predicted movement, position at the initial epoch, velocity at the epoch, drag coefficient, predicted position and predicted velocity. The data used for training and evaluation were sourced from satellite simulation environments. First, they trained an XGBoost model to choose the most appropriate combination of features. Then, based on that parameters, the PCA–XGBoost model was trained,in which PCA was used to reduce the dimensionality of the data. The target variable was the true error prediction, which consists of six elements, three axis for position and another three for velocity, such that a total of six models are trained. They concluded that this proposal improves satellite prediction in their simulation enviroment and claimed that its capability of generalization is good for all six components. Despite the previous approaches, the application of machine learning techniques to satellite propagation has several unexplored areas that present opportunities for significant advancements. The impact of some decisions regarding the development of models is unclear. For example, when Liu et al. use a model for each positional coordinate, it is unknown to what degree the separation of a model into separate coordinates contributes towards improving results. Other aspects suffer 3https://www.n2yo.com/satellite/?s=27845 (accessed 15 October 2024) 4https://www.n2yo.com/satellite/?s=20774 (accessed 15 October 2024) 5https://www.n2yo.com/satellite/?s=23975 (accessed 15 October 2024) from a similar lack of insights, such as the consequences of evaluating a refinement model across different time windows (short, medium, and long term) or even training models that focus on a single time window. But perhaps the most remarkable gap is the absence of models that incorporate features related to the state of celestial bodies. Features used by existing models are based on the description of the state of the satellite at the time of inference, namely its position and velocity, and features that are derived from these and the application of the SGDP4 model such as orbital elements. However, the accurate propagation of a satellite orbit relies on the position and speed of all surrounding celestial bodies, which are not included in existing models. These are also more challenging to properly exploit, since their position changes during propagation time, and therefore the way in which they impact the satellite position. Related to satellite trayectory predictions, we encounter time series processing. The literature commonly addresses the classification of multivariate time series (Xiao et al., 2023). In these cases, models learn discriminative representations to assign a label to an entire sequence of temporal data. However, our objective diverges from this approach. Rather than classifying time series into discrete labels, our goal is to predict continuous values at a specific future time step. Another dominant paradigm in the field of machine learning for time series is sequence prediction, as seen in traffic forecasting models (Wu, Pan, Long, Jiang, & Zhang, 2019) or meteorology (Pathak et al., 2022). In these cases, the goal is to infer the future values of the entire time series based on previous observations. A prominent example of such models is the use of Transformers (Vaswani, 2017), which have proven highly effective in capturing long-term dependencies within temporal sequences. However, in our case, we do not aim to predict the entire future evolution of the sequence but rather only the state of the satellite at a specific instant. In summary, while our study shares fundamental aspects with existing research on time series, such as leveraging historical data to infer future values, our objective is different. We neither seek to classify temporal sequences into fixed categories nor to model their entire future trajectory. Instead, our focus is on precise state estimation at a particular future instant. 3. Our proposal We propose a hybrid model that combines the SGDP4 propagator with a neural network to improve satellite propagation predictions while leveraging novel features. The model operates in two steps: 1. Initial Prediction: Let 𝐩SGDP4(𝑡) denote the position prediction of the satellite at time 𝑡 generated by the SGDP4 propagator. 2. Error Correction: The initial prediction 𝐩SGDP4(𝑡) is used as input for a neural network , which forecasts the error vector 𝐞(𝑡) made by the SGDP4 propagator. The corrected position prediction 𝐩corrected(𝑡) is then given by: 𝐩corrected(𝑡) = 𝐩SGDP4(𝑡) + 𝐞(𝑡)(1) By integrating the traditional SGDP4 propagator with the advanced capabilities of a neural network, our hybrid model aims to enhance the accuracy of satellite position predictions. This approach leverages the strengths of both methods: the established reliability of SGDP4 and the adaptive learning potential of neural networks. Neural networks allow us to feed our model a variety of features 𝐟 to study their effect on the outcome, and assess which variants offer the best results. In the following sections, we describe the aims and scopes of our work in order to outline its limitations, the specific goals of our experiments, the neural network architectures we implemented, and the design of the datasets we used in our evaluation. 3.1. Aims and scopes Our research is focused on the evaluation of neural network models when applied to the prediction of satellite trajectories, in combination Expert Systems With Applications 281 (2025) 127453 3 J. Calderón et al. Fig. 2. Architecture for the Densely Connected Network. with the traditional propagator SGDP4. Our aim is to compare the impact on performance of different strategies to assess their predictive potential. Particularly, we aim to study how the addition of features related to the state of celestial bodies impacts performance, since we consider that the advance of the application of machine learning in this field greatly relies on the development of novel ways to model relevant information as features. Our goal is not to replace traditional models and their physical simulations. Instead, we aim to refine and build upon the foundation provided by traditional models, leveraging their robustness and accuracy to introduce progressive improvements and insights into which strategies are most effective. We focus on the creation of models for a single satellite, as opposed to creating a model that, after being trained, can be used to infer the trajectory of any satellite. The latter would require modelling detailed information about the physical characteristics of each satellite, which is generally not available, in order to properly make predictions about different satellites (e.g. the section of the satellite used to measure the force of radiation used by the HPOP model). Our experiments correspond to a use case in which the user is able to gather data about the precise trajectory of a satellite in order to train a model and enable more precise prediction of its trajectory in the future. By focusing models on a single satellite, the model is implicitly adjusted to the specific conditions of that satellite, eliminating the need to know exactly all the parameters that influence its trajectory, thus simplifying the model development process. We do not intend to propose a definitive and closed solution for trajectory propagation, but rather to validate the effectiveness of certain techniques in a controlled context. The first features set used in our experimentation corresponds to the features used by the techniques in the State of the Art, serving as a baseline representing the information typically used in existing proposals, namely the position and velocity of the satellite. The conclusions derived from our results are of independent interest to any technique. Likewise, we do not seek an exhaustive evaluation of complex neural network architectures or extensive combinations of hyperparameters. Instead, we focus on using simple and effective architectures that are sufficient to exploit the information of interest. 3.2. Goals Our objective is to assess the effect of considering different architectures with different input data on neural network performance for satellite propagation. In particular, our experiments aim to assess the impact of the following variables: 1. Base Features: The base features used for the prediction of the trajectory, including the position 𝐩(𝑡) and velocity 𝐯(𝑡) of the satellite, as well as the SGDP4 prediction 𝐩SGDP4(𝑡). 2. Celestial Bodies at Initial Epoch: The inclusion of features related to celestial bodies at the initial epoch 𝑡0, and the inclusion of celestial bodies beyond the most influential ones. 3. Sequence of Celestial Bodies: The inclusion of features related to celestial bodies in a sequence of positions from the initial epoch 𝑡0 to the inference epoch 𝑡. Specifically, we feed the models a sequence of fixed length 10, corresponding to evenly distributed time points between 𝑡0 and 𝑡. 4. Separate Models for Coordinates: The creation of a separate model for each predicted coordinate, as opposed to a model that predicts all coordinates simultaneously, replacing  with 𝑥, 𝑦, and 𝑧 for the 𝑥, 𝑦, and 𝑧 coordinates, respectively. 5. Time Frame Specific Models: The creation of a model for single time frames, as opposed to a generic model applicable to any time frame, replacing  with 0.5ℎ, 1ℎ, 2ℎ and so on. 3.3. Network architecture Next, we describe the neural network architectures used in our models. In particular, we implemented three different networks: a simple one for processing static features, a second one with GRU (Gated Recurrent Unit) layers (Chung, Gulcehre, Cho, & Bengio, 2014) and a third architecture combining a Transformer encoder with GRU layers to process the sequences of features related to celestial bodies. The first model is a densely connected network (DCN), also known as a fully connected network, which primarily uses dense layers. These layers connect every input feature to every output feature with a unique weight for each connection. They are particularly useful for tasks in which the features do not have an inherent order or spatial relationships between them, as is the case with most of the features we use, such as the prediction provided by the SGDP4 model, or the satellite position and velocity at the initial epoch. This model, depicted in Fig. 2, contains 7 intermediary dense layers of sizes 128, 128, 256, 256, 128, 128, 64. Each layer is followed by a ReLU activation layer, and a skip connection is added after every layer with the same output size as the previous one. Finally, the output layer has size 3 corresponding to the three predicted coordinates, or size 1 for the cases in which an independent model is trained for each coordinate. The use of residual skip connections was essential, specially in models that received sequential data. This change in the architecture led to big improvements in the prediction performance, reducing the error in 80% approximately. The second model adds GRU layers. This model uses as its main component one of said layers, which are designed to consume sequences of data and produce an output for each element in the sequence. The output corresponding to the last element of the sequence will be influenced by the entire sequence, making GRU layers ideal to process sequential data (Cho et al., 2014). GRU networks have been applied to all kinds of problems in which the input contains some kind of sequence, such as electrocardiogram classification (Lynn, Pan, & Kim, 2019) or short-term power load forecasting (Zheng et al., 2018). Compared to other recurrent neural networks, GRU networks are simpler and faster to train thanks to their inner simplicity, having potential for faster convergence in cases where the input sequence is not very long. In our model, the input sequence has a fixed length of Expert Systems With Applications 281 (2025) 127453 4 J. Calderón et al. Fig. 3. Architecture for the GRU network. 10, corresponding to the 10 evenly distributed time points between the initial and inference epochs. Informal tests confirmed that the performance of GRU layers was slightly better than that of popular LSTM layers. As shown in Fig. 3, our GRU network handles two types of inputs: sequential and non-sequential data. For the sequential data, we define an input layer with a shape that matches the sequence length (10) and the number of sequential features. This input is processed by a GRU layer with 64 units using the tanh activation function. For the nonsequential data, we define another input layer with a shape matching the number of non-sequential features. This input is processed by a dense layer with 64 units. The outputs from the GRU layer and the dense layer are concatenated to form a combined input. This combined input is then passed through a series of layers equivalent to layers from the first model, starting at layers of size 258. The third model incorporates a Transformer encoder, a mechanism known for its success in handling sequential data across various domains. Transformers rely on self-attention mechanisms to capture relationships between all elements in a sequence, regardless of their distance, overcoming the limitations of recurrence-based models like GRU or LSTM, which struggle with long-range dependencies. Transformer encoders, specifically, focus on encoding input sequences into a rich, context-aware representation by attending to every element in the sequence and weighing their relevance (Vaswani, 2017). Our Transformer-Based model, shown in Fig. 4, handles two types of inputs: sequential and non-sequential data. For the sequential data, we define an encoder to process sequences of length 10. First, we concatenate positional encoding, which vary according to the position in the sequence, and celestial body encoding, which are trainable and vary according to the celestial body. After that, we incorporate multihead attention for processing the sequence with the encodings. This attention mechanism uses 3 heads with an embedding dimension of 16, helping the network capture the spatial and contextual relationships in the data. The output of the attention mechanism is followed by a skip connection with the same output size as its input and a normalization layer. After that comes two dense layers with 13 and 27 neurons and ReLU activations, respectively. The output of this second layer is followed by a skip connection with the input of the first dense layer and a normalization layer. After the encoder, the architecture aligns with that of a GRU network. However, the number of units in the GRU layer has been increased to 256. 3.4. Computational complexity 3.4.1. Densely connected network Let us denote by 𝐿 the total number of dense layers and by 𝑛𝑖 the number of units in layer 𝑖, with 𝑛0 being the number of input features. The core operation in each layer is the matrix multiplication between the input of size 𝑛𝑖−1 and the layer’s trainable weights of size 𝑛𝑖−1 ×𝑛𝑖. Hence, the forward pass for layer 𝑖 costs (𝑛𝑖𝑛𝑖−1). Summing across all dense layers yields: (𝐿 ∑ 𝑖=1 𝑛𝑖𝑛𝑖−1). Residual or skip connections primarily add element-wise operations, which remain (𝑛𝑖) per skip link and do not alter the total complexity. The memory footprint is determined by the trainable parameters inside these dense layers, which total ∑𝐿 𝑖=1(𝑛𝑖𝑛𝑖−1 +𝑛𝑖) including biases. Therefore, the dense architecture is well-suited when the input dimensionality is moderate, as is the case. 3.4.2. GRU-based model In this architecture, let 𝑇 be the sequence length, 𝑑 the dimension of each time step’s input, and ℎ the number of hidden units per GRU layer. A GRU layer processes each time step sequentially, mixing the current input with the previous hidden state. Each step requires multiple matrix multiplications of size 𝑑×ℎ and ℎ×ℎ, resulting in a per-step complexity of ((𝑑+ℎ)ℎ). Over 𝑇 time steps, this yields: (𝑇(𝑑+ℎ)ℎ). The following dense layers add an additional cost of (∑𝐿 𝑖=1 𝑛𝑖𝑛𝑖−1). Thus, the computational complexity remains proportional to 𝑇, 𝑑, and ℎ, and can be expressed as: (𝑇(𝑑+ℎ)ℎ+ 𝐿 ∑ 𝑖=1 𝑛𝑖𝑛𝑖−1). In terms of parameter count, a single GRU layer has on the order of 3ℎ(ℎ+𝑑) trainable parameters (including biases), because each gate in the GRU cell has its own set of weights. Overall, for short sequences or smaller hidden dimensions, concatenating the GRU output with non-sequential dense components adds moderate overhead. In our implementation, 𝑑= 3, 𝑇= 10, and ℎ= 64, which should result in reasonable training times. 3.4.3. Transformer-based model In this architecture, a classical Transformer encoder is applied before the GRU layer. Let 𝐸 be the dimension of each time step’s input after the addition of positional and celestial body embeddings, 𝐻 the number of attention heads, and 𝑑 the desired size for the query, key, and value projections. The addition of embeddings has a complexity of (𝑇 𝐸). The Transformer encoder consists of multi-head self-attention mechanisms and feed-forward neural networks. The self-attention mechanism involves projecting the input features into query, key, and value Expert Systems With Applications 281 (2025) 127453 5 J. Calderón et al. Fig. 4. Architecture for the Transformer-Based Network. vectors, each with a complexity of (𝑇 𝐸𝑑). The self-attention mechanism itself has a complexity of (𝑇2𝑑) per head, and with 𝐻 heads, this becomes (𝐻𝑇 2𝑑). The feed-forward network adds a complexity of (𝑇 𝑑2). The total complexity of the Transformer encoder layer is: (𝑇 𝐸𝑑 +𝐻𝑇 2𝑑+𝑇 𝑑2). After the Transformer encoder, the GRU layer processes the sequence with a complexity of (𝑇(𝑑+ℎ)ℎ). The following dense layers add an additional cost of (∑𝐿 𝑖=1 𝑛𝑖𝑛𝑖−1). Therefore, the overall computational complexity of the Transformerbased model is: (𝑇(𝐸𝑑 +𝐻𝑇 𝑑 +𝑑2+ (𝑑+ℎ)ℎ) + 𝐿 ∑ 𝑖=1 𝑛𝑖𝑛𝑖−1) The Transformer encoder has on the order of 𝐻𝐸𝑑 parameters for the attention mechanism and 2𝑑2 parameters for the feed-forward network per layer. The GRU layer has 3ℎ(ℎ+𝑑) trainable parameters, including biases. In our implementation, 𝑑= 16, 𝑇= 10, ℎ= 256, 𝐸= 3 + 16 + 8 = 27, and 𝐻= 3. These values were chosen according to informal tests in which we found that they provided the best performance. 3.5. Dataset design For the creation of the dataset we used data from the RINEX repository of the Institute of Statistics and Cartography of Andalusia (IECA) (Instituto de Estadística y Cartografía de Andalucía, 2024), specifically from the Córdoba station (CRDB). The data consisted of RINEX (Receiver Independent Exchange Format) files from GPS and GLONASS satellites covering the period from 2019 to 2022. RINEX files are a standard data format used for storing and exchanging GNSS (Global Navigation Satellite System) data. Developed to facilitate the exchange of GNSS data between different receiver manufacturers and processing software, RINEX files provide positions and velocities of satellites. For the treatment of the data and the application of orbital propagators we used asteRisk (Ayala et al., 2023), an R package for computation of satellite positions. This package provides tools for precise orbital mechanics calculations, ensuring accurate determination of satellite trajectories. The Kosmos 2514 satellite from the GLONASS constellation provided the most extensive data (number of messages detailing the position and velocity at a given epoch). Using this information, we created an auxiliary dataset, which contained positions, velocities and dates for each available instant. We selected a list of prediction time windows for our study, including 30 min, 1 h, 2 h, 5 h, 10 h, 24 h, 48 h, 120 h, 240 h, and 720 h. These intervals represented the various times at which we aimed to make predictions and train and test our model, covering a wide range of magnitudes ranging from half an hour to 30 days. To create the final dataset, we started by selecting random dates from the auxiliary dataset. For each chosen date and each prediction time, we generated a new row in the dataset. These new rows include: the selected initial epoch, the true position and velocity of the satellite at that epoch, the true position and velocity at the prediction epoch, and the position and velocity predicted by the SGDP4 method for that prediction epoch. For our SGDP4 prediction, we first needed to convert the positions and velocities from RINEX files into the TLE (Two-Line Element Set) format Vallado and Cefola (2012), which is required by the SGDP4 model. The TLE format is a standardized format used to represent orbital elements of satellites. In addition to the orbital elements, we also had to provide the 𝐵∗ drag coefficient. The 𝐵∗ coefficient is a parameter that encapsulates the effects of atmospheric drag on the satellite orbit. However, calculating the 𝐵∗ drag coefficient requires specific data about the shape of the satellite and its mass that is not available. Given this limitation and the fact that the 𝐵∗ value is typically very small, we decided to set it to 0 for all of our examples. This simplification allowed us to proceed with the predictions using the available data. During this process, there were some cases in which there were no satellite data available for certain prediction epochs. To address this Expert Systems With Applications 281 (2025) 127453 6 J. Calderón et al. Fig. 5. Dataset samples collection. Black dots represent epochs at which positional information is available. Arrows represent samples with different desired prediction times. If there is no data for a desired prediction time, the nearest point can be used instead. Table 2 Number of samples in the dataset for each time band. Prediction time Samples 0.5 h 3929 1 h 3557 2 h 3231 5 h 2399 10 h 3405 24 h 3890 48 h 3918 120 h 3914 240 h 3917 720 h 3839 issue, we implemented a threshold to select a time close to the intended prediction time. For instance, if data were not available for a prediction time of 5 h, we allowed the selection of data from a nearby point, such as 5.5 h. This approach ensured that the prediction epoch was near the intended time windows. This is exemplified in Fig. 5. After creating the dataset, we incorporated data from various celestial bodies within our solar system. This dataset includes the positions and velocities of the Sun, the Moon, Mercury, Venus, Mars, Jupiter, Saturn, Uranus, Neptune, and Pluto, all referenced relative to Earth. Additionally, we gathered information on lunar libration, which describes the oscillating motion of the Moon as observed from Earth. These celestial bodies were chosen by virtue of being the major gravitational forces that should influence the trajectory of a satellite. Note that other celestial bodies, such as asteroids or comets, were not included in our dataset. This decision was made due to the negligible effect that the gravitational influence exerted by such objects has in the trajectory of artificial satellites in the timescales considered in this study. We were able to acquire these data through the JPL DE440 model (Park, Folkner, Williams, & Boggs, 2021), developed by NASA. The JPL DE440 is a highly advanced model that enables quick and accurate calculations of the positions and velocities of celestial bodies, as well as other relevant astronomical data. We added celestial data to our dataset in two different versions. In the first version, we incorporated the data at the initial epoch of each row. In the second version, we included the celestial data in sequences of length 10, evenly distributed between the initial and prediction epochs. This captures the progression of celestial conditions over the entire prediction period, during which the planets exert gravitational force on the satellite. The final dataset has a total of 36000 samples, distributed as shown in Table 2. 4. Experimental results For training all the models in our experiments, the hold-out technique was used for data splitting. The data was divided equitably, with 50% assigned to the training set and the remaining 50% to the test set. The Adam optimizer was selected to dynamically adjust the learning rates and the Mean Squared Error (MSE) loss function was used to evaluate the model during training. Additionally, a callback was implemented to automatically save the best model based on the Mean Absolute Error (MAE) metric, ensuring the selection of the most accurate model during the training process. Mathematically, the MSE loss function is defined as: MSE =1 𝑛 𝑛 ∑ 𝑖=1 (𝐲𝑖− 𝐲𝑖)2(2) where 𝐲𝑖 is the true value and  𝐲𝑖 is the predicted value. The MAE metric is defined as: MAE =1 𝑛 𝑛 ∑ 𝑖=1 |𝐲𝑖− 𝐲𝑖|(3) To obtain a reliable measure of model performance, each model underwent 10 training sessions, each with a different random seed, ranging from seed 667 to 676. The mean error was calculated across these 10 sessions, taking into account the SGDP4 correction in all models. All our experiments were conducted on a computer equipped with an Intel Core i9-9900K CPU, 64 GB of DDR4 RAM and an Nvidia RTX 3080-Ti GPU. Their code is available on https://github.com/DEALUS/satellite-propagation-refinement When we indicate that a model receives inputs like position, displacement, or velocity, we refer to the three coordinates or components of the corresponding variable (𝑋, 𝑌, and 𝑍). This approach is similarly applied to other magnitudes, such as lunar libration. Dates are provided in UNIX format for consistency and ease of processing. All models in the experiment receive the following common input features: the initial epoch 𝑡0 , the position 𝐩(𝑡0) of the satellite at that instant, and the epoch 𝑡 on which the prediction is to be made. Thus, these inputs are referred to as the base features. We compare the different models against the SGDP4 baseline, since our purpose is to improve the results obtained by said model. The full results of our experiments can be found at Appendix, in which the average of each model we developed is included. Next, we present the results of our experiments in a progressive manner, starting with the simplest models and gradually introducing more complex models with additional features. This allows us to identify the impact of different features and additions, serving as an ablation study to understand the importance of each component. The experimentation began with an experiment on the DCN (Densely Connected Network) models, evaluating the impact of incorporating velocity 𝐯(𝑡) as part of the input data and whether it was more effective to use the prediction of the SGDP4 model as an absolute position 𝐩SGDP4(𝑡) or as the relative movement 𝛥𝐩SGDP4(𝑡) = 𝐩SGDP4(𝑡) − 𝐩(𝑡0). The training parameters configuration was 100 epochs and a batch size of 128. We also tested the performance of three existing techniques in the State of the Art, particularly those by Peng and Bai (2017), Peng and Bai (2018), and Zhai et al. (2022). Fig. 6 shows the results of this experiment. It is from the 48-hour predictions onwards that some models start to outperform SGDP4, with the improvement being very clear in the 720-hour (30-day) time window. This represents a significant advance for long-term predictions, which should be easier for the model, since at that point predictions by the SGDP4 model deviate largely. At smaller time windows, it is hard for the model to generalize, resulting in a distorted prediction with a much higher error. The inclusion of the SGDP4 prediction as an absolute position 𝐩SGDP4(𝑡) results in minor improvement at some time windows (2 h, and 24 h onwards), while the inclusion of the initial velocity 𝐯(𝑡0) seems to be irrelevant. However, in subsequent experiments using additional features, the use of the predicted movement 𝛥𝐩SGDP4(𝑡) resulted in a Expert Systems With Applications 281 (2025) 127453 7 J. Calderón et al. Fig. 6. Spatial error in base models and state-of-the-art proposals. ‘‘mov’’ denotes the inclusion as a feature of the SGDP4 predicted position as a differential of position (𝛥𝐩SGDP4(𝑡)). ‘‘pos’’ denotes the inclusion of the former as an absolute position (𝐩SGDP4(𝑡)). ‘‘vel’’ denotes the inclusion of the velocity of the satellite at the initial epoch (𝐯(𝑡0)). ‘‘svm’’ denotes the technique by Peng and Bai (2017). ‘‘ann’’ denotes the technique by Peng and Bai (2018).‘‘pca + xgboost’’ denotes the technique by Zhai et al. (2022). Fig. 7. Spatial error in celestial data models. ‘‘smj’’ denotes the inclusion as features of the position of the Sun, Moon, and Jupiter at the initial epoch. ‘‘svmmj’’ denotes the inclusion of the position of the Sun, Venus, Moon, Mars, and Jupyter. ‘‘planets’’ denotes the inclusion of the position of all Solar System planets (including Pluto). ‘‘lib’’ denotes the inclusion of the lunar libration. slightly better average. Therefore, from now on we report on the results obtained with the ‘‘mov’’ features. Out of the techniques in the literature, Peng and Bai (2017) obtains the best results, but still achieves worse performance than our models. Next, we perform the same experiments, but adding as features different combinations of the position of the bodies in the solar system and the lunar libration, maintaining the training configuration. Fig. 7 shows the results of this experiment. The inclusion of these features leads to a significant improvement of the results. Now, we can observe some improvement starting at the 24 h time windows, as well as better overall results. However, short-term predictions are still worse than the baseline. Fig. 8. Spatial error in sequential celestial data models. ‘‘gru(...)’’ denotes the inclusion of a set of features as a sequence that is fed to GRU layers of the network. Fig. 9. Spatial error in separate coordinate models. ‘‘(XYZ)’’ denotes the use of three independent models for the prediction of the three target coordinates. It is noteworthy that the inclusion of features about additional celestial bodies does not result in better results. On the contrary, the simplest set of features (‘‘smj + mov’’) yields slightly better results at shorter time windows. Subsequently, Residual GRU was implemented, to include the same features while expanding those related to celestial bodies to a sequence of ten points. The number of epochs was adjusted to 300 in order to account for the increased complexity of the network. Fig. 8 shows the results of this experiment. While the features were the same, their inclusion as an expanded sequence results in significant improvements. In the short-term, said improvement allows the distortion introduced by the neural network to be relatively minimal, while in the long-term there is a substantial improvement. The SGDP4 model is now improved at the 10 h time window. Interestingly enough, the use of features involving additional celestial models (‘‘svmmj’’) now results in a slightly better median. In the next experiment, we test a variant of the models in which coordinate separation was performed by dividing the predictions into Expert Systems With Applications 281 (2025) 127453 8 J. Calderón et al. Fig. 10. Example of evolution of spatial error for a single starting position. three distinct components, resulting in the creation of three models instead of one: 𝑥, 𝑦, and 𝑧. This process was carried out for the models that showed the best performance, which were all of the Residual GRU architecture. At first glance, the models may perform similarly. However, the loss from one coordinate might negatively impact the others. This can lead to worse overall performance, as errors in predicting one coordinate can affect the predictions of the others. Therefore, we test if training a separate model for each coordinate, despite being more costly, improves the prediction results. By having individual models, each one can focus on the specific characteristics and patterns of its coordinate without interference from others. Fig. 9 shows the results of this experiment. The improvement in predictions is present at every prediction time. This indicates that separating coordinates has a positive impact on both short-term and long-term prediction accuracy. Note that this addition leads to the first models that consistently do not worsen the prediction of the SGDP4 model at short-term, but achieve an equivalent or better median at every time window. Fig. 10 shows an example in which the satellite of an orbit is propagated from a random single point, displaying the spatial error in time of the SGDP4 model and the correction with the best neural network refinement. The model is, however, still not of much use for predictions below the 10 h time window. It is our hypothesis that the general training leads the model to focus on the longer prediction times that cause a higher loss. Therefore, we explore training models for each time window that may lead to better performance compared to a general approach. Based on the success of the incorporation of a GRU layer we implemented a Transformer-Based architecture, incorporating an encoder and maintaining the same features as in the last model. The number of epochs was adjusted to 200. Fig. 11 shows the results of this experiment. Although the features remained unchanged, the incorporation of the encoder significantly increased both the training and the execution time. However, this modification did not enhance the predictive performance compared to the previous GRU model. The architecture showing the best accuracy was gru(svmmj) + mov + vel (XYZ). Based on this architecture, 10 models were trained corresponding to each time window. In this approach, for a specific prediction time, the nearest specialized model is used, such as applying the 10-hour window model for an 8-hour prediction. Fig. 12 shows the results of this experiment. The difference of performance is dramatic at the short-term, where it has achieved better results across all time windows. A lower improvement can be observed at the 2 h time window, which may be caused by some cyclical astronomic factor. In the medium and long term, the new model also outperforms the SGDP4 baseline. However, it is important to note that it does not surpass its equivalent general version, probably because of the presence of more training data when not restricting the model to a single time window. Regarding computational efficiency, the training of the best model, gru(svmmj) + mov + vel (XYZ), took approximately 27 min, and the specialized models took around 2.3 min each. Considering that the model actually requires the training of three separate models, this is a reasonable time that would allow the user to train with a bigger dataset or to explore model hyperparameters if more complexity is needed for the model. If the user is interested in further reducing the training time, apart from traditional generic techniques, a lower number of points for the sequence of celestial bodies could be used. 5. Conclusions and future work We have presented a thorough study on the application of neural networks to satellite propagation, specifically to refine the predictions provided by the SGDP4 model. While existing literature contains a variety of machine learning techniques applied to this domain, our research distinguishes itself through a number of novel contributions that allow the reduction of the error in the propagation by an order of magnitude while. An exhaustive evaluation of different variants has enabled us to properly study the contribution of several design choices, such as the creation of distinct models for each coordinate. Fig. 11. Spatial error in Transformer-Based models. Expert Systems With Applications 281 (2025) 127453 9