Deep-learning Reconstruction of Sunspot Vector Magnetic Fields for Forecasting Solar Storms

The SHARP features considered (listed in Figure 2 and Table 2) are correlated among each other and are divided into four groups based on mutual Pearson correlations (Dhuri et al. 2019): (i) features that depend on the area of ARs, i.e., extensive features that include AR area, total unsigned flux, total unsigned vertical current, total unsigned current helicity, total free energy density, and total Lorentz force; (ii) features that depend on the electric current in ARs, i.e., absolute net current helicity and sum of net current per polarity; (iii) features that depend only on the nonpotential energy in ARs, i.e., mean free energy density and area with shear >45°; and finally (iv) Schrijver R_value (Schrijver 2007), viz., the sum of flux on the polarity inversion line. Overall, we develop four different CNNs (Figure 2) to estimate SHARP features from these four respective groups. For each CNN, the output layer comprises K neurons to estimate K SHARP features corresponding to each of the four groups.

Extensive features are strongly correlated with AR total unsigned flux, which depends only on the radial component of the magnetic field. The radial component is traditionally estimated from AR LOS magnetic field using a potential field approximation (Leka et al. 2017). Using a CNN, we directly estimate these extensive features without first requiring to estimate the radial magnetic field. The R_value depends only on LOS magnetic field and can be directly calculated using GONG LOS magnetograms. However, to match HMI SHARP R_value, GONG LOS magnetic fields require a cross-calibration. CNN models are expected to implicitly learn the cross-instrument calibration during training (Munoz-Jaramillo et al. 2022), and estimated SHARP values are also expected to be automatically cross-calibrated.

Unlike the extensive features, an accurate estimation of SHARP depending on electric current and mean free energy requires explicit knowledge of the full vector magnetic fields. Such features are important for understanding triggers of solar storms and are typically estimated assuming magnetic field models, e.g., linear and nonlinear force-free models (Régnier & Priest 2007). Here we provide a purely data-driven estimation of these features using a CNN. In order to assess the performance of the CNN, we use linear regression (LR) models as a baseline. We develop two separate LR models, one each for features that depend on electric current and free energy, respectively. As input, the LR models have extensive features and R_value. The first LR model produces absolute net current helicity and the sum of net current per polarity as the output, while the second one produces mean free energy density and area with shear >45° as the output. Figure 2 shows a schematic of the CNN models and the baseline models.

We use Pearson and Spearman correlations for measuring the performance of the CNN and baseline models. Pearson correlation measures a linear correlation between the true and estimated values of the vector magnetic field features. Spearman correlation is a rank correlation that captures the monotonic relationship between the true and estimated values in addition to the linear relationship measured by the Pearson correlation. Pearson and Spearman correlations for the CNN-estimated vector magnetic field features are listed in Tables 2 and 3, respectively. For the baseline models, these correlations are listed in Table 4.

From Table 2, the Pearson correlations of CNN-estimated vector field features are higher for HMI than GONG and thus appear to be dependent on the spatial resolution of LOS magnetograms. For HMI, the CNN-estimated extensive features yield a Pearson correlation of ∼95% for the validation data and ∼90% for the test data. For GONG data, the corresponding correlation is ∼90%. The Pearson correlations of CNN-estimated values of extensive features are not a perfect ∼100%, since the SHARP calculation does not consider all pixels, but rather only taking into account those for which disambiguation of the azimuthal component of the magnetic field is reliable (Bobra et al. 2014). From Table 3, the Spearman correlations for the extensive features are only slightly lower than the corresponding Pearson correlations, implying that the ranking of the estimated features is generally consistent with the true ranking.

The Pearson and Spearman correlation values for features that depend on the nonpotential energy are significantly high (>90%) across the validation and test data sets. These correlation values are also ∼10%–20% higher compared to the LR baseline model (Table 4). These features, namely, mean free energy density and area with shear >45°, explicitly depend on the full vector magnetic field.

For features that depend on electric current—i.e., absolute net current helicity and sum of net current per polarity—the CNN does not perform better than the baseline model. While the Pearson correlations for the validation are 20% higher compared to the baseline of 70%, Spearman correlations are approximately equal (up to the error bars) at 60%. In addition, the CNN fails to generalize to the test data, with low Pearson and Spearman correlation scores of 60% each.

Figure 3 shows scatter plot visualizations of the correlation between true and CNN-estimated SHARP features for HMI and GONG from the 10 validation sets. For HMI, the true and CNN-derived values mostly match relatively closely, except at only very small values (<200 G² − m⁻¹) of absolute net current helicity, where the CNN estimates are significantly larger. For the HMI test data and GONG data, the CNN estimation of absolute net current helicity for large values (>1000 G² − m⁻¹) is consistently on the lower side (∼500 G² − m⁻¹). Figure 4 explicitly shows mean absolute errors in the CNN estimation as a function of the true values for HMI and GONG. Mean absolute errors in CNN-estimated values from GONG magnetograms show higher dependence on true values compared to HMI and increase significantly with increasing true values of the respective features, particularly for the validation data. For total unsigned flux, mean absolute errors of CNN-estimated features of both HMI and GONG are significantly higher for the extreme values (15–20) × 10²² Mx. For the HMI test data and GONG data, mean absolute errors are more than 12 times higher at large magnitudes (>1000 G² − m⁻¹) of absolute net current helicity compared to the HMI validation data. The average relative errors for GONG and HMI are comparable at ≈80% ± 10%, 900% ± 100%, and 25% ± 2% for total unsigned flux, absolute net current helicity, and mean free energy density, respectively. The high average relative errors imply that the CNN estimates are far off from true values, particularly for SHARP features with low true values. SHARP features from ARs that produce at least one major flare (M5 or greater) show a significant drop in average relative errors, at approximately 30% ± 15%, 300% ± 80%, and 16% ± 2%, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For understanding AR magnetic field dynamics and improving forecasting of solar storms, it is important that temporal variations of the CNN-estimated SHARP is faithful to the true SHARP. We measure trends in the time evolution of SHARP features of an AR by fitting the observed and the CNN-estimated values with smooth spline curves and calculate numerical time derivatives. Table 5 lists Pearson correlations between time derivatives of splines, fitted to the true and the CNN-estimated values of total unsigned flux, absolute net current helicity, and mean free energy density. We find that the Pearson correlations are high, ∼80%, for HMI, with the exception of absolute net current helicity values from the test data. For GONG, only the Pearson correlations for mean free energy density are high enough (∼70%) to suggest that the corresponding trends are captured reasonably accurately in the CNN-estimated features. These discrepancies between trends of the CNN-estimated GONG and true values appear to be a consequence of the lower resolution of GONG magnetograms.

A comparison of the time evolution of true and CNN-estimated features obtained from HMI and GONG for individual ARs that produce at least one major flare (M5 or greater) is shown in Figure 5. The true and CNN-estimated values of total unsigned flux, absolute net current helicity, and mean free energy density are in agreement, particularly for HMI, capturing evolution of these features before and after flares. Disagreements between the true and CNN-estimated features occur only at the extreme values of these features. For example, for the X9.3 flare in NOAA Active Region 12673 in 2017 September, which was the largest flare in cycle 24, the CNN accurately estimates the rise of the total unsigned flux and also mean free energy density prior to the flare. The absolute net current helicity rises to unusually high values prior to the X9.3 flare (30% more than the maximum values encountered in the training data), and therefore the corresponding CNN estimates are inaccurate. More examples of comparisons of the time evolution between true and CNN-estimated features from ARs that produce at least one M5 or greater flare are included in Appendix A.

Zoom In Zoom Out Reset image size

The CNN estimation of SHARP features on flaring ARs is thus useful for understanding AR magnetic field evolution leading to particularly violent solar storms in the past. The Halloween storms of 2003 October produced extreme flares from NOAA AR 10486 of magnitudes X17.0, X10.0, and the largest recorded flare X28.0 (Pulkkinen et al. 2005). The magnetic field evolution leading to these extreme flares was characterized by rotation of a major positive polarity of the delta sunspot as shown in the top panel of Figure 6 (Zhang et al. 2008). Without the knowledge of vector magnetic fields, free energy and current helicity during these storms are previously modeled based on the magnetic virial theorem (Metcalf et al. 2005; Régnier & Priest 2007), linear/nonlinear force-free field extrapolation (Régnier & Priest 2007), and a Minimum Current Corona model (Kazachenko et al. 2010). We obtain model-free and purely data-driven CNN estimates of total unsigned flux, absolute net current helicity, and mean free energy density during these storms using LOS magnetograms. However, the HMI observations are not available for this period. We therefore use the CNN trained with HMI magnetograms to process LOS observations from MDI during the Halloween storms to estimate time evolution of total unsigned flux, absolute net current helicity, and mean free energy density. Flare X28.0 is excluded because it occurred outside 45° of the central meridian. In particular, the CNN-estimated absolute net current helicity of NOAA AR 10486 rises continuously by 25% between the X1.2 flare and the X17.0 flare corresponding to the observed sunspot rotation. A similar gradual rise of a modeled helicity flux by 50% between the X1.2 and X17.0 flares has been reported (Kazachenko et al. 2010), caused primarily by helicity injection from the rotation of the sunspot. The CNN estimates show that the absolute net current helicity stays high leading to the X10.0 flare and falls thereafter. The CNN-estimated mean free energy density also rises leading to the X17.0 and X10.0 flares. Note that these CNN-estimated values from MDI magnetograms are not expected to be corrected for the instrument cross-calibration between the MDI and HMI since the CNN is trained with only HMI magnetograms. Table 8 in Appendix B lists the Pearson and Spearman correlations between the true values and the CNN-estimated values using MDI LOS magnetograms, during the overlap period of MDI and HMI. These correlation values are significantly lower compared to those estimated from HMI magnetograms (Tables 2 and 3). Therefore, a rigorous estimation first requires standardization of MDI and HMI magnetograms (e.g., with other approaches such as super-resolution; Munoz-Jaramillo et al. 2022). We also used the GONG magnetograms to estimate the vector field features during the storms using the CNN trained with GONG (see Appendix C). The values of the vector field features estimated using GONG magnetograms are in the extreme range, as expected during the storms. However, the sensitivity of these estimated values to pre- and post-flare magnetic field variations is lower compared to the features estimated from MDI.

Zoom In Zoom Out Reset image size

The SHARP features have been extensively used for building flare forecasting models using ML (Bobra & Couvidat 2015; Bobra & Ilonidis 2016; Nishizuka et al. 2017; Chen et al. 2019; Dhuri et al. 2019; Ahmadzadeh et al. 2021). In order to assess the utility of CNN-estimated SHARP for flare forecasting tasks, we compare their flare forecasting performance to the true SHARP. We set up the problem of forecasting M-/X-class flares with 24 hr warning, similarly to Bobra & Couvidat (2015). We use two approaches for the comparison. First, we build Linear Discriminant Analysis (LDA) classification models using one SHARP feature at a time. This allows for direct comparison of the true and the CNN-estimated values of each SHARP feature for flare forecasting. Second, we use all SHARP features together to train a support vector machine (SVM) for flare forecasting. We measure the flare forecasting performance using accuracy, recall, and the True Skill Statistics (TSS) score (Peirce 1884). However, only the latter two are robust to the class imbalance prevalent in the flare forecasting problem (Bobra & Couvidat 2015; Ahmadzadeh et al. 2021) and therefore reliable for comparison. Our definitions of positive and negative classes are identical to the operational approach described in Bobra & Couvidat (2015). In addition, we use the 10 times repeated-holdout validation described in Section 3. Unlike Bobra & Couvidat (2015), we explicitly ensure that the samples from a given AR are not mixed in training and validation sets (Ahmadzadeh et al. 2021). In addition, as mentioned in Section 2, we only consider ARs with maximum area >25 Mm². Both the LDA and SVM are implemented using the scikit-learn library in Python.

Table 6 lists performance metrics for the classification of M-/X-class flares using the LDA of one SHARP feature at a time. The accuracy, recall, and TSS values obtained using each of the CNN-estimated features from HMI and GONG magnetograms are consistent with those of the true SHARP features up to the validation error bars. We note that Schrijver's R_value (Schrijver 2007) gives the highest TSS values for flare forecasting using individual features.

Table 7 lists the performance metrics for the SVM classification of M-/X-class flares using all SHARP features together. TSS (∼68%) and recall (∼86%) values obtained using an SVM trained with the CNN-estimated features from HMI are consistent with those obtained using the true SHARP. TSS (∼62%) and recall (∼80%) values from an SVM trained with the CNN-estimated features from GONG are slightly lower. For a comparison, we list TSS (∼76%) and recall (∼83%) from Bobra & Couvidat (2015) that are higher. The systematically lower TSS of the SVM in forecasting flares when using true SHARP values here as compared with Bobra & Couvidat (2015) is due to exclusion of observations from ARs with maximum area <25 Mm² (all nonflaring) and the explicit restriction that samples from an AR are part of either training or validation sets. Largely consistent performance metrics for flare forecasting with the CNN-estimated SHARP imply that, high relative errors notwithstanding, the CNN-estimated features can be useful for building space weather forecasting tools. This is a consequence of (true) SHARP feature values varying over several orders of magnitudes and thus being significantly different for flaring and nonflaring ARs for forecasting of flares (Dhuri et al. 2019). Accuracy of the CNN-estimated SHARP features may be improved by significantly increasing the resolution of LOS magnetograms from, e.g., GONG, using techniques such as super-resolution (Munoz-Jaramillo et al. 2022). Our method is thus suitable for reconstructing vector field features from historical LOS magnetograms, ultimately useful for reliable space weather forecasting.

CNNs and, in general, deep learning are extremely efficient at identifying correlations in the data. In this case, the CNN builds a useful model of AR vector magnetic fields from the observed LOS magnetograms. In particular, the CNN-estimated SHARP features may be reliably used to study energy buildup and time evolution of magnetic fields in flaring ARs. Yet it is very challenging to open up the trained network and understand the CNN to uncover the information absorbed. Nevertheless, weights learned by the CNN can shed some light on its working. There are also attribution methods to quantify the contribution of different parts of the input image to the CNN's output. Here we analyze the weights of the CNN and obtain attribution maps for input magnetograms to interpret the trained CNN.

The CNN architecture (Figure 1) comprises a fully convolutional network for processing LOS magnetograms and a fully connected layer of neurons for processing information about the location of ARs on the solar disk. The penultimate concatenation layer comprises a global-average-pooling layer, a global-max-pooling layer that processes the LOS magnetograms, and a fully connected layer that processes the location of the ARs. The global-average-pooling neurons are sensitive to the entire spatial extent of LOS magnetograms, while the global-max-pooling neurons are sensitive to spatially local patterns. The fully connected neurons are sensitive to AR coordinates on the disk. Figure 7 illustrates the distribution of the top weights of each of the three components in the penultimate layers as their contribution to the output of the CNN that estimates total unsigned flux, absolute net current helicity, and mean free energy density. Neurons associated with global average pooling contribute dominantly to the total unsigned flux and mean free energy, implying that their estimation depends on the consideration of the entire LOS magnetograms. For absolute net current helicity, key contributors are neurons from the global-max-pooling layer, and its estimation is sensitive to spatially local patterns from LOS magnetograms. Without the global-max-pooling layer, absolute net current helicity and related CNN-estimated SHARP features show ∼30% less Pearson correlation with the true values. Weights from neurons related to AR location on the solar disk are ∼0, and thus the CNN estimation does not strongly depend on the AR location. Indeed, the CNN may be trained equally well without the additional input of the AR location. This may be a consequence of considering AR patches only within ±45°, where the projection effects are not significant.

Zoom In Zoom Out Reset image size

While there are many attribution methods, gradient-based methods such as saliency maps (Simonyan et al. 2013), grad-CAMs (Selvaraju et al. 2017), integrated gradients (IG; Sundararajan et al. 2017), etc., are favored over perturbation-based methods such as occlusion masks (Zeiler & Fergus 2014) because of computational efficiency and higher-resolution attribution maps. IG attribution maps are of the same resolution as the input magnetograms and are thus superior to grad-CAMs obtained from the CNN feature maps. In addition, unlike saliency maps, IG attribution maps are calculated using a reference input image that facilitates assigning a cause for the attribution, e.g., by comparing the magnetic field evolution (Sun et al. 2022). Thus, here we use IG attribution maps to identify pixels, and hence the magnetic field features in the input, that are important for the CNN output. The IG attribution map for a given input image is calculated by integrating gradients in the CNN output along the path from a reference image. Formally,

$\begin{eqnarray}&&{L}^{f}\left(x,{x}_{0}\right)=\left(x-{x}_{0}\right)\times {\int }_{\alpha =0}^{1}\displaystyle \frac{\partial {Y}^{f}\left({x}_{0}+\alpha \times \left(x-{x}_{0}\right)\right)}{\partial x}d\alpha ,\end{eqnarray} \tag{ 1 }$

where x₀ is the reference image and Y^f is the CNN output for SHARP vector field feature f.

Figure 8 shows contour plots of typical IG attribution maps for a few example magnetograms from flaring ARs (bottom rows). The red/blue contours include regions of net positive/negative contribution toward the CNN output. The IG attribution maps for the three SHARP features—total unsigned flux, absolute net current helicity, and mean free energy density—are shown separately along with the reference magnetograms (top rows) used. In general, increasing/decreasing positive polarity flux corresponds to net positive/negative attribution. For total unsigned flux, almost all magnetic field regions, even relatively smaller regions with weaker magnetic fields, constitute a positive/negative attribution. In contrast, for absolute net current helicity and mean free energy density, only relatively larger and stronger magnetic field regions constitute a positive/negative attribution. For mean free energy, positive/negative attribution regions typically correspond to the uniformly increasing/decreasing positive flux. In the case of absolute net current helicity, attributions correspond to regions with "mixed" magnetic fields of the positive–negative polarities closely located. The appearance of a spurious magnetic field polarity inversion line (PIL) is a known artifact in the LOS magnetograms whenever the magnetic field inclination relative to the LOS exceeds 90° (Leka et al. 2017). We find that in many cases (e.g., HARPs 407, 3291, and 3311) when the PIL artifact exists for magnetic fields within penumbrae, it wrongly constitutes an important attribution. The misattribution results from the failure of the CNN to learn the PIL artifact (Sun et al. 2022) and, as a consequence, limits the accuracy of the reconstructed the vector field features. Additional examples of IG attribution maps are shown in Appendix D.

Zoom In Zoom Out Reset image size