Predicting Coronal Mass Ejections Using SDO/HMI Vector Magnetic Data Products and Recurrent Neural Networks

Following Bobra & Ilonidis (2016), we intend to use past observations of a flaring AR to predict its future CME productivity. Specifically, we want to solve the following binary classification problem: will an AR that produces an M- or X-class flare within the next T hours also produce a CME associated with the flare? We consider T ranging from 12 to 60 in 12 hr intervals. These prediction times are commonly used by researchers (Ahmed et al. 2013; Bobra & Ilonidis 2016; Inceoglu et al. 2018).

Our data set contains data samples collected within the T hours prior to the peak time of an M- or X-class flare regardless of whether the flare produces a CME. Those data samples collected within the T hours prior to the peak time of an M- or X-class flare that produces a CME belong to the positive class; the other data samples belong to the negative class. Data samples collected in years 2010–2014 are used for training, and those in years 2015–2019 are used for testing. Thus, the training set and test set are disjoint, and hence our proposed RNNs will make predictions on ARs that they have never seen before. Table 1 summarizes the numbers of positive and negative samples for different T values used for training and testing respectively.

If a data sample is missing at some time point or if there are insufficient data samples within the T hours prior to the peak time of an M- or X-class flare, we adopt a zero-padding strategy by adding synthetic data samples with zeros for all feature values to yield a complete, gapless time-series data set. This zero-padding method is used after normalizing the feature values. Therefore, the zero-padding method does not affect the normalization procedure. For a given time point t and an AR that produces an M- or X-class flare within the next T hours of t, the proposed RNNs predict whether or not the flare will initiate a CME.

We employ two kinds of RNNs: a GRU network (Cho et al. 2014; Goodfellow et al. 2016) and an LSTM network (Hochreiter & Schmidhuber 1997; Goodfellow et al. 2016). A GRU contains three interactive parts—a memory content, an update gate, and a reset gate—as illustrated in Figure 1, where boldface is used for matrix notations. The new memory content h_t is updated by the previous memory content h_t−1 and the candidate memory content ${\tilde{{\boldsymbol{h}}}}_{t}$ as follows (Cho et al. 2014; Goodfellow et al. 2016):

$$\begin{eqnarray}&&{{\boldsymbol{h}}}_{t}={{\boldsymbol{z}}}_{t}\odot {{\boldsymbol{h}}}_{t-1}+(1-{{\boldsymbol{z}}}_{t})\odot {\tilde{{\boldsymbol{h}}}}_{t},\end{eqnarray} \tag{ 1 }$$

where the update gate z_t that determines how much of the past information from previous time steps needs to be passed to the future is calculated as follows (Cho et al. 2014; Goodfellow et al. 2016):

$$\begin{eqnarray}&&{{\boldsymbol{z}}}_{t}=\sigma ({{\boldsymbol{W}}}_{z}\cdot [{{\boldsymbol{h}}}_{t-1},{{\boldsymbol{x}}}_{t}]+{{\boldsymbol{B}}}_{z}),\end{eqnarray} \tag{ 2 }$$

and the reset gate r_t that determines how much of the past information to forget is computed as follows (Cho et al. 2014; Goodfellow et al. 2016):

$$\begin{eqnarray}&&{{\boldsymbol{r}}}_{t}=\sigma ({{\boldsymbol{W}}}_{r}\cdot [{{\boldsymbol{h}}}_{t-1},{{\boldsymbol{x}}}_{t}]+{{\boldsymbol{B}}}_{r}).\end{eqnarray} \tag{ 3 }$$

Here x_t represents the input vector at time step t. The candidate memory content ${\tilde{{\boldsymbol{h}}}}_{t}$ is computed as follows (Cho et al. 2014; Goodfellow et al. 2016):

$$\begin{eqnarray}&&{\tilde{{\boldsymbol{h}}}}_{t}=\tanh ({{\boldsymbol{W}}}_{h}\cdot [{{\boldsymbol{r}}}_{t}\odot {{\boldsymbol{h}}}_{t-1},{{\boldsymbol{x}}}_{t}]+{{\boldsymbol{B}}}_{h}).\end{eqnarray} \tag{ 4 }$$

In the above equations, W and B contain weights and biases respectively, which need to be learned during training; [.] denotes the concatenation of two vectors; σ(·) is the sigmoid function, i.e., $\sigma ({\boldsymbol{z}})=\tfrac{1}{1+{e}^{-{\boldsymbol{z}}}};$ tanh(·) is the hyperbolic tangent function, i.e., $\tanh ({\boldsymbol{z}})=\tfrac{{e}^{{\boldsymbol{z}}}-{e}^{-{\boldsymbol{z}}}}{{e}^{{\boldsymbol{z}}}+{e}^{-{\boldsymbol{z}}}};$ ⊙ denotes the Hadamard product (element-wise multiplication).

Zoom In Zoom Out Reset image size

Our GRU network contains a GRU layer with m GRUs (in the study presented here, m is set to 20). We add an attention layer with m neurons above the GRU layer to focus on information in relevant time steps as done in Liu et al. (2019). We then add a fully connected layer with 100 neurons on top of the attention layer. Finally, the output layer with one neuron, which is activated by the sigmoid function, produces predicted values. Our LSTM network is similar to the GRU network except that the GRU layer is replaced by an LSTM layer with m LSTM cells. This is reminiscent of the LSTM network presented in Liu et al. (2019).

Let x_t represent the data sample collected at time point t. During training, for each time point t, we take m consecutive data samples ${{\boldsymbol{x}}}_{t-m+1},{{\boldsymbol{x}}}_{t-m+2},...,{{\boldsymbol{x}}}_{t-1},{{\boldsymbol{x}}}_{t}$ from the training set and use the m consecutive data samples to train the proposed RNNs including the GRU and LSTM networks. The label of these m consecutive data samples is defined to be the label of the last data sample x_t. Thus, if x_t belongs to the positive class, then the input sequence ${{\boldsymbol{x}}}_{t-m+1},{{\boldsymbol{x}}}_{t-m+2},...,{{\boldsymbol{x}}}_{t-1},{{\boldsymbol{x}}}_{t}$ is defined as positive; otherwise the sequence is defined as negative. Because the data samples are collected continuously at a cadence of 12 minutes and missing values are filled up by our zero-padding strategy, the input sequence spans m/5 hr.

Also during training, we use a weighted cross-entropy cost function for optimizing model parameters, where the cost function is computed as follows (Goodfellow et al. 2016):

$$\begin{eqnarray}&&J=\sum _{n=1}^{N}{\omega }_{0}{y}_{n}{\rm{log}}({\hat{y}}_{n})+{\omega }_{1}(1-{y}_{n}){\rm{log}}(1-{\hat{y}}_{n}).\end{eqnarray} \tag{ 5 }$$

Here, N is the total number of sequences each having m consecutive data samples in the training set, ω₀ and ω₁ are the weights of the positive and negative classes respectively, which are derived from the ratio of the sizes of the positive and negative classes with more weight given to the minority class.⁷ We use y_n and ${\hat{y}}_{n}$ to denote the observed probability (which is equal to 1 if the nth sequence belongs to the positive class) and the estimated probability of the nth sequence, respectively.

The proposed RNN methods are implemented in Python, TensorFlow, and Keras. A mini-batch strategy (Goodfellow et al. 2016) is used to achieve faster convergence during backpropagation. The optimizer used is RMSprop, which is a method for gradient descent, where the learning rate is set to 0.001. The batch size is set to 256 and the number of epochs is set to 20. The length of each input sequence, m, is set to 20, meaning that every time 20 consecutive data samples are used as input to our RNNs. The hyperparameter values, the optimizer, and the cost function in Equation (5) are chosen to maximize the TSS scores to be defined in section 4.1.

During testing, to predict whether a given AR that produces an M- or X-class flare within the next T hours of a time point t will also produce a CME associated with the flare, we take x_t and its preceding m − 1 data samples, and then feed the m consecutive test data samples ${{\boldsymbol{x}}}_{t-m+1},{{\boldsymbol{x}}}_{t-m+2},...,{{\boldsymbol{x}}}_{t-1},{{\boldsymbol{x}}}_{t}$ into the trained RNNs. The output of the RNNs, i.e., the predicted result, is a scalar number with a value of 1 or 0, indicating that x_t is positive (i.e., the AR will also produce a CME associated with the flare) or x_t is negative (i.e., the AR will not produce a CME associated with the flare). This value is determined by comparing the probability calculated by the sigmoid function in the output layer of the RNNs with a threshold. If the probability is greater than or equal to the threshold, then x_t is predicted to be positive; otherwise x_t is predicted to be negative. It should be pointed out that the way we use the m consecutive test data samples ${{\boldsymbol{x}}}_{t-m+1},{{\boldsymbol{x}}}_{t-m+2},...,{{\boldsymbol{x}}}_{t-1},{{\boldsymbol{x}}}_{t}$ to predict whether a given AR that produces an M- or X-class flare within the next T hours of a time point t will also produce a CME associated with the flare is different from the previously published machine learning methods (Bobra & Ilonidis 2016), which used only the test data sample x_t to make the prediction.

Given an AR that produces an M- or X-class flare within the next T hours of a time point t and a data sample x_t observed at time point t, we define x_t to be a true positive (TP) if our RNNs predict that x_t is positive, and x_t is indeed positive, i.e., the M- or X-class flare produces, or is associated with, a CME. We define x_t to be a false positive (FP) if our RNNs predict that x_t is positive while x_t is actually negative, i.e., the M- or X-class flare does not produce, or is not associated with, a CME. We say x_t is a true negative (TN) if our RNNs predict x_t to be negative and x_t is indeed negative; x_t is a false negative (FN) if our RNNs predict x_t to be negative while x_t is actually positive. We also use TP (FP, TN, FN, respectively) to represent the total number of true positives (false positives, true negatives, false negatives, respectively) produced by our RNNs.

The performance metrics used in this study include the following:

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

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

$$\begin{eqnarray}&&{\rm{Accuracy}}\,{\rm{(ACC)}}=\displaystyle \frac{\mathrm{TP}\ +\ \mathrm{TN}}{\mathrm{TP}\ +\ \mathrm{FP}\ +\ \mathrm{TN}\ +\ \mathrm{FN}},\end{eqnarray} \tag{ 8 }$$

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

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

The HSS (Heidke 1926) is used to measure the fractional improvement of our prediction over the random prediction (Florios et al. 2018). The TSS score is the recall minus the false alarm rate (Bloomfield et al. 2012). We also use the area under the curve (AUC) in a receiver operating characteristic (ROC) curve (Marzban 2004), which represents the degree of separability, indicating how well a method is capable of distinguishing between two classes with the ideal value of one. These performance metrics are commonly used when dealing with binary classification problems such as the one at hand as defined at the beginning of Section 3.1. In general, the larger HSS, TSS, and AUC score a binary classification method has, the better performance the method achieves.

To gain a better understanding of the behavior of the proposed RNNs, we adopt the following cross-validation methodology. We partition the training (test, respectively) set into 10 equal-sized folds. For every two training (test, respectively) folds i and j, $i\ne j$, fold i and fold j are disjoint; furthermore, fold i and fold j contain approximately the same number of positive training (test, respectively) data samples and approximately the same number of negative training (test, respectively) data samples. In the ith run, $1\leqslant i\leqslant 10$, all training data samples except those in training fold i are used to train a model, and the trained model is used to make predictions on all test data samples except those in test fold i. We calculate the performance metric values based on the predictions made in the ith run. There are 10 runs. The means and standard deviations over the 10 runs are calculated and recorded.

We conduct a series of experiments to analyze the importance of each of the 18 features studied here using the cross-validation methodology described above and the feature assessment method introduced in Section 4.3 of Liu et al. (2019). Each time only one feature is used to make predictions. The probability threshold used by our RNNs is set to maximize the TSS score in each run. There are 10 runs and the corresponding mean TSS score is calculated and recorded. There are 18 features, so 18 mean individual TSS scores are recorded. These 18 mean individual TSS scores are sorted in descending order, and the 18 corresponding features are ranked from the most important (with the highest mean individual TSS score) to the least important (with the lowest) accordingly. Table 2 presents the 18 features ranked by our GRU and LSTM networks respectively for different T values where T ranges from 12 to 60 in 12 hr intervals. It can be seen from the table that MEANPOT, SHRGT45, ABSNJZH, and SAVNCPP are consistently ranked in the top 10 list by both LSTM and GRU networks. In particular, MEANPOTS plays the most important role in CME prediction, and is ranked as the top one in all cases. Other features such as TOTPOT, USFLUX, MEANJZH, MEANALP, and MEANSHR also show strong predictive power and are ranked in the top 10 list in most cases. Compared to the top 10 lists of Bobra & Ilonidis (2016), who used a different feature ranking method for T = 24, 48 (see Table 1 of Bobra & Ilonidis 2016), these lists overlap to some extent (seven features simultaneously occur in the top 10 lists given by Bobra & Ilonidis 2016 and our methods).

Next, mean cumulative TSS scores are calculated according to the ranked features. Specially, the mean cumulative TSS score of the top k, 1 ≤ k ≤ 18, most important features is equal to the mean TSS score of using the top k most important features altogether for CME prediction. We calculate the mean cumulative TSS scores for T = 12, 24, 36, 48, and 60 hr. It is found that using all the 18 features together does not yield the highest mean cumulative TSS scores. In fact, using the top 16, 12, 9, 14, and 5 features for T =  12, 24, 36, 48, and 60 hr respectively yields the highest mean cumulative TSS scores, achieving the best performance for our GRU network. Using the top 15, 12, 8, 15, and 6 features for T = 12, 24, 36, 48, and 60 hr respectively yields the highest mean cumulative TSS scores for our LSTM network. This happens probably because low ranked features are noisy features, and using them may deteriorate the performance of the methods. In subsequent experiments, we use the best features for each method.

We compare our proposed RNNs with three closely related machine learning methods including an MLP (Inceoglu et al. 2018), an SVM (Bobra & Ilonidis 2016), and the random forest algorithm (RF) (Liu et al. 2017). MLP and SVM have been previously used for CME prediction (Bobra & Ilonidis 2016; Inceoglu et al. 2018). RF has been used in flare prediction with good performance (Liu et al. 2017, 2019; Florios et al. 2018). These three machine learning methods are inherently probabilistic forecasting models (Florios et al. 2018) in the sense that each of them predicts a probability. Since we are dealing with a binary classification problem, as defined at the beginning of Section 3.1, we convert each probabilistic forecasting model into a binary classification model (Jonas et al. 2018; Nishizuka et al. 2018) by comparing the predicted probability with a threshold. If the predicted probability is greater than or equal to the threshold, then the model predicts that a flare will produce, or is associated with, a CME; otherwise, the model predicts that the flare will not produce a CME.

MLP consists of an input layer, an output layer, and two hidden layers both with 200 neurons. SVM uses the radial basis function (RBF) kernel. RF has two parameters: mtry (number of features randomly selected to split a node) and ntree (number of trees to grow in the forest). We vary the values of ntree ∈ {300, 500, 1000} and mtry ∈ [2, 8], and set ntree to 500 and mtry to 3. The hyperparameter and parameter values used by these three related machine learning methods are chosen to maximize their TSS scores. As in the proposed RNNs, we use data samples collected in years 2010–2014 for training and data samples in years 2015–2019 for testing. To deal with the imbalanced data sets described in Table 1, we give more weight to the minority class during training, as done for the RNNs. The same cross-validation methodology as described in Section 4.1 is adopted to evaluate the performance of the three related machine learning methods.

Table 3 presents the confusion matrix, which lists mean TP, FP, FN, TN (with standard deviations enclosed in parentheses) of the five machine learning methods for different T values, where T ranges from 12 to 60 hr in 12 hr intervals. The probability thresholds used by the machine learning methods are set to maximize their TSS scores. Table 4 presents the mean performance metric values (with standard deviations enclosed in parentheses) of the five machine learning methods. The best performance metric values are highlighted in boldface. Figure 2 shows the ROC curves for the five machine learning methods. It can be seen from Table 4 and Figure 2 that our GRU and LSTM networks perform better than the three related machine learning methods in terms of HSS, TSS, AUC, and ROC curves. There is no clear distinction between the GRU and LSTM networks. These results indicate that both of the proposed RNNs are suitable for solving the binary classification problem at hand.

Zoom In Zoom Out Reset image size

The three related machine learning methods are inherently probabilistic forecasting models. Our proposed RNNs can be easily converted from a binary classification model to a probabilistic forecasting model as follows. Instead of comparing the probability, calculated by the sigmoid function in the output layer of the RNNs, with a threshold, the RNNs simply output the probability. For a given time point t and an AR that will produce an M- or X-class flare within the next T hours of t, this output now represents a probabilistic estimate of how likely it is that the flare will initiate a CME.

We use the Brier Score (BS) (Brier 1950) and the Brier Skill Score (BSS) (Wilks 2010) to quantitatively assess the performance of a probabilistic forecasting model, where

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

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

Here N is the total number of sequences each having m consecutive data samples in the test set, y_n and ${\hat{y}}_{n}$ denote the observed probability and the estimated probability of the nth sequence respectively as defined in Equation (5), and $\overline{y}=\tfrac{1}{N}{\sum }_{n=1}^{N}{y}_{n}$ is the average value of all the observed probabilities. The values of BS range from 0 to 1, with the perfect score being 0. The values of BSS range from minus infinity to 1, with the perfect score being 1.

Table 5 presents the mean BS and BSS scores and standard deviations of the five machine learning methods for different T values where T ranges from 12 to 60 hr in 12 hr intervals. It can be seen from the table that the proposed GRU and LSTM networks are comparable, and again outperform the three related machine learning methods in terms of BS and BSS. These results in Tables 4 and 5 suggest that our RNNs work well when used as either binary classification models or probabilistic forecasting models.