In this study, we used two kinds of data from May 2010 to October 2022, including the 3-hour Kp index from the National Oceanic and Atmospheric Administration (NOAA) and the 193Å wavelength images measured at the SDO/AIA (O’Dwyer et al., 2010; Lemen et al., 2012; Pesnell et al., 2012). They were divided into the training subset (including data from May 2010 to December 2019) and the testing subset (including data from January 2020 to October 2022). A prediction model was trained by the training subset and tested by the independent testing subset.

It should be mentioned that we only focus on the background solar wind (including the coronal hole high-speed streams and co-rotating interaction regions) in this study. Therefore, we removed the time periods when the geomagnetic conditions are affected by interplanetary coronal mass ejections (ICMEs), according to the near-Earth ICMEs list (Cane and Richardson, 2003; Richardson and Cane, 2010). For each day, we selected two samples according to the timepoints at 0:00 UTC and 12:00 UTC. Thus, we obtained 6,042 samples from the training subset and 1853 samples from the testing subset.

As shown in the flow chart in Figure 1, we prepared the 54-day Kp index time series before the current timepoint (T), and the AIA 193Å image was measured at T as inputs and we took the 3-day Kp index time series after T as outputs. Figure 2 shows how to divide the Kp time series into the previous 54-day Kp input and the next 3-day Kp output for the current T sample. Take the sample at 00:00 UTC on 22 Oct 2019 as an example, we took the Kp time series (432 points) from 00:00 UTC on August 29 to 00:00 UTC on October 22 as the input and the Kp time series (24 points) from 00:00 UTC on October 22 to 00:00 UTC on October 25 as the output. Figure 3 shows how to prepare the AIA 193Å image measured at T as input for the current T sample. Take the sample at 00:00 UTC on 22 Oct 2019 as an example, the observed image is shown in the top left panel. Then, a slice covering the [S40, N40] and [E40, W40] areas was cut from the observed image and shown in the top right panel. A median value was calculated based on the previous monthly slices. After all values above the median value in the slice at current T were replaced by 0, we obtained a slice as shown in the bottom left panel referring to the properties of coronal holes. The slice in the bottom left panel was then resized into a smaller size of 32 × 32 pixels and used as one of the inputs.

To compare the difference between the training and testing subsets, we calculated the maximum value of the 3-day Kp output and took it as a representative for each sample. Then, we draw the sample distribution histogram in Figure 4. The blue and red bars represent the training and testing subsets, respectively. It was found that the maximum Kp of the samples in the two subsets had similar distribution properties. There were 1,614 (26.7%) and 390 (21.0%) storm samples (maximum Kp ≥ 5) in the training and testing subsets, respectively.

Two samples at timepoints T _i and T _j from the training subset were called a pair (where the timepoint T _i was more than 3 days farther from T _j ). For each sample at timepoint T _i , we picked up more than 5,000 other samples from the training subset to form sample pairs with it. All the samples that were within a 2-month interval of the current timepoint were excluded. As a result, we obtained more than 25 million sample pairs from the training subset.

We assume that if the inputs of two samples in a pair have high similarity (that is, their previous and current geomagnetic conditions are similar), their geomagnetic conditions in the next 3 days should be similar too. In this case, the two samples in the pair are similar to each other so that one sample can be a proper representative of the other one. In this study, the similarity parameters of the pair inputs were calculated and used to develop a model for the 3-day Kp forecast based on this assumption.

Figure 1 shows that based on the 54-day Kp time-series inputs, we derived 11 similarity parameters for each pair. The first two parameters were the mean absolute error (MAE _Kp ) and the root mean square error (RMSE _Kp ). They were calculated by the following formulas, where x and y represent the inputs of a pair at timepoints T _i and T _j , respectively: M A E x , y = 1 N ∑ i = 1 N x i − y i , (1) R M S E x , y = 1 N ∑ i = 1 N x i − y i 2 . (2)

The next two parameters were the maximum absolute error (MaxDiff _Kp ) and the sum of the absolute error for the storm (Kp ≥ 5) period (SDiff _Kp ). Also, four spatial distances were calculated, namely, the Euclidean distance (D (Euclidean) _Kp ), cosine distance (D (Cosine) _Kp ), correlation (D (Correlation) _Kp ), and Hamming distance (D (Hamming) _Kp ), by the following formulas and were taken as similarity parameters. The SciPy package ⁶ was used to calculate these spatial distances. D E u c l i d e a n = ∑ i = 1 N x i − y i 2 (3) D C o s i n e = 1 − x ⋅ y x 2 y 2 , (4) D C o r r e l a t i o n = 1 − x − x ̄ ⋅ y − y ̄ x − x ̄ 2 y − y ̄ 2 . (5)

We also derived three features (Feature _KPCA(linear) , Feature _KPCA(rbf) , Feature _KPCA(cosine) ) from each Kp time-series input using the kernel principle component analysis (KPCA) algorithm with the “linear,” “rbf,” and “cosine” kernel using the Scikit-learn package ⁷ . Through the KPCA algorithm, a Kp time-series input can be represented by a number (the feature). So, for each pair, we obtained two numbers (Feature _{x, KPCA(kernel)} and Feature _{y, KPCA(kernel)}) representing the Kp time-series pair inputs (x and y). Then, we calculated the absolute error of the two numbers and took it as the similarity parameter for the pair. In this way, we obtain three similarity parameters, namely, F (linear) _Kp , F (rbf) _Kp , and F (cosine) _Kp , by the following formula: F k e r n e l = F e a t u r e x , K P C A k e r n e l − F e a t u r e y , K P C A k e r n e l , k e r n e l ∈ l i n e a r , r g b , c o s i n e . (6)

Figure 1 shows that based on the coronal holes image inputs, we derived another seven similarity parameters for each pair. The first two parameters were the mean absolute error (MAE _CH ) and root mean square error (RMSE _CH ), and were calculated by Formulas 1 and 2. The next two were the widely used image similarity parameters, such as structural index similarity (SSIM _CH ) and peak signal-to-noise ratio (PSNR _CH ), which were calculated using the Scikit-image package ⁸ . Then, in a similar way, we calculated the absolute error of the KPCA features derived by the image pair inputs using Formula 6 and took them as the similarity parameters, such as F (linear) _CH , F (rbf) _CH , and F (cosine) _CH .

Then, we obtained 18 similarity parameters from the pair inputs in the training subsets. Except for the correlation (D (Correlation) _Kp ), structural index similarity (SSIM _CH ), and peak signal-to-noise ratio (PSNR _CH ), the larger the other 15 parameters, the higher similarity the pair inputs reach. Hence, it is the opposite for the correlation (D (Correlation) _Kp ), structural index similarity (SSIM _CH ), and peak signal-to-noise ratio (PSNR _CH ). Then, we standardized the 18 similarity parameters independently by computing the relevant statistics on the pairs in the training subsets using the Scikit-learn package.

As we established the pair inputs and calculated the similarity parameters, the pair outputs (3-day Kp time series) should be analyzed and labeled as a number (1 or 0). Here, we adopted a batch labeling method according to the following principles:

1) For each sample at timepoint T _i , if the next 3-day geomagnetic conditions reached the storm level (maximum value of Kp ≥ 5), it is a storm sample. Otherwise, it is a non-storm sample.

Finally, we prepared our dataset into a standard form consisting of inputs (18 similarity parameters) and output (binary labels) that is suitable for a binary classification task.

Feature selection is a widely used approach in the machine learning community to select the best features from the dataset in which the model can perform better with less training time. In this study, we selected the best features by comparing the Pearson’s correlation of the 18 similarity parameters with the binary labels of the pairs using the SelectKBest from the Scikit-learn package. The correlations were standardized to scores ranging from 0 to 1, as shown in Figure 5.

It was found that there were three best features that had a significantly higher score than others. The top three best features were the Euclidean distance of the 54-day Kp time series (D (Euclidean) _Kp ), the Hamming distance of the 54-day Kp time series (D (Hamming) _Kp ), and the differences between the KPCA features of the current coronal hole images for each pair (F (linear) _CH ). They were selected from the dataset and used for model development, as shown in the flowchart in Figure 1.

After the three best features had been selected by feature selection, we determined to develop a classification model for pairs in the training subset to predict whether a pair of two samples is a good match or not. If the pair is a good match, their geomagnetic conditions are similar to each other so that one sample can be considered as a forecast of the other.

The scalable tree boosting system (XGBoost) is a widely used machine learning algorithm for binary classification tasks. It implements gradient boosting which performs additive optimization in functional space and incorporates a regularized model to prevent over-fitting (Chen and Guestrin, 2016).

We applied the XGBoost algorithm to develop a classification model using the Scikit-learn package ⁹ to predict whether a pair of two samples is a good match or not. If the pair is a good match, their geomagnetic conditions are similar to each other so that one sample can be considered as a forecast of the other. Considering the “f1-weighted” (it is a balanced score of a standard binary classification task) as metrics, the best hyper-parameters of the XGBoost model are n_estimators = 0.01 and max _depth = 20.

As shown in Figure 1, to predict the 3-day Kp time series at T _i , we give a forecast of the 3-day Kp time series following the steps:

1) We pick up a sample at T _j to establish a pair and calculate the similarity parameters of the pair.

We developed a prediction model using the training subset and applied the model on the testing subset from January 2020 to October 2022 and evaluated its performance. Both the geomagnetic storm prediction (maximum Kp ≥ 5) and non-storm (maximum Kp < 5) prediction are important in space weather forecasts. Thus, we will evaluate the model for the geomagnetic storm prediction (maximum Kp ≥ 5) as well as the non-storm (maximum Kp < 5) prediction tasks.

For a binary classification task like the next 3-day geomagnetic storm (maximum Kp ≥ 5) prediction, the confusion matrix is shown in Table 1. The model is evaluated by its ability to predict the next 3-day geomagnetic storm, the first day (day 1) geomagnetic storm, the second day (day 2) geomagnetic storm, and the third day (day 3) geomagnetic storm. We also evaluate the model’s ability to predict the next 3-day non-storm conditions, the first day (day 1) of non-storm conditions, the second day (day 2) of non-storm conditions, and the third day (day 3) of non-storm conditions.

There are three metrics for binary classification used in the study, including precision, recall, and F1-score. They are calculated by the following formulas: r e c a l l = T P T P + F N , (7) p r e c i s i o n = T P T P + F P , (8) F 1 = 2 × p r e c i s i o n × r e c a l l p r e c i s i o n + r e c a l l , (9) where the true positive (TP), false positive (FP), false negative (FN), and true negative (TN) are calculated following the confusion matrix, as shown in Table 1 for storm prediction and Table 2 for non-storm prediction. For example, for storm prediction, the higher the recall, the more storm samples have been correctly predicted. The higher the precision, the fewer false alarms have been made. An F1-score is a balance metric for recall and precision. A larger F1-score shows a better ability to classify a sample into the correct category.

The three metrics, recall, precision, and F1-score for the geomagnetic storm prediction task (considering the storm samples as positive samples) were listed as “For storm category (maximum Kp ≥ 5),” as shown in Table 3. The metrics for the non-storm prediction task (considering the non-storm samples as positive samples) were listed as “For non-storm category (maximum Kp < 5),” as shown in Table 3. We found that

1) Our model reaches an F1-score of 0.96, a recall of 0.99, and a precision of 0.93 for the next 3-day period of non-storm category prediction and an F1-score of 0.82, a recall of 0.70, and a precision of 0.98 for the next 3-day period of storm category prediction.

2) However, both for the non-storm category and storm category, the three metrics (F1-score, recall, and precision) by our model are lower at the first day, second day, and third day period prediction than that in the next 3-day period prediction. It indicates that it is difficult to accurately predict the storm periods in a shorter time period (less than 3 days).

We also conducted an error analysis of the model forecasts with the observations in the testing dataset. For each sample (containing 24 points of the 3-hour Kp values) in the testing subset, we compared the forecasts with the observations by four statistics, including the mean absolute error (MAE), the root mean square error (RMSE), the maximum of the absolute error (MaxDiff), and the sum of the absolute error for the storm (Kp ≥ 5) period (SDiff). For 1853 samples in the testing subset, we calculated the average and standard deviation of the four statistics by the following formulas: A v e r a g e x = 1 N ∑ i = 1 N x i , (10) S t d D e v x = 1 N − 1 ∑ i = 1 N x i − x ̄ . (11)

The eight statistics of the 3-day period, namely, day 1, day 2, and day 3 predictions, for all samples in the testing subset by our model are shown in Table 3. It was found that

1) The average and standard deviation of the mean error for the next 3-day Kp prediction is 0.03 and 0.65, respectively. The average and standard deviation of the mean absolute error for the next 3-day Kp prediction is 1.06 and 0.32, respectively.

Moreover, we compared our model results with the daily product of the 3-day Kp index forecasts provided by SWPC from 30 Nov 2020 to 27 Oct 2022. After we had removed the time-periods when the geomagnetic conditions are affected by interplanetary coronal mass ejections (ICMEs), there were 50 storm samples and 548 non-storm samples. Evaluation metrics for SWPC’s 3-day Kp index forecast products and results obtained by our model are shown in Table 4. It was found that

1) For the next 3-day prediction, our model provides a positive mean error average, while SWPC provides a negative mean error average. It indicates that statistically, our results are usually higher than the observations and SWPC’s products are usually lower than the observations.