2.1. Dataset description and input parameters

The importance of geomagnetic and ionospheric indexes is emphazised in the introduction section. They provide valuable information about the complex interactions between the Earth’s magnetic field and the ionosphere and aid in improving the accuracy of ionospheric predictions. The geomagnetic and ionospheric information was downloaded from NASA https://omniweb.gsfc.nasa.gov/ow.html (accessed on 15 January 2023) (https://ccmc.gsfc.nasa.gov/) and Data Analysis Center for Geomagnetism and Space Magnetism https://wdc.kugi.kyoto-u.ac.jp/dstdir/ (accessed on 15 January 2023) (https://wdc.kugi.kyoto-u.ac.jp/dstdir/).

The input parameters of the ionospheric foF2 forecasting model, related to the ionospheric variability in time, space, solar, and geomagnetic activity, are described as follows:

1. 1. To account for the influence of the number of days in a year on ionospheric $ foF2 $ variations (Chen et al., Reference Chen, Liu and Chen2000), (Myles et al., Reference Myles, Feudale, Liu, Woody and Brown2004), a day number ( $ DN $ ) is converted into sine and cosine components as follows:

(2.1) $$ DNS=\sin \left(\frac{2\pi \times DAY}{365}\right) $$

(2.2) $$ DNC=\cos \left(\frac{2\pi \times DAY}{365}\right) $$

1. 2. The universal time altered the ionospheric storm occurrence, and adjustments are made accordingly for different times (Wintoft & Cander, Reference Wintoft and Cander2000).

(2.3) $$ UT S=\sin \left(\frac{2\pi \times UT}{24}\right) $$

(2.4) $$ UT C=\cos \left(\frac{2\pi \times UT}{24}\right) $$

1. 3. The solar zenith angle is transformed using the provided formulas to calculate the sine ( $ CHIS $ ) and cosine ( $ CHIC $ ) components:

(2.5) $$ CHI S=\sin \left(\frac{2\pi \times CHI}{360}\right) $$

(2.6) $$ CHI C=\cos \left(\frac{2\pi \times CHI}{360}\right) $$

1. 4. The $ ap $ index uses a time-weighted accumulation series calculated from the geomagnetic planetary index $ ap $ (Wrenn, Reference Wrenn1987). This study uses a constant factor of attenuation multiplier, called $ \tau $ , ranging from 0 to 1. In this article, $ \tau $ is set to 0.8. The initial value of the magnetic index is $ {ap}_0 $ , while $ {ap}_{-1} $ , $ {ap}_{-2} $ , and $ {ap}_{-11} $ represents the values 3 hours before, 6 hours before, and 33 hours before, respectively.

(2.7) $$ ap\left(\tau \right)=\left(1-\tau \right)\sum \limits_{i=0}^{11}{\left(\tau \right)}^i\times {ap}_{-i} $$

1. 5. Based on previous research, there is a significant correlation between the geomagnetic index $ Dst $ , $ Kp $ and $ foF2 $ (Kutiev & Muhtarov, Reference Kutiev and Muhtarov2001). As a result, we also consider $ Dst $ and $ Kp $ as input parameters similar to $ ap $ and define them as follows:

(2.8) $$ Dst\left(\tau \right)=\left(1-\tau \right)\sum \limits_{i=0}^{11}{\left(\tau \right)}^i\times {Dst}_{-i} $$

1. 6. The first differences in the $ Kp $ values are taken to achieve stationarity of the $ Kp $ index. This involves computing the difference between each $ Kp $ value $ \left( Kp(t)\right) $ and its previous value $ \left( Kp\left(t-1\right)\right) $ , denoted as $ \Delta Kp(t)= Kp(t)- Kp\left(t-1\right) $ .

(2.9) $$ \Delta Kp(t)={\beta}_0+{\beta}_1\Delta Kp\left(t-1\right)+{\beta}_2\Delta Kp\left(t-2\right)+\dots +{\beta}_{18}\Delta Kp\left(t-18\right)+\varepsilon (t) $$

Here, $ \Delta Kp(t) $ refers to the first differences of the $ Kp $ index at time $ t $ and $ \Delta Kp\left(t-1\right) $ to $ \Delta Kp\left(t-18\right) $ represents the lagged first differences of the $ Kp $ index up to 18-time steps. The coefficients $ {\beta}_0 $ to $ {\beta}_{18} $ in the model capture the impact of these lagged differences on the current value $ \Delta Kp(t) $ , and $ \varepsilon (t) $ represents the error term. To estimate the coefficients $ \varepsilon (t) $ , we apply the least squares method using a matrix ( $ X $ ) containing the lagged first differences and a vector ( $ Y $ ) with the first differences of the $ Kp $ index. The formula to estimate the coefficients is $ \beta ={\left({X}^{\prime }X\right)}^{\left(-1\right)}{X}^{\prime }Y $ . In this work, set $ T $ = 18 hours (Wang et al., Reference Wang, Shi, Wang, Zherebtsov and Pirog2008).

1. 7. The characteristics of the ionospheric $ F $ layer are affected by the wind in the thermosphere. According to the vertical ion drift equation (Oyeyemi et al., Reference Oyeyemi, Poole and McKinnell2005), this is the case.

(2.10) $$ W=U\times \sin \left(\theta -D\right)\times \cos I\times \sin I $$

$ W $ represents the vertical ion drift velocity, $ U $ represents the horizontal wind velocity, and $ \theta $ is the geographic azimuth angle. Additionally, we denote magnetic declination and inclination as $ D $ and $ I $ , respectively. To decompose $ D $ into sine and cosine components, we use $ DS $ and $ DC $ . Similarly, $ I $ is transformed into its sine component, $ IS $ .

(2.11) $$ DS=\sin \left(\frac{2\pi \times D}{360}\right) $$

(2.12) $$ DC=\cos \left(\frac{2\pi \times D}{360}\right) $$

(2.13) $$ IS=\sin \left(\frac{2\pi \times I}{360}\right) $$

1. 8. According to studies on the $ F2 $ layer, $ foF2 $ values are significantly influenced by previous hour values (Oyeyemi et al., Reference Oyeyemi, Poole and McKinnell2005). Our model uses the relative deviation formula to account for this dependence.

(2.14) $$ {\Delta}_1(h)={f}_0{F}_0(h)+{f}_0{F}_2\left(h-1\right) $$

(2.15) $$ \Delta R(h)=\frac{\Delta_1(h)}{f_0{F}_2(h)} $$

1. 9. The sunspot number, also known as $ Rz $ , is widely used as an input in various ionosphere models, such as the International Radio Consultative Committee $ foF2 $ model (Wrenn, Reference Wrenn1987). This is because $ Rz $ can effectively map the ionospheric $ foF2 $ response to changes in solar output. Therefore, we have included $ Rz $ as an input factor in our models.

This research focuses on ionospheric $ foF2 $ critical frequency forecasting at mid and high latitudes, using hourly time-series data from 5 stations worldwide (Table 1 and Figure 1). The data is collected for 2012, 2015, and 2009, 2019, as these years represent the maximum and minimum solar activity, correspondingly (Barta et al., Reference Barta, Sátori, Berényi, Kis and Williams2019). Solar activity variations of 11 years are essential for decadal variations in the solar-terrestrial environment (Gnevyshev, Reference Gnevyshev1977; Ohmura, Reference Ohmura2009). This article investigates the critical frequency variations based on the 24th solar cycle from 2008 to 2019.

Table 1. The geographical locations and coordinates of the selected stations

Figure 1. A map of the stations’ spatial distribution, marked by red dots.

2.2. The developed LSTM model

Working with time-series data requires understanding a system’s dynamics, such as its periodic cycles, how the data changes over time, its regular variations, and its sudden changes. LSTM networks are recurrent neural networks designed to avoid the long-term dependency problem by learning order vulnerability in sequence from time-series data (Hochreiter & Schmidhuber, Reference Hochreiter and Schmidhuber1997).

The LSTM network consists of the following stages: the forget gate, cell state, input gate, and output gate (Figure 2). One of the essential properties of the LSTM is its ability to memorize and recognize information that enters the network and to discard information that is not required by the network to learn and make predictions (Yu et al., Reference Yu, Si, Hu and Zhang2019). The forget gate determines whether or not information can pass through the network’s layers. It expects two types of input: the information from the previous layers and information from the current layer.

Figure 2. The LSTM block diagram, the repeating module in LSTM, with four interacting layers.

Furthermore, LSTMs have a chain structure with four neural network layers interacting uniquely (Gonzalez & Yu, Reference Gonzalez and Yu2018). The LSTM algorithm relies on the cell state, represented by the horizontal line running across the top of the diagram, and is crucial in information transmission (Gers et al., Reference Gers, Schmidhuber and Cummins2000). The LSTM can remove or add information to the cell state, which is carefully controlled by gate structures. Gates allow information to pass through if desired (Yu et al., Reference Yu, Si, Hu and Zhang2019). The “forget gate layer” is a sigmoid layer that removes information that is no longer relevant or useful for further predictions. It examines $ {h}_{t-1} $ (a hidden state at the timestamp $ t-1 $ ) and $ {x}_t $ (the input vector at the timestamp $ t $ ) and returns a number between 0 and 1 for each number in the cell state $ {C}_{\left(t-1\right)} $ . Number 1 indicates “completely keep,” while number 0 indicates “completely remove” where $ {h}_t $ represents a hidden state at the current timestamp $ t $ . The updated cell from the cell state is passed to the $ \tanh $ , an activation function, which is then multiplied by the output state’s sigmoid function. After calculating the hidden state at the timestamp $ t $ , the value is returned to the recurrent unit and combined with the input at the timestamp $ t+1 $ . The same procedure is repeated for $ t+2,t+3,\dots, t+n $ timestamps until the desired number $ n $ of timestamps is reached.

(2.16) $$ {f}_t=\sigma \left({x}_t\times {u}_t+{h}_{t-1}\times {w}_t\right) $$

(2.17) $$ {\hat{c}}_t=\tanh \left({x}_t\times {u}_c+{h}_{t-1}\times {w}_c\right) $$

(2.18) $$ {i}_t=\sigma \left({x}_t\times {u}_i+{h}_{t-1}\times {w}_i\right) $$

(2.19) $$ {o}_t=\sigma \left({x}_t\times {u}_o+{h}_{t-1}\times {w}_o\right) $$

(2.20) $$ {c}_t={f}_t\times {u}_{t-1}+{i}_t\times {\hat{c}}_t $$

(2.21) $$ {h}_t={o}_t\times \tanh \left({c}_t\right) $$

Where: $ {x}_t $ is the input vector, $ {h}_{t-1} $ is the previous cell output, $ {C}_{t-1} $ is the previous cell memory, $ {h}_t $ is the current cell output, $ {c}_t $ is the current cell memory and $ w $ , $ u $ are weight vector for forget gate ( $ f $ ), candidate ( $ c $ ), $ i/p $ is the gate ( $ i $ ) and $ o/p $ is the gate ( $ o $ ).

2.3. Network configuration

Figure 3 depicts the input parameters of the LSTM model. The developed LSTM model has 14 input parameters fed into the LSTM neural network, as described in Section 2.1. The proposed model has 14 input variables connected to 36 LSTM units, 12 units of dense layer, and one output layer. The objective of the LSTM model is to extract significant and representative features from the historical data. Various experiments are conducted using different hyperparameters and LSTM units to obtain an optimal architecture, taking inspiration from previous research studies (Reimers & Gurevych, Reference Reimers and Gurevych2017). A developed configuration for the LSTM model is derived after a meticulous evaluation process. The activation function is set to $ Relu $ to handle non-linear relationships effectively (Yadav et al., Reference Yadav, Jha and Sharan2020). A batch size of 15 is found to balance efficiency and capture meaningful patterns. Through iterative fine-tuning of parameters and learning, the model is trained for 200 epochs. This LSTM configuration has demonstrated superior performance and is anticipated to produce accurate predictions for the given task.

Figure 3. A block diagram of the proposed LSTM model with input parameters.

2.4. Output parameters

The current study focuses primarily on the storm-time $ foF2 $ forecast. As output parameters, 1, 2, 6, 12, 18, and 24 hours ahead prediction results are demonstrated for years with maximum and minimum solar activity.

2.5. The machine learning process

Figure 4 presents a comprehensive overview of the methodology utilized in constructing a machine-learning model. The process entails sequential steps, including pre-processing, normalization, model training, model testing, and model evaluation (El Naqa & Murphy, Reference El Naqa and Murphy2015). The article uses the temporal segmentation technique to address seasonal changes in time-series data, which involves dividing the data into segments corresponding to different months or recurring temporal patterns. Segmenting data by several months and time ensures that each season is equally represented, reducing the risk of skewed model performance. Additionally, splitting the data optimizes model performance by aligning with the network’s strengths, mitigating biases, enhancing generalization, and facilitating comprehensive evaluation (Lovrić et al., Reference Lovrić, Milanović and Stamenković2014).

Figure 4. A generalized diagram depicting the process of developing a machine learning model.

The dataset is divided into three subsets: 80% for training, 10% for validation, and 10% for testing. However, a non-shuffled approach is adopted to preserve temporal information and mitigate potential seasonal bias errors. It is important to note that a non-shuffled approach preserves temporal information. To improve the accuracy of our methodology, we use a refined approach. Firstly, the data is divided into four separate segments, each representing a three-month period. These segments are then assigned to training, validation, and test subsets based on predetermined percentages. This process is repeated every three months, ensuring that the training, validation, and test datasets are evenly distributed. The reason for using this approach is based on the robust capabilities of LSTM networks in dealing with sequential data. This method capitalizes on LSTM’s ability to capture temporal patterns, thus enhancing model performance. Essentially, the sequential partitioning method not only utilizes the power of LSTM but also effectively eliminates any potential seasonal distortions. This meticulous approach results in improved overall quality in model development and outcomes.

Data pre-processing is an important step that involves cleaning, transforming, normalizing, and organizing raw data into a format suitable for analysis (Károly et al., Reference Károly, Galambos, Kuti and Rudas2020; Rinnan et al., Reference Rinnan, Berg, Thygesen, Bro and Engelsen2009; Rinnan et al., Reference Rinnan, Berg, Thygesen, Bro and Engelsen2020). Pre-processing’s primary function is handling missing data and inconsistencies in raw data. After clean data is performed, selected parameters are used as model inputs; parameter preparation is described in Section 2.1.

The time-series data, which is associated with ionospheric deviations, is inputted into the network. The input sequence is gradually shifted to predict the ionospheric behavior at various future time points. This involves incrementally shifting the data from 1 hour up to 24 hours for each testing instance. The result of this shifting procedure is a set of input-output pairs, where the input consists of past observations up to a particular time, and the output is the prediction for a specific time point in the future. The generalized process can be described as follows:

(2.22) $$ foF{2}_{\left(t+i\right)}=f\left({UTS}_t,{UTC}_t,{HIS}_t,{HIC}_t,{DNS}_t,{DNC}_t,{DS T}_t,{ap}_t,{Kp}_t,{DS}_t,{DC}_t,{IS}_t,{Rz}_t,\Delta {R}_t\right) $$

where $ t $ indicates time and $ i $ depicts prediction hour (up to 24 hours).

The mean absolute percentage error ( $ MAPE $ ), root mean square error ( $ RMSE $ ) and mean absolute error ( $ MAE $ ) were chosen to measure the predictive accuracy of the developed model. (Jakowski et al., Reference Jakowski, Stankov, Schlueter and Klaehn2017). The $ foF2 $ forecasting performance was compared to the other three algorithms. The $ MAPE $ , $ RMSE $ , and $ MAE $ are defined as follows:

(2.23) $$ MAPE=\frac{1}{n}\sum \limits_i^n\left|\frac{\ foF{2}_{act}- foF{2}_{pred}}{foF{2}_{act}}\times 100\right| $$

(2.24) $$ RMSE=\sqrt{\sum \limits_i^n\frac{{\left(\ foF{2}_{act}- foF{2}_{pred}\right)}^2}{N}} $$

(2.25) $$ MAE=\sum \limits_i^n\left|\frac{\ foF{2}_{act}- foF{2}_{pred}}{N}\right| $$

where n is the number of summation iterations, $ foF{2}_{act} $ is the actual value, $ foF{2}_{pred} $ is the predicted value. A low accuracy rate value indicates that the model’s predictions are close to the actual values, whereas a high accuracy rate value indicates that the model’s predictions are far off.

2.6. Compared models

To evaluate the developed model performance, LSTM was compared to linear regression, decision trees, and MLP algorithms. Linear regression is a widely used method for data analysis (Montgomery et al., Reference Montgomery, Peck and Vining2021). It is good at finding the linear relationship in data sets by establishing the relationship between independent and dependent variables when performing a regression task (finding the best-fitting line). It denotes the existence of a linear relationship between dependent and independent variables. In a linear regression model, the dependent variable is modeled as a linear combination of the independent variables plus an error term (Weisberg, Reference Weisberg2005). The equation is written as follows:

(2.26) $$ y={\beta}_0+{\beta}_1\times {x}_1+{\beta}_2\times {x}_2+\cdots +{\beta}_n\times {x}_n+\varepsilon $$

The dependent variable is y, and the independent variables are $ {x}_1,{x}_2,\cdots, {x}_n $ . The coefficients ( $ {\beta}_0,{\beta}_1,{\beta}_2,\cdots, {\beta}_n $ ) describe the strength of each independent variable’s relationship with the dependent variable. $ \varepsilon $ is the error, representing the difference between the predicted and actual value of the dependent variable.

Decision Tree (DT), an effective tool that employs a tree-like model of decisions and their potential outcomes, was implemented for $ foF2 $ critical frequency forecasting up to 24 hours in advance (Iban & Şentürk, Reference Iban and Şentürk2022). A series of branches represent possible decisions and outcomes (Myles et al., Reference Myles, Feudale, Liu, Woody and Brown2004). A decision tree model’s prediction is based on the probabilities or scores associated with the leaf nodes where the data ends up (Lan et al., Reference Lan, Zhang, Jiang, Yang and Zha2018).

MLP was developed for $ foF2 $ ionospheric disturbances forecasting, composed of multiple layers of interconnected “neurons.” MLP’s ability to have many hidden layers between the input and output increases complexity and density. The MLP models have been used to forecast air quality, daily solar radiation, and solar activity (Elizondo & McClendon, Reference Elizondo and McClendon1994). A comprehensive description of the developed Decision Tree and MLP models, along with their corresponding block diagrams and hyperparameters, is provided in Supplementary Material.