Long Short‐Term Memory Neural Network for Ionospheric Total Electron Content Forecasting Over China

An increasing number of terrestrial‐ and space‐based radio‐communication systems are influenced by the ionospheric space weather, making the ionospheric state increasingly important to forecast. In this study, a novel extended encoder‐decoder long short‐term memory extended (ED‐LSTME) neural network, which can predict ionospheric total electron content (TEC) is proposed. Useful inherent features were automatically extracted from the historical TEC by LSTM layers, and the performance of the proposed model was enhanced by considering solar flux and geomagnetic activity data. The proposed ED‐LSTME model was validated using 15‐min TEC values from GPS measurements over one solar cycle (from January 2006 to July 2018) collected at 15 GPS stations in China. Different assessment experiments were conducted in different geographical locations and seasons as well as under varying geomagnetic activities, to comprehensively evaluate the model's performance. These comparative experiments were conducted using an ED‐LSTM, a traditional LSTM, a deep neural network, autoregressive integrated moving average, and the 2016 International Reference Ionosphere models. The results indicated that the ED‐LSTME model is superior to the other statistical models, with R2 and root mean square error values of 0.89 and 12.09 TECU, respectively. In addition, TEC was adequately predicted under different ionospheric conditions, and satisfactory results were obtained even under geomagnetically disturbed conditions. These results suggest that the prediction performance could be significantly improved by utilizing auxiliary data. These observations confirm that the proposed model outperforms several state‐of‐the‐art models in making predictions at different times and under diverse conditions.

The results of previous studies suggest that NN models can reflect TEC variations well. A genetic algorithm-based neural network (GA-NN) TEC-prediction model developed for use over China was proposed by Song et al. (2018) and Huang et al. (2015). Its forecasting performance was significantly better than those of the back propagation-based NN (BP-NN) and IRI 2012 models. Okoh et al. (2016) integrated IRI's critical plasma frequency parameter as an additional neuron in the input layer of a regional global navigation satellite system (GNSS)-TEC NN model to monitor Nigeria. The potential extrapolation abilities and the limitations of ANNs were investigated by Habarulema and colleagues, who proposed a regional TEC prediction model for South Africa by determining the relationship between multiple inputs and TEC (Habarulema et al., 2007(Habarulema et al., , 2011. In addition, many hybrid methods involving NN components with high spatiotemporal resolution data have been proposed to model and predict the ionospheric TEC; these include the wavelet NN model (Ghaffari & Voosoghi, 2016), empirical orthogonal function NN model (Uwamahoro & Habarulema, 2015), adaptive NN model with an in situ learning algorithm (Acharya et al., 2011), and Gaussian mixture model-improved NN or radial basis function (RBF) (Huang & Yuan, 2014).
However, the ionospheric TEC is intrinsically complex. Some complex factors that affect forecasting effectiveness such as different geographical locations, seasons, and geomagnetic activities make it difficult to achieve accurate TEC prediction, which requires us to build a more powerful and complex model network. Finally, the data used in most previous studies did not cover a full solar cycle. However, complete solar cycle data coverage is necessary, given the considerable impact of solar activities on the ionospheric dynamics.
Deep learning is considered to be a function of second-generation NNs (Hinton & Salakhutdinov, 2006) and may be employed to better model spatiotemporal variations in ionospheric TEC (Orus Perez, 2019). Hochreiter and Schmidhuber (1997) developed a special deep-learning architecture known as a long short-term memory NN (LSTM NN). LSTM NNs are capable of long-term series learning and are not affected by a vanishing gradient. However, to date, LSTM NNs have been rarely used for ionospheric TEC prediction; their application for a few GPS stations (Ruwali et al., 2020;Srivani et al., 2019) has been attempted only a few times and only to generate a predicted TEC map (Cherrier et al., 2017;Kaselimi et al., 2020); Chen (2019) proposes an improved deep learning algorithm to fill in the TEC map, which shows satisfactory ionospheric peak structures at different times and under different geomagnetic conditions. Therefore, this study aimed to develop an encoder-decoder LSTM extended (ED-LSTME) model to reflect the spatiotemporal relationships between GPS stations and to predict ionospheric TEC. This model was developed and evaluated based on geomagnetic indexes, and measurements were recorded at 15 GPS stations over one solar cycle in China.

Vertical TEC Derivation
All GPS stations provide pseudorange and carrier phase measurements at two L-band frequencies. The difference between the code and carrier phase measurements of the two frequencies was calculated to obtain the pseudorange TEC (STEC a ) and phase TEC (STEC r ) along the path from a satellite to a receiver (Mannucci et al., 1998). The formula of the pseudorange TEC is as follows: where A = 40.3 m 3 /s 2 ; f 1 and f 2 are GPS signal frequencies; P 1 and P 2 indicate the respective recorded pseudo ranges; c is the speed of light; and b s,1 −b s,2 and b r,1 −b r,2 are the Differential Code Biases for the satellite and receiver, respectively. The Differential Code Bias must be estimated by eliminating the differences between the ionospheric delays of the corresponding observations (Z. Li et al., 2012;Yuan et al., 2015). The TEC phase is expressed as follows: where  1 and  2 are the carrier phases;  1 and  2 indicate the wavelengths; and     -Pajares et al., 2011-Pajares et al., , 2012.
The slant TEC was converted into the vertical TEC (VTEC) to allow the convenient examination of the distribution of TEC in different regions. This conversion was achieved using a mapping function at different ionospheric pierce points (IPPs), which are the points of intersection between the line of sight and the ionospheric shell. The charged particles of the ionosphere are assumed to be concentrated in a thin shell concentric with the earth, and this thin layer is located in the ionosphere. The shell's height above the ground varies based on the time of day or night, geographical location, sun zenith angle, and other factors; typical values range from 350-480 km and can be rapidly calculated. When the elevation angle exceeds 30°, the calculation yields highly accurate results and is suitable for application in most areas of the world. The VTEC is obtained via the equation below (Afraimovich et al., 2001;Hernández-Pajares et al., 2011).
where R E represents the average radius of the Earth; h m represents the F2 layer peak height; E 0 indicates the elevation angle of the satellite station (Xiong et al., 2014(Xiong et al., , 2016. In this study, we used dense GNSS observation station data to fix the satellite Differential Code Biases (DCB), assuming that the vertical TEC over the grid point within a certain time and space were equal and that the vertical TEC over the grid and the hardware delay of the GNSS system were solved simultaneously through the observation equation (Choi et al., 2011;Ma & Maruyama, 2003), The specific equation used is as follows: Among them, B Si is the hardware delay between different frequency signals of the ith satellite, and B rj is the hardware delay between different frequency signals of the jth receiver. The vertical TEC is based on the technique of Ma and Maruyama (2003) involving mapping and all data obtained from the Crustal Movement Observation Network of China were used. In the above approach, we used a zero-mean condition for the separation of the receiver-dependent part of the bias. In practice, the singular value decomposition (SVD) method was used to estimate DCB, a more detailed description of the method is given in Choi et al. (2011).

Data Collection and Preprocessing
In this study, GPS TEC data were obtained from the Crustal Movement Observation Network of China (CMONOC), which consists of 260 GNSS stations covering mainland China. As geographical location affects TEC estimation and data from evenly distributed stations can better reflect the geographic variability in the model performance, this study focuses on 15 stations evenly distributed throughout China ( Figure 1). Table 1 lists the geographic locations of the 15 GPS stations that were used to construct and validate the model.
The data include three main parts, including (1) vertical TEC values with a temporal resolution of 15 min were obtained from dual-frequency (1575.42 and 1227.6 MHz) observations (Mannucci et al., 1998); (2) The geomagnetic index a p was used as an indicator of the overall geomagnetic activity and magnetic storms. The magnetic effect of the K p index was utilized to measure solar particle radiation. Bartels introduced the threehour-range K p index in 1949 (Bartels & Veldkamp, 1949); the data can be downloaded from http://wdc.kugi. kyoto-u.ac.jp.
(3) Finally, the F10.7 index, the solar radio flux at 10.7 cm (2800 MHz), was used owing to its excellence and consistency as an indicator of solar activities.
XIONG ET AL.

10.1029/2020SW002706
4 of 20 We calculated the mean and standard deviation of the TEC time series for each station and then calculated appropriate thresholds to identify and exclude outliers, which were considered as values that were more than three standard deviations from the mean. To fill in the removed and missing data, we opted for an approach drawn from statistics and control theory called Kalman smoothing, which is available in the imputeTS package (Moritz, 2016) in R. The Kalman filter takes a forward pass through the data up to the current time point and can be applied in real-time. Kalman smoothing adds a backward pass through the data, thereby using all the data. After the data preprocessing and integration, a total of 2,632,410 records was collected for model development (175,494 records for each station).
This study covered a period of 12.5 years (January 2006 to April 2018). The training sets were used as input parameters, and TEC values from January 2006 to December 2016 were used, representing the period of one solar cycle. The model performance was verified using data from the same stations from January 2017 to April 2018 as validation datasets.

Time Series Sample Modeling
Preprocessing was conducted after data acquisition, including sample modeling, which is an essential step for time-series data preprocessing. The data can be converted using a machine-learning or deep-learning model. For time series problems, the sliding window is a representative modeling technique and was used in this study (Zivot & Wang, 2003).
Based on the assumption that samples and data were available for all periods, sliding windows were used for back tests to check the prediction ability of several series models. Back testing involves the following steps: 1. Select the size of the sliding window, m (i.e., the number of consecutive observations per sliding window), which depends on the sample size, T, and periodicity of the data 2. Select the prediction horizon, h, which depends on the application and the periodicity of the data. c) calculate the prediction error for each prediction; where • e nj represents the prediction error of window n for the j−step−ahead prediction • y is the response, and • ˆn j y represents the j−step−ahead prediction of window subsample n 5. Calculate the root mean square error (RMSE) values according to the prediction error for each stepahead prediction type 6. Compare the RMSEs of these models. The model with the lowest RMSE is considered to yield the best prediction performance 3. Methodology

ED-LSTME Model
The LSTM, a type of recurrent NN, was specifically designed to prevent the attenuation or explosion of model output when a given input is circulated through a feedback loop (Hochreiter & Schmidhuber, 1997). It can effectively model time dependencies and has been utilized in many fields include meteorology (Hu & Chen, 2018;Qing & Niu, 2018), network traffic (Zhao et al., 2017), and air-pollution forecasting (Yang et al., 2020).
However, the lengths of the input and output sequences in the LSTM model may differ, which could lead to vanishing and explosion gradients. The ED-LSTM has been established as an effective means to address these prediction problems (Cho et al., 2014). The key advantage of this model is that it can be easily constructed with little more than a list of inputs and outputs. The ED-LSTM consists of two models; the first is used to read the input sequence and encode it into a fixed-length vector, and the second is used to decode the fixed-length vector and output the predicted sequence.
In this study, an encoder-decoder LSTM NN was used to predict the regional TEC. In general, encoder-decoder LSTM models are used to map input and output sequences of arbitrary lengths (Sutskever et al., 2014). The encoder component can be obtained by applying one or more LSTM layers, and the model output is a vector of a fixed size (i.e., the internal representation of the input sequence) defined by the number of XIONG ET AL.  Table 1 List of GPS Stations With Geographic Coordinates memory cells in the layer. The decoder component can also be implemented using one or more LSTM layers and converts the internal representation into the correct output sequence. Thus, the fixed-sized output of the encoder model is read using the decoder.
In the ED-LSTME model, we propose (shown in Figure 3), representative features were extracted from historical TEC data and auxiliary input data via LSTM layers. The dense layer in a TimeDistributed wrapper (Chollet, 2017) was utilized as the network output, and the RepeatVector layer (Chollet, 2017) served as an adapter for both the network's encoding and decoding sections (the RepeatVector can be deployed to repeat the fixed-length vector once at each time step in the output sequence). We implemented the ED-LSTME model using the Keras (https://github.com/fchollet/keras available from GitHub) (Chollet, 2015) XIONG ET AL.
10.1029/2020SW002706 7 of 20 Previous studies have shown that solar and geomagnetic activities affect TEC variability (Blagoveshchensky et al., 2018;Kumar & Singh, 2010;Purohit et al., 2015;Verkhoglyadova et al., 2013) and that the use of auxiliary solar and geomagnetic activity index data potentially improves prediction performance. Specifically, time-delayed historical data and current auxiliary data of the solar flux index K p and geomagnetic activity index a p from all selected GPS stations were concatenated and stacked to compose an input sequence for the LSTM layers ("Main Inputs" in Figure 3).
The encoder LSTM layers process the input sequence and generate a summary of the past input sequence through the unit-state vector. After N recursive updates, the encoder LSTM layers convert the overall input sequence into the final unit-state vector. The encoder then passes the vector to the decoder LSTM, which uses it as the initial unit state for sequence generation. The decoder LSTM then recursively generates the output sequence. Subsequently, further representations of the merged features are obtained by utilizing one or more fully joined layers ("FCs" in Figure 3). Finally, the prediction output is produced based on the fully joined layer.

Parameter Optimization
Better results can be obtained by adjusting the ED-LSTME model parameters (i.e., the number of neurons in every LSTM layer where y i represents the observed value of the ith case, and  i y is the predicted value. A smaller RMSE indicates better performance. Hyperparametric optimization trials (Data Set S1 in the Supplementary Materials) were conducted. After several rounds of grid searching, an optimal result was yielded with an RMSE value of 12.09 TECU when using the following parameters:   

Model Comparison
To assess its performance, the performance of the proposed model was compared with those of other algorithms such as the classic ED-LSTM, ordinary LSTM (Hochreiter & Schmidhuber, 1997), deep neural network (DNN) (Goodfellow et al., 2016;LeCun et al., 2015), autoregressive integrated moving average (ARIMA) (Hyndman & Athanasopoulos, 2018;Makridakis & Hibon, 1997), and IRI 2016 models (Bilitza, 2018;Rawer et al., 1978). Hyperparametric optimization, which is used in all the methods except for IRI 2016, allow for the identification of the optimum parameters. Based on this optimization, relevant technology can be selected with high confidence, ensuring the use of a robust method. For the classic ED-LSTM, ordinary LSTM, and DNN models, the hyperparameters of the neuron number in every LSTM and FC layer, training epoch number, and batch size were selected for model optimization. The optimized ARIMA model was obtained by applying the Hyndman-Khandakar algorithm (Hyndman & Khandakar, 2008) for automatic modeling based on the combination of unit root tests, the minimization of the Akaike information criterion (AICc), and maximum likelihood estimation (MLE). Hyperparametric optimization trials conducted for these models are described in Datasets S1 of the Supplementary Material. In addition, the degrees of freedom (the number of parameters in the model) for the proposed ED-LSTME model; the encoder-decoder LSTM model; the LSTM, DNN, and ARIMA models; and the IRI 2016 model were 140116, 138916, 226992, 235292, 8, and 28, respectively.
To assess the performance of the various models, four different statistical indicators were used in this study: RMSE, the coefficient of determination (R 2 ), mean absolute error (MAE), and the correlation coefficient (ρ), given by Equations 6-9, respectively. MAE refers to the absolute error between the true value and the predicted value. RMSE is the square of the difference between the true value and the predicted value. When MAE and RMSE are used together, the degree of dispersion of the sampling error can be seen. The correlation coefficient is used to describe the degree of linear correlation between two variables. The R-square is generally used in regression models to evaluate the degree of agreement between the predicted value and the actual value.
where n is the number of cases, y i is the observed value of the ith case,

Model Performance
The performance of the various models is detailed in  As an intelligent algorithm, the encoder-decoder LSTM model offered better TEC prediction, with R 2 and RMSE values of 0.5462 and 24.1536 TECU, respectively. The ED-LSTME model was superior to the standard ED-LSTM owing to its consideration of space environment information. The ED-LSTME R 2 value increased by 0.34 over that of the standard ED-LSTM (from 0.5462 to 0.8862) and the RMSE value decreased by 12.06 TECU (from 24.1536 to 12.0954 TECU). These findings indicate that the proposed ED-LSTME model delivered the best performance among the considered models, followed by that of the encoder-decoder LSTM, LSTM, DNN, ARIMA, and IRI 2016 models.
While conventional models (e.g., ARIMA) have delivered reasonable results under geomagnetically quiet conditions (Xiaohong et al., 2014;Zhang et al., 2013), they did not deliver the expected performance on a larger scale, with R 2 and RMSEs values of 0.0689 and 34.5988 TECU, respectively. However, the R 2 and RMSEs values ranged from 0.2784 to 0.8862 and from 12.0954 to 30.4596 TECU, respectively. Therefore, the more advanced models (ED-LSTME, encoder-decoder LSTM, LSTM, and DNN) demonstrated better TEC estimation performance.
As the superiority of the ED-LSTME model was confirmed under the study scenario, the performance of the model was further evaluated. The model's efficiency was investigated using Taylor diagrams (Elvidge et al., 2014;Taylor, 2001) (Figure 4). The diagram shows that the ED-LSTME model performed significantly better than did the other models, as ED-LSTME has a stronger correlation to the expected results and a smaller standard deviation than those of the other models. The error standard deviation of IRI 2016 was larger than those of the more advanced models (ED-LSTME, encoder-decoder LSTM, LSTM, and DNN) and ARIMA (albeit only slightly). The standard deviation values for encoder-decoder LSTM, IRI 2016, and DNN were similar and were superior to those of ARIMA and LSTM. Compared with those of the other methods, the standard deviation of the TEC values predicted using the ED-LSTME model displayed very little bias. Thus, based on a consideration of the collective results, the ED-LSTME model significantly outperformed the other models.
In addition, from the International GNSS Service (https://www.igs.org), the global ionospheric TEC levels from January 2017 to April 2018 were relatively low, with an average of around 12 TECU. To prove our proposed model has better generalization ability during higher ionospheric activity, we select the relatively high value in the test set (TEC values greater than 20 TECU) and their corresponding predicted values, and conducted a comparative study. From Figure S3, it can be concluded that while the prediction performance of other models is worse than before, the ED-LSTME not only performs best in the error standard deviation and correlation, but also is the best model in bias. These results show that the ED-LSTME has good generalization ability.

Variation of Model Performance With Geographical Location
Geographical location affects TEC estimation. Therefore, models were established for each station to determine the geographic variability in the model performance. Table 1 provides the geographical coordinates of the stations, and Tables S1 and S2 present the overall evaluation results for the ED-LSTME, encoder-decoder LSTM, LSTM, DNN, ARIMA, and IRI 2016 models for each station.  the RMSEs of all six models were lower in the middle and high latitudes. This can possibly be explained by the equatorial ionospheric anomaly (EIA) in lower latitudes, which may lead to larger RMSEs as suggested by Song et al. (2018).
In addition, the R 2 (upper left), MAE (upper right), RMSE (lower left), and ρ (lower right) of the six models are compared in the bar chart in Figure 6. Notably, the ED-LSTME model had the highest prediction efficiency and outperformed the other five models across latitudes.

Seasonal Variation in Model Performance
Previous studies indicated TEC prediction model performance varies seasonally (Mukesh et al., 2020;Ruwali et al., 2020;Song et al., 2018;Tebabal et al., 2018Tebabal et al., , 2019. Therefore, all considered models were operated in each season to investigate their performance throughout the solar year. The seasons were delineated as spring (March-May), summer (June-August), autumn (September-November), and winter (December-February), for which the number of data records was 138,240, 139,776, 141,312, and 141,312 respectively.
XIONG ET AL.
10.1029/2020SW002706 Figure 4. Taylor diagram of model bias and standard deviation of errors. The azimuth indicates the correlation, the radial distance represents the standard deviation, and the semicircles centered on the "Expected" label represent the error standard deviation. The color scale represents the degree of bias. Each quantity was standardized to draw multiple parameters. The appropriate "factors" in the upper right of the chart can be utilized to modify the original values. The chart presents the performances of the ED-LSTME, encoder-decoder LSTM, LSTM, DNN, ARIMA, and IRI 2016 models. ARIMA, autoregressive integrated moving average; DNN, deep neural network; ED-LSTME, encoder-decoder long short-term memory extended; IRI, International Reference Ionosphere; LSTM, long short-term memory. Figure 7 presents the seasonal variations in model performance using the testing data. Tables S3 and S4 list the R 2 , MAE, RMSE, and correlation coefficient (ρ) values for the predicted TEC and seasonal variations in these values, respectively.
Tables S3 and S4 indicate that the ED-LSTME model delivered the best performance in each season (R 2 = 0.8209, 0.8425, 0.8769, and 0.9066 for the testing data from spring, summer, autumn, and winter, respectively). It was followed by the encoder-decoder LSTM, LSTM, DNN, and ARIMA models (R 2 values of 0.7332, 0.7963, 0.0659, and 0.0204, respectively); the IRI 2016 model delivered the worst performance. This indicates that relatively advanced models achieved significantly better TEC-estimation accuracy than did simpler models based on their consideration of more space environment parameters.
XIONG ET AL.

10.1029/2020SW002706
12 of 20 In addition, despite outperforming the conventional methods, the performance of the DNN model was slightly inferior to those of the LSTM-based models in each season; this suggests that the LSTM models could learn (or "remember") larger temporal dependence and to perform well in time-series forecasting. Overall, prediction error was largest in winter, although the winter TEC is smaller than those in spring and summer. This may reflect the effects of solar activities.
The performances of most models on the seasonal scale were considerably better than those on the annual scale, thus highlighting the effect of seasons on TEC estimation. The R 2 values obtained using conventional models (DNN, ARIMA, and IRI 2016) at the seasonal scale (especially in spring and summer) were significantly improved over those obtained at the annual scale, thus verifying that conventional models are more appropriate for seasonal observations.

Variation in Model Performance Owing to Solar and Geomagnetic Activity
To further investigate the prediction abilities of the proposed models under varying ionospheric conditions at diverse locations and times, the TEC values obtained on quiet days (K p < 3.0 or a p < 56) and under disturbed conditions (K p > 3.0 or a p > 56) were compared. The geomagnetic activity indices K p and a p were utilized to determine the quietest and most disturbed days of every month, as it is known that ionospheric variability significantly increases during geomagnetic activity.
XIONG ET AL.

10.1029/2020SW002706
13 of 20 Tables S5 and S6 and Figure 8 present the specific errors of the models. The ED-LSTME model yielded the best prediction results based on the comparison of errors for quiet days and disturbed conditions. Under disturbed conditions, the RMSE increased by 5.04/1.89 TECU (K p /a p ), and the R 2 value decreased by 0.11/0.08 TECU (K p /a p ). The other five models exhibited diverse changes. Notably, only the ED-LSTME model could accurately capture complex trends during storms and other geomagnetic disturbances, likely because this model considered solar and geomagnetic activity values as input data. In order to further verify whether the performance can be improved by adding auxiliary input, extended LSTM (LSTME) and extended DNN (DNNE) models have been developed by adding solar and geomagnetic activity values as auxiliary inputs.
Taylor diagrams in Figure S2 shows that both the error standard deviation and correlation of LSTME and DNNE performed better than that of LSTM and DNN, respectively. The biases are also less in the LSTME and DNNE models. Overall, the ED-LSTME is still the best performing model and shows a strong capability in TEC forecasting across the solar and geomagnetic activity.
Further, to better verify the effectiveness of the proposed ED-LSTME model, we selected four quiet days (days with a K p index under 3.0) in spring (March 11, 2017), summer (August 16, 2017), autumn (September 10, 2017), and winter (December 10, 2017) as test cases. The TEC results obtained using the proposed ED-LSTME model, encoder-decoder LSTM, LSTM, ARIMA, and IRI 2016 models, NeQuick model (Nava et al., 2008), and GPS observations at 15 stations are compared in Figure 9. ED-LSTME, encoder-decoder long short-term memory extended; IRI, International Reference Ionosphere; LSTM, long short-term memory; MAE, mean absolute error; RMSE, root mean square error; TEC, total electron content. Figure 9 demonstrates that the TEC estimated by the ED-LSTME model agrees with the observed TEC, implying that the model is notably superior to other methods with respect to capturing diurnal variations in TEC, while NeQuick and IRI 2016 models only captured the overall trends of the actual data. At low-latitude stations such as XIAM and XIAG, model performance was weaker due to the EIA.
To study the prediction performance of the models under conditions of geomagnetic disturbance, we selected four intense geomagnetic storms (with K p > 3.0) in 2017 and analyzed their effects on the predicted TEC. These storm events (K p index above 3.0) occurred on March 2 (spring), August 19 (summer), September 8 (autumn), and December 5 (winter). Figure S1 shows that the ED-LSTME model could predict the TEC disturbances associated with these geomagnetic storms for all stations. The ED-LSTME model's performance was particularly strong at WUSH, YANC, and KMIN stations and was typically better in spring than in winter.
The uncertainty of the proposed ED-LSTME model is shown in Figure 9 and Figure S1. We calculated the standard deviation between the results obtained using the proposed ED-LSTME model and GPS observations and considered the result to represent the uncertainty of the proposed ED-LSTME model. As parameter tuning of the proposed deep-learning model is time-consuming and challenging, we cannot guarantee that optimized parameters were obtained for the models trained at each GPS station; however, most cases were covered through the grid search method employed in our study. Still, this introduced additional uncertainty to the TEC forecasts of the proposed model. XIONG ET AL.

Discussion
The results imply that all the tested deep-learning models had high learning levels and prediction abilities, with potential variations owing to their diverse learning methods (Tien et al., 2020). The deep-learning models were more flexible than conventional models with respect to nonlinear learning, particularly of large datasets, although their implementation was more complex than that of traditional machine-learning models. In conventional machine-learning models, most input features require expert classification to reduce the complexity of the data and make the patterns more accessible to the learning algorithms. In contrast, the deep-learning models incrementally learned to extract high-level features from the data, negating the need for expert feature classification.
Although the results of this study confirmed that LSTM-based deep-learning models can be used for shortterm predictions (e.g., using one-day historical input TEC data to forecast 3 h of TEC data), future studies should focus on the capability of deep-learning models to estimate TEC in the medium and long terms.
Preprocessing methods, such as the combination of wavelet decomposition with deep-learning models, could also be explored for TEC forecasting. Decomposition methods can be utilized to remove noise from the data, improving the model accuracy. Furthermore, the hybrid CNN-LSTM model, deep reinforcement learning, generative adversarial network, and other deep-learning models should also be explored. Moreover, the storm-time ionospheric prediction is usually related to seasonal variation. Tang  that the LSTM-based ionospheric storm-time prediction model shows a significant dependence on seasonal variation, which provides a new dimension for our future research.

Conclusions
In this study, an ED-LSTME model was proposed to predict ionospheric TEC based on historical TEC, solar flux, and geomagnetic-activity data, and it demonstrated good performance in modeling time series with long time dependencies and in identifying optimum time lags and prediction horizons. To evaluate the performance of the proposed model, 15-min TECs from GPS measurements over one solar cycle were collected from 15 GPS stations in China. The same datasets were used to compare six different models: ED-LSTME, encoder-decoder LSTM, LSTM, DNN, ARIMA, and IRI 2016. The models were trained using optimal hyperparameters. Regarding R 2 , RMSE, MAE, and correlation coefficient, the ED-LSTME outperforms the other algorithms, and compared with conventional models, such as ARIMA and IRI 2016, deep learning-based models delivered better prediction performance. In comparison with the DNN model, the proposed ED-LSTME, Encoder-Decoder LSTM, and conventional LSTM identify spatiotemporal correlations more effectively and deliver better prediction performances. The performances were significantly improved with the use of auxiliary data on solar flux and geomagnetic activity.
To more accurately assess the prediction performances of the models, geographical location, seasonal variation, and geomagnetic-activity variations were also considered as factors affecting TEC estimation. The RMSEs obtained for the six models were larger at low latitudes, which might have been due to the EIA. Data sets for each season were established to determine seasonal variations in model performances. The results revealed that the proposed model had good predictive power in each season and that seasonal observations were more suitable than annual observations. Under both quiet and disturbed conditions, the ED-LSTME model effectively predicted short-term changes in TEC during solar flux and geomagnetic events.
Long-term predictions are more difficult than short-term predictions, as they require more pertinent historical input data. Hence, the investigation of the performance of different deep-learning models to make medium-and long-term TEC predictions is recommended for future studies, as is the operational implementation of the proposed model.