Solar flares are one of the solar magnetic phenomena that suddenly bursts out large amounts of electromagnetic energy from the solar atmosphere. When they combine with other solar activities like coronal mass ejections (CMEs), they may potentially provoke devastating geomagnetic storms, which can cause a massive scale of electrical blackouts, disrupt power grids and radio communications, and even damage sensitive electronic devices in space. Flare phenomena usually occur when the intense magnetic fields of the Sun tangle up and become unstable. The relationship between the photospheric magnetic field and flaring (or eruptive) phenomena still remains questioned. Therefore, researchers have put efforts on predicting the flare event prior to its occurrence to avoid crucial costs and destruction.
Over the past two decades, many machine learning-based algorithms attempted to remedy the problem [1] by issuing binary predictions on major flares (i.e., M-and X-class) instead of considering where the relationship could be derived from. However, most of the predictions are generated by applying point-in-time values within the maps of magnetic field strength (i.e., magnetograms) meaning that only one individual value in multiple physical parameters observations is representative for each solar flare event [2] [3]. Such approaches could be limited in predicting power because it does not take the temporal evolution aspect of these parameters into account. Within the recent proposed classification algorithms, Deep Neural Networks (DNNs) have achieved remarkable performance in the field of computer vision [4], image recognition [5], and natural language processing (NLP) [6] [7] [8]. As one of the modern machine learning architectures, this approach has the ability to automatically transform input data into a relatively more abstract representation and use the composite representation to perform feature extraction or classification during the learning process. This benefits especially for prediction in high dimensional space since it requires less human cognitive process.
Many researchers focus on improving the deep learning model by improving its generalization ability, which is the performance difference of a model between training and testing phase. This means that models with poor generalization ability are easier to cause overfitting than others. To understand or predict solar flares better and issue accurate predictions, there is a need for more data since the high performance of deep learning architecture relies heavily on data volume. However, one major issue for data-driven solar flare forecasting is the limited amounts of high quality data, especially with the flaring instances. Within the most recent records, less than 20% of the active regions produce large flaring (e.g., ≥M) instances [9], which lead to insufficient features derived from these instances for training and validating complex prediction models.
Data augmentation is a way to synthetically modify the raw data for the models to recognize generalizable features and hidden patterns and therefore improve the overall performance. The techniques implemented for data augmentation aim to improve the quantity of the datasets so that better models can be trained. With augmented data, models can cover those unexplored input spaces, and therefore improve the generalization ability. Inspired by this idea, our work is further extended to seek the feasibility of applying data augmentation techniques to optimize the performance of different deep learning models. Four most common time series data augmentation techniques are selected and five deep learning architectures are refined to monitor the performance compared to the baseline model that is trained without applying any augmentation techniques. By this, we conducted a multivariate time series classification by utilizing the time series data from a well-known solar flare prediction dataset [10] and synthetically modified the time series instances with various data augmentation techniques.
The rest of the paper is organized as follows. In Section II, we discuss some related work on well-known flare prediction with a focus on deep learning and utilization of data augmentation. In Section III, we explain our research methodology on data collection and preparation as well as our training and evaluation processes. In Section IV, we present our experimental evaluation results for our models. In Section V and VI, we provide a summary of our findings and discuss some future work avenues.
In the field of solar flare forecasting, there are various studies on training deep learning models on magnetic field data products, either the original magnetic field rasters or the derived metadata products. One major source for such prediction is the Michelson Doppler Imager (MDI) observation program from Solar and Heliospheric Observatory (SOHO) [11], which has been superseded in 2011 for the Helioseisemic and Magnetic Imager (HMI) onboard the Solar Dynamics Observatory (SDO) [12]. As for working with the magnetic field rasters, the study proposed by [13] transformed the Xray flux time series data from GOES at one-minute cadence into 64×64 Markov Transition Field images and provided the Convolutional Neural Network (CNNs) with these image data to predict flare occurrences. FlareNet [14], a framework designed for experimentation with solar flare prediction, utilizes high dynamic image data directly from SDO where each pixel contains a high dynamic range of flux that can potentially affect gradient updates and encourage overfitting. The problem is addressed by treating each timestamp of temporally adjacent images captured in the same spatial locations as additional image channels and appending side channel information vectors to the fully connected layers to avoid memorization. Sumi et al. [15] also uses soft X-ray data with longer wavelength at one-minute cadences from Geostationary Operational Environmental Satellite (GOES) to predict Xray flux in continuous target space and evaluates with three DNN architectures. In their research, the hierarchical dense residual network has clearly better performance compared to the baselines and other DNN architectures.
Deep Flare Net-Reliable (DeFN-R) [16], a probabilistic flare forecasting model composed of multilayer perceptrons, uses 79 extracted features from each active region (i.e, AR) observed by the multiwavelength emissions and annotated with flare occurrence labels (i.e, X-, M-, and C-classes). This model forecasts the maximum flare classes of the next 24 hours as well as the occurrence probability of each event. LSTM DNN Flare Net (LDFN), in [17], also uses a derived dataset where its experimental dataset contains multivariate time series (MVTS) of magnetic field metadata. This approach combines Long Short-Term Memory (LSTM) [18] architecture with multi-scale skip connections. The overall interpretability of this architecture is improved by introducing human cognitive processes based on physical definitions that were extracted and combined with similar ones using deep learning models. Hui et al. [2] proposed an ordinal logistic regression model by using three predictive parameters derived from the MDI magnetograms, with categories of the maximum magnitude of flares, to predict the probability of a given AR in the next twenty-fours hours. As in [19], Huang et al. generates a meta dataset that obtained multiple instruments (SDO/HMI and SOHO/MDI) magnetic field measurements and evaluates with a four layered DNN.
As solar flares are relatively less frequent events, the scarcity of the major solar flares often result in an extreme classimbalance ratio where the number of flaring samples contained (usually the intensity levels above M-class in NOAA/GOES flare classification) are much higher than the number of nonflaring samples contained (with intensity levels of C, B, and flare-quiet regions). Even with the available data, there still exists the problem of limited high quality data instances for models like deep neural networks, which relies heavily on large scale data. To alleviate these two problems, oversampling and undersampling are the typical approaches to remedy the class imbalance issue [9] [20] [21]. However, a big portion of data is left out when performing undersampling, and a large amount of data is replicated when oversampling is applied [21]. Data augmentation, on the other hand, tackles this classimbalance issue by modifying existing instances to generate potentially more generalizable features. In this way, more data is introduced to the models and the generalization capabilities of these classifiers can be improved by expanding the model's decision boundaries [22]. Many original proposals of the CNN architectures take advantage of data augmentation. AlexNet [23], for example, archived a benchmark record on the Ima-geNet Large Scale Visual Recognition Challenge (ILSVRC) dataset by using cropping, mirroring, and color augmentation. Other network architectures such as DenseNet [24], Visual Geometry Group network [25], or Residual Networks [26] also uses augmentation.
Our work implements various data augmentation techniques for solar flare prediction via five deep learning models. This not only provide a multivariate schema for such deep learning models but also intend to improve their performance using augmented data instances. We hope this work affords researchers and practitioners in related fields more exposure to multivariate time series analysis for space weather forecasting.
As solar flares are mainly electromagnetic phenomena, they are usually predicted by employing the adequate magnetic field parameters. Inspired by Bobra and Couvidat that applied a large set of SDO/HMI Vector magnetic field data for predicting M-and X-class flares [27], Space Weather Analytics for Solar Flares (SWAN-SF) [10] is a comprehensive multivariate time series dataset consists of physical parameters derived from solar photospheric magnetic fields in the Space weather HMI Active Region Patches (SHARPs) data series.
There are five partitions from the SWAN-SF dataset, where every partition contains a similar number of instances labeled with X-and M-class flares. All data in the partitions are time segmented, which means that all data points are temporally non-overlapping time series slices and none of the active regions cross in time or space. In this study, we used Partition 1, 2, 3, and 5 for training purposes and excluded Partition 4 for evaluating the model performance since this partition has a relatively more representative class imbalance ratio compared to the overall dataset. Every reported flare in each partition is associated with a designated active region identifier and related time series data of the magnetic field parameters extracted from this active region. As flare intensity can be classified logarithmically by their peak X-ray flux, we assigned level X and M flare as the flaring classes while level C, B, and flare-quiet regions (labeled as N) as the non-flaring classes. Additionally, five parameters employed in this study are USFLUX (Total unsigned flux), TOTUSJH (Total unsigned current helicity), ABSNJZH (Absolute value of the net current helicity), SAVNCPP (Sum of the modulus of the net current per polarity) and TOTBSQ (Total magnitude of Lorentz force), which are shown to be most representative of the solar activity in [27]. This dataset will be utilized as our benchmark dataset throughout this study for consistency between our results.
In order to be compatible with the five selected physical parameters, we choose the commonly used machine learning classifiers -Support Vector Machine (i.e., SVM) [28] and five well-known deep learning classifiers that are able to classify multivariate time series instances.
We adopt a SVM classifier, which makes use of summary statistical parameter values derived from our benchmark SWAN-SF dataset. From each original physical parameter employed above, a total of 43 summary statistics data is generated based on every individual time series. This results in a frame of data with the same amount of flaring versus non-flaring instances in each partition.
One Direction CNN (ODCNN) is designed as one of the most common CNNs where convolutional layers are applied on a layer-by-layer architecture to the given time series [29]. In this framework, we create the model based on a sequential classifier meaning that the layers in this classifier are linearly stacked and do not communicate with layers other than the previous and next one. Therefore, the model flexibility is limited as it only focuses from the outputs of previous layer. Each layer contains weights that passed to the posterior layer and is followed by an activation function that allows the model to take nonlinear relationships into account. We set up most standard 5 filters with length of 3 and a Rectified Linear Unit (ReLU) as the activation function [30].
Encoder for time series classification, originally proposed in [31], is a Fully Connected Convolutional Neural Networks (FCNs) that uses dense connections between layers. We generate 128, 256, and 512 filters in the first three convolutional layers with respective lengths of 5, 11, and 21. Each of the convolution layers is operated by an instance normalization [32] and the result from such normalization is then fed into the Parametric Rectified Linear Unit (PReLU) [33] activation function. Furthermore, a dropout operation is applied after each activation function as well as a max pooling with length of 2. The main difference of such architecture from a traditional FCN is replacing the last pooling layer with an attention mechanism that allows the network to learn which parts of the time series are more important for the classification task. This mechanism is implemented by multiplying the input time series with another time series (that has gone through a softmax function). It is worth mentioning that both of the two time series multiplied have the same length, where only the second time series is acting as the weight for the first one. In this way, the network can learn the importance of each timestamp in its corresponding weights.
Residual Network (ResNet) is another architecture originally proposed by [34], which builds shortcut residual connection between consecutive convolutional layers. This architecture reduces the vanishing gradient effect [35] by infusing a linear shortcut between the input and output of each residual block so that the gradient can flow through these shortcuts and make the training smoother. In our case, three residual blocks are used where each residual block contains three convolutions followed by a ReLU activation function. Each of our convolutions consists of 64 filters and the filter lengths are set to 8, 5, and 3, respectively. The output from all three convolutions will be added back to the input of each block and combined as a whole for the final fully-connected layer.
Multi Channel Deep Convolutional Neural Network (MCCNN) was proposed by [36] where the convolutions of the model are parallelized and applied individually to each channel (i.e., dimension) of the input time series. Each channel of the input series will pass two convolutional stages of 8 filters. Each filter's length is 5 and followed by a ReLU activation function as well as a max pooling operation with length of 2. The output from these channels will be concatenated together and fed into a fully-connected layer of 732 neurons. At last, a softmax classifier is generated where the number of neurons matches the number of classes.
Inception [37] is designed to tackle the problem of stacking large numbers of convolutional layers by implementing multiple kernels with different sizes that operate on the same stage. In our case, we set the kernel sizes to 3, 5, 8, 11, 17 and generate a module for saving these kernels. Each of these different convolution layers works differently. For example, convolution size of 1x1 allows the network to learn the patterns across the input series while convolution size of 3x3 (and 5x5) allows the network to learn spatial patterns across all dimensions of the input including height, width, and depth. Similar to the MCCNN mentioned above, this network also ends with a max pooling operation as well as a ReLU activation function.
In the second part of this research, we explore the feasibility of using time series augmentation techniques to improve the performance of deep learning-based multivariate time series classifiers. We will use four well-known augmentation techniques [38] described as follows.
Jittering is a type of noise injection, in which a matrix of random values are drawn from a Gaussian distribution and added onto the original data. A general form for using this technique is provided below:
where ϵ is the Gaussian noise added to the time series with length of T. In this case, the standard deviation α of the selected Gaussian distribution is a hyperparameter that needs to be determined before training the model. This method increases the generalization of the networks [39] [40] by assuming that the testing data are similar to the training data but only with a difference of a factor of noise. It has been mostly used for testing the robustness of a given model against noise and ideally improve model performance.
On the other end of the spectrum, we also used smoothing, which works by applying a weighted average to the original data. This method often reduces noise and can benefit more for models that are strongly impacted by noisy data points.
Scaling is a method based on enlarging (or shortening) the global magnitude with a random scalar. The equation of this method is defined as:
where α is the random value drawn from the Gaussian distribution. Similar to jittering, scaling also introduces noise to the data. However, instead of adding individual noise directly on each data point, scaling multiplies the entire time series with a scaling parameter α (which is also drawn from a Gaussian distribution). In such cases, the scaling parameter α is the hyperparameter of the model.
Magnitude Warping is another method that aims at changing the magnitude of each time series. The formulation of this method is defined as:
where α 1 to α T is a series interpolated by a cubic spline S(u).
The only difference of such method is that each time series magnitude is warped by a smooth curve in range of 0 to 1. This is based on the idea that increasing or decreasing random regions in the time series can add small fluctuations into the dataset. Moreover, instead of many other transformationbased methods that use only one hyperparameter, this method assumes random transformation is practical and uses two predefined hyperparameters instead.
IV. EXPERIMENTAL EVALUATION As for the model evaluation, we implemented multiple metrics as well as some essential skill score measurements. In these measurements, TP (true positives), TN (true negatives), FP (false positives), and FN (false negatives) are used in standard settings. Positive stands for occurrence of a large flare (i.e., X-and M-class) while negative stands for comparatively smaller flare and flare-quite region (i.e., C-, B-class, and flare-quite regions). Note here TP represents when the model predicts correctly with the positive class while TN is where the model predicts correctly with the negative class. In space weather forecasting, some classical measures like accuracy and precision should not be considered as a key score to be relied upon because such measure scores can be significantly impacted by the imbalanced class distribution. In the case of substantially larger amount of non-flaring classes than the flaring classes, a classifier can always predict a flare will not occur and still reach a high accuracy score. Similarly, by introducing more negative samples to the dataset while maintaining the same amount of positive samples, the number of FP (false positives) may potentially increase, which will result in a lower precision score.
Two main skill score measurements we used for model evaluation in the experiments are True Skill Statistic score (TSS) and Heidke Skill Score (HSS). True Skill Statistic score (shown in Eq. 4) compares the difference between the probability of detection (Recall) and the probability of false detection (the ratio of inaccurate predictions of nonflaring/negative class over all the actual negative class).
Heidke Skill Score (shown in Eq. 5) measures the improvement of the forecast over a class-imbalance-aware random forecast where P = T P + F N and N = F P + T N . HSS values range from -1 to +1 where the value of 1 indicates the ideal forecasting and the value of -1 indicates an all incorrect forecasting. The closer the value to 0, the model has less power to distinguish between labels, meaning it has no skill over a class-imbalance-aware random prediction.
Three experiments are designed in this study discovering the performance of time series classifiers as well as how such deep learning models are affected after applying standard and synthetic augmented data instances.
In the first experiment, we evaluated our baseline SVM model (trained with the summary statistical data) and the other five selected time series deep learning models (trained with time series data).
From the results shown in Fig 1, models that are trained with deep neural networks have better performance than the baseline SVM model. It demonstrates that TSS increased from 12% to an average of 80% for ODCNN, Encoder, and MCCNN. Although we do not see a significant improvements in HSS, there is still a up to 5% improvement in Encoder, MCCNN, and Inception models. This strongly suggests that the overall performance of the baseline can be improved using deep learning techniques.
As for the second experiment, we modified the training samples with each augmentation techniques and trained the five aforementioned deep learning models with these generated samples. This is designed to evaluate the augmentation techniques itself with the original data instances. The total amount of instances are obtained the same as original. Then, we compared these results within our previous models, which are trained without infusing any augmentation instances.
As presented in Fig 2, models that are trained with augmented data perform generally better than the original. However, Inception, ResNet, and Encoder do not have a strong ability to distinguish between flaring and non-flaring classes when they are exposed to additional noise (in the case of jittering experiments) in the dataset. As for comparison, smoothing works the opposite way by reducing potentially noisy data points and achieve a 4% improvement in TSS. It shows a slightly better performance than the baseline as regards to these models. On the other hand, for models that are not affected much by noise (e.g., MCCNN and ODCNN), the augmentation method of jittering actually helps to improve the model performance more than smoothing. Besides jittering, the other two methods of scaling and magnitude warping have similar performance in the results for Encoder and MCCNN as both of these methods scale the range of the data.
In the third experiment, we generated synthetic augmented instances by stacking our original data with our modified data. The oversampled instances will be fed as training samples into models for training and then evaluated based on the measurements.
As the results shown in Fig. 3, the oversampled synthetic data has the ability to improve certain performance of the model (like ResNet). Both TSS and HSS scores for ResNet have a steady increase especially by using smoothing and scaling. However, the augmentation techniques also do not show much effect for models like Encoder, ODCNN, and MCCNN. The average improvement for these three models are less than 5% compared to our baseline model. It is worth noticing that the performance of Inception is not consistent. The model was boosted around 10% in both measuring scores with scaling and magnitude warping, but decreased about 5% with the use of smoothing.
From the experiments conducted above, we can observe that the overall performance of our deep learning classifiers have achieved better than the commonly used machine learning classifier. Even with the least improvement from ResNet model, there is still around a 40% increase in TSS. Such improvement might be due to the improvement on True Negative Rate where the deep learning models has better ability to specify the True Negatives correctly.
From the experiment of comparing original and synthetic modified data, we can observe that only applying jittering on the original data does not improve the performance for models that are strongly influenced by the noise. On the other hand, applying smoothing raises up the performance for these noise sensitive models. It can also be obeserved that by training these noise sensitive models (i.e., Encoder, ResNet, Inception) with synthetic oversampled dataset, jittering actually provides comparable results as the baseline models especially in HSS measurements.
In this work, we have observed the performance of utilizing deep learning models onto multivariate time series data for solar flare forecasting. We also explored the feasibility of applying various data augmentation techniques on these models. Our results revealed that the model performance slightly improves when applying with augmentation. However, individual augmentation technique coupled with deep learning models showed fluctuations, making some baseline models perform better than their data augmented counterparts. There is still a high number of false positives exists in our evaluation phase, which highly impacts our results of TSS and HSS scores. This is a common problem in flare forecasting, which poses itself as a tradeoff between false positives and false negatives [20]. The possibility of combining different techniques to improve the performances also remains. Our work can be extended by adding more parameters to the classifiers or implementing different models that are more robust.