A COMPARISON OF FLARE FORECASTING METHODS. I. RESULTS FROM THE "ALL-CLEAR" WORKSHOP

The data set provided for analysis contains sub-areas extracted from the full-disk magnetogram and continuum images from the SOHO/MDI. These extracted magnetogram and intensity image files, presented in FITS format, are taken close to noon each day, specifically daily image #0008 from the M_96m magnetic field data series, which was generally obtained between 12:45 and 12:55 UT, and image #0002 from the Ic_6h continuum intensity series, generally obtained before 13:00 UT.

To extract regions for a given day, the list of daily active-region coordinates was used, as provided by the National Oceanic and Atmospheric Administration (NOAA), and available through the National Center for Environmental Information (NCEI).¹⁷ The coordinates were rotated to the continuum image/magnetic field time using differential rotation and the synodic apparent solar rotation rate. A box was centered on the active region coordinates whose size reflects the NOAA listed size of the active region in micro-hemispheres (but not adjusted for any evolution between the issuance and time of the magnetogram or continuum image). A minimum active region size of 100 μH was chosen, corresponding to a minimum box size of 125'' × 125'' at disk center. The extracted box size was scaled according to the location on the solar disk to reflect the intrinsic reported area and to roughly preserve the area on the Sun contained within each box regardless of observing angle. This procedure most noticeably decreases the horizontal size near the solar east/west limbs, although the vertical size is impacted according to the region's latitude. The specifics of the scalings and minimum (and maximum) sizes were chosen empirically for ease of processing, and do not necessarily reflect any solar physics beyond the reported distribution of active region sizes reported by NOAA.

During times of high activity, this simple approach to isolating active regions using rectangular arrays often creates two boxes which overlap by a significant amount. Such an overlap, especially when it includes strong-field areas from another active region, introduces a double-counting bias into the flare-prediction statistics. To avoid this, regions were merged when an overlapping criterion based on the geometric mean of the regions' respective areas and the total was met. If two or more of these boxes overlapped such that ${ \mathcal A }/\sqrt{(\mathrm{area}\ {\mathrm{box}}_{1}\times \mathrm{area}\ {\mathrm{box}}_{2})}\gt 0.95$, where ${ \mathcal A }$ is the area of overlap (in image-grid pixels) of the two boxes, then the two boxes were combined into one region or "merged cluster." Clustering was restricted to occur between regions in the same hemisphere. No restrictions were imposed to limit the clustering, and in some cases more than two (at most six) regions were clustered together. In practice, in fewer than 10 cases, the clustering algorithm was over-ridden by hand in order to prevent full-Sun clusters, in which case the manual clustering was done in such a way as to separate clusters along areas of minimum overlap. The clustering was performed in the image-plane, although as mentioned above the box sizes took account of the location on the observed disk. When boxes merged, a new rectangle was drawn around them. The area not originally included in a single component active region's box was zeroed out. An example of the boxes for active regions and a two-region cluster for 2002 January 3 is shown in Figure 1. JPEG images of all regions and clusters similar to Figure 1 (top) are available at the workshop website. Note that a morphological analysis method based on morphological erosion and dilation has been used as a robust way to group or reject neighboring active regions (Zhang et al. 2010) although it is not implemented here.

Zoom In Zoom Out Reset image size

The final bounds of each extracted magnetogram file are the starting point for extracting an accompanying continuum file. The box was shifted to adjust for the time difference between the acquisition time of the magnetogram and that of the continuum file. If the time difference was greater than four hours, continuum extracted files were not generated for that day. Also, if the MDI magnetogram 0008 file was unavailable and the magnetogram nearest noon was obtained more than 96 minutes from noon, then neither the continuum nor the magnetogram extracted files were made for that day.

Data were provided as FITS files, with headers derived from the original but modified to include all relevant pointing information for the extracted area and updated ephemeris information. Additionally, the NOAA number, the NOAA-reported area (in μH) and number of spots, and the Hale and McIntosh classifications (from the most recent NOAA report) are included. The number of regions is included in the header, which is >1 only for clusters. In the case of clusters, all of the classification data listed above are included separately for each NOAA region in the cluster.

There was no additional stretching or re-projection performed on the data; the images were presented in plane-of-the-sky image coordinates. No pre-selection was made for an observing-angle limit, so many of the boxes are close to the solar limb. Similarly, no selection conditions were imposed for active-region size, morphology or flaring history, beyond the fact that a NOAA active region number was required. The result is 12,965 data points (daily extracted magnetograms) between 2000–2005.

Event lists were constructed from flares recorded in the NCEI archives. Three definitions of event were considered:

• 1.  
  at least one C1.0 or greater flare within 24 hr after the observation,
• 2.  
  at least one M1.0 or greater flare within 12 hr after the observation, and
• 3.  
  at least one M5.0 or greater flare within 12 hr after the observation.

Table 1 shows the flaring sample size for each event definition.¹⁸ No distinction is made between one and multiple flares above the specified threshold: a region was considered a member of the flaring sample whether one or ten flares occurred that satisfied the event definition. This flaring sample size defines the climatological rate of one or more flares occurring for each definition, which in turn forms the baseline forecast against which results are compared.

A subtlety arises because not all of the observations occur at exactly the same time of day, whereas the event definitions use fixed time intervals (i.e., 12 or 24 hr). As such it is possible for a particular flare to be double-counted. For example, if magnetograms were obtained at 12:51 UT on day #1 and then 12:48 UT on day#2, and a flare occurred at 12:50 UT on day #2 then both days would be part of the C1.0+, 24 hr flaring sample. This situation only arises for the 24 hr event interval and is likely to be extremely rare. However, it means that not all the events are strictly independent.

The NOAA active region number associated with each event from NCEI is used to determine the source of the flare. When no NOAA active region number is assigned, a flare is assumed not to have come from any visible active region, although it is possible that the flare came from an active region but no observations were available to determine its source. For large-magnitude flares, the vast majority of events (≈85% for M-class and larger flares, ≈93% for M5.0 and larger, and ≈93% for X-class flares) are associated with an active region, but a substantial fraction of small flares are not (≈38% for all C-class flares).

The request was made for every method to provide a forecast for each and every data set provided. Many methods, as alluded to in the descriptions in Appendix A, had restrictions on where it was believed they would perform reliably, and so each method did not provide a forecast in every case. The resulting variation in the sample sizes, as shown in the summary tables, is not fully accounted for in the skill scores reported (Section 3), meaning that direct comparisons between methods with different sample sets is not reliable. Even the H&KSS is not fully comparable among data sets if the samples of events and non-events are drawn from different populations, for example all regions versus only those regions with strong magnetic neutral lines.

To account for the different samples, three additional data sets are considered for performance comparison. The first is all data (AD), with an unskilled forecast (the climatology, or event rate) used if a method did not provide a forecast for that particular target. In this way, forecasts are produced for all data for all methods. However, this approach penalizes methods that produce forecasts for only a limited subset of data.

The second approach is to extract the largest subset from AD for which all methods provided forecasts. Two such maximum common data sets were constructed, one (MCD#1) for all methods except the event-statistics, while the second (MCD#2) included the additional event-statistics restrictions imposed by requiring at least ten prior events.²¹ One method (MSFC/Falconer, Appendix A.4) did not return results for C1.0+, 24 hr, so strictly speaking these should be null sets. However, for the purposes of this paper, that method was removed for constructing the MCDs for C1.0+, 24 hr. A drawback to the MCD approach arises for the methods that were trained using larger data samples, i.e., samples which included regions that were not part of the MCD. For methods that trained on AD, for example, many additional regions were used for training, while for methods with the most restrictive assumptions, almost all the regions used for training are included in MCD#1, hence the impact of using the maximum common data sets varies from method to method. In the case of MCD#1, where the primary restriction is on the distance of regions from the disk center, this may not be a large handicap since the inherent properties of the regions are not expected to change based on their location on the disk, although the noise in the data and the magnitude of projection effects do change. However, for MCD#2, the requirement of a minimum number of prior events means the samples are drawn from very different populations. For any method, training on a sample from one population then forecasting on a sample drawn from a different population adversely affects the performance of the method.

The magnitude of the changes in the populations from which the samples are drawn can be roughly seen in the changes in the event rates shown in Table 3. Between AD and MCD#1, the event rates change by no more than about 10%; between AD and MCD#2, the event rates change by up to an order of magnitude for the smaller event sizes such that MCD#2 for C1.0+, 24 hr is the only category with more events than non-events. The similarity in the event rates of AD and MCD#1 is consistent with the hypothesis that they are drawn from similar underlying populations. However, based on the changes in the event rates, it is fairly certain that MCD#2 is drawn from different populations than AD and MCD#1. Despite this difference, measures of the performance of the methods are presented for all three subsets to illustrate the magnitude of the effect and the challenge in making meaningful comparisons of how well the methods perform.

The Appleman skill scores and BSSs for each method are listed for each event definition in Tables 4–6, separately for each of the three direct-comparison data sets: AD, MCD #1 (without the event-statistics restrictions) and #2 (with the further restrictions from event statistics). As described in Section 3, the probability thresholds for generating the binary, categorical classifications used for calculating the ApSS were chosen to maximize the ApSS computed for each event definition using each method's optimal data set (see Appendix A).²² When a method produces more than one forecast (e.g., the generalized correlation dimension, Appendix A.7, which produces a separate forecast for each q value), the one with the highest BSS is presented. Using a different skill score to select which forecast is presented generally results in the same forecast being selected, so the results are not sensitive to this choice.

Table 4.  Performance on All Data with Reference Forecast

Parameter/	Statistical	C1.0+, 24 hr		M1.0+, 12 hr		M5.0+, 12 hr
Method	Method	ApSS	BSS	ApSS	BSS	ApSS	BSS
Beff	Bayesian	0.12	0.06	0.00	0.03	0.00	0.02
ASAP	Machine	0.25	0.30	0.01	−0.01	0.00	−0.84
BBSO	Machine	0.08	0.10	0.03	0.06	0.00	−0.01
WLSG2	Curve fitting	N/A	N/A	0.04	0.06	0.00	0.02
NWRA MAG 2-VAR	NPDA	0.24	0.32	0.04	0.13	0.00	0.06
$\mathrm{log}({ \mathcal R })$	NPDA	0.17	0.22	0.01	0.10	0.02	0.04
GCD	NPDA	0.02	0.07	0.00	0.03	0.00	0.02
NWRA MCT 2-VAR	NPDA	0.23	0.28	0.05	0.14	0.00	0.06
SMART2	CCNN	0.24	−0.12	0.01	−4.31	0.00	−11.2
Event Statistics, 10 prior	Bayesian	0.13	0.04	0.01	0.10	0.01	0.00
McIntosh	Poisson	0.15	0.07	0.00	−0.06	N/A	N/A

Note. An entry of N/A indicates that the method did not provide forecasts for this event definition.

Download table as:  ASCIITypeset image

Table 5.  Performance on Maximum Common Dataset #1

Parameter/	Statistical	C1.0+, 24 hr		M1.0+, 12 hr		M5.0+, 12 hr
Method	Method	ApSS	BSS	ApSS	BSS	ApSS	BSS
Beff	Bayesian	0.23	0.06	0.00	0.12	0.00	0.04
ASAP	Machine	0.29	0.32	0.07	0.05	0.00	−0.81
BBSO	Machine	0.24	0.30	0.12	0.17	0.00	−0.07
WLSG2	Curve fitting	N/A	N/A	0.14	0.24	0.00	0.10
NWRA MAG 2-VAR	NPDA	0.30	0.38	0.08	0.16	0.00	0.07
$\mathrm{log}({ \mathcal R })$	NPDA	0.29	0.38	0.07	0.21	0.00	0.08
GCD	NPDA	0.05	0.13	0.00	0.07	0.00	0.03
NWRA MCT 2-VAR	NPDA	0.29	0.37	0.09	0.21	0.04	0.08
SMART2	CCNN	0.27	−0.22	0.03	−4.46	0.00	−12.49
Event Statistics, 10 prior	Bayesian	N/A	N/A	N/A	N/A	N/A	N/A
McIntosh	Poisson	0.12	−0.03	0.00	−0.05	N/A	N/A

Note. An entry of N/A indicates that the method did not provide forecasts for this event definition.

Download table as:  ASCIITypeset image

Table 6.  Performance on Maximum Common Dataset #2

Parameter/	Statistical	C1.0+, 24 hr		M1.0+, 12 hr		M5.0+, 12 hr
Method	Method	ApSS	BSS	ApSS	BSS	ApSS	BSS
Beff	Bayesian	0.12	0.13	0.00	0.08	0.00	0.01
ASAP	Machine	0.22	0.22	0.14	0.09	0.00	−0.72
BBSO	Machine	0.23	0.17	0.17	0.11	0.00	−0.13
WLSG2	Curve fitting	N/A	N/A	0.19	0.18	0.00	0.08
NWRA MAG 2-VAR	NPDA	0.38	0.29	0.11	0.08	0.00	0.04
$\mathrm{log}({ \mathcal R })$	NPDA	0.23	0.26	0.10	0.13	0.00	0.05
GCD	NPDA	−0.47	−0.37	0.00	−0.10	0.00	−0.02
NWRA MCT 2-VAR	NPDA	0.23	0.25	0.11	0.10	0.05	0.04
SMART2	CCNN	0.15	0.18	0.03	−0.15	0.00	−1.47
Event Statistics, 10 prior	Bayesian	0.05	−0.21	0.06	0.13	0.00	−0.03
McIntosh	Poisson	0.02	−0.09	0.00	0.01	N/A	N/A

Note. An entry of N/A indicates that the method did not provide forecasts for this event definition.

Download table as:  ASCIITypeset image

Recalling that skill scores are normalized to unity, none of the methods achieves a particularly high skill score. No method for any event definition achieves an ApSS or BSS value greater than 0.4, and for the large event magnitudes, the highest skill score values are close to 0.2. Thus there is considerable room for improvement in flare forecasting. The skill scores for some methods are much lower than might be expected from prior published results. This is likely a combination of the data set provided here not being optimal for any of the algorithms, and variations in performance based on the particular time interval and event definitions being considered.

In each category of event definition and for most data sets, at least three methods perform comparably given a typical uncertainty in the skill score, so there is no single method that is clearly better than the others for flare prediction in general. For a specific event definition, some methods achieve significantly higher skill scores. There is a tendency for the machine learning algorithms to produce the best categorical forecasts, as evidenced by some of the highest ApSS values in Tables 4–6, and for NPDA to produce the best probability forecasts, as evidenced by some of the highest BSS values.

For rare events, most of the machine learning methods (ASAP, SMART2/CCNN, and to a lesser degree BBSO) produce negative BSS values, even when the value of the ApSS for the method is positive. This is likely a result of the training of the machine learning algorithms, which were generally optimized on one or more of the categorical skill scores with a probability threshold of 0.5. The maximum H&KSS is obtained for a probability threshold that is much smaller than 0.5 when the event rate is low (Bloomfield et al. 2012). When a machine learning algorithm is trained to maximize the H&KSS with a threshold of 0.5, it compensates for the threshold being higher than optimal by overpredicting (see the reliability plots in Appendix A). That is, by imposing a threshold of 0.5 and systematically overpredicting, a similar classification table is produced as when a lower threshold is chosen and the method is not overpredicting. The former results in higher categorical skill score values, but reduces the value of the BSS because the latter does not use a threshold but is sensitive to overprediction. In contrast, discriminant analysis is designed to produce the best probabilistic forecasts and so it tends to have high reliability (does not overpredict or underpredict). This results in good BSS and ApSS values.

Skill scores from the different methods for the three data sets and different event definitions are shown in Figure 3. Several trends are seen in the results. From the left panels, it can be seen that forecasting methods perform better on smaller magnitude events, whether evaluated based on the Brier or the Appleman skill score, with the most notable exceptions being ones for which the BSS is negative for C1.0+, 24 hr, including the event statistics method. The event statistics method uses the small events to forecast the large events, and thus is not as well suited to forecasting smaller events. The other exceptions are for the MCD#2, and thus are likely a result of a mismatch between the training and the forecasting data sets. The overall trend for most methods is likely due to the smaller sample sizes and lower event rates for the M1.0+, 12 hr and M5.0+, 12 hr categories, and holds for AD and both MCD sets. The smaller sample sizes make it more difficult to train forecasting algorithms, and the lower event rates result in smaller prior probabilities for an event to occur.

Zoom In Zoom Out Reset image size

The right panels of Figure 3 show that almost all methods achieve higher skill scores on MCD#1 than on AD, although the improvement is modest. For methods that did not provide a forecast for every region, this is simply a result of the methods making better predictions than climatology. For methods providing a forecast for every region, it suggests that restricting the forecasts to close to disk center improves the quality of the forecasts, although the effect is not large. However, methods trained on AD and then applied to MCD#1 may show a more substantial improvement if trained on MCD#1. One of the main restrictions in many methods is the distance from disk center due to projection effects. Thus it is likely that improved results would be achieved by using vector magnetograms.

Most methods also achieve higher skill scores on MCD#2 than on AD, but a considerable fraction have lower skill scores, and there is more variability in the changes between AD and MCD#2, as measured by the standard deviation of the change in skill score, than in the changes between AD and MCD#1. This supports the hypothesis that AD and MCD#1 draw from similar populations, while MCD#2 draws from significantly different populations. The ranking of methods changes somewhat between AD and MCD#1, but more significantly between AD and MCD#2. This highlights the danger of using a subset of data to compare methods, particularly when the subset is drawn from a different population than the set used for training.

As discussed, different skill scores emphasize different aspects of performance. This is demonstrated by the results for the B_eff and the BBSO methods shown in Table 5, for C1.0+, 24 hr using MCD#1. The forecasts using these methods result in essentially the same values of the ApSS of 0.23 and 0.24, respectively. However, the corresponding BSS values of 0.06 and 0.30 are quite different. This is graphically illustrated in comparing Figure 4 with Figure 2, right. The B_eff method (Figure 4) systematically overpredicts for forecast probabilities less than about 0.6, but slightly underpredicts for larger forecast probabilities. Thus the probabilistic forecasts result in a small BSS, but by making a categorical forecast of an event for any region with a probability greater than 0.55, the method produces a much higher ApSS. The SMART/CCNN method produces a similar systematic under- and overprediction (see Appendix A.9) while most methods (e.g., the BBSO method shown in Figure 2) have little systematic over- or underprediction, so the Brier and Appleman skill scores are similar.

Zoom In Zoom Out Reset image size

Inspection of reliability plots for a single method (Figure 5) shows a phenomenon common to most methods: the maximum forecast probability typically decreases with increasing event size, so most methods only produce low-probability forecasts for, say, M5.0+, 12 hr event definitions. This explains the small values of the ApSS for larger event definitions as very few or no regions have high enough forecast probabilities to be considered a predicted event in a categorical forecast. It also suggests that all-clear forecasts have more promise than general forecasts. However, attention to the possibility of missed events would be critical from an operations point of view.

Zoom In Zoom Out Reset image size

All groups were given the same data, and many computed the same parameters. However, implementations differ significantly, so values for the same parameter are substantially different in some instances. In other cases, two parameters computed using completely different algorithms lead to parameter values that are extremely well correlated.

The total unsigned magnetic flux of a region, $\sum | {B}_{z}|$, is often considered a standard candle for forecasting. Larger regions have long being associated with greater propensity for greater-sized events, and the total flux is a direct measure of region size, hence it provides a standard for flare-forecast performance. Four groups calculated the total magnetic flux for this exercise, and provided the value for each target region. By necessity, since the MDI data provide only the line-of-sight component of the field vector, approximations were made, which varied between groups, and one group (NWRA) calculated the flux in two ways, using different approximations. There are also different thresholds used to mitigate the influence of noise, and different observing-angle limits beyond which some groups do not calculate this parameter. How large are the effects of these different assumptions and approximations on the inferred value of the total unsigned flux, and what is the impact on flare forecasting?

The distributions of the total flux parameter for the four groups are shown in Figure 6, estimated from the MCD#1, which includes the same regions. There are considerable differences among the distributions. In general, the NWRA values of the flux are larger than the other groups, although the SMART values have similar peak values, but with a tail to lower values than seen in the NWRA flux distribution. The NWRA distributions are also more sharply peaked than the other groups' distributions.

Zoom In Zoom Out Reset image size

To understand these differences, Figure 7 shows the values for each group plotted versus the NWRA values, for all regions for which the flux was computed by that group. The NWRA method is used as the reference because it has a value for every region in the data set. The values for regions in the MCD#1 (black) generally show much less scatter than when all available regions are considered. However, there are systematic offsets in values for NWRA versus BBSO and MSFC even for the MCD#1. The MCD#1 only includes regions relatively close to the disk center, so much of the scatter may be a result of how projection effects are accounted for. This is particularly apparent for the SMART values.

Zoom In Zoom Out Reset image size

How much influence do the variations in the total flux resulting from different implementations have on forecasting flares? To separate out the effects of the statistical method, NPDA was applied to all of the total-flux related parameters. The forecast performance solely as a single-variable parameter with NPDA is summarized in Table 7. Because different approximations may not influence how well the determination of total flux works toward the limb, entries are included using both AD files (without restriction), and using only the MCD.

Despite differences in the inferred values of the total unsigned flux, the resulting skill scores for the MCD#1 are effectively the same for all the implementations. The AD forecasting results use a climatology forecast where a method did not provide a total flux measurement (because it was beyond their particular limits, for example). This difference can be seen in the variation between the results for AD and MCD#1, in particular for the BBSO implementation, which had the most restrictive condition on the distance from the disk center for computing the total unsigned flux. Overall, the different implementations result in significantly different skill scores for AD.

There is also evidence that different approximations applied to the ${B}_{{\rm{los}}}$ to retrieve an estimate of ${B}_{{\rm{z}}}$, the radial field, can make a difference. The NWRA potential field implementation performs nearly as well on AD as on the MCD#1, and better than the other implementations on AD. This suggests that, for flare forecasting, the potential field approximation is better than the μ-correction when only measurements of the line-of-sight component are available.

Details of the implementation are also important for other parameters. Figure 8 illustrates the impact of implementation on the measure ${ \mathcal R }$ of the strong gradient polarity inversion lines proposed by Schrijver (2007), and on the measure B_eff of the connectivity of the coronal magnetic field proposed by Georgoulis & Rust (2007) as computed by different groups. The difference is more pronounced in the connectivity measure. Although the same mathematical formula for B_eff is used, the differences are due to the distinct approaches followed for partitioning a magnetogram to determine the point sources, and for inferring the connectivity matrix for a given set of sources, as noted in Appendix A.8.

Zoom In Zoom Out Reset image size

In contrast to the differences in a parameter as implemented by different groups, parameters proposed by different researchers can also be strongly correlated. For example, the parameter ${ \mathcal R }$ proposed by Schrijver and the parameter WL_SG2 proposed by Falconer are both measures of strong gradient polarity inversion lines. The methods by which they are calculated are quite different, but the linear correlation coefficient between the two is r = 0.95 (Figure 9, left). Perhaps even more surprising is that WL_SG2 is also strongly correlated with the parameter B_eff of Georgoulis & Rust (2007), which is a measure of the connectivity of the coronal magnetic field (Figure 9, right).

Zoom In Zoom Out Reset image size

Two conclusions can be drawn from this exercise. First, implementation details can greatly influence the resulting parameter values, although this makes surprisingly little difference in the forecasting ability of the quantities considered. Second, there may be a limited amount of information available for flare forecasting from only the line-of-sight magnetic field without additional modeling. Even if additional modeling is used, such as when the coronal connectivity is used to determine B_eff, solar active regions that are small tend to be simpler, larger regions tend to be more complex, and differentiating those with imminent flare potential remains difficult.

The analysis presented in Georgoulis & Rust (2007) describes the coronal magnetic connectivity using the effective connected magnetic field strength B_eff. The B_eff parameter is calculated following inference of a connectivity matrix in the magnetic-flux distribution of the target active region (employing ${B}_{{\rm{z}}}\approx {B}_{{\rm{los}}}/\cos \theta$, where θ is the angle from the disk center). An additional multiplicative correction factor of 1.56 (plage) and 1.45 (sunspots) was applied to correct for systematic insensitivity (Berger & Lites 2003).

The resulting flare forecast relies on an event probability using Bayesian inference and Laplace's rule of succession (Jaynes & Bretthorst 2003, p. 154). In particular, the steps used to produce a forecast are:

• 1.  
  Partition the photospheric vertical magnetic field into (e.g.,) N₊ positive-polarity and N₋ negative-polarity magnetic-flux concentrations (Figure 10, left), generally following Barnes et al. (2005). Particulars include: 50 G noise threshold, minimum partition flux of 10²⁰ Mx, minimum partition size of 40 pixels, and a slight smoothing (mean-neighborhood by a factor 2) applied prior to partitioning, only in order to draw smoother partition outlines.
• 2.  
  Determine a connectivity matrix, ψ_ij (Figure 10, right), whose elements are the magnetic flux connecting the pairs of sources (i, j), by using simulated annealing to minimize the functional

  $$\begin{eqnarray*}&&F=\displaystyle \sum _{i,j}\left(\displaystyle \frac{| {{\boldsymbol{x}}}_{i}-{{\boldsymbol{x}}}_{j}| }{{R}_{{\rm{\max }}}}+\displaystyle \frac{| {{\rm{\Phi }}}_{i}^{\prime }+{{\rm{\Phi }}}_{j}^{\prime }| }{| {{\rm{\Phi }}}_{i}^{\prime }| +| {{\rm{\Phi }}}_{j}^{\prime }| }\right),\end{eqnarray*}$$

  where ${{\boldsymbol{x}}}_{i}$, ${{\boldsymbol{x}}}_{j}$ are the vector positions of the flux-weighted centroids of two opposite-polarity partitions i and j ($i\equiv \{1,\,\ldots ,\,{N}_{+}\}$, $j\equiv \{1,\,\ldots ,\,{N}_{-}\}$) with respective flux content ${{\rm{\Phi }}}_{i}^{\prime }$ and ${{\rm{\Phi }}}_{j}^{\prime }$ and R_max is a constant, maximum distance within the studied magnetogram, typically its diagonal length. The implementation of R_max in the functional F is a refinement over the initial approach of Georgoulis & Rust (2007). Another refinement is the introduction of a mirror "flux ring" at a distance well outside of the studied magnetogram, that makes the flux distribution exactly balanced prior to annealing (Georgoulis et al. 2012). This step consists of introducing a ring of mirror flux (as much positive-/negative-polarity flux as the negative-/positive-polarity flux of the active region) at large distances from the region, typically three times larger than the largest dimension of the immediate region. The ring of flux participates in the connectivity process via the simulated annealing. Magnetic connections between active-region flux patches and the flux ring are considered open, that is, closing beyond the confines of the active region. These connections are not taken into account in the calculation of B_eff.
• 3.  
  Define the effective connected magnetic field as

  $$\begin{eqnarray*}&&{B}_{{\rm{eff}}}=\displaystyle \sum _{i=1}^{{N}_{+}}\displaystyle \sum _{j=1}^{{N}_{-}}\displaystyle \frac{{\psi }_{{ij}}}{| {{\boldsymbol{x}}}_{i}-{{\boldsymbol{x}}}_{j}{| }^{2}}.\end{eqnarray*}$$
• 4.  
  Construct the predictive conditional probability of an event above a certain size according to Laplace's rule of succession:

  $$\begin{eqnarray*}&&{P}_{{\rm{flare}}}^{{\rm{th}}}=\displaystyle \frac{{N}_{{\rm{mag(flare)}}}({B}_{{\rm{eff}}}\gt {B}_{{\rm{eff}}}^{{\rm{th}}})+1}{{N}_{{\rm{mag}}}({B}_{{\rm{eff}}}\gt {B}_{{\rm{eff}}}^{{\rm{th}}})+2},\end{eqnarray*}$$

  where ${N}_{{\rm{mag(flare)}}}$ is the number of event-producing magnetic structures (magnetograms) with B_eff greater than a given threshold ${B}_{{\rm{eff}}}^{{\rm{th}}}$ (for a particular event definition), and N_mag is the total number of magnetograms with B_eff greater than the same threshold. Increasing B_eff-thresholds are successively selected, and the resulting curve of ${P}_{{\rm{flare}}}^{{\rm{th}}}$ as a function of ${B}_{{\rm{eff}}}^{{\rm{th}}}$ (for the targeted event definition) is then fitted by a sigmoidal curve; this curve returns the flaring probability for an incoming B_eff measure. The lowest threshold used for a particular event definition is the minimum B_eff found for all magnetograms for which an event (as defined) was recorded. Probabilities for magnetograms with B_eff less than this cutoff value are set to zero.

Zoom In Zoom Out Reset image size

This approach does not use magnetic-field extrapolation. Instead, the coronal model relies on the minimum value of the functional F which finds the shortest connections between opposite-polarity flux concentrations. The B_eff values were calculated for all but one data set. However, following Georgoulis & Rust (2007), a limit of ±41° from the disk center is imposed to minimize projection-effect artifacts for the measures of this method's optimal performance shown in Table 8, thus the values differ slightly from those presented in Tables 4–6 as a different subset of data is used here. All values of B_eff are used in the later AD comparisons (Section 5.2).

The skill-score results are unusual in that there is a large discrepancy in the performance based on which skill score is used to evaluate the method, with the HSS and H&KSS skill score values being much larger than the ApSS or BSS for C1.0+, 24 hr. The reliability plots (Figure 11, top) show a systematic over-prediction for lower probabilities and underprediction for higher probabilities. As is typical for most methods, the maximum probability forecast decreases as the event threshold increases, while the maximum H&KSS value and the ROC curve (Figure 11, bottom) improve with increasing event-threshold magnitude.

Zoom In Zoom Out Reset image size

A real-time automated flare prediction tool has been developed at the University of Bradford/Centre for Visual Computing (Colak & Qahwaji 2008, 2009). The Automated Solar Activity Prediction (ASAP)²⁴ system uses the following steps to predict the likelihood of a solar flare:

• 1.  
  A feature-recognition system generates McIntosh classifications (McIntosh 1990) for active regions from MDI white-light images.
• 2.  
  A machine learning system is trained using these classifications and flare event databases from NCEI.
• 3.  
  New data are then used to generate real-time predictions.

The system relies upon both MDI magnetic and white-light data, and hence is unable to make a prediction for those data for which the white-light data are unavailable. Generally, ASAP generates a McIntosh classification, as determined by step #1 above. This is difficult for the active region clusters made up of multiple active regions in the database we presented. Thus, the recorded McIntosh classifications included in the file headers were used for the predictions (i.e., step #1 was essentially skipped), and forecasts were made for the entire database. However, the performance may have been less than optimal due to this peculiarity.

The results for ASAP (Table 9) compared to the average show higher values for HSS and H&KSS skill but lower for ApSS and BSS, a common result for machine-learning based forecasting algorithms which are typically trained to produce the largest HSS or H&KSS. The reliability plots (Figure 12, top) show a systematic over-prediction for the larger event thresholds, and ROC plots (Figure 12, bottom) show lower probability of detection for high false alarm rates than most other methods.

Zoom In Zoom Out Reset image size

Another approach that uses a machine learning technique as the statistical forecasting method has been developed at the New Jersey Institute of Technology (Yuan et al. 2010, 2011). The steps in this method are:

• 1.  
  Compute three parameters describing an active region:

  • i.  
    total unsigned magnetic flux, computed using only the pixels for which the absolute value of the field is greater than the mean value plus three times the standard deviation of all field in the area under consideration;
  • ii.  
    the length of the strong-gradient neutral line (above 50 G Mm⁻¹);
  • iii.  
    the total magnetic energy dissipation E_diss = $\int 4$[(∂${B}_{{\rm{z}}}$/∂x)² + (∂${B}_{{\rm{z}}}$/∂y)²] + 2(∂${B}_{{\rm{z}}}$/∂x + ∂${B}_{{\rm{z}}}$/∂y)²dA, following Abramenko et al. (2003).
• 2.  
  Use ordinal logistic regression and support vector machines (SVMs) to make predictions.

An SVM is a supervised learning method used for classification (Boser et al. 1992), whose principle is to minimize the structural risk (Vapnik 2000). An SVM tries to find a plane in an n-dimensional space that separates input data into two classes. The larger the distance from the plane to the two different classes of data points in the n-dimensional space, the smaller the classification error (Cortes & Vapnik 1995). The results presented here used the open-source SVM implementation called LIBSVM (Chang & Lin 2011).

The summary of results is given in Table 10. The method used and made forecasts only on extracted data that had a single NOAA active region within ±40° of the disk center, reducing the sample to less than half that provided. This restriction also presents an example of a method whose requirements are not well met by the data used for this workshop, and thus whose skill scores may be penalized as a result. When calculating the flux, only pixels for which the absolute value of the field is greater than the mean value plus three times the standard deviation of all field in the area under consideration are included. The summary plots (Figure 13) show a weak trend toward overprediction at higher probabilities, but with larger error bars due to the smaller sample sizes. The ROC plots are fairly typical.

Zoom In Zoom Out Reset image size

Falconer et al. (2008) presented a prediction method based on parameterizing active region magnetic morphology that was first applied to CME prediction (Falconer et al. 2002, 2003) but more recently extended to flares. The modern version has been implemented in a code called MAG-4 (Falconer et al. 2011), running at the Community Coordinated Modeling Center and elsewhere. The results presented here are from a predecessor of MAG-4.

Two parameters are calculated. A measure of the free magnetic energy is estimated based on the presence of strong gradient neutral lines:

$$\begin{eqnarray}&&{{WL}}_{\mathrm{SG}2}={\int }_{{NL}}{({\boldsymbol{\nabla }}{B}_{{\rm{los}}})}^{2}\,{dl}.\end{eqnarray} \tag{ 12 }$$

An example of strong gradient neutral lines is shown in Figure 14. The measure WL_SG2 is supplemented with the total unsigned magnetic flux (for pixels above 100 G) in the target active region. No correction to the ${B}_{{\rm{los}}}$ data is performed. The utility of WL_SG2 was shown in Falconer et al. (2008) where parameterizations using vector magnetogram data were favorably compared to the WL_SG2 proxy, which can be calculated from line-of-sight data. In this study, the parameterizations for both WL_SG2 and total magnetic flux Φ_tot were calculated for all extracted data sets, but the method is generally restricted to regions within 30° of the disk center. Table 11, summarizing the method's performance, is restricted to this subset of predictions. It is expected that line-of-sight data beyond 30° do not provide reliable estimates of WL_SG2.

Zoom In Zoom Out Reset image size

With these parameterizations, a least-squares power-law fit to the event rates as a function of WL_SG2 and total flux was constructed for each event definition using only the data within the 30° limit. These event rates were converted into probabilities as a function of the length of time t of the forecast interval, assuming that the event rate is constant over the larger forecast interval (in this case 2000–2005). The conversion from event rate to probability assumes Poisson statistics (Moon et al. 2001; Wheatland 2001): the probability of an event is $P(t)=100\times (1-{e}^{-\lambda t})$ where λ is the flaring rate for the particular event definition.

The rate-fitting algorithm is best for larger flares, and so no forecasts were made for the C1.0+, 24 hr events. It has been found that combining this strong gradient neutral line with secondary measures (e.g., total magnetic flux, and more recently with prior flare history, Falconer et al. 2012) is likely to give more accurate forecasts, but these secondary measures are not included for the forecasts presented here. The results for the workshop data set are shown in Table 11, and only include the regions with higher confidence. The values of the BSS are among the best, although the small sample sizes are reflected in the larger error bars in the reliability plots (Figure 15, top). The ROC curves show that the probability of detection remains near one to relatively small values of the false alarm rate, hence this is one of the better methods for issuing all-clear forecasts.

Zoom In Zoom Out Reset image size

The NWRA moment analysis parameterizes the observed magnetic field and its spatial derivatives using the first four moments (mean, standard deviation, skew and kurtosis), plus totals and net values when appropriate (Leka & Barnes 2003a). The moments of the observed field strength describe its distribution, while an approximation of the total flux indicates how large the region is based on both its area and field strength. The higher-order moments can be sensitive to the presence of small areas of complex field that are missed in the lower-order moments and summations (Leka & Barnes 2003a).

The default boundary is the observed line-of-sight data. A correction is applied to these data: ${B}_{{\rm{z}}}(x,y)\approx {B}_{{\rm{los}}}(x,y)/\cos \theta (x,y)$, meaning that each observed point is divided by the cosine of its observing angle, θ, rather than an average observing angle or the observing angle of the center of the field of view. This correction is especially important at large observing angles. Additionally, a magnetic-field boundary that more closely approximates the radial field is prepared by computing the potential field that matches the line-of-sight observed component (e.g., Alissandrakis 1981) and using the resulting radial field as a boundary for some parameters. This approach has the effect of mitigating the appearance of the magnetic neutral lines that are solely a manifestation of projection effects. For the potential-field approximation of the radial field, no additional correction using the observing angle is needed. For both boundaries, however, a limit of $\cos \theta \gt 0.1$ is imposed: parameters are not computed for observing angles ≳85°.

The parameters considered are the subset of those presented in Leka & Barnes (2003a) which can be computed using solely the line-of-sight component of the field. Over 50 parameters are considered, describing the photospheric field distribution using three basic categories:

• 1.  
  The distribution of the approximated radial (vertical) component of the magnetic field.
• 2.  
  The horizontal gradient of that distribution $| {\boldsymbol{\nabla }}{B}_{{\rm{z}}}|$.
• 3.  
  The character of inferred magnetic neutral lines.

The last category includes a variation on the ${ \mathcal R }$ parameter described in Appendix A.6.

The forecast is produced with discriminant analysis (DA, e.g., Kendall et al. 1983; Silverman 1986). This statistical approach classifies new measurements as belonging to one of two populations by dividing parameter space into two regions based on where the probability density of one population (e.g., flaring regions) exceeds the other (flare-quiet regions). A set of new measurements, i.e., a new active region, that falls on the appropriate side of the division is then predicted to flare (see Figure 16). Using Bayes's theorem, the probability that a new measurement belongs to a given population can be estimated from the probability density estimates (Barnes et al. 2007). Given an accurate representation of the probability density of a parameter, DA will maximize the overall accuracy of predictions (see Leka & Barnes 2003b, 2007, for examples).

Zoom In Zoom Out Reset image size

The probability density is typically either assumed to be Gaussian, which results in a linear discriminant function, or it is estimated nonparametrically. For the results presented here, a nonparametric density estimate was made using the Epanechnikov kernel with the smoothing parameter chosen optimal for a Gaussian distribution (Silverman 1986). One strength of DA is that multiple variables can be considered simultaneously. Two-variable combinations are employed here, as noted in the tables. Cross-validation is employed by default in order to remove bias in the skill scores.

The results for the variable combination with the highest BSS for each event definition are given in Table 12. Other variable combinations may produce higher values of other skill scores, but only the results for the variable combination with the highest BSS are listed here. The reliability plots for Schrijver's implementation of ${ \mathcal R }$ and one-variable DA show a slight tendency for over-prediction for the M1.0+, 12 hr threshold at higher probabilities, but less of such a trend for M5.0+, 12 hr (Figure 17, top).

Zoom In Zoom Out Reset image size

Table 12.  Optimal Performance Results: NWRA Field Parameterizations, Non-parametric Discriminant Analysis, Two-variable Combinations

Event	Sample	Event	RC	HSS	ApSS	H&KSS	BSS
Definition	Size	Rate	(Threshold)	(Threshold)	(Threshold)	(Threshold)
C1.0+, 24 hra	12965	0.201	0.85 (0.48)	0.50 (0.35)	0.24 (0.48)	0.56 (0.22)	0.32
M1.0+, 12 hrb	"	0.031	0.97 (0.42)	0.29 (0.12)	0.04 (0.42)	0.58 (0.03)	0.13
M5.0+, 12 hrc	"	0.007	0.99 (0.47)	0.20 (0.07)	0.00 (0.47)	0.72 (0.01)	0.06

Notes.

^aVariable combination: $[\sigma ({\boldsymbol{\nabla }}| {B}_{{\rm{los}}}| ),\ \mathrm{log}({{ \mathcal R }}_{{\rm{nwra}}}^{{\rm{pot}}})$]. ^bVariable combination: $[\varsigma ({\boldsymbol{\nabla }}{B}_{{\rm{los}}}),\ \sigma ({B}_{{\rm{z}}}^{{\rm{pot}}})$]. ^cVariable combination: $[\varsigma ({\boldsymbol{\nabla }}{B}_{{\rm{los}}}),\ \sigma ({B}_{{\rm{z}}}^{{\rm{pot}}})$].

Download table as:  ASCIITypeset image

In Schrijver (2007), a parameter ${ \mathcal R }$ was proposed as a proxy for the emergence of current-carrying magnetic flux, i.e., magnetic systems with significant free magnetic energy which would be carried through the photosphere and into the solar corona, enabling solar flares. The ${ \mathcal R }$ (and, specifically, $\mathrm{log}({ \mathcal R }$)) parameter was proposed and demonstrated as useful for forecasting solar flares in Schrijver (2007) by determining an empirical threshold above which flare activity was significant.

The parameter ${ \mathcal R }$ was computed from the line-of-sight magnetic field maps using the following steps:

• 1.  
  Dilate bitmaps of the magnetograms where the positive or negative flux density exceeds a threshold (150 Mx cm⁻²).
• 2.  
  Define high-gradient polarity-separation lines as areas where the bitmaps overlap.
• 3.  
  Convolve the resulting high-gradient polarity-separation line bitmap with a Gaussian to obtain a "weighting map."
• 4.  
  Obtain ${ \mathcal R }$ by multiplying the weighting map by the unsigned line of sight field and computing the total.

For the results here, the ${ \mathcal R }$ parameter was calculated as part of the NWRA magnetic field analysis (Appendix A.5), but is also included in the parameterizations by other groups (e.g., SMART, see Appendix A.9), with slightly different implementations (see Section 5.3). Within the NWRA magnetic parameterization, ${ \mathcal R }$ is calculated first so as to replicate the parameter in Schrijver (2007): ${B}_{{\rm{los}}}$ is used directly, and a fixed width of 10 MDI pixels is used for the Gaussian used in the convolution, to identify an area within roughly 15 Mm of a magnetic neutral line. The targets are restricted to those within 45° of disk center, which limits the sample size. An image of the boundary of the inferred high-gradient neutral lines from this replication of Schrijver (2007) is shown in Figure 18.

Zoom In Zoom Out Reset image size

The second implementation within the NWRA magnetic parameterization uses the ${B}_{{\rm{z}}}$ potential-field boundary as described in Appendix A.5, and varies the width of the convolution function such that it effectively preserves the target 15 Mm physical distance on the Sun and performs all calculations in a helioplanar coordinate system. A comparison of this implementation with the original is also shown in Figure 18.

For the results here, predictions using ${ \mathcal R }$ were made using one-variable NPDA with cross-validation (Appendix A.5), and the resulting skill scores are shown in Table 13. The NWRA implementation is not included here, but it is included in one high-performing two-variable combination in Table 12. The reliability plots (Figure 19, top) show a slight underprediction at the lower probabilities for C1.0+, 24 hr, and little trending otherwise with the exception of points in the highest predicted probability bins for both M1.0+, 12 hr and M5.0+, 12 hr. The ROC curves (Figure 19, bottom) show that the probability of detection for the M1.0+, 12 hr and M5.0+, 12 hr event definitions remain close to one for relatively small values of the false alarm rate, thus this is another method that may be well suited to issuing all-clear forecasts.

Zoom In Zoom Out Reset image size

The morphology of a flux concentration (active region) can be described by the fractal dimension and related quantities. This approach is useful for understanding the underlying influence of turbulence on solar structures (McAteer et al. 2010).

The generalized correlation (akin to a fractal) dimension D_BC of an active region is computed using images of the magnetic field distribution. Images were processed in the same manner as McAteer et al. (2005b), which includes applying multiplicative corrections to stronger/weaker areas and an observing-angle correction ${B}_{{\rm{z}}}\approx {B}_{{\rm{los}}}/\cos \theta$. The images were also subjected to a thresholding algorithm in the same manner as McAteer et al. (2005b) and a box counting algorithm (with cancellation) was applied to the resulting binary images. Only regions within θ < 60° were considered. The generalized correlation dimensions were calculated for q values of 0.1, 0.5, 1.5, 2.0, 8.0 (q referring to the "q-moment," which governs the influence of strong versus weak areas in the D_BC measure, see McAteer et al. 2010).

The fractal dimension is a parameterization of the active region, and does not include a statistical flare prediction method per se. In McAteer et al. (2005a), the fractal dimension was shown to be related to a region's flare productivity: lower limits of D_BC = 1.2(1.25) were found to be required for M-class (X-class) flares (with q = 8). Extensions of this work involving the multifractal spectrum (Conlon et al. 2008, 2010) were not computed for this study. Georgoulis (2012) show that one multifractal algorithm is resolution-dependent. However, McAteer (2015) and McAteer et al. (2016) show that a multifractal analysis that uses the wavelet transform modulus maxima approach, performed over a series of images, may produce a useful measure of the energy stored in the coronal magnetic field.

For the results here, predictions were made using one-variable NPDA with cross-validation (Appendix A.5) with these fractal-related parameterizations. The best-performing parameterization is shown in Table 14, and is the q = 8.0 variable as expected from prior work, and thus is subsequently the only variable presented here. The performance as measured by the BSS is comparatively higher than for the H&KSS as expected from the NPDA. The ROC curves (Figure 20) show the usual improvement with increasing event magnitude, but overall show worse performance than most other methods.

Zoom In Zoom Out Reset image size

A corona with a very complex magnetic topology is one which will more readily allow initiation of fast magnetic reconnection, and hence an energetic event. Indeed, Barnes & Leka (2006) demonstrated that parameters derived using a magnetic charge topology coronal model (MCT; Barnes et al. 2005, and reference therein) out-performed those describing the photospheric field in distinguishing flare-imminent from flare-quiet times for seven active regions (see Table 5 of Barnes & Leka 2006.) To implement the MCT model, the following steps are followed:

• 1.  
  Partition the photospheric field into magnetic field concentrations and represent each partition as a single point source, located at the flux-weighted center of the partition, ${{\boldsymbol{x}}}_{i}$, with magnitude equal to the flux in the partition, Φ_i.
• 2.  
  Represent the coronal magnetic field as the potential field associated with this collection of point sources.
• 3.  
  Define the magnetic connectivity matrix as the amount of magnetic flux in each coronal connection, ψ_ij, and estimate it by tracing field lines.
• 4.  
  Calculate the location of magnetic null points (places where the magnetic field vanishes) using the method described in Barnes (2007).

An example of the results is shown in Figure 21. Compared to the partitioning used to compute B_eff (Appendix A.1), this implementation represents the plage areas with fewer sources. Despite this, the resulting connectivity matrix typically has more connections per source.

Zoom In Zoom Out Reset image size

Almost 50 parameters were calculated based on moments and totals of the distributions of properties of the sources and the connectivity matrix. These properties characterize quantities such as the number, orientation, and flux in the connections (Barnes et al. 2005; Barnes & Leka 2006, 2008), and include ${\phi }_{{\rm{tot}}}=\sum {\psi }_{{ij}}/| {{\boldsymbol{x}}}_{i}-{{\boldsymbol{x}}}_{j}|$, the total of each connection's distance-weighted connection flux, as well as moments of the distribution of ${\psi }_{{ij}}/| {{\boldsymbol{x}}}_{i}-{{\boldsymbol{x}}}_{j}|$, and the moments of the distribution of ${\psi }_{{ij}}/| {{\boldsymbol{x}}}_{i}-{{\boldsymbol{x}}}_{j}{| }^{2}$. The total of this quantity, denoted ${\phi }_{2,\mathrm{tot}}$, is essentially the same as B_eff (Appendix A.1; Georgoulis & Rust 2007, but see the discussion in Barnes & Leka 2008) except in how the connectivity matrix is inferred. NPDA with cross-validation is used to make a prediction (Appendix A.5, Barnes & Leka 2006), with the resulting skill scores for the variable combinations that resulted in the best BSS values shown in Table 15; the M5.0+, 12 hr BSS is one of the better quoted in this study, but still not overly impressive. The reliability plots (Figure 22, top) show a slight tendency for underprediction at high probabilities for the M1.0+, 12 hr and M5.0+, 12 hr categories. The ROC curves (Figure 22, bottom) are fairly typical.

Zoom In Zoom Out Reset image size

Table 15.  Optimal Performance Results: Magnetic Charge Topology, Two-variable Combinations, Non-parametric Discriminant Analysis

Event	Sample	Event	RC	HSS	ApSS	H&KSS	BSS
Definition	Size	Rate	(Threshold)	(Threshold)	(Threshold)	(Threshold)
C1.0+, 24 hra	12965	0.201	0.85 (0.55)	0.48 (0.29)	0.23 (0.55)	0.52 (0.17)	0.28
M1.0+, 12 hrb	"	0.031	0.97 (0.50)	0.31 (0.19)	0.05 (0.50)	0.61 (0.03)	0.14
M5.0+, 12 hrb	"	0.007	0.99 (0.33)	0.20 (0.08)	0.00 (0.33)	0.70 (0.01)	0.06

Notes.

^aVariable combination: $[{\phi }_{2,\mathrm{tot}},\sigma ({\phi }_{2})$]. ^bVariable combination: $[{\phi }_{{\rm{tot}}},\overline{| {{\boldsymbol{x}}}_{i}^{{\rm{null}}}-{{\boldsymbol{x}}}_{j}^{{\rm{null}}}| }$]. ^cVariable combination: $[{\phi }_{{\rm{tot}}},\overline{| {{\boldsymbol{x}}}_{i}^{{\rm{null}}}-{{\boldsymbol{x}}}_{j}^{{\rm{null}}}| }$].

Download table as:  ASCIITypeset image

The Solar Monitor Active Region Tracker (SMART) method is incorporated into the well-known SolarMonitor.org resource, as the SMART2 code package, which performs a combination of detecting, tracking, and characterizing active regions (Higgins et al. 2011). For the present comparison, the provided patches are used rather than the full-disk MDI data, but all subsequent analysis proceeds as described in Higgins et al. (2011). This includes smoothing and thresholding to differentiate plage from spot regions, and an observing angle correction similar to other methods, but no additional multiplicative factor correction.

Since only one magnetogram per day was provided, the parameters that characterize temporal variations were not calculated. Twenty parameters were calculated including:

• 1.  
  Area of the final dilated mask, and relevant totals, net, min/max and moments of the flux within the final dilated mask, and the area of the mask.
• 2.  
  Maximum, mean, and median of the field spatial gradients.
• 3.  
  Length of the line of sight polarity separation lines and strong-gradient neutral lines.
• 4.  
  Amount of flux near strong-gradient polarity separation lines and non-potentiality gauges following Schrijver (2007) (Appendix A.6) and Falconer et al. (2008) (Appendix A.4), with two different thresholds applied.

These parameters are used to make forecasts using the cascade correlation neural networks method (CCNN; Qahwaji et al. 2008) implemented in Ahmed et al. (2013). For the machine-learning task, all active regions with a finite value of the parameters were considered, providing 8137 forecasts, and the six years are alternately rotated to enable training on five years of data at a time, in a jackknife manner (e.g., Efron & Gong 1983). This training is performed separately for each flare event definition, using all parameters simultaneously. The CCNN is trained to optimize H&KSS. The skill results shown in Table 16 confirm this, with the highest H&KSS and HSS scores of any method, especially for the larger thresholds—but with the worst BSS skill of any method, as well.

The reliability plots (Figure 23, top) show that the method is substantially overpredicting events in almost all cases. This is likely a consequence of optimizing on the H&KSS: by overpredicting and using a threshold of 0.5 to convert to a categorical forecast, a similar classification table is achieved to making accurate probabilistic forecasts but using a lower threshold for converting to categorical forecasts. This appears to be supported by the ROC plots (Figure 23, bottom), in which the largest H&KSS is obtained for a threshold value close to 0.5, as compared to most methods for which the threshold for the maximum H&KSS is much lower. The ROC curves also show a relatively rapid decrease in the probability of detection as the false alarm rate decreases, indicating that this method may not be well suited to issuing all-clear forecasts.

Zoom In Zoom Out Reset image size

Wheatland (2004) presented a Bayesian method to predict flaring probability for different flare sizes using only the flaring history of observed active regions. Since it uses only flare history as input, the results from this approach serve as a baseline for comparison with methods using magnetogram and/or white-light image data.

The event statistics method assumes that solar flares (the events) obey a power-law frequency-size distribution:

$$\begin{eqnarray}&&N(S)={\lambda }_{1}(\gamma -1){S}_{1}^{\gamma -1}{S}^{-\gamma },\end{eqnarray} \tag{ 13 }$$

where N(S) denotes the number of events per unit time and per unit size S, where λ₁ is the mean rate of events above size S₁, and where γ is the power-law index. The method also assumes that events occur randomly in time, on short timescales following a Poisson process with a constant mean rate, so that the event waiting-time distribution above size S is

$$\begin{eqnarray}&&P(\tau )=\lambda \exp (-\lambda \tau ),\end{eqnarray} \tag{ 14 }$$

where τ denotes the waiting time for events above size S, and where λ is the corresponding mean rate of events. Given a past history of events above a small size S₁, the method infers the current mean rate λ₁ of events subject to the Poisson assumption, and then uses the power-law distribution to infer probabilities for occurrence of larger events within a given time.

Wheatland (2005) presented an implementation of the event statistics method for whole-Sun prediction of GOES soft X-ray events. Figure 24 shows a GOES light curve, from which GOES event lists are routinely produced. In the application to GOES data the event size S is taken to be the peak 1–8 Å GOES flux in the GOES X-ray event lists, and the times of events are identified with the corresponding tabulated peak times. The whole-Sun method presented in Wheatland (2005) is used in this study, and the method is also applied (for the first time) to events in individual active regions.

Zoom In Zoom Out Reset image size

Three applications of the method were run for these tests: active region forecasts for which a minimum of five prior events was required for a prediction, active region forecasts for which ten prior events were required, and a full-disk prediction. Since none of the other methods produced a full-disk flare probability, those results are included for this study just for completeness. This method has no explicit restriction on the position of an active region on the disk, although the minimum number of events means that predictions are not made for regions that have just rotated into view. The sample sizes for the region forecasts with ten prior events are the smallest of any method, indicating that a substantial fraction of regions do not produce even ten events.

A summary of skill results is given in Table 17. The results for C1.0+, 24 hr are arguably the worst of any method, but the approach is biased against regions that produce only a few small events. The values of the skill scores generally increase when more prior events are included, so the five prior event case has most skill score values lower than the ten prior event case, which in turn has most skill score values lower than the full disk case. It is likely that the larger number of prior events simply allows for a better estimate of the power-law index and the mean rate of events above a given size. The accuracy of the forecast is expected to scale as ${N}_{1}^{-1/2}$, where N₁ is the number of small events observed and used to infer the rate λ₁ (Wheatland 2004).

The summary plots (Figure 25) are shown only for the region forecasts requiring ten prior events. The reliability plots show a systematic overprediction for all event definitions, and the ROC plots do not noticeably improve with increasing threshold, in contrast to most other methods. The peak H&KSS score is somewhat low and remarkably consistent across event definitions.

Zoom In Zoom Out Reset image size

Another method that was applied to the active region patch data set which did not use the magnetogram data was one based on historical flare rates from McIntosh active region classifications. This method is the same as that presented in Bloomfield et al. (2012), where the occurrence of GOES X-ray flares from individual McIntosh classifications were collated over 1969–1976 and 1988–1996. These average 24 hr flaring rates, μ₂₄, lead to a Poisson probability of one or more flares occurring in any 24 hr interval from Gallagher et al. (2002),

$$\begin{eqnarray}&&{P}_{{\mu }_{24}}(N\geqslant 1)=1-\exp (-{\mu }_{24}),\end{eqnarray} \tag{ 15 }$$

or the probability of one or more flares in any 12 hr interval from,

$$\begin{eqnarray}&&{P}_{{\mu }_{12}}(N\geqslant 1)=1-\exp (-{\mu }_{24}/2).\end{eqnarray} \tag{ 16 }$$

The nature of the statistics collated in Bloomfield et al. (2012) places a limitation on the forecasts outlined in Section 2.2 that are able to be studied by this method. Forecasts of at least one C1.0 or greater flare within 24 hr (C1.0+, 24 hr) are directly achieved by the Poisson flare probabilities for "Above GOES C1.0" published in Table 2 of Bloomfield et al. (2012), while forecasts of at least one M1.0 or greater flare within 12 hr (M1.0+, 12 hr) were calculated by combining the 24 hr M- and X- class flare rates published in the same Table (μ₂₄) and applying Equation (16) above. Forecasts of at least one M5.0 flare or greater within 12 hr (M5.0+, 12 hr) were not capable of being achieved by this method because the flares that were collated from 1969–1976 (Kildahl 1980) were identified only by their GOES class and not the complete class and magnitude.

For each magnetogram patch the observation date and NOAA region number(s) contained within it were used to cross-reference the McIntosh class of that active region in that day's NOAA Solar Region Summary file. It should be noted that, even for the C1.0+, 24 hr and M1.0+, 12 hr forecasts, flare probabilities were not able to be issued for some patches because the magnetogram was recorded before the active region had received a NOAA number designation. For the case of magnetogram patches containing multiple active regions, the reported flare probability was the largest probability from any of the regions within that patch.

The summary plots (Figure 26) for the C1.0+, 24 hr and M1.0+, 12 hr thresholds show an overprediction tendency in the reliability plots, with marginal improvement with increased threshold in the ROC plots. The Poisson method as implemented here is an example of a method whose requirements are not well met by the data used for this workshop, and the method is likely penalized as a result, as evidenced by the relatively low skill score values in Table 18.

Zoom In Zoom Out Reset image size