Forecasting Global Ionospheric TEC Using Deep Learning Approach

Global ionospheric total electron content (TEC) maps are widely utilized in research regarding ionospheric physics and the associated space weather impacts, so there is a great interest in the community in short‐term ionosphere TEC forecasting. In this study, the long short‐term memory (LSTM) neural network (NN) is applied to forecast the 256 spherical harmonic (SH) coefficients that are traditionally used to construct global ionospheric maps (GIM). Multiple input data, including historical time series of the SH coefficients, solar extreme ultraviolet (EUV) flux, disturbance storm time (Dst) index, and hour of the day, are used in the developed LSTM NN model. Different combinations of the above parameters have been used in constructing the LSTM NN model, and it is found that the model using all four parameters performs the best. Then the best performing LSTM model is used to forecast the SH coefficients, and the global hourly TEC maps are reproduced using the 256 predicted SH coefficients. A comprehensive evaluation is carried out with respect to the CODE GIM TEC. Results show that the first/second hour TEC root mean square error (RMSE) is 1.27/2.20 TECU during storm time and 0.86/1.51 TECU during quiet time, so the developed model performs well during both quiet and storm times. Moreover, typical ionospheric structures, such as equatorial ionization anomaly (EIA) and storm‐enhanced density (SED), are well reproduced in the predicted TEC maps during storm time. The developed model also shows competitive performance in predicting global TEC when compared to the persistence model and two empirical models (IRI‐2016 and NeQuick‐2).


Introduction
The ionosphere is the ionized part of the Earth's upper atmosphere where sufficient plasma exists to affect the propagation of radio waves. It is arguably the most important region in terms of space weather impacts of modern technologic society, because irregularities in the ionosphere frequently degrade and disrupt satellite navigation and communication. The total electron content (TEC) refers to the total number of electrons integrated along a tube of 1-m 2 cross section (unit: TECU, 1 TECU = 10 16 electrons/m 2 ). TEC is a significant ionospheric parameter, which is widely applied to the satellite ionospheric delay correction and scientific research area about space weather impacts on the ionosphere and thermosphere (Liu et al., 2020). As a consequence, there is a great interest in the community in the TEC forecasting.
The ionospheric TEC along the ray path between the satellite and receiver can be derived by using dual frequency global navigation satellite system (GNSS) measurements, because of the dispersive characteristics of the ionosphere. For instance, the International GNSS Service (IGS) ionospheric working group, consisting of seven Ionospheric Associate Analysis Centers (IAACs), has started generating reliable global TEC maps based on worldwide GNSS data since 1998, which are widely used in the ionosphere monitoring (Hernández-Pajares et al., 2009;Li et al., 2015;Schaer, 1999;Yao et al., 2018).
In terms of the short-term GNSS-based TEC forecasting, many approaches have been developed in recent decades, such as least squares collocation, discrete cosine transform (DCT), linear regression, adaptive autoregressive model, and neural networks (NN). Based on the least squares collocation method, Schaer (1999) predicted the global TEC parameters in the next 2 days by extrapolating the spherical harmonic (SH) coefficients used to fit the TEC map in the previous 30 days. Tulunay et al. (2006) presented a NN technique to forecast the TEC maps over Europe, which is trained based on the Levenberg-Marquardt algorithm, and results show that the developed model learned the shape of the inherent nonlinearities during space weather events. García-Rigo et al. (2011) developed a global TEC prediction model by applying the Discrete Cosine Transform (DCT) to the TEC maps and then used a linear regression module to forecast the time evolution of the DCT coefficients, but physical information, including the solar and geomagnetic activity, was not included in the prediction model. Habarulema et al. (2011) developed an ionospheric TEC prediction model in South Africa using artificial NN and evaluated its capability by using independent ionosonde data and observed TEC values. They found that the developed model performed well during quiet periods while the accuracy decreased under storm conditions. Wang et al. (2018) proposed an adaptive autoregressive model to predict the SH coefficients used in TEC map fitting, and then the global TEC maps can be generated based on the predicted SH coefficients. Results show that the predicted TEC has better performance during the low solar activity periods than the middle-to-high solar activity periods.
Nowadays, deep learning techniques, such as the convolutional neural network (CNN) and recurrent neural network (RNN), are very promising and have been increasingly applied to the short-term ionosphere forecasting. For instance, Sun et al. (2017) developed a single-station TEC forecasting model at Beijing station in China based on the long short-term memory (LSTM) NN using 5-day sequence of TEC, AP, and F10.7 as inputs. Results showed that their model outperformed the conventional multilayer perceptron network, though the prediction performance under disturbed conditions still needs further improvement. Similarly, Srivani et al. (2019) also developed a LSTM forecasting model over the Bengaluru station in India using similar inputs as Sun et al. (2017) but for different network configuration parameters (i.e., the input time span is 24 hr). They found that the model is able to capture the TEC variations even during the active solar and geomagnetic activities. Boulch et al. (2018) presented three convolutional RNN architectures for global ionospheric TEC prediction, and their results showed that the prediction performance was comparable with state-of-the-art methods, such as the autoregressive (AR), autoregressive moving average (ARMA) model, and the radial basis function (RBF) NN.
However, most of the above-mentioned deep learning approaches for the TEC forecast are either for a single-station forecasting model or have not considered the physical parameters representing the solar and geomagnetic disturbances. In this paper, we aim at developing a novel deep learning model to forecast the SH coefficients used in constructing the global TEC map by using time series of the SH coefficients and multiple physical parameters representing solar and geomagnetic activity levels. The model performance is then evaluated by comparing the prediction with independent TEC data under both geomagnetic quiet and storm conditions.

Global Ionosphere Maps Reconstruction
Based on deferential delay measurements from globally distributed GNSS tracking sites, the Centre for Orbit Determination in Europe (CODE) develops hourly global ionosphere maps (GIM) in the solar-geomagnetic reference frame using SH expansion (Schaer, 1999), where β and s are the geomagnetic latitude and Sun-fixed longitude of the ionospheric pierce point (IPP), respectively. The IPP is the intersection point between the satellite-receiver ray path and the single-layer height. The Sun-fixed longitude is defined as s = 15(t − 12) + λ, where λ is the geographic longitude of the IPP and t is the universal time. TEC(β, s) is the vertical TEC (VTEC, hereafter referred as TEC for simplification) at the location of the IPP. e P nm denotes the normalized associated Legendre function with degree n and order m. e A nm and e B nm are the hourly SH coefficients to be determined, which are then used for constructing the global TEC map. n max is the maximum degree of the SH expansion (n max = 15), which means that there are 256 ((n max + 1) 2 = 256) SH coefficients for 1-hr interval.
In general, we are not interested in the individual SH coefficients but rather in the TEC value, which is represented as a function of geographic coordinates (β, λ) and time t, namely, TEC(β, λ, t). For global TEC map 10.1029/2020SW002501 Space Weather reconstruction, once obtaining the 256 SH coefficients ( e A nm and e B nm ) at time t, we are able to generate GIM using Formula 1 with a spatial resolution of 2.5°in latitude and 5°in longitude. These resolutions are chosen to be consistent with the TEC maps developed at CODE for model validation and performance evaluation purpose. These 256 SH coefficients are estimated routinely with 1-hr interval and are available from the CODE analysis center (http://ftp.aiub.unibe.ch/CODE/) on a daily basis since October 2014. We predict the time series of these SH coefficients in the next 2 hr using the CODE SH coefficients in the past and then construct the GIM.

SH Coefficients Forecast Using the LSTM NN
One challenge of this study is to design an appropriate NN architecture to extract temporal trend and features from previous states and to predict the SH coefficients from historical sequences and related external drivers, such as solar and geomagnetic activities.

Data Preparation
In this study, the most important data set is the time series of the SH coefficients, which is freely available from the CODE analysis center. Moreover, considering that the solar and geomagnetic activities play significant roles in changing the ionospheric TEC, realistic solar irradiance and geomagnetic activity index are also needed. Here, the solar extreme ultraviolet (EUV) flux from Flare Irradiance Spectral Model (FISM) and the Disturbance Storm Time (Dst) index representing the ring current strength and storm phases are used in this study. One of the reasons for choosing Dst instead of other geomagnetic indices, such as the 3-hr Kp index, is because the time resolution of Dst index is 1 hr, which is consistent with the SH coefficients provided by the CODE center. All of the above data sets are collected from 19 October 2014 to 31 December 2016 (see Figure S1 in the supporting information). It should be noted that FISM estimates the solar irradiance at wavelengths from 0.1 to 190 nm at 1-nm resolution every minute (Chamberlin et al., 2007(Chamberlin et al., , 2008, meaning that there are 190 wave bands for the solar EUV flux. Figure 1a shows the Pearson correlation between the 190 FISM solar EUV fluxes. We recognize that these 190 features are highly correlated with each other, and not all of these parameters are needed as the input data of the network, because too many redundant features might result in overfitting problems in the machine learning algorithm. The principal component analysis (PCA) is a very powerful procedure to transform many correlated parameters into a smaller number of uncorrelated parameters via dimension reduction, so here we extracted several principal components from the 190 FISM EUV fluxes using the PCA method. Figure 1b shows the first three principal components from the 190 original features. We found that the first principal component could account for more than 99.6% of the variability of the 190 features, so in this study only the first principal component extracted from 190 FISM features is used to represent the solar EUV flux.
In addition, it is well known that ionosphere has a diurnal variation, so we introduce time (hour of day, HD) as an input to enable the network learn the diurnal variations of the SH coefficient sequences. To allow

10.1029/2020SW002501
Space Weather continuity across the midnight boundary, the HD is defined as two quadrature components (Poole & McKinnell, 2000), where HD s and HD c are the sine and cosine components of the HD, respectively.

NN Architecture
The LSTM NN, a special type of the RNN, is capable of learning time dependence in sequential prediction through the forget gate weights and is developed to avoid gradient exploding and vanishing problems existing in conventional RNN (Chen et al., 2019;Gers et al., 1999;Hochreiter & Schmidhuber, 1997). In this study, we try to forecast the SH coefficients separately 2 hr ahead by using LSTM NN first and then use the predicted SH coefficients to construct the GIM by using the SH expansion function.
The overall architecture for a single SH coefficient prediction is shown in Figure 2, which includes the input data in the input layer, two LSTM layers, and one dense layer in the hidden layer, as well as the predicted SH coefficients in the output layer.
In the input layer, the input time span is set to be 24 hr, as the autocorrelation for each SH coefficient at 24 hr is the highest due to the ionospheric diurnal variation. The final input candidates after the data preparation procedure are the SH coefficients, HDc, HDs, the first component of the FISM solar irradiance, and the Dst index of the previous 24 hr.
In the hidden layer, two stacked LSTM layers are used to extract representative features from historical input data and then are connected to the fully connected (FC) hidden layer. The neuron node number in each hidden layer is chosen to be 48 (n = 48) throughout the paper, which is determined by training the NN from Node 1 to 60 (see details in in section 3.1 and Figure 3). In the output layer, the output parameters are the individual SH coefficients for the next 2 hr, which could be extended further in the future.
The entire data set is split into two groups, namely, Group 1 (Training set: 1 January 2015 to 26 May 2016) and Group 2 (Testing set: 19 October to 31 December 2014 and 27 May to 31 December 2016) (see Figure  S1 in the supporting information). To find the optimum NN parameters, fivefold cross-validation strategy is utilized in the Training set of Group 1. More specifically, Group 1 is divided into fivefolds, where fourfolds are used for training the NN, and the remaining onefold is used for validation. This process is repeated for five times, and the final performance is the averaged results of these five times. Finally, the NN parameters showing best performance for both the fourfold training and onefold validation will be used as the final NN parameters. Once the optimum NN parameters are determined by the fivefold cross validation, we test their performance in TEC prediction using Group 2 (Testing set), which is not involved in training the NN model.

Evaluation Metrics for the Deep Learning-Based Model
One of the core tasks in building deep learning model is to evaluate its effectiveness and reliability. In this paper, several important evaluation metrics, such as mean absolute error (MAE), mean square error (MSE), and root mean square error (RMSE), as well as mean error (ME), are used to qualify how well the proposed deep learning algorithm in predicting the SH coefficients and the resulting TEC results. 10.1029/2020SW002501

Space Weather
where y j and b y j are observations and predictions, respectively, j = 1, 2, 3, …, N, and N is the total number of data samples. Considering the various ranges of the input and output features, each feature in this study is normalized to [0,1] before training the NN model, which is a commonly used method in machine learning. The outputs or predictions can be recovered to the actual dynamic range using an inverse normalization process. In our study, if the skill scores are calculated based on the normalized observations or predictions, we will call them as normalized skill scores.

Construction of GIM Using the Predicted SH Coefficients
As mentioned in section 2.1, global ionospheric TEC maps can be represented as a function of geomagnetic coordinates and time when the 256 SH coefficients are plugged into Formula 1. Thus, the hourly GIM can be produced when the 256 SH coefficients are predicted successfully from the proposed NN model in section 2.2. The spatial resolution of the resulting GIM is 2.5°in latitude and 5°in longitude, which is consistent with the global TEC products released by the CODE ionosphere analysis center.

Determination of the NN Parameters
Neural nodes and inputs are two types of key parameters in the NN architecture, so in this part we will try to determine them appropriately.
Here, we take the first SH coefficient prediction results as an example to identify the neural nodes, and similar procedures are performed for the rest of the 255 SH coefficients. We train the proposed NN by changing the node number from 1 to 60 in the step of 1 and evaluate the Training/Testing performance in terms of MAE and MSE. Figures 3a-3d show the MAE/RMSE error distributions for various neural node numbers on Train/Test data sets, and Figure 3e shows the learning curve variations with different training epochs. It can be seen that the SH prediction achieved the smallest error when the neural node number is 48, except in Figure 3d where that optimal number is 49. Moreover, the learning curves in Figure 3e show overall very good fitting, because the normalized MAE errors of both Train and Test sets decrease gradually with the training epochs and are overlapped with each other. Therefore, the ideal neural node number in hidden layer is set to 48 throughout this paper.
Another important task is to determine the most appropriate input parameters for the input layer. As mentioned in section 2.2, the candidates in the input layer are the SH coefficient, HDc, HDs, the first principal component of the FISM solar irradiance, and the Dst index. In this study, we try different input combinations and evaluate the corresponding prediction performance in terms of the normalized MAE/RMSE by using fivefold cross-validation strategy. Table 1 presents the normalized MAE/RMSE for different input parameter combinations, including I (SH coefficient + first component of FISM solar irradiance + Dst + HDc + HDs), II (SH coefficient + first component of FISM solar irradiance + HDc + HDs), III (SH coefficient + Dst + HDc + HDs), and IV (SH coefficient + HDc + HDs). It is found that the MAE/RMSE obtained from both the training and validation sets of Group 1 reaches minimum when the input Combination I is used in the input layer. This suggests that the combination of the solar and geomagnetic activity information is able to improve the SH coefficient prediction results to some degree. Moreover, similar prediction performances are achieved when the input Combination I is applied to next four SH coefficients (see Tables S1-S4 in the supporting information), which are also important ones for TEC determination (see Figure S2 in the supporting information). Therefore, the input Combination I is used

Space Weather
in the input layer of the proposed NN architecture, although the improvement in predicting the SH coefficient is not significant compared to the other input combinations.

10.1029/2020SW002501
Space Weather Figure 5. Global TEC maps with 6-hr interval under both quiet (11 October, Figure 5a) and storm (13 October, Figure 5b) days in 2016. The predicted TEC maps from our proposed NN and corresponding ones from CODE GIM, along with their differences, are, respectively, given in the left, middle, and right panels of each figure. Full sets of global TEC maps with 1-hr interval for both first and second hours are available in Movies S1-S12 of the supporting information.

The Performance of the Predicted Global TEC Maps
After predicting the 256 SH coefficients separately, we are able to produce global TEC maps according to Formula 1. To test the accuracy of the proposed algorithm, the resulting TEC maps obtained from the predicted SH coefficients are evaluated by comparing with the CODE GIM products. It should be noted that the 256 SH coefficients used for generating the CODE GIM are the Test data mentioned in section 2, and they are not involved in training the LSTM NN. Figure 5 shows the resulting predicted TEC maps and the corresponding CODE GIM for the first hour, along with their difference, for randomly selected geomagnetic quiet and storm days. Compared to the CODE TEC maps, the typical equatorial ionization anomaly (EIA) crest regions can also be clearly observed from the predicted TEC maps, and the global residuals between them are very small, mostly within 1 TECU. These indicate that the resulting TEC prediction is satisfactory. Residuals may reach to 5 TECU near the EIA region. Moreover, the storm time TEC maps also agree well with CODE TEC products, since complicated ionosphere features with large TEC magnitudes, that is, storm-enhanced density (SED) structure (e.g., Foster, 1993;Foster et al., 2005;Zou et al., 2013Zou et al., , 2014, are well reproduced in the predicted TEC map (see Figure 5b). This result suggests that the developed deep learning approach is also reliable even for the storm time TEC forecasting. Figure 6 shows the forecasted TEC residual histogram during both quiet (blue) and storm (red) days for the Testing set. It is clear that the first/second hour TEC residual distributions for both quiet and storm time are concentrated around zero with very small RMSE errors, though the first/second hour RMSE for the storm time (1.27/2.20 TECU) is slightly larger compared to those during quiet time (0.86/1.51 TECU). Therefore, the comparison demonstrates that the proposed deep learning algorithm performs well in predicting the global TEC for both storm and quiet times. It is also evident that the TEC residual distribution for the first hour is more concentrated around zero than that for the second hour, and the averaged RMSE increases slightly from 1.27/0.86 to 2.20/1.51 TECU. This suggests that our proposed deep learning approach performs better in predicting the first hour TEC than the second hour, but the degradation with time is minor. Figure 6. Statistical evaluation of global predicted TEC with respect to the CODE GIM products during both quiet (blue curves) and storm days (red curves) for the testing set. Left and right panels represent the TEC prediction results for the first and second hours, respectively.

Space Weather
We also compared the proposed method's performance with respect to the persistence model. The persistence implies that the predicted values remain unchanged between current time t and future time t + δt (Kleissl, 2013), namely, where δt = 1, 2 hr. TEC(t) represents the TEC observation at current time t and the TEC(t + δt) is the predicted TEC of the first/second hour based on TEC(t). It is found that the TEC forecasting errors of both the first and second hours from the proposed NN model are much smaller than that from the persistence model, and the proposed NN model shows consistent improvements of above 60% when compared to the persistence model. This suggests that the proposed NN model significantly outperforms the persistence model.
To better demonstrate the merit of the proposed NN model, NeQuick-2 and IRI-2016 models are also involved in comparison. The NeQuick-2 is the latest version of the NeQuick ionosphere electron density model, and its outputs (e.g., electron density, TEC) depend on solar activity (monthly mean sunspot number R12 or 10.7-cm solar radio flux F10.7), location, and time (Nava et al., 2008). The IRI-2016 model, as the most widely used empirical climatological model of the ionosphere, can also provide electron density and TEC parameters in the ionospheric altitude range for a given location, time, and date, and it is driven by solar (F10.7, R12), ionospheric, and magnetic indices (Bilitza et al., 2017;Liu et al., 2019). In our study, the daily F10.7 indices are used for driving both IRI-2016 and NeQuick-2 models, while other optional input parameters are set to default ones. Also, the storm mode of the IRI-2016 model is turned on in this study. Figure 7 shows the hourly TEC RMSE from the NN, IRI-2016 and NeQuick-2 models with respect to the CODE GIM over the entire Testing set. The first period, that is, 19 October to 31 December 2014, is shown on the left, and the second period, that is, 27 May to 31 December 2016, is shown on the right. The first and second hours RMSE errors are shown in the top and bottom rows, respectively, and the hourly Dst index, which represents the condition of geomagnetic activities, is also shown in the top panel. It can be seen that both the first/second hour TEC RMSE errors from the proposed NN are smaller than those from the IRI-2016 and NeQuick-2 models. The NN prediction errors over the entire Testing set are very low and stable with averaged RMSE within 1.06/1.84 TECU for the first/second hour, while the RMSE errors from the IRI-2016 and NeQuick-2 models are around 9.21/5.5 TECU, respectively. Therefore, the proposed NN approach is competitive in predicting the global TEC when compared to the traditional IRI-2016 and NeQuick-2 models. One can also see that the time series of the NN RMSE errors remain small during storm periods despite relative increase comparing with quiet times.  Comparing to those models with no external solar or geomagnetic activity drivers concatenated into the LSTM layer, our developed LSTM NN model is able to improve the forecast of the SH coefficients by including the solar and magnetic indices. Therefore, the optimal input data used in the input layer of the LSTM NN architecture are time series of the SH coefficients, solar EUV flux and the Dst index, and time of day. The advantage of this input combination can also be confirmed by the good SH prediction performance in Figure 4, because the Pearson correlation values between the predicted SH coefficients and observations from the Train/Test are generally higher than 0.99, and the predicted errors are also extremely small. We found that the SH predicted error for the first hour is slightly smaller than those for the second hour, suggesting that the model has better performance in predicting the SH coefficients for the first hour than the second hour.
After using the developed LSTM model to forecast all SH coefficients, the global TEC maps are generated and used for evaluation. The global residuals between the predicted TEC and the CODE GIM TEC for both geomagnetic quiet and storm days are very small, mostly within 1 TECU and rarely reaching 5 TECU. Prominent structures in the ionosphere, such as EIA and SED during storm time, have been reproduced. A comprehensive evaluation of the predicted TEC has been carried out with respect to the independent CODE GIM products, and results show that the first/second hour TEC residual distributions for both quiet and storm times are concentrated around zero with very small RMSE errors: the first/second hour TEC RMSE are 0.86/1.51 and 1.27/2.20 TECU during the quiet and storm times, respectively. These indicate that the developed LSTM model performs well during both quiet and storm times. Moreover, the proposed approach also shows competitive TEC forecast performance in terms of RMSE when compared to the persistence model and the traditional IRI-2016 and NeQuick-2 models. Both the IRI-2016 and NeQuick-2 models are valuable and widely used community models. The conclusion regarding the comparison of both the IRI-2016 and NeQuick-2 models is based on the current version of the model trained and validated using relatively limited data sets. In the future, we will incorporate more data, in particular larger geomagnetic storms, for training and validation and perform a more thorough evaluation.
Currently, the proposed deep learning approach has shown satisfactory short-term TEC prediction performance based on experimental results obtained, and we will extend its prediction horizon to longer term (beyond 2 hr) in further study. At this point, real-time SH coefficients from the CODE center are not available and there is usually a several days' delay for the SH coefficients to become available. In the future, if the SH coefficients become available in real time or near real time, the developed deep learning model can be readily used in forecasting the coefficients and TEC maps.

Data Availability Statement
The SH coefficient time series and GIM maps are freely available from the CODE analysis center (http://ftp. aiub.unibe.ch/CODE/). The solar EUV data are obtained from the FISM. The geomagnetic data are available from the NASA Goddard Space Flight Center (https://spdf.gsfc.nasa.gov/index.html). The IRI-2016 and NeQuick-2 models are publicly available at http://irimodel.org and https://t-ict4d.ictp.it/nequick2/ nequick-2-web-model, respectively.