Predicting Solar Energetic Particles Using SDO/HMI Vector Magnetic Data Products and a Bidirectional LSTM Network

As mentioned in Section 1, we aim to solve the following two binary prediction problems. [FC_S problem] Given a data sample x_t at time point t in an AR where the AR will produce an M- or X-class flare within the next T hours of t and the flare initiates a CME, we predict whether x_t is positive or negative. Predicting x_t to be positive means that the AR will produce an SEP event associated with the flare/CME. Predicting x_t to be negative means that the AR will not produce an SEP event associated with the flare/CME. [F_S problem] Given a data sample x_t at time point t in an AR where the AR will produce an M- or X-class flare within the next T hours of t regardless of whether or not the flare initiates a CME, we predict whether x_t is positive or negative. Predicting x_t to be positive means that the AR will produce an SEP event associated with the flare. Predicting x_t to be negative means that the AR will not produce an SEP event associated with the flare. For both of the two binary prediction problems, we consider T ranging from 12 to 72 in 12 hr intervals as frequently considered in the literature (Ahmed et al. 2013; Bobra & Ilonidis 2016; Inceoglu et al. 2018; Liu et al. 2020).

In solving the two binary prediction problems, we first show how to collect and construct positive and negative data samples used in our study. Figure 1(a) (Figure 1(b), respectively) illustrates how to construct positive (negative, respectively) data samples for the FC_S problem where T = 24 hr. Refer to the FC_S database described in Section 2, which indicates whether a flaring AR that already produces an M- or X-class flare/CME will initiate an SEP event associated with the flare/CME. For the flaring AR, we collect data samples that are within the T = 24 hr prior to the peak time of the flare.

• 1.  
  If the flare/CME are associated with an SEP event, the data samples belong to the positive class and are colored (labeled) by blue as shown in Figure 1(a). Thus, for each blue (positive) data sample, there is an M- or X-class flare that is within the next 24 hr of the occurrence time of the data sample, the flare initiates a CME, and the flare/CME are associated with an SEP event.
• 2.  
  If the flare/CME are not associated with an SEP event, the data samples belong to the negative class and are colored (labeled) by green as shown in Figure 1(b). Thus, for each green (negative) data sample, there is an M- or X-class flare that is within the next 24 hr of the occurrence time of the data sample, the flare initiates a CME, but the flare/CME are not associated with an SEP event.

Constructing positive and negative data samples for the F_S problem is done similarly, and its description is omitted.

Zoom In Zoom Out Reset image size

Table 1 shows the numbers of positive and negative data samples constructed for the FC_S and F_S problems respectively. Consider the FC_S problem. The positive and negative data samples are constructed based on the 31 records tagged by "yes" and 97 records tagged by "no" in the FC_S database described in Section 2. When T = 24 hr and the cadence is 12 minutes, one would expect the total number of positive data samples to be 24 hr × 60 minutes hr⁻¹ × (1/12 minutes) × 31 = 3720, and the total number of negative data samples to be 24 hr × 60 minutes hr⁻¹ × (1/12 minutes) × 97 = 11,640. However, the total number of positive (negative, respectively) data samples is 2017 (5522, respectively). This happens because we removed many data samples of low quality as described in Section 2. If a gap occurs in the middle of a time series due to the removal, we use a zero-padding strategy as done in Liu et al. (2020) to create a synthetic data sample to fill the gap. The synthetic data sample has zero values for all the 18 SHARP parameters. The synthetic data sample is added after the normalization of the SHARP parameter values, and hence the synthetic data sample does not affect the normalization procedure.

After explaining how to construct the positive and negative data samples, we now show how to solve the binary prediction problems. Consider again the FC_S problem where T = 24 hr. Here we want to predict whether a given test data sample x_t at time point t is positive (blue) or negative (green) given that there will be an M- or X-class flare within the next 24 hr of t, and the flare initiates a CME. If there is an SEP event associated with the flare/CME, and we predict x_t to be positive (blue), then this is a correct prediction as illustrated in Figure 1(c). If there is an SEP event associated with the flare/CME, but we predict x_t to be negative (green), then this is a wrong prediction as illustrated in Figure 1(e). On the other hand, if there is no SEP event associated with the flare/CME, and we predict x_t to be negative (green), then this is a correct prediction as illustrated in Figure 1(d). If there is no SEP event associated with the flare/CME, but we predict x_t to be positive (blue), then this is a wrong prediction as illustrated in Figure 1(f). The F_S problem is solved similarly. In the following subsection, we describe how to train our model and use the trained model to make predictions.

We consider one of the recurrent neural networks (RNNs) that is called LSTM (Hochreiter & Schmidhuber 1997; Goodfellow et al. 2016) to build our model. LSTM has shown good results in solar eruption prediction (Liu et al. 2019, 2020). We create a model using biLSTM. Generally, a bidirectional RNN (Schuster & Paliwal 1997) functions by duplicating the initial recurrent layer in the network to obtain two layers so that one layer uses the input as is, and the other duplicated layer uses the input in a reverse order. This design allows biLSTM to discover additional patterns that cannot be found by LSTM with only one recurrent layer (Siami-Namini et al. 2019). In addition, the data used in our study is time series, and biLSTM has shown an improvement over LSTM for general time series forecasting (Althelaya et al. 2018; Kang et al. 2020). As our experimental results show later, biLSTM also outperforms LSTM in SEP prediction.

Figure 2 presents the architecture of our neural network, which accepts as input a data sequence with m consecutive data samples. (In the study presented here, m is set to 10.) The neural network consists of a biLSTM layer configured with 400 neurons. In addition, the neural network contains an attention layer motivated by Goodfellow et al. (2016) to direct the network to focus on important information and characteristics of input data samples. The attention layer is designed to map and capture the alignment between the input and output by calculating a weighted sum for input data sequences. Specifically, the attention context vector for the output ${\hat{y}}_{i}$, denoted CV _i , is calculated as follows:

$$\begin{eqnarray}&&{{\boldsymbol{CV}}}_{i}=\sum _{j=1}^{m}{{\boldsymbol{W}}}_{i,j}{{\boldsymbol{H}}}_{j},\end{eqnarray} \tag{ 2 }$$

where m is the input sequence length, H _j is the hidden state corresponding to the input data sample x_j , and W contains weights applied to the hidden state. W is computed by a softmax function as follows:

$$\begin{eqnarray}&&{{\boldsymbol{W}}}_{i,j}=\displaystyle \frac{{e}^{{S}_{i,j}}}{\sum _{k=1}^{m}{e}^{{S}_{i,k}}}.\end{eqnarray} \tag{ 3 }$$

Here S_i,j is a score function calculated as follows:

$$\begin{eqnarray}&&{S}_{i,j}={\boldsymbol{V}}\times {\rm{\tanh }}({{\boldsymbol{W}}}^{{\prime} }({{\boldsymbol{S}}}_{i},{{\boldsymbol{H}}}_{j})),\end{eqnarray} \tag{ 4 }$$

where ${\rm{\tanh }}(\cdot )$ is the hyperbolic tangent function, S _i is the output state corresponding to the output ${\hat{y}}_{i}$, V and ${{\boldsymbol{W}}}^{{\prime} }$ are weight matrices learned by the neural network. The attention layer passes its resulting vector to a fully connected layer.

Zoom In Zoom Out Reset image size

During training, our biLSTM network takes as input overlapping data sequences where each data sequence contains m = 10 consecutive data samples. The label (color) of a training data sequence is defined to be the label (color) of the last (i.e., tenth) data sample in the data sequence while the labels (colors) of the other nine data samples in the data sequence are ignored. Thus, if the tenth data sample is positive (blue), then the training data sequence is positive; if the tenth data sample is negative (green), then the training data sequence is negative. We feed one training data sequence at a time to our biLSTM network when training the model. Figure 3(a) illustrates three positive data sequences used to train our biLSTM model. Figure 3(b) illustrates three negative data sequences used to train our biLSTM model.

Zoom In Zoom Out Reset image size

The loss function used in our biLSTM model is the weighted binary cross-entropy (WBCE; Goodfellow et al. 2016; Liu et al. 2020). Let N denote the total number of data sequences each having m consecutive data samples in the training set. Let w₀ denote the weight for the positive class (i.e., minority class), and let w₁ denote the weight for the negative class (i.e., majority class). The weights are calculated based on the ratio of majority and minority class sizes with more weight assigned to the minority class. Let y_i denote the observed probability of the ith data sequence; y_i is 1 if the ith data sequence is positive and 0 if the ith data sequence is negative. Let ${\hat{y}}_{i}$ denote the predicted probability of the ith data sequence. The WBCE, calculated as follows, is suitable for imbalanced data sets such as those tackled here where the negative class has more data samples than the positive class; see Table 1.

$$\begin{eqnarray}&&{\rm{WBCE}}=\sum _{i=1}^{N}{w}_{0}{y}_{i}\mathrm{log}({\hat{y}}_{i})+{w}_{1}(1-{y}_{i})\mathrm{log}(1-{\hat{y}}_{i}).\end{eqnarray} \tag{ 5 }$$

We configure the network to use a fraction (1/10) of the training set as the internal validation subset. We employ progressive learning with early stopping and adopt the strategy of saving the highest-performing model during the iterative learning process. The performance of a model is measured by the WBCE on the internal validation subset where the smaller the WBCE is, the better performance the model has. In each iteration, the process checks the performance of the models in the current and previous iterations to decide which model to use for the next iteration. If the model in the current iteration has a better performance, the process copies its weights as starting weights for the next iteration; otherwise, it copies the weights of the model in the previous iteration as starting weights for the next iteration. This progressive process improves the weights of the network's hidden layers, and as a result, the overall performance of the network is also improved. In addition, during the iterations, if the performance of the network degrades, the process stops and selects the highest-performing model it identifies within the iterations.

During testing/prediction, we are given a test data sample x_t , and our biLSTM model will predict the label (color) of x_t , i.e., predict whether x_t is positive or negative. We pack the m − 1 data samples preceding x_t , namely x_t−m+1, x_t−m+2, ..., x_t−1, along with x_t into a test data sequence with m consecutive data samples and feed this test data sequence to our biLSTM model as shown in the input layer in Figure 2. Figure 3(c) illustrates a test data sequence where m is 10. The output layer of our biLSTM model calculates a probability between 0 and 1 for the test data sequence. We compare the probability with a threshold, which is set to 0.5. If the probability is greater than or equal to the threshold, our biLSTM model outputs 1 indicating the test data sequence, more precisely the test data sample x_t , is positive; otherwise our model outputs 0 indicating the test data sequence, more precisely x_t , is negative.

We conducted a series of experiments to evaluate the performance of the proposed method and compare it with related ML methods. For the data sample x_t at time point t, we define the following:

• 1.  
  x_t to be true positive (TP) if our model (network) predicts that x_t is positive and x_t is indeed positive, i.e., an SEP event will be produced with respect to x_t ;
• 2.  
  x_t to be false positive (FP) if our model predicts that x_t is positive while x_t is actually negative, i.e., no SEP event will be produced with respect to x_t ;
• 3.  
  x_t to be true negative (TN) if our model predicts x_t is negative and x_t is indeed negative;
• 4.  
  x_t to be false negative (FN) if our model predicts x_t is negative while x_t is actually positive.

We also use TP (FP, TN, FN, respectively) to denote the total number of TPs (FPs, TNs, FNs, respectively) produced by a method.

The following performance metrics are used in our study:

$$\begin{eqnarray}&&\mathrm{Recall}=\displaystyle \frac{\mathrm{TP}}{\mathrm{TP}\ +\ \mathrm{FN}},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\mathrm{Precision}=\displaystyle \frac{\mathrm{TP}}{\mathrm{TP}\ +\ \mathrm{FP}},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{Balanced}\ \mathrm{Accuracy}\ (\mathrm{BACC})=\displaystyle \frac{1}{2}\\ \quad \quad \times \left(\displaystyle \frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}+\displaystyle \frac{\mathrm{TN}}{\mathrm{TN}+\mathrm{FP}}\right),\end{array}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{Heidke}\ \mathrm{Skill}\ \mathrm{Score}\ (\mathrm{HSS})\\ =\,\displaystyle \frac{2\times (\mathrm{TP}\times \mathrm{TN}\ -\ \mathrm{FP}\times \mathrm{FN})}{(\mathrm{TP}\ +\ \mathrm{FN})\times (\mathrm{FN}+\mathrm{TN})+(\mathrm{TP}\ +\ \mathrm{FP})\times (\mathrm{FP}+\mathrm{TN})},\end{array}\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\mathrm{True}\ \mathrm{Skill}\ \mathrm{Statistics}\ (\mathrm{TSS})=\displaystyle \frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}-\displaystyle \frac{\mathrm{FP}}{\mathrm{FP}+\mathrm{TN}}.\end{eqnarray} \tag{ 10 }$$

BACC (García et al. 2009) is an accuracy measure mainly for imbalanced data sets. HSS (Heidke 1926) and TSS (Bloomfield et al. 2012) are commonly used for flare, CME, and SEP predictions (Bloomfield et al. 2012; Florios et al. 2018; Inceoglu et al. 2018; Liu et al. 2019, 2020). HSS ranges from −∞ to +1. The higher HSS value a method has, the better performance the method achieves. TSS ranges from −1 to +1. Like HSS, the higher TSS value a method has, the better performance the method achieves. In addition, we use the weighted area under the curve (WAUC; Bekkar et al. 2013) in our study. The area under the curve (AUC) in a receiver operating characteristic curve (Marzban 2004) indicates how well a method is capable of distinguishing between two classes in binary prediction with the ideal value of one. When calculating the AUC, we do not distinguish between the accuracy on the minority class (positive class) and the accuracy on the majority class (negative class). In contrast, when calculating the WAUC, which is an extension of the AUC and mainly for imbalanced data sets like those tackled here, the accuracy on the minority class has a larger contribution to the overall performance of a model than the accuracy on the majority class. As a consequence, we assign more weight to the accuracy on the minority class where the weight is defined to be the ratio of the sizes of the minority and majority classes. All the metrics mentioned above are calculated using the confusion matrices obtained from the cross-validation (CV) scheme. With CV, we train a model using a subset of data, called the training set, and test the model using another subset of data, called the test set, where the training set and test set are disjointed. We consider six years, namely 2011–2015 and 2017, as mentioned in Section 2. Data samples from each year in turn are used for testing in a run, and data samples from all the other five years together are used for training in the run. There are six years, and hence there are six runs in total. For each performance metric, the mean and standard deviation over the six runs are calculated and recorded.

We first assessed the importance of the 18 SHARP parameters described in Section 2 to understand which parameters are the most important ones with the greatest predictive power by utilizing a parameter ranking method, called Stability Selection (Meinshausen & Bühlmann 2010). This method is based on the Least Absolute Shrinkage and Selection Operator algorithm (Tibshirani 1996). Table 2 presents the rankings of the parameters with respect to T = 12, 24, 36, 48, 60, and 72 for the FC_S and F_S problems respectively. The parameter ranked first is the most important one while the parameter ranked 18th is the least important one. ABSNJZH is ranked consistently high for the FC_S problem while SAVNCPP and TOTUSJH are ranked high for the F_S problem. AREA_ACR, TOTUSJZ, and USFLUX are ranked consistently low for both of the FC_S and F_S problems.

We then used the recursive parameter elimination algorithm (Butcher & Smith 2020) in combination with our biLSTM model to select a set of parameters that achieves the best performance where the performance is measured by TSS. The parameter elimination algorithm is an interactive procedure. It selects parameters by recursively considering smaller and smaller sets of parameters where the least important parameters are successively pruned from the current set of parameters. Figure 4 presents the parameter selection results for the FC_S and F_S problems respectively. It can be seen from the figure that using the top 15 most important parameters achieves the best performance for both of the FC_S and F_S problems. When using the top k, 1 ≤ k ≤ 14, most important parameters, the less parameters we use, the worse performance our model achieves. Using the top-ranked, most important parameter alone would yield a lower TSS than using all the top 15 most important parameters together. In subsequent experiments, we used the top 15 most important parameters for our biLSTM model. That is, we removed the three least important parameters AREA_ACR, TOTUSJZ, and USFLUX from data samples, and each data sample contained only the top 15 most important SHARP parameters.

Zoom In Zoom Out Reset image size

Next, we compared our biLSTM network with four related ML methods, including MLP, SVMs, RF, and long short-term memory (LSTM; Liu et al. 2019). These four methods are commonly used to predict solar flares, CMEs, and SEPs (Bobra & Ilonidis 2016; Liu et al. 2017; Florios et al. 2018; Inceoglu et al. 2018; Chen et al. 2019; Liu et al. 2019, 2020; Wang et al. 2020; Abduallah et al. 2021).

MLP (Rosenblatt 1958; Arias del Campo et al. 2021) is a feed-forward artificial neural network (Braspenning et al. 1995) that consists of an input layer, an output layer, and one or more hidden layers. The number of hidden layers is set to 3 with 200 neurons in each hidden layer. SVM (Cristianini & Ricci 2008) is trained with the Radial Basis Function kernel, and the cache size is set to 20,000 to speed the training process. RF (Breiman et al. 1984) is an ensemble algorithm that has two hyperparameters for performance tuning: m (the number of SHARP magnetic parameters randomly selected and used to split a node in a tree of the forest) and n (the number of trees to grow). We set m to 2 and n to 500. The implementation of LSTM follows that described in Liu et al. (2020). The hyperparameters not specified here are set to their default values provided by the scikit-learn library in Python (Pedregosa et al. 2011).

As done in biLSTM, we used the recursive parameter elimination algorithm (Butcher & Smith 2020) to identify and select the best parameters for the four related ML methods based on the importance rankings of the 18 SHARP parameters in Table 2 for the FC_S and F_S problems respectively. Our experiments showed that, like biLSTM, using the top 15 most important parameters achieved the best performance for the four related ML methods. Consequently, we used the top 15 most important parameters for the four ML methods in the experimental study.

Figures 5 and 6 present the confusion matrices of the five ML methods (RF, MLP, SVM, LSTM, biLSTM) for the FC_S and F_S problems respectively. For each T, T = 12, 24, 36, 48, 60, 72, and each ML method, the figures show the minimum, average, maximum (displayed from top to bottom) TP, FN, TN, FP, respectively, from the six runs based on our CV scheme. For example, refer to T = 12 and biLSTM in Figure 5. The minimum (maximum, respectively) TP obtained by biLSTM from the six runs is 87 (246, respectively); the average TP over the six runs is 139. It can be seen from Figures 5 and 6 that the average TN values are much larger than the average FP values for both of the FC_S and F_S problems. This happens because there are many negative training data samples in our data sets (see Table 1). As a consequence, the ML methods gain sufficient knowledge about the negative data samples and hence can detect them relatively easily. For the FC_S problem, the average TP values (TN values, respectively) are consistently larger than the average FN values (FP values, respectively), indicating that the ML methods can solve the FC_S problem reasonably well. For the F_S problem, the average TP values are close to, or even smaller than, the average FN values in many cases, suggesting that the ML methods have difficulty in detecting positive data samples. This is understandable given that there are much fewer positive training data samples than negative training data samples for the F_S problem (see Table 1).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Tables 3 and 4 compare the performance of the five ML methods for the FC_S and F_S problems respectively. The tables present the mean performance metric values averaged over the six runs based on our CV scheme with standard deviations enclosed in parentheses. Best average metric values are highlighted in boldface. It can be seen from Tables 3 and 4 that our biLSTM network outperforms the four related ML methods in terms of BACC, HSS, TSS, and WAUC. Furthermore, the five ML methods generally perform better in solving the FC_S problem than in solving the F_S problem. This result indicates that one can predict SEP events more accurately when ARs will produce both flares and associated CMEs. Using flare information alone to predict SEP events is harder and would produce less reliable prediction results.

The five ML methods (RF, MLP, SVM, LSTM, biLSTM) studied here are inherently probabilistic forecasting models in the sense that they calculate a probability between 0 and 1. We compare the probability with a threshold, which is set to 0.5, to determine the output produced by each ML method. The output is either 1 or 0 (see Figure 2), and hence each method is essentially a binary prediction model. In addition to comparing the methods used as binary prediction models, we also compare the methods used as probabilistic forecasting models, where the output produced by each model is interpreted as follows. [FC_S problem] Given a data sample x_t at time point t in an AR where the AR will produce an M- or X-class flare within the next T hours of t and the flare initiates a CME, based on the SHARP parameters in x_t and its preceding m − 1 data samples x_t−m+1, x_t−m+2, ..., x_t−1, we calculate and output a probabilistic estimate of how likely it is that the AR will produce an SEP event associated with the flare and CME. [F_S problem] Given a data sample x_t at time point t in an AR where the AR will produce an M- or X-class flare within the next T hours of t regardless of whether or not the flare initiates a CME, based on the SHARP parameters in x_t and its preceding m − 1 data samples x_t−m+1, x_t−m+2, ..., x_t−1, we calculate and output a probabilistic estimate of how likely it is that the AR will produce an SEP event associated with the flare.

The distribution and behavior of the predicted probabilistic values may not match the expected distribution of observed probabilities in the training data. One can adjust the distribution of the predicted probabilities to better match the expected distribution observed in the training data through calibration. Here, we adopt isotonic regression (Kruskal 1964; Sager & Thisted 1982) to adjust the probabilities. Isotonic regression works by fitting a free-form line to a sequence of data points such that the fitted line is nondecreasing (or nonincreasing) everywhere, and lies as close to the data points as possible. Calibrated models often produce more accurate results. We add a suffix "+C" to each model to denote the calibrated version of the model.

To quantitatively assess the performance of a probabilistic forecasting model, we adopt the Brier Score (BS; Wilks 2010) and Brier Skill Score (BSS; Wilks 2010), defined as follows:

$$\begin{eqnarray}&&\mathrm{BS}=\displaystyle \frac{1}{N}\sum _{i=1}^{N}{({y}_{i}-{\hat{y}}_{i})}^{2},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&\mathrm{BSS}=1-\displaystyle \frac{{\rm{BS}}}{\tfrac{1}{N}\sum _{i=1}^{N}{({y}_{i}-\bar{y})}^{2}}.\end{eqnarray} \tag{ 12 }$$

Here, N is the total number of data sequences each having m consecutive data samples in the test set (see Figure 2 where a test data sequence with m consecutive data samples is fed to our biLSTM model); y_i denotes the observed probability, and ${\hat{y}}_{i}$ denotes the predicted probability of the ith test data sequence respectively; $\bar{y}=\tfrac{1}{N}{\sum }_{i=1}^{N}{y}_{i}$ denotes the mean of all the observed probabilities. The BS values range from 0 to 1, with 0 being a perfect score, whereas the BSS values range from −∞ to 1, with 1 being a perfect score.

Table 5 compares the performance of the five ML methods used as probabilistic forecasting models for the FC_S and F_S problems respectively. The table presents the mean BS and BSS values averaged over the six runs based on our CV scheme with standard deviations enclosed in parentheses. Best BS and BSS values are highlighted in boldface. It can be seen from Table 5 that the probabilistic forecasting models generally perform better in solving the FC_S problem than in solving the F_S problem, suggesting that F_S is a harder problem, and hence the forecasting results for the F_S problem would be less reliable. These findings are consistent with those in Tables 3 and 4 where the ML methods are used as binary prediction models. Furthermore, the calibrated version of a model is better than the model without calibration. Overall, biLSTM+C performs the best among all the models in terms of both BS and BSS.