Prediction of ionospheric total electron content data using NARX neural network model

ABSTRACT


INTRODUCTION
The signals from the global positioning system (GPS) satellites orbiting round the earth propagate through the ionosphere.Ionosphere is the ionized region of the earth's upper atmosphere ranging from about 60-1,000 km which contains significant amount of free electrons.The ionosphere is known for its dispersive nature and thus delays the electromagnetic signals and affects precise measurements in satellite navigation and communication.This influence is directly proportional to the density of electrons between the satellite transmitter and the receiver which quantifies the total electron content (TEC).Ionospheric TEC exhibits complex spatial and temporal fluctuations along with certain solar and geomagnetic related ionospheric disturbances [1].The daily diurnal variations, seasonal variations, 27-day variations, 11-year solar cycle variations, geographic latitude, and some events like geomagnetic storms which generate sudden disturbances in the ionosphere are some obvious variations of ionospheric TEC [2].Prediction of ionospheric TEC can thus help in improving the performance of GPS navigational systems.
There are a number of empirical models like international reference ionosphere (IRI) [3], NeQuick [4], Klobuchar [5] which are available for reconstructing the state of the ionosphere.These models effectively predict the state of the ionosphere during quiet times but are not very accurate during periods of Bulletin of Electr Eng & Inf ISSN: 2302-9285  Prediction of ionospheric total electron content data using NARX neural network model (Nayana N. Shenvi) 549 high ionospheric activity.Several time-series models developed for prediction of TEC include fourier series [6], discrete cosine transform [7], and autoregressive distributed lag [8].However, these models are not effective in predicting the non-stationary fluctuations present in the ionosphere.During last few decades, neural network techniques have been applied to predict the ionospheric TEC as well as model the various parameters in the ionospheric study [9]- [11].These models consider the solar and geomagnetic factors in improving the prediction accuracy.Neural network models are found to outperform statistical methods.However, these models do not consider the temporal feature of TEC [12].The TEC data exhibits non-stationarity which implies that the behavior of the ionosphere is extremely complex and nonlinear [13].nonlinear autoregressive with exogenous inputs (NARX) neural network has been widely used for modeling non-linear dynamic systems [14].Hence, in this paper we are proposing a NARX neural network model fed with time of the day, solar and geomagnetic indices as exogenous parameters to predict the TEC.Literature suggests number of applications modeled by NARX neural network.NARX model was tested on chaotic laser time series and traffic time series and it outperformed standard neural network architectures like time delay neural network (TDNN) and Elman architectures [14].NARX gave better results over autoregressive integrated moving average (ARIMA) for wind speed prediction [15].NARX was successfully tested for prediction of rainfall and showed improved results compared to linear algorithms like support vector machines [16].NARX was used for traffic prediction and had superior prediction accuracy to RNN models [17].The ability to predict symmetric horizontal (SYM-H) storm time index using NARX showed better accuracy compared to back propagation and recurrent Elman networks [18].NARX has also been used for TEC prediction using time, day of the year and season as exogenous inputs [1].
In this study, we discuss prediction of TEC using NARX neural network model using exogenous parameters such as solar index F10.7 and geomagnetic indices such as disturbed storm (Dst), Kp, and Ap along with local time.We have studied TEC prediction at different latitudinal regions namely equatorial (0 0 -10 0 N), low-latitude (10 0 -30 0 N), mid-latitude (30 0 -60 0 N), and high-latitude (60 0 -80 0 N), during both the minimum solar activity year (2008) and maximum solar activity year (2014).We compare the results of our NARX model (henceforth called NARX1) with the NARX model developed for TEC prediction by Guoyan et al. [1] (henceforth called NARX2).

METHOD 2.1. Data and Input parameters
The ionospheric TEC data is downloaded from IONOLAB (www.ionolab.org)[19].To explore latitudinal influence on the prediction of ionospheric TEC, we selected four locations in the northern region, each situated at different latitudes.The selected stations are qaq1 (Greenland and Denmark), baie (Baie-Comeau and Canada), mas1 (Maspalomas and Spain), and bogt (Bogota and Columbia).The latitude and longitude co-ordinates of these stations is shown in Table 1 The exogenous parameters for the NARX model were selected based on the factors influencing TEC.We calculated the Pearson correlation coefficient for various exogenous parameters and selected the ones that correlate well with TEC data.Table 2 shows the values of Pearson correlation coefficient for the various exogenous parameters with TEC data.From the table it is clear that Ap and Kp show high positive correlation, Dst shows high negative correlation and F10.7 shows moderate correlation with TEC.Time of the day was selected since TEC exhibits diurnal characteristics.Thus, we selected solar index F10.7 and geomagnetic indices Dst, Kp and Ap along with time of the day.The exogenous paramters were obtained from the space physics data facility at NASA's Goddard Space Flight Center, accessible through their website at https://omniweb.gsfc.nasa.gov/form/dx1.html[20].The TEC data along with the solar and geomagnetic indices were sampled every one hour.We used the data from 1 st July, 2008 to 26 th July, 2008 (total samples=26 days*24= 624) for training the model and then tested the prediction capability of the model for unseen data from 27 th July, 2008 to 31 st July, 2008 (total samples= 5 days*24=120).The same period was selected for training and testing the model for solar maximum year 2014.The various exogenous parameters used are as follows: -Local time: TEC exhibits diurnal variation because the electron density progressively rises to a maximum value at noon before falling to a minimum value at midnight.We compute sin ( ) ℎ and cos ( ) ℎ where local time, ℎ ranges from 0-23.The sin and cos functions help in normalizing the time input in the range (-1,1) [1].-Solar activity: solar activity significantly influences TEC values, with lower TEC during solar minimum years and higher TEC during solar maximum years [21].F10.7, which correlates directly with the number of sunspots and solar radiation in the ultraviolet and visible spectrum, serves as a valuable predictor of solar activity [22].-Geomagnetic activity (Dst, Kp, and Ap index): the Dst, Kp, and Ap represent geomagnetic activity.
-Dst index: the Dst is a metric utilized to gauge the intensity of geomagnetic storms.It quantifies the deviation of the Earth's horizontal magnetic field component from its typical variation during calm days [23].Under tranquil conditions, the Dst ranges from 0 to -50 nT, but during the occurrence of a powerful storm, it can surpass -200 nT.-Kp index: the planetary index Kp is the average standardized K-index obtained from 13 geomagnetic observatories between 44 0 and 60 0 in either hemisphere [21].This parameter serves as a measure of geomagnetic activity.The higher the Kp index, more severe the geomagnetic activity, and greater the impact on the TEC.-Ap index: the Ap index, also known as the planetary index, is a numerical scale that ranges from 0 to 400, with higher values signifying more intense geomagnetic activity.This index serves as a valuable tool for monitoring the effects of space weather on the earth's magnetic field and assessing potential impacts on critical systems such as power grids, satellites, and radio communications [21].

Architecture of NARX1 model
A NARX neural network is a dynamic model that incorporates recurrent feedback connections from the output layer to the input layer with time-delayed units.This unique architecture makes NARX networks an efficient tool for modeling and validation purposes, as they demonstrate better generalization and faster convergence compared to other conventional neural network models [24].This can model nonlinear timeseries systems and can be mathematically represented by:

551
The NARX neural network, like a multilayer perceptron (MLP), comprises an input layer, a hidden layer, and an output layer in its architecture.However, the key difference is that the NARX network includes recurrent feedback connections from the output layer to the input layer with time-delayed units, enabling it to model dynamic systems and time-dependent data more effectively.The past exogenous inputs along with the past predicted TEC values are fed at the input of the model.The feedback connections allow the network to take into account the previous behavior of the system and to capture any patterns that may exist between the input and output variables over time.Figure 2 depicts the framework of NARX neural network.The input layer takes in the exogenous inputs along with past predicted TEC values.The output of the intermediate hidden layer   () is given by: where, each hidden layer node is activated by the function F 1 .The input neuron x(t-p) is connected to the i th hidden layer node with a connection weight of w ip and the output feedback neuron y(t-q) is also connected to the i th hidden layer node with a connection weight of   .Additionally, each hidden layer node is associated with a bias term a i .These connection weights and bias terms are learned during the training process and they determine the strength and direction of the connections between the neurons in the network [25].The predicted TEC value at the output y ̂ (t+1) is given by: where w i is the connection weight between the i th hidden layer node and output layer, b i is the bias of the network,  is the number of hidden layer nodes.Therefore, the predicted TEC value at time t+1 is a function of the weighted sum of the inputs passed through a nonlinear activation function and is determined by the weights and biases of the neural network.In the open-loop architecture, as depicted in Figure 3(a), the NARX neural network predicts the output TEC y ̂ (t+1) using historical values of both the exogenous data and the actual past TEC values.This configuration enables the model to utilize accurate information from the actual TEC time series, contributing to more precise inputs [26].In the closed-loop configuration, as depicted in Figure 3 y ̂(t+1)=f (y(t),y(t-1),…,y(t-n y ), x(t), x(t-1),…,x(t-n x )) (4) and for the closed-loop architecture; y ̂(t+1)=f (y ̂(t),y ̂(t-1),…,y ̂(t-n y ),x(t), x(t-1),…,x(t-n x )) (5) where f (.) is the mapping function of the neural network, y ̂(t+1) is the predicted TEC output of the NARX using either the past true TEC values y(t),y(t-1),…,y(t-n y ) or the past predicted TEC values y ̂(t), y ̂(t-1),…,y ̂(t-n y ) and the past values of the exogenous time series x(t),x(t-1)…,x(t-n x ) and   and   are the input and output delays.

Implementation of NARX1 model in MATLAB
We have simulated the NARX model using the neural net toolbox in MATLAB.The simulated model is shown in Figures 4 and 5.The number of hidden layer neurons was set to 20.The delays n x and n y was set to two.These hyperparameters were finalized after multiple simulations.The model was run with different number of network feedback delays and then finalized to two where the prediction model performed well.The neurons in the hidden layer determine the accuracy of the model.Less neurons in the hidden layer fail to establish relationship of the data but increasing the number of neurons leads to more computational time [27].Initially we trained the model in open loop configuration as shown Figure 4.The advantage of doing this is that the NARX model uses the previous values of the actual TEC time series and hence it receives more precise input during the training phase.Then we converted the model into closed loop for making predictions into the future as shown Figure 5.We perform recursive multi-step forecasting in which we predict one-step TEC value and use this predicted value in the next step to make next prediction.This process is continued till we forecast TEC values for the next 5 days.Training stops when error starts increasing indicating that there is no further improvement in generalization.The output provides the predicted TEC values.The NARX model was developed using a hyperbolic tangent function in the hidden layer and a linear transfer function in the output layer.We calculate the root mean square error (RMSE), mean absolute error (MAE), correlation coefficient (r), and symmetric mean absolute percentage error (sMAPE) using the following formulas to evaluate the performance of the model.
where y i ̂ is the predicted TEC by the model and y i is the actual TEC value,  ̅  is the mean of the actual TEC values, y ̅i is the mean of the predicted TEC values and  is the frequency of observations.An ideal prediction model is indicated by only one non-zero value in the autocorrelation function, which occurs at zero lag, signifying a perfect match between the model's predictions and the actual data [28].The predicted errors are below 95% confidence limit (red dotted line) and thus are uncorrelated with each other.This also shows that the network has been trained satisfactorily and can be used for prediction.The error histogram plot in MATLAB is a graphical representation of the distribution of errors in a dataset.The plot shows the frequency of occurrence of errors, binned into intervals, with the x-axis representing the difference between predicted and actual TEC values and the y-axis representing the instances of those errors.Figure 7 shows a well-behaved error histogram plot which has a roughly symmetric distribution around the mean error value, with very few or no outliers.The plot clearly indicates that maximum samples from the training, validation and test dataset have a very small error equal to 0.01736.The regression plot denotes an equation between the predicted TEC (output) and true (target) value of TEC. Figure 8 shows the regression line plot for baie station during training, validation and test phases during the solar maximum year 2014.Table 3 shows the slope and intercept values for the regression plot shown in Figure 8.The regression equation relating the predicted and the true TEC is shown on the y-axis.The coefficient of the equation shows proportionality between predicted TEC and the true TEC and is very close to 1.The constant term represents an error which should be added to the scaled target to make it as close as possible to the predicted output and is very close to 0 (see Figure 8).Overall, high values of R were achieved, specifically 99.75% for the training set, 99.04% for the validation set, 99.16% for the test data set and overall 99.54%.R values greater that 95% suggest excellent fit between true and predicted values [27].Figures 6-8 clearly indicate that NARX network has been trained efficiently and can be effectively used for prediction of TEC data.The validity of our model's results could not be meaningfully compared to existing literature due to variations in the methodologies used in each study.These differences include the use of different sets of stations, varying amounts of training data, and use of different solar minimum and active years.To accurately compare our results and validate our model, we simulated the NARX2 model using the identical architecture and parameters as described by Guoyan et al. [1] using time of the day, day of the year and season as exogenous inputs.

For the solar quiet year 2008
We check the prediction capability of our model by testing it with unseen data from 26 th July to 31 st July, 2008.The plots in Figure 9 show the observed and the predicted TEC at Figure 9

Solar maximum year 2014
We also check the prediction capability of our model by testing it with unseen data from 26 th July to 31 st July, 2014 for the solar maximum year 2014 (see Figure 10).The plots in Figure 10 show the observed and the predicted TEC at Figure 10

Comparison with NARX2 model
To validate our model, we compared our results with NARX2.NARX1 model was implemented using time, solar index (F10.7),and geomagnetic indices (Dst, Kp, and Ap).The literature used NARX2 model for TEC prediction using time of the day, day of the year and season as exogenous inputs.This model was implemented for China region for low and mid latitude regions for the solar maximum year 2011 and quiet year 2017.To ensure a fair comparison, we implemented NARX2 for our regions and also for the same solar minimum (2008) and active years (2014).Using the same model architecture (NARX), we wanted to see the improvement in prediction accuracy using solar and geomagnetic indices as exogenous parameters at all latitudinal regions and for both the solar minimum and active years.The RMSE, MAE, correlation coefficient and sMAPE values obtained for each of the station for both the solar minimum and active years using both NARX1 and NARX2 model are plotted in Figures 11(a As anticipated, the RMSE increases as latitude decreases, implying better prediction performance in higher and mid-latitude regions compared to lower and equatorial regions.The challenging aspect of modeling and predicting ionospheric TEC in low-latitude regions is attributed to the presence of the equatorial ionization anomaly [29].This anomaly introduces complex and variable ionospheric behavior, making accurate modeling and prediction more difficult in these specific areas.The RMSE and MAE values for solar maximum year are higher compared to solar minimum year.This is to increased solar radiation, ionization, and disturbances in the ionosphere during solar maximum year and is consistent with the literature.The RMSE, MAE, correlation coefficient and sMAPE metrics show that

557
NARX1 model can predict TEC with better accuracy even for solar maximum year.According to our study findings, the NARX1 model consistently demonstrated better predictive performance for TEC values compared to the NARX2 model.

CONCLUSION
This study examined the use of NARX neural network for prediction of TEC data during the solar minimum and active years using time of the day, solar, and geomagnetic indices as exogenous parameters.The NARX model's ability to capture nonlinear dynamics, incorporate autoregressive terms, and include exogenous inputs make it a powerful tool for predicting TEC data.The error analysis based on RMSE, MAE, correlation coefficient and sMAPE indicate that use of solar index (F10.7)and geomagnetic indices (Dst, Kp, and Ap) as exogenous data helps in effective prediction of TEC.These parameters deviate from their normal values with changes in the ionospheric TEC.Any fluctuations in TEC due to disturbance in the ionosphere is quickly captured by these exogenous parameters and are effective in TEC prediction.But, NARX2 uses time of the day, day of the year and season as exogenous parameters.These parameters cannot capture the sudden ionospheric disturbances and are only useful in TEC prediction during quiet stable ionospheric conditions.Further, the TEC data, the solar and the geomagnetic data are all used in raw form without normalizing indicating the suitability of this technique for real time TEC prediction.The prediction of TEC at different latitudinal regions is also consistent with the literature and shows increasing RMSE with decreasing latitudes.The error autocorrelation function, error histogram and regression line plots also indicate the best fit of the model for TEC prediction.The RMSE between 0.8212 TECU to 4.0277 TECU for the solar minimum year and between 1.8952 TECU to 6.9822 TECU for the solar maximum year indicate the suitability of this model and correct selection of the exogenous parameters for ionospheric TEC prediction.All the evaluation metrics clearly showed that the NARX1 model outperformed the NARX2 model in accurately predicting TEC values.

Figure 1 .
Figure 1.NARX model with tapped delay line at the input and output,  −1 represents unit time delay NARX neural network can be used in open-loop mode and closedloop mode.The open-loop architecture and the closed-loop architecture of NARX neural network model is shown in Figures 3(a) and (b) respectively.

Figure 2 Figure 3 .
Figure 2. NARX neural network structure (b), the NARX model feeds back the predicted TEC output to the input of the feedforward neural network.Consequently, the  ISSN: 2302-9285 Bulletin of Electr Eng & Inf, Vol. 13, No. 1, February 2024: 548-558 552 prediction of y ̂ (t+1) is based on the previous values of the exogenous data and the previously predicted TEC values.This closed-loop architecture allows the model to exploit its own predictions to make further estimations, potentially capturing dynamic dependencies and time-dependent patterns effectively.The mathematical expression for the NARX open-loop architecture is given by:

Figure 4 .
Figure 4. NARX open loop architecture simulated in MATLAB during the training phase

1 .
Stability of the modelAfter training the model we check the stability of the model by plotting the error autocorrelation function, error histogram, and regression coefficient.
Figure 6 shows the error autocorrelation plot for qaq1 station during the solar maximum year 2014.The error autocorrelation function is a valuable tool for validating the performance of a model.It helps describe the relationship between prediction errors over time, revealing how errors are correlated at different time lags.

Figure 6 .
Figure 6.Error autocorrelation plot for qaq1 station during solar active year 2014

Figure 9 .
Figure 9. Observed and predicted TEC using NARX model for the period from 26 th July 2008 to 31 st July 2008 for; (a) qaq1, (b) baie, (c) mas1, and (d) bogt stations during the solar minimum year 2008

Figure 10 .
Figure 10.Observed and predicted TEC using NARX model for the period from 26 th July 2014 to 31 st July 2014 for; (a) qaq1, (b) baie, (c) mas1, and (d) bogt stations during the solar maximum year 2014

Table 1 . Selected stations in the northern hemisphere
.

Table 3 .
Slope and intercept values for the regression plot in Figure8