RELIABLE PROBABILITY FORECAST OF SOLAR FLARES: DEEP FLARE NET-RELIABLE (DeFN-R)

We used a solar flare prediction model employing DNNs named DeFN to investigate how to realize reliable probabilistic forecasts (Nishizuka et al., 2018). We performed probabilistic forecasts for two categories: (i) $\geq$M-class and $<$M-class/non-flare events and (ii) $\geq$C-class and $<$C-class/non-flare events. We used the database of observation data taken by Geostationary Operational Environmental Satellite (GOES) and Solar Dynamic Observatory (SDO; Pesnell et al., 2012) from June 2010 to December 2015. During this period, 90% of the observed flares, such as 26 X-class, 383 M-class, and 4054 C-class flares, occurred on the solar disk. Since the number of X-class flares was too small to evaluate the prediction results, X-class flares were not predicted in this paper.

Occurrence probabilities of flares are predicted by the following steps. (i) First, we download observation data from the data archives of SDO and GOES, such as the light curves of the soft X-ray intensity, line-of-sight magnetograms, vector magnetograms, and images obtained by the 1600  Å and 131  Å filters. (ii) Second, active regions (ARs) are detected from the full-disk line-of-sight magnetograms and numbered by tracking their time evolution. (iii) For each AR, physics-based features are extracted from multiwavelength observation images, and the feature database is annotated with flare labels if an X-, M-, or C-class flare occurs within 24 h after obtaining an image. (iv) The flare occurrence probabilities in the following 24 h is predicted by supervised machine learning using a DNN for each region with a 1 h cadence. Because we investigated reliable outputs by tuning DeFN, we rename it DeFN-Reliable (DeFN-R) in this paper.

First, we made the database from 3$\times$10⁵ solar images of the full-disk line-of-sight magnetograms obtained by the Helioseismic and Magnetic Imager (HMI; Scherrer et al., 2012; Schou et al., 2012; Hoeksema et al., 2014) on board SDO, from which we automatically detected ARs with a reduced cadence of 1 h (see details of the detection rules described in Nishizuka et al., 2017). We neglected ARs on the limb and numbered the same ARs by tracking their time evolution. To extract various kinds of features, the frame coordinates of the detected ARs were applied to other observation images with different wavelengths, such as vector magnetograms taken by HMI/SDO and 1600  Å and 131  Å filter images by the Atmospheric Imaging Assembly (AIA; Lemen et al., 2012) on board SDO.

Next we calculated 79 features from each AR observed by multiwavelength emissions. The features are parameters based on physics that were considered to be effective for solar flare prediction in previous papers and in daily forecasting operations. The extracted features are the same as used by Nishizuka et al. (2018); for example, line-of-sight/vector magnetograms in the photosphere, brightening in the corona at 131  Å (T$\geq$10⁷ K), and the 131  Å and X-ray emissivity 1 and 2 h before an image. The sunspot area, the magnetic flux, the number of magnetic neutral lines, and the Lorentz force were calculated from the line-of-sight/vector magnetograms. Preflare brightening in the bottom chromosphere was detected in the UV–continuum taken by the 1600  Å filter, and the coronal heating in the flaring region over 10⁷ K was taken by the 131  Å filter, as well as the integrated X-ray emission in the range of 1-8  Å observed by GOES.

Furthermore, we standardized the feature database before the input into the DNN (e.g., Bishop, 2006). The database of 2010–2015 was chronologically split into two: the dataset in 2010–2014 for training and the one in 2015 for validation and testing. These chronological datasets for training and testing are more challenging for predicting flares than the randomly shuffled and divided datasets (e.g., Bobra & Couvidat, 2015; Muranushi et al., 2015; Nishizuka et al., 2017, 2018). The optimization and evaluations were repeated several times, using subsets of the training and testing datasets.

DeFN-R is based on algorithms of deep neural networks and composed of multilayer perceptrons. The input is the database of 79 features standardized in the preprocess, i.e., 79-dimensional vectors of features. The output is the forecast probability of each class of flares, $p(y)$. Here, $y$ = (0, 1) is the class of $\geq$M-class flare events and $y$ = (1, 0) is the one of $<$M-class or non-flare events. We calculate the two probabilities of $\geq$M-class flare events and $<$M-class or non-flare events.

The architecture of DeFN-R is fundamentally the same as that of DeFN. DeFN-R consists of several layers. At each layer of the neural networks, the input is converted to the output with a linear combination and an activation function. As the activation function, we used the rectified linear units (ReLU; Nair & Hinton, 2010) for the first layer to the penultimate layer and a softmax function for the last layer,

which gives the outputs of probabilistic forecasts. We determine the parameters in Table 1.

Figure 1(a) shows the architecture of DeFN (see Nishizuka et al., 2018, for simplified representations), while Figures 1(b) and 1(c) show the architectures of DeFN-R, which were retuned to increase BSS for $\geq$M-class and $\geq$C-class flare predictions (see Appendix A.1). DeFN-R includes simple skip connections (see also the residual network in He et al., 2015) and batch normalization represented by the notation BN (Ioffe & Szegedy, 2015) to stabilize the training and improve the precision of the model. The number of nodes was investigated in the range of 79–200, and the batch size was in the range of 50–200. The architecture of DeFN-R with 5–12 layers was surveyed by attaching or detaching skip connections and dropouts. Finally, we selected the network architecture that provides representative results in the pilot experiment, as shown in Table 2.

For classification problems, models are usually trained by optimizing parameters to minimize the cross entropy. The cross entropy between $p(y_{k}^{*})$ and $p(y_{k})$ is determined by the following equation:

where we omitted $n$ representing the $n$th sample for simplicity. $p(y_{k}^{*})$ is the initial probability of correct labels $y_{k}^{*}$, i.e., 1 or 0, while $p(y_{k})$ is forecast probability. The components of $y_{k}^{*}$ are 1 or 0, and thus, $p(y_{k}^{*})$ = $y_{k}^{*}$.

In DeFN, since the flare occurrence rate is imbalanced, the following summation of the weighted cross entropy was used as the loss function instead:

Here, $w_{k}$ is the weight of each class, which was set to the inverse ratio of the class occurrence in DeFN, e.g., [1, 50] for $\geq$M-class flare events and [1, 12] for $\geq$C-class flare events. However, this results in reduced reliability. Thus, in this paper, we set the weight to be constant, i.e., $w_{k}$ = [1, 1].

The parameters used here are summarized in Tables 1 and 2. As the method for stochastic optimization, we used adaptive moment estimation (Adam; Kingma & Jimmy, 2014), which is an extension of AdaGrad, RMSprop, and AdaDelta. We used the recommended values for Adam’s hyperparameters, which are given in Table 2, because overfit hyperparameters only work with certain architectures.

Figures 2(a) and 3(a) respectively show reliability diagrams of probabilistic forecast results for $\geq$M-class flare events obtained by DeFN and DeFN-R. The blue line graph depicts the conditional expectation values of the outcome. The diagonal dotted line indicates perfect reliability, on which 95% consistency bars are attached (Brocker & Smith, 2007; Jolliffe & Stephenson, 2012). The 95% consistency bars show the ranges within which 95% of the conditional expectation values of the outcome given the probability would fall if it were assumed that the original data were sampled from a perfectly reliable probabilistic forecast system, namely, the null hypothesis that the probabilistic forecast is perfectly reliable is rejected with the 95% significant level if the sampled data located outside the consistency bars. The horizontal dotted line is the climatological event rate, and the line between the perfect reliability and the climatological event rate corresponds to BSS = 0. Red histograms show the number of probabilistic forecasts within bins.

Plots located above the diagonal dotted line indicate under-forecasts, while plots located below the line indicate over-forecasts. In Figure 2(a), plots are located below BSS=0 line. This means that DeFN lacks reliability. On the other hand, in Figure 3(a), the reliability diagram is markedly improved, achieving BSS = 0.298. Reliability diagrams for $\geq$C-class flare events are shown in Figures 4(a) and 5(a), and a similar improvement is shown, although the reliability for $\geq$C-class flare prediction is higher (BSS = 0.412). Larger BSS does not always correspond to a better reliability diagram. When BSS is improved, a reliability diagram sometimes gets worse. Since the number of forecasts is larger for a smaller forecast probability, it is more efficient to improve a small forecast probability; thus, the model learns to fit a lower forecast probability preferentially, and large BSS does not always coincide with a better reliability plot.

The output of the probabilistic forecast for a dataset is denoted as $p(y)$. The skill scores for deterministic forecasting require a probability threshold ($P_{th}$), above which an event is predicted to occur. When we set the probability threshold, the output is determined depending on whether $p(y)$ is greater than $P_{th}$ or not. In DeFN or machine learning algorithms trained for deterministic forecasting, the threshold is set to 50% by default, which effectively optimizes TSS, because the errors of forecast probabilities in both categories are treated equally. In probabilistic forecast models such as DeFN-R, the threshold can be freely selected. In these models, categorical skill scores are calculated from contingency tables.

Figures 2(b) and 3(b) show the ROC curve, which tracks the performance of TSS as a function of $P_{th}$, because the two axes, probability of detection (POD) and probability of false detection (POFD), are the first and second terms of TSS, respectively (= POD – POFD) (see Appendix A.2). The ROC curve shows the discrimination performance of correctly recognizing the flare occurrence. Note that the data with the maximum TSS are mostly away from the diagnostic line, but not the data closest to the point (0, 1).

Figures 2(b) and 3(b) are almost the same. For these figures, the ROC score, defined as the area under the ROC curve (AUC) (Mason & Graham, 1999), is 0.96 and 0.93, and the ROC skill score (ROCSS), defined as ROCSS = 2AUC-1, is 0.91 and 0.86, respectively. Therefore, DeFN-R improved the reliability diagram while maintaining the ROC curve. Namely, DeFN-R achieved both high discrimination performance and high reliability. Similar results are also shown in Figures 4(b) and 5(b) for $\geq$C-class flare probabilistic forecasts, where AUC is 0.89 and 0.89, and ROCSS is 0.77 and 0.78, respectively. In other words, the ROC curve and the variability of TSS in section 4.3 are not significantly affected by small changes of a loss function and an architecture, while a reliability diagram and BSS considerably change.

TSS is a function of threshold probability ($P_{th}$), which is shown in Figures 2(c) and 3(c) for $\geq$M-class flare predictions. We took a threshold step of 0.02, resulting in 49 data, and a grid covering the region [0, 1] excluding the edges. In Figure 2(c), TSS continuously changes with the threshold probability and peaks at around 0.50 (50%) for DeFN. On the other hand, in Figure 3(c), TSS peaks at around 0.032 (3.2%), which is the climatological event rate of $\geq$M-class flare events in the datasets. In more detail, the maximum values of TSS are 0.80 and 0.77 at $P_{th}$ = 0.52 and 0.04 in Figures 2(c) and 3(c), respectively. This is confirmed by similar ROC curves of DeFN and DeFN-R. Similar findings are obtained for $\geq$C-class flare predictions, as shown in Figures 4(c) and 5(c). In Figure 4(c), TSS peaks at around 0.50 (50%), while it peaks near the climatological event rate of 0.196 (19.6%) in Figure 5(c). The maximum values of TSS are 0.63 and 0.64 at $P_{th}$ = 0.50 and 0.22 in Figures 4(c) and 5(c), respectively.

As a result, DeFN-R for $\geq$M-class and $\geq$C-class flare predictions shows that the maximum TSS occurs when $P_{th}$ is approximately the climatological event rate. This is consistent with previous works by Bloomfield et al. (2012), Barnes et al. (2016), and Leka et al. (2018). It has been mathematically verified that a reliable probabilistic forecasting model has the peak of TSS at the probability threshold of the climatological event rate (Kubo, 2019). In our case, DeFN-R is reliable, so the maximum TSS occurred at the climatological event rate. However, since the maximum TSS did not occur at this rate, DeFN lacks reliability. This is related to the fact that, in DeFN, the threshold is set to 50% to effectively optimize TSS and that the prior distribution is normalized by the weight of the cross entropy, $w_{k}$.

We also show the variability of the Heidke skill score (HSS) as a function of threshold probability in Figures 2(c), 3(c), 4(c), and 5(c) for comparison. The maximum values of HSS are 0.51 and 0.48 at $P_{th}$ = 0.90 and 0.24 for $\geq$M-class flare prediction in Figures 2(c) and 3(c), respectively. On the other hand, it is 0.57 in both Figures 4(c) and 5(c) at $P_{th}$ = 0.64 and 0.42, respectively, for $\geq$C-class flare prediction. The variability of HSS is also discussed in Bobra & Couvidat (2015) and Florios et al. (2018).