A Causal Long Short-Term Memory Sequence to Sequence Model for TEC Prediction Using GNSS Observations

: The necessity of predicting the spatio-temporal phenomenon of ionospheric variability is closely related to the requirement of many users to be able to obtain high accuracy positioning with low cost equipment. The Precise Point Positioning (PPP) technique is highly accepted by the scientiﬁc community as a means for providing high level of position accuracy from a single receiver. However, its main drawback is the long convergence time to achieve centimeter-level accuracy in positioning. Hereby, we propose a deep learning-based approach for ionospheric modeling. This method exploits the advantages of Long Short-Term Memory (LSTM) Recurrent Neural Networks (RNN) for timeseries modeling and predicts the total electron content per satellite from a speciﬁc station by making use of a causal, supervised deep learning method. The scope of the proposed method is to compare and evaluate the between-satellites ionospheric delay estimation, and to aggregate the Total Electron Content (TEC) outcomes per-satellite into a single solution over the station, thus constructing regional TEC models, in an attempt to replace Global Ionospheric Maps (GIM) data. The evaluation of our proposed recurrent method for the prediction of vertical total electron content (VTEC) values is compared against the traditional Autoregressive (AR) and the Autoregressive Moving Average (ARMA) methods, per satellite. The proposed model achieves error lower than 1.5 TECU which is slightly better than the accuracy of the current GIM products which is currently about 2.0–3.0 TECU.


Introduction
In the past few years, among the Global Navigation Satellite Systems (GNSS) community, the Precise Point Positioning (PPP) technique has generated strong interest and is being used for a broad number of applications today, mainly due to its capability to provide low cost and very precise positions using a single receiver, with errors as small as a few centimeters under suitable conditions. Compared with other conventional GNSS positioning methods, it is recognized that PPP still has some disadvantages, with the major weakness being that it takes considerable time (up to tens of minutes) to converge to full (cm-level) positioning accuracy. Recently, various approaches based on combining single-or multi-constellation, single-frequency GNSS measurements with global ionospheric map (GIM) data have been proposed as one way to accelerate the PPP convergence time in real-or near-time applications. GIM data provide the necessary external total electron content (TEC) information required to overcome most of the error imposed by the ionosphere. In order to ensure single frequency users. Additionally, the International Reference Ionosphere (IRI) model, is a global climatological model for the terrestrial ionosphere which is recommended for international use by the Committee On Space Research (COSPAR) and the International Union of Radio Science (URSI) [11].
Apart from the above physical and empirical models which are being used widely in many applications to obtain TEC measurements, information about the ionosphere can also be extracted from GNSS observations through the use of linear combinations between multi-frequency observables [12][13][14] from single-and/or multi-constellations. Based on such computational procedures, the International GNSS Service (IGS) provides on regular basis Global Ionosperic Maps (GIM) in the IONospheric EXchange (IONEX) format.
The GIM models are constructed by the Ionosphere Associate Analysis Centers (IAACS) using TEC data derived from a few hundred GNSS stations worldwide. The overall effort is coordinated by the IGS Working Group on Ionosphere which was established in 1998 with the major task to provide routine products based on a combination of ionosphere maps regularly produced by four IGS Ionosphere Associate Analysis Centers (IAACs): the Centre for Orbit Determination in Europe (CODE), Jet Propulsion Laboratory (JPL), European Space Agency (ESA), and Polytechnic University of Catalonia (UPC). These centers use different approaches and techniques to compute global ionospheric models [15] and through averaging and interpolation techniques to produce Global Ionospheric Maps (GIMs) as part of IGS precise GNSS products. These maps are produced every 2 h, providing TEC values on a (latitude, longitude) grid spacing 2.5 • × 5 • , with the accuracy of TEC values published in IONEX format varying between 2-8 TECU. The Global Ionospheric Maps can provide better results than the Klobuchar model [16], although the resolution is generally not sufficient to model ionospheric delay at a regional (and even more so on local) scale.
It is worth mentioning that recently, in an effort to complement the ionospheric data provided by the GIM maps, relevant studies have turned their attention to obtaining improved global ionospheric parameters based on data fusion between multi-GNSS, Satellite Altimetry and Very Long Baseline Interferometry (VBLI) measurements [17]. This is possible because, VLBI allows probing the ionosphere in an absolute sense, whereas satellite altimetry can fill unprobed gaps of ionosphere sensing especially over the oceans (i.e., in areas lacking GNSS data).
The common characteristic of all of the aforementioned models are global representations of TEC values, with their accuracy being limited by various factors such as, for instance: the adoption of a simplified model in the case of the Klobuchar model; considering that the NeQuick model provides typical median condition of the ionosphere (so, it cannot handle successfully extreme ionospheric conditions) [18], or the non-uniform accuracy of IGS TEC data in the global scope for the GIM maps [19]. Furthermore, most of these models are stationary, difficult to be adapted to the irregular and complex distribution of TEC in the ionosphere, and additionally, to learn complex abstractions. Nowadays, many significant research efforts are spent exploting suitable TEC forecasting methods, including Support Vector Machine (SVM) [20], the nonlinear radial basic function (RBF) neural network [21,22], artificial neural network [18,19], and LSTM [23]. In particular, because of their capability of learning from noisy and nonlinear relationships, Neural Networks are considered an interesting approach in the domain of ionosphere TEC prediction [24]. Recurrent neural networks (RNNs) can handle with computational efficiency large amounts of TEC observations and can be adapted in order to learn the temporal dependencies from context [25]. LSTMs learn long term correlations in a sequence [8], they can model complex sequences with various features, as it is shown in [8], they can model complex sequences with various features, as it is shown in [26], where solar radio flux at 10.7 cm and magnetic activity index are taken into consideration to provide more accurate results. However, the above mentioned studies still focus on TEC prediction based on estimates derived from TEC models, such as GIM maps, and they don't try to predict ionosphere delays separately from every different visible GNSS satellite, as in our approach.

Our Contribution
In this work, we create regional TEC models with autoregressive character, using deep learning methods. Our intention, as a next step, is to extend the applicability of this approach towards the implementation of this model in near-real time PPP-RTK processing in order to use it as part of the integer ambiguity resolution (IAR)-enabled PPP owing to the use of predicted ionospheric delays in addition to other corrections. Until now, CODE provides GIM maps in a grid of (2.5 • × 5.0 • ), in 2-h temporal resolution, leading to a spatio-temporal sparse model. This means that the corrections applied improve the accuracy of the provided solution for positioning, but fail to remove the total amount of noise caused by the ionosphere. Stations or roving users near the reference area in which regional LSTM models have been created, could use corrections from the regional model and not a global one.
In order to accurately create the regional TEC model, RNN is employed as the formalism for solving the ionosphere variations prediction problem at hand. RNN is a class of neural network that allows previous outputs to be used as inputs, has hidden states and also, uses weights that are shared across time. Long Short-Term Memory network (LSTM) is a special variant of the typical RNN structure that deals with the vanishing gradient problem encountered by traditional RNNs. Adopting RNN method for ionosphere variations prediction is a suitable approach, as this is a time-sequence prediction problem. Recently, encouraging results have been obtained by employing LSTM-based configurations in domains such as language modelling, speech recognition and action recognition [27][28][29]. As already mentioned, ionosphere modeling is a problem of high complexity, as TEC values are space-and timevarying. Our model exploits LSTM [30] model's non-linearity, which is learned from the data. On the contrary, traditional time series prediction models, such as AR and ARMA [6], are particularly simple structures which are essentially linear update models plus some noise thrown in.
Our model handles TEC values at a micro-scale level and examines how consistent the extracted TEC values per observed satellite are, and whether the remaining biases affect them. Here, we clarify that DCBs have been removed using P1-P2.dcb and .bsx products. Therefore, our analysis examines the ability of a GNSS data-based method to provide accurate estimates of ionosphere variations in the likely presence of unaccounted noisy conditions. Then, an ensemble solution at station-level is extracted with proper combination of the information about TEC variations that provided from the analysis of each satellite results.
The contribution of this study is to explore the potential of LSTM in time series analysis for TEC modeling and prediction. Compared with traditional time series prediction models, LSTM can obtain better prediction results, by taking advantage of the following features of the method:

•
All-in-one modeling: Usually, existing models for timeseries prediction pertain to univariate timeseries, thus predicting TEC timeseries per satellite implies using a number of models equal to the number of satellites being tracked. On the contrary, our model uses at the same time all the satellites, thus forming a generic model. • Sequential modeling: Most existing models use a fixed data duration as input, whereas the TEC sequences per satellite, are not all of the same length, as this depends on the time duration that a satellite is visible from a ground station. LSTMs can successfully deal with data sequences of varying lengths. • Long term modeling: LSTMs with their internal memory, remember previous information that reflects the past behavior of the ionosphere and find patterns across time for accurate prediction of the next guess-estimates, making them ideal for time series forecasting.
The remainder of this paper is structured as follows. Section 2 deals with the formulation of the TEC prediction problem and describes the proposed LSTM model architecture. Section 3 describes the adopted processing strategy in order to extract meaningful features and the necessary TEC values that are used for training of the supervised algorithms. In Section 4, experimental results and analyses are provided, together with various comparisons, as for instance: (i) between traditional AR and ARMA methods, as well as, (ii) with well known TEC models. Conclusions are given in the Section 5.

Problem Formulation
We use the undifferenced and uncombined observations in dual-frequency PPP processing to extract STEC values. The code P and phase φ observations in a given frequency band f i between a receiver r and a GNSS satellite s, are written as [31]: where P s denote pseudorange and carrier phase observables, respectively; ρ s r is the geometric distance between the receiver to the satellite; c is the speed of light in vacuum; dt r and dt s are the receiver and satellite clock offsets, respectively; d ( f i) r is the frequency-dependent receiver uncalibrated code delay (UCD) while d s f i is the frequency dependent satellite UCD(in seconds); T is troposphere delay; I ( f i) is the line of sight (LOS) (slant) ionospheric delay on the frequency f i; δ ( f i) r and δ s f i are the frequency-dependent receiver and satellite uncalibrated phase delay, respectively (in seconds); N s ( f i) is the phase ambiguity; s are the sum of measurement noise and multi-path error for pseudorange and carrier phase observations. Equation (1) for dual-frequency GPS observations is: since the code biases are commonly generated and distributed as differential code biases (DCBs): Given γ 2 = f 2 1 / f 2 2 , we have [32]: Based on Equation (3), the term I 1 is grouped with DCBs, thus [13]: The difference between the uncombined PPP (UPPP) model and the traditional ionosphere free (IF) PPP model, lies in the way ionosphere delays are treated. The UPPP model estimates the ionosphere delays as unknown parameters, while the IF model eliminates the ionospheric delays by combination of dual frequency observations. However, when estimating ionospheric total electron content (TEC) using pseudorange and phase GNSS data, careful consideration should be given to the DCBs, as 1ns in a DCB causes an ∼ 2.9 TECU error in TEC estimation. After correcting the DCBs, then, based on the single layer model assumption, the STEC can be converted into the vertical total electron content VTEC as follows: where MF I is the mapping function: and R e is the mean radius of the Earth; θ is the elevation angle of the satellite; h is the height of the ionospheric layer (typically taken at 350 km). The TEC values can be easily estimated from GNSS data through the combination of VTEC values from individual satellites (distinguished by their respective PRNs). The typical ionospheric prediction model assumes that the ionosphere is a thin shell above the Earth, located near the mean altitude of maximum TEC (approximately 350 km).
The intersection between a line of sight and this shell is called the Ionospheric Pierce Point (IPP). At different epochs, at every IPP point, each individual satellite provides a different VTEC measurement ( Figure 1). A commonly followed approach is to model these VTEC values from GNSS measurements using mathematical functions such as bi-linear models [33][34][35][36], to represent the spatial and temporal variability of the VTEC. In all these models, the vtec parameter is related to time t as well as the coordinates (φ, λ) of each IPP point. Thus, a generalized form of [33], is: The parameters of the function f (·) in equation [11] were estimated in previous works for each station, using a Kalman filtering approach. In our approach, we are using a deep learning recurrent neural network. In the single layer model, the latitude IPP and the longitude IPP is related to azimuth and elevation angles [37], thus these two variables have been used as input features to our model. IPP points are also related to receiver's coordinates, which we don't include in our model, as for each network, these coordinates are constant, so they will not improve the model. However, in order to maintain the information regarding the line-of-sight direction, we include as features satellites coordinates, aiming at: (i) distinguishing each different satellite, (ii) providing time information indirectly, as at every epoch satellite coordinates change. Given a satellite's position vector (x SatX , x SatY , x SatZ ), and its azimuth x Az and elevation x El angles, we infer its direction whose intersection with the ionosphere shell is the IPP point ( Figure 1). Then, for each satellite, at every epoch t, the vertical TEC is functionally expressed as: where f (·) is a nonlinear relationship modeled by a learning process (Section 2.2) and e(t) is the noise factor. The task in hand is to compute the timewise vetrical TEC value from every visible satellite, denoted by vtec s (t), in an effort to provide accurate TEC measurements, under the assumption that at a close by-region, for a given time epoch, VTEC estimates from individual satellites should be similar. Then, the various satellite-specific VTEC estimates are combined to a single spatio-temporal solution, which is compared to the solution of the GIM maps.

Unidirectional Long-Term Recurrent Regression
Our proposed deep learning model is based on timeseries modelling. The inputs during training are not independent but they follow a time dependent order. In contrast to other deep learning schemes, mostly applied in classification problems, in which the common procedure is to randomize the entire input data, here we form our data in a time-dependent sequence to train our network. Thus, the proposed network is specially designed for time-series modelling and infers successfully the time-related information. Assuming L hidden neurons and one linear output layer, the estimate vtec s every epoch, is given by [38] vtec s (t) = u T (t) · v (13) tanh(·) that refers to the hyperbolic tangent, so each element of the vector u(t) ranges between −1 and 1. The weights w i , i = 1, . . . , L, connect the input x(t) with the i-th hidden neuron. Similarly, v are the weights that connect the hidden neurons with the output neuron. The term u(t) gathers the outputs of all L hidden neurons u i with values ranging between −1 and 1.
Since the ionosphere delay follows a causal relationship, the value of a state depends on its previous values, forming a short-term recurrent neural network structure. Therefore, An LSTM network is a variation of the classic RNN model, enriched with slowly changing weights to address the vanishing gradient problem that RNN suffer from. LSTM networks have internal mechanisms called gates that can regulate the flow of information, maintain the worth-remembering pieces of information and forget the unnecessary ones. LSTMs are of similar structure with the recurrent regression model, but each node in the hidden layer is replaced by a memory cell, instead of a single neuron [5]. It processes the available data passing on information as it propagates forward.
The memory cell ( Figure 2) contains three different components: (i) the forget gate, (ii) the input gate, and (iii) the output gate. Each component applies a non-linear relation to the inner product between the input vectors and respective weights (estimated through a training process). Some of the components have the sigmoid function, expressed as σ(·), while others have the hyperbolic tangent function, tanh(·).
Adding a "forget" mechanism. The model adopts a separate forgetting/remembering mechanism; this part of the LSTM cell memory is responsible for deciding which part of the information should be omitted.
Adding a "save" mechanism. Primarily, with the advent of a new input, the memory cell mechanism forgets any worthless information from previous terms, to subsequently learn the remaining part of the input entering the model.
Focusing long-term memory into working memory. Finally, forget gates can learn to reset the internal state of the memory cell when the stored information is no longer needed.
The forget gate f (t) separates the worth-remembering information from the unnecessary one, by keeping the latter out of the memory cell [39]. Similarly, the input node g(t) activates appropriately the respective state (true or false output from the tanh activation). The input gate i(t) regulates whether the respective hidden state is "significant enough" for the accurate estimation of vertical TEC. The output gate o(t) regulates whether the response of the current memory cell is "significant enough" to contribute to the next cell. It should be clarified that the nomenclature "Long term" used for LSTM networks, is a deep learning related terminology and is associated with the model's ability to learn from previous historic past values, and should not be confused with ionosphere's long term temporal variations that implies an analysis of many years. In our approach we are trying to model daily variations of the ionosphere and not long term ones. Figure 3 shows the adopted strategy and the basic structure of the proposed LSTM model. The network configuration setup consists of: (i) the input layer, (ii) two stacked unidirectional LSTM layers and (iii) the regression layer. In our case, we have a sequence-to-sequence model, in which data are organized as sequences of input and output pairs. This implies that LSTM model infers the relative time information during training, as it is not trained for every pair (input features, output) at every epoch t individually, but is trained to recognize a whole sequence/set of these pairs during the day, for a specific satellite. Thus, the input layer receives the current data in a sequence-like form, of a time window with varying duration. The network includes two bidirectional LSTM layers: the first layer consists of 60 filters with a kernel size 1 × 5, while the second layer consists of 72 filters of the same kernel size. The regression layer consists of one fully connected hidden layer and the final output regression layer. The output is a sequence of VTEC predicted values for each satellite, with size that covers the whole duration where the satellite is visible from the receiver.  In order to evaluate the ionospheric delay effects on GPS positioning, days between 10 and 13 December 2018 (day of year (DOY), from 344 to 347), where 80% of them has selected for training and the remaining 20% as validation set. Furthermore, 10 days have been selected as test set (14 December and 19 to 27 December 2018). Thus, we feed our network with the sequence-pairs ( Figure 3):

PPP Processing Strategy
This section describes the pre-processing strategy to extract the necessary inputs and outputs for our LSTM model, using observational data from a small group of permanent GNSS stations in the IGS network. The evaluation of the LSTM method for vectical TEC values prediction has been conducted on various ground stations of the global network of International GNSS Service (IGS). The time resolution is selected to be 30s. Table 1 shows the geographical distribution of the six selected IGS stations [40], which are spread across central Europe and are located in close proximity to each other. We used the GAMP software [31], a secondary development software based on RTKLIB [41] for multi-constellation and multi-frequency Precise positioning. GAMP allows the use of undifferenced and uncombined observations in dual-frequency PPP processing to extract STEC values. Using the uncombined PPP (UPPP) model allows estimating the ionospheric effects as unknown parameters, without the need to impose ionospheric-free constraints imposed by the traditional PPP model which amplifies the measurement noise when dual frequency observations are combined to cancel the ionospheric effects. The necessary RINEX observation and navigation files, along with precise orbit and clock information, IGS ANTEX (igs14.atx) and SINEX files, as well as ocean tide loading coefficients and DCBs are entered into the GAMP software for static PPP processing. Wuhan University's satellite orbits and clock offsets were provided as input for GAMP software. Differential code biases (DCBs), are essential in many navigation and non-navigation applications (such as ionospheric analysis). As new signals are currently provided by the modernized GNSS systems, the need for a comprehensive multi-GNSS DCB product arises. DCB products, as part of the IGS Multi-GNSS Experiment (MGEX), are provided by the Chinese Academy of Sciences (CAS) in Wuhan (Table 2).

Performance Comparison among Different Methods
We compare the proposed LSTM method against traditional stochastic models for timeseries modeling and forecasting, such as the autoregressive AR model and the autoregressive moving average ARMA model. Our LSTM model has been trained and deployed using Python with Tensorflow and Keras libraries. We trained the model using the adaptive moment estimation optimization algorithm (adam) with a learning rate of 10 −4 . Model weights and coefficients are updated using a mini-batch size of 28 samples at each training iteration. The maximum number of epochs for training is selected to be 600. Deep learning performance is improved through data balance and normalization of the data used as input to the LSTM network to the range [0, 1]. We have individual networks, one per station. AR and ARMA models were also implemented in Python using the stats models library. Table 3 shows performance metrics between LSTM model, as well as ARMA and AR models. For each satellite s i in every station r, a sequence of VTEC values is produced using the TEC model, for a given time period T. For each individual satellite, we compute the absolute difference between the ground truth vtec  Table 3. Performance metrics (mean, minimum (Min), maximum (Max) and root mean squared error (rmse)) for VTEC (TECU) prediction for six selected stations (bor1, ganp, graz, leij, pots and wtzz) among the proposed solution (LSTM) and two other commonly used methods (ARMA, AR). For every station r, we have also computed two additional metrics for every individual PRN (s i -satellite): the mean absolute error MAE s i r (i.e., the average sum of all absolute errors) and the root mean squared error RMSE s i r . However, Table 3 summarizes an overall accuracy for each station, using the following performance metrics: (i) mae r is the average value of the combined MAE s i r for all s i satellite estimations per station, (ii) min r is the average of all s i minimum difference values per station, (iii) in accordance to min r metric, max r is the average of all s i maximum differences values per station, and (iv) rmse r is the average of s i root mean squared error of all PRNs (satellites) per station. LSTM, ARMA ans AR are causal models, thus they are appropriate for ionospheric time series modeling. However, our proposed unidirectional LSTM model performs best (Table 3) mainly due to its capability to effectively identify structure and pattern of TEC data, such as non-linearity and complexity in time series prediction and infer TEC behavior based on previous states.

Receiver
LSTM model is undoubtedly, more complicated and difficult to train and in most cases does not exceed the performance of a simple ARMA model. Table 4 shows the computational time required for model training by each method. Evidently, in terms of computational load, the AR model requires the minimum time (1 s) for training, while LSTM is the most time-consuming model (>600 s), but this is to be anticipated mainly due to its complexity and its large number of trainable parameters. However, by inspecting the various statistical metrics, it can also be seen that the LSTM model achieves the best performance in terms of its agreement with the conventional GNSS-derived VTEC data.  Figure 4 illustrates the mae r metric for a testing period of 10 days, for the 6 different stations. MAE values ranging from 0.9 to 1.5 TECU. The 'ganp' station shows the worse values with low accuracy, while 'bor1' station achieves the highest accuracy among the other stations. Regarding the test period, the results are quite similar, however we observe a slightly increasing trend in MAE error, which is to be expected, as the predictions are carried out further away from the training period of the LSTM model.

Comparison between Different Satellites
In order to gain a further insight into the LSTM's model performance, we evaluated the VTEC prediction performance for each satellite. Figure 5 shows the estimated VTEC values for each station, per satellite on December 14, 2018. As noted, VTEC values show significant variation during the day, depending on the electron density in the ionosphere. The ionospheric delay changes slowly through a daily cycle, exhibiting its minimum values (2-4 TECU) between midnight and early morning, and reaching its peak (6-8 TECU) during the daylight hours, in the mid-latitudes.  As a next step, we evaluated the LSTM model's performance using the root mean squared error (RMSE) metric of VTEC values for each visible satellite, in a bar graph format. In general, RMSE TEC values are ranging between 0.5 and 2 TECU for each reference station ( Figure 6). For station 'bor1', the highest RMSE values are observed for G07, G26 and G11 satellites in view. For station 'ganp', G18 and G11 show the highest error values. Station 'graz' exhibits significant variation between the satellites' RMSE values, with the satellites G15, G13 and G09 showing the maximum RMSE TEC values. Station 'leij' RMSE VTEC values reach the maximum for satellites G18 and G05. Satellite G18 has also the maximum RMSE for the other two remaining stations 'pots' and 'wtzz'.
In particular, Figure 7 shows the predicted LSTM-aided vertical TEC values for each one of the GPS satellites (when they are visible from 'bor1' station), along with the vertical TEC values being derived from GAMP software (ground truth data). As noted, LSTM model can adequately predict VTEC's model performance for every satellite separately. The satellites G07, G11 and G26 have the worst performance, as expected, based on results from Figure 6. On the contrary, the satellites G05, G16 and G22 have the best performance in VTEC prediction.

Comparison among Different Ionosphere Models
To further analyze the LSTM model's performance, Figure 8 shows the mean VTEC values at every station, as obtained from NeQuick, IRI2001, IRI01-cor and GIM TEC estimates compared to the GPS-based TEC and LSTM TEC predictions, during the day. The VTEC maxima is appeared for all stations between 8:00 and 12:00 a.m. In most cases, LSTM predicted values are similar to GPS TEC derived values, however, it is noted that our model underestimates VTEC values, showing lower values than those of GPS TEC. GIM-aided TEC values are also close to LSTM derived TEC values.
NeQuick, IRI2001 and IRI01-cor values show greater variability during the day with higher maximum and lower minimum than those of GPS TEC and GIM values.  Figure 9 shows the mean difference values having as base GPS TEC estimates. As noted, in most cases LSTM predictions are closer to GPS TEC estimates, as they have small difference, except for station 'graz', where GIM model's estimates seem to have a slightly better performance.

Conclusions
We propose a supervised, unidirectional regression LSTM model for efficient prediction of vertical TEC values. Typical and widely used AR and ARMA models are linear models, suitable for TEC prediction, as ionosphere exhibits periodic changes. Thus, the AR and ARMA models' autoregressive character offers an advantage in VTEC prediction. However, the additional inclusion of selected features (such as the satellite vector, and the satellite elevation and azimuth) in the supervised LSTM algorithm, as well as, the fact that LSTM networks are equipped with memory cell and consequently they are able to preserve the worth-remembering information content in the input STEC data, gives superiority and an extra boost to our proposed method. Non-linearity and long-term prediction are additional advantages of our proposed LSTM method. As next steps, we intend to direct our studies into investigating the performance of this LSTM model under typical non-quiet conditions, including data from periods with intense ionospheric activity (e.g., periods with ionosphere anomalies and disturbances), and explore the possibility of using, with appropriate expansions, our LSTM model as a tool capable of adapting its performance to detect/catch extreme ionosphere conditions. To this end, alternative recurrent neural network methods, such as bidirectional LSTMs, that introduce non-causality, should be explored or/and check the performance of advanced Convolutional Neural Network (CNN) -based methods and how they can be optimally adapted in timeseries prediction problems for further benefit of the variant PPP techniques. In addition, we intend to explore models' behavior in datasets consisting of multiple GNSS constellations (such as GLONAS, Galileo, etc.