ORIGINAL RESEARCH article
Multi-Variate LSTM Prediction of Alaska Magnetometer Chain Utilizing a Coupled Model Approach
- 1Department of Physics and Geophysical Institute, University of Alaska, Fairbanks, AK, United States
- 2NASA Goddard Space Flight Center, Greenbelt, MD, United States
- 3Department of Physics and Astronomy and Space Science Center, University of New Hampshire, Durham, NH, United States
- 4Department of Electrical and Computer Engineering, University of New Hampshire, Durham, NH, United States
- 5ASTRA, Boulder, CO, United States
- 6Department of Electronics and Electrical Engineering, University of Bath, Bath, United Kingdom
During periods of rapidly changing geomagnetic conditions electric fields form within the Earth’s surface and induce currents known as geomagnetically induced currents (GICs), which interact with unprotected electrical systems our society relies on. In this study, we train multi-variate Long-Short Term Memory neural networks to predict magnitude of north-south component of the geomagnetic field (|BN|) at multiple ground magnetometer stations across Alaska provided by the SuperMAG database with a future goal of predicting geomagnetic field disturbances. Each neural network is driven by solar wind and interplanetary magnetic field inputs from the NASA OMNI database spanning from 2000–2015 and is fine tuned for each station to maximize the effectiveness in predicting |BN|. The neural networks are then compared against multivariate linear regression models driven with the same inputs at each station using Heidke skill scores with thresholds at the 50, 75, 85, and 99 percentiles for |BN|. The neural network models show significant increases over the linear regression models for |BN| thresholds. We also calculate the Heidke skill scores for d|BN|/dt by deriving d|BN|/dt from |BN| predictions. However, neural network models do not show clear outperformance compared to the linear regression models. To retain the sign information and thus predict BN instead of |BN|, a secondary so-called polarity model is utilized. The polarity model is run in tandem with the neural networks predicting geomagnetic field in a coupled model approach and results in a high correlation between predicted and observed values for all stations. We find this model a promising starting point for a machine learned geomagnetic field model to be expanded upon through increased output time history and fast turnaround times.
Geomagnetically induced currents (GICs) are produced when the solar wind interacts with the Earth’s magnetic field, driving disturbances that map to the Earth’s surface (Oliveira and Ngwira, 2017). Electrically conductive materials, like the Earth’s crust, in the presence of these disturbances experience electric fields proportional to that of the changing geomagnetic field which are known as geomagnetically induced electric fields (Pirjola, 2000) and drive currents, GICs, as a response. These GICs, if strong enough, can disrupt and damage sensitive electrical devices on the ground that are not designed to handle these currents. GICs have been known to cause power outages, transformer damage, and pipeline corrosion on the ground which impacts our technology and fossil fuel dependent economy; such an event was the cause for a 9 h power grid blackout in Quebec, Canada on 13 March 1989, where strong GICs overloaded and damaged a transformer of the Hydro-Quebec electric company.
In response to the damage done by these events the science community has put focus on the prediction of GICs. While GICs can happen in any part of the globe, there is a higher occurrence of these events in higher magnetic latitude regions. Large geomagnetic field disturbances are often observed in these regions due to geomagnetic field lines at the surface connected to dynamic regions of the magnetosphere (e.g., polar cusps and the magnetotail). During times of high geomagnetic activity, strong geomagnetic field disturbances may propagate lower into middle magnetic latitudes, leading to a majority of GIC studies being focused on high and mid magnetic latitudes (Pirjola, 2005; Pulkkinen et al., 2005; Pirjol et al., 2007; Fiori et al., 2014; Blake et al., 2016). Some studies have also shown geomagnetic field disturbances in the low magnetic latitudes formed from oblique pressure shocks (Carter et al., 2015; Zhang et al., 2016; Oliveira et al., 2018). The widespread nature of geomagnetic field disturbances in response to fluctuating solar wind parameters can also be seen in Supplementary Movie S1 of the supplemental material, where a global response to a geomagnetic storm is observed, with the strongest field fluctuations located at high magnetic latitudes in the midnight sector, reaffirming the focus on high latitude GICs.
Various mechanisms are known to cause large geomagnetic field disturbances on the ground. The study by Carter et al. (2015) has shown interplanetary shocks being a viable creation mechanism for large geomagnetic field disturbances at high magnetic latitudes and the magnetic equator. Recent studies from Heyns et al. (2021) and Rogers et al. (2020) analyzed GICs as a function of geomagnetic pulsations, indicating that ULF waves can drive GICs for extended periods at high and mid latitudes. The studies by Rodger et al. (2017) and Dimmock et al. (2019) identified extreme geomagnetic storm activity and sudden geomagnetic storm commencement as drivers of GICs in the New Zealand and Fennoscandia regions. However, the studies by Ngwira et al. (2015) and Dimmock et al. (2020) have shown the timing of GICs in relation to geomagnetic storm activity can vary based on location. Local dB/dt is a function of ionospheric and magnetospheric currents, generally associated with geomagnetic storm activity and local conductivity gradients within the Earth’s crust.
To understand these localized peaks, high resolution physics-based models have been utilized to determine these fluctuations (Welling, 2019), however these models, while important in the progress of our understanding, are computationally expensive and time consuming, which make them inefficient for real-time predictions of GICs. The need for efficient and computationally inexpensive models has led to the utilization of machine learned neural networks, such as the ones done by Wintoft et al. (2015), Lotz and Cilliers (2015), and Keesee et al. (2020), to capture these disturbances from upstream drivers, mainly the solar wind parameters.
Machine learned algorithms are efficient due to their fast computation times after training and range in complexity defined by the user, allowing for versatile solutions. The bulk of computational requirements needed for most neural networks are during training, with the finished models being significantly lightweight at runtime. The lightweight models this study aims to produce are focused on the Alaska region, which is a high magnetic latitude area susceptible to pipeline corrosion from GICs (Gummow and Eng, 2002; Pirjola et al., 2003; Khanal et al., 2019; Liu et al., 2019). The models produced for this region make use of real-time magnetometer data that are available locally, allowing for an enhancement in model performance.
The first step in forecasting GICs in Alaska starts with a geomagnetic field prediction model utilizing 16 years of SuperMAG geomagnetic field observations across Alaska and NASA OMNIweb solar wind and interplanetary magnetic field (IMF) conditions. The data is used as input to three different model types: a multi-variate LSTM model, a multi-variate linear regression (MLR) model, and a coupled set of LSTM models. In the following section the data sources and different model types will be explained. Sections 3, 4 will cover the model results and a discussion on the effectiveness of these models, their shortcomings, and routes of improvement, followed by a summary of the work presented.
Data and Models
SuperMAG and OMNIweb Data
The study presented utilizes the SuperMAG magnetometer database and the NASA OMNIweb solar wind database from 01/01/2000 to 12/31/2015. The data selected from OMNIweb database provides solar wind plasma and interplanetary magnetic field (IMF) propagated from solar wind monitors at Lagrangian point 1 to a subsolar bow shock location. The data provided from OMNIweb are solar wind density, flow speed, dynamic pressure, temperature, IMF magnitude, and IMF BZ in geocentric solar magnetic (GSM) coordinates with a 1 min resolution. The parameters were selected from a combination of known indicators of geomagnetic storms and substorms utilized in similar studies such as Lotz and Cilliers (2015) and Keesee et al. (2020). The use of derived parameters (i.e., parameters calculated from other data values, such as dynamic pressure) was limited to avoid data redundancy within the neural network which can result in poor performance within these models. Due to a non-continuous data coverage, a linear interpolation was applied to fill in gaps of 10 min or less, increasing the coverage from 70 to 80% during the 16 years segment.
The SuperMAG database hosted by Johns Hopkins University (Gjerloev, 2012) collects data from world wide magnetometer stations and provides baseline removed geomagnetic field data with a consistent coordinate system. From SuperMAG we selected four Alaska stations at geographic coordinates: Fort Yukon (FYU) at 66.56° N 214.87° E, College (CMO) at 64.87° N 212.15° E, Poker Flat (PKR) at 65.12° N 212.57° E, and Kaktovik (KAV) at 70.14° N 216.35° E; and utilized the north-south component of the observed field with 1 min resolution. Geomagnetic field information from SuperMAG is in their local magnetic coordinate system where the Z-component is kept static and the other two components are rotated to maximize North-South and minimize East-West with respect to a slowly varying declination angle. The magnetic field datasets utilized have had standard yearly and daily baseline removal as detailed by Gjerloev (Gjerloev, 2012). The data from CMO was used as an initial testing set to determine the best performing model configuration, and then the configuration was re-applied to the other three stations. The stations were chosen for their locations, which are roughly on a magnetic meridional line, making them perpendicular to the typical auroral oval in Alaska. This property of the stations chosen makes them suitable for ionospheric current predictions above Alaska.
A basic recurrent neural network (RNN) takes incoming data, computes the output, and sends the output as an input parameter to be used with the next incoming dataset (Brownlee, 2017). This means that each output is directly reliant on the previous known output and the newest data, making RNN able to predict time dependent sequences where the features act on a small time scale. The downside to RNN is that features acting on a long time scale are not properly accounted for, since the newest information is always the most relevant in the prediction process. These neural networks also rely on continuous datasets, treating a known gap in data as a standard time step in progression. Further, these networks are prone to vanishing or exploding error gradients during trainings which can result in excessive computation for minimal gains in performance.
For datasets where the input features have long or variable time scale implications RNN are not suitable, as the network will forget the information. To combat this, advanced RNNs, such as LSTM, are developed and utilized to retain information on a longer time scale (Hochreiter and Schmidhuber, 1997). LSTM achieves this through the use of two internal features known as gates and the cell state. Within the LSTM kernel three gates are utilized, these gates determine which information is added and removed from the cell state, and determine which parts of the incoming data are relevant to the current prediction. Unlike the standard RNN, where the output prediction is fed to the next iteration, the cell state is passed to the next LSTM kernel, and is used in unison with the incoming dataset to make a prediction. Due to the activation of the gates, a cell state may not be updated every sequence, allowing it to retain information from previous datasets. This model type has been successfully applied for predicting magnetometer data (Keesee et al., 2020).
Linear regression is a widely used approach for many empirical models (Verbeek, 2017). This study developed a MLR model for each station for comparison against the state-of-the-art machine learning models. The MLR models developed utilize the same inputs and outputs as the LSTM and are applied with the same 16 years of SuperMAG and OMNIweb data to fit the following equation:
where xi (t − n) is an input variable from n minutes prior, C is an offset variable, and αin is a scaling factor for variable xi (t − n). This model type is applied to each station and compared to their respective LSTM models. We chose MLR for comparison due to its ease of implementation and wide use for empirical studies.
Model Development and Results
Geomagnetic Field Prediction Model Using CMO Dataset
The initial phase of models produced was aimed at determining the length of time history of each feature to provide to the LSTM model. Four models were trained off a smaller 2009–2014 dataset utilizing 1-, 30-, 60-, and 120- min time histories as input for predicting |BN| at the next minute and tested for performance against the 2015 test set. Each model takes the same basic inputs of IMF magnitude, its Z-component in GSM coordinates, solar wind density, speed, flow (ram) pressure, temperature, and magnetic local time (MLT) of the station. The MLT has been converted into sin (MLT) and cos (MLT) to maintain a cyclic dependence of this variable, which is important for preserving nighttime and daytime dependencies and transitions. While the LSTM kernel is adept at remembering via the implemented cell state, the previous 1-, 30-, 60-, and 120 min of each variable were used as input in a 2D array of shape [m x n], where m is the number of features and n is the amount of time history, as seen in Figure 1. For consistency, each model utilized a single LSTM layer with 32 hidden neurons, a rectified linear unit (ReLU) activation layer, a 68.75-25-6.25 training-validation-testing set split, adaptive moment estimation (ADAM) optimizer, a learning rate of 0.001, mean squared error loss, cube root normalization, 360 batches during training before updating weights, and a single unit dense layer to pass the final output. One may think to stack multiple LSTM networks to achieve better performance, however, in testing, the use of 2 or more LSTM layers degraded performance rather than enhancing it. Finally, each of the models was trained for up to 100 epochs, however the training implemented early stopping, occurring at around 20 epochs on average, based on the validation loss statistic and saved only the best model during the training to avoid overfitting the model to the dataset. From this testing we are able to determine the best amount of time history to train with, aiming to supply immediate information to LSTM with the known previous solar wind data. One minute time history was chosen as a starting choice because it is the standard LSTM setup. The other 30-, 60-, and 120-min time histories are selected in consideration of variable propagation of solar wind and IMF information from the subsolar bow shock to geomagnetic activity (Connor et al., 2014; Maggiolo et al., 2017).
FIGURE 1. Subfigure (A) showcases the setup of the LSTM |BN| model utilized in the study. The inputs are structured in a 2D array where each column is one feature of the solar wind and IMF inputs detailed in Section 2.2 while the rows refer to each of the features at a previous time point from t−1 to t−n, where n is the amount of time history selected. The input array is provided directly to the LSTM model which sends the output to a single dense layer to provide the final output. In subfigure (B) the polarity model is displayed in a similar manner to the LSTM model of subfigure (A). The input array is a 1D array sent to an LSTM model which sends the final prediction to a single dense layer. The output of the polarity model is decoded and multiplied against the LSTM |BN| model to obtain the coupled output.
Figures 2A–G compares geomagnetic field predictions of the four LSTM models described for the 07-Sep-2015 geomagnetic storm. Figures 2A–E show IMF in GSM coordinates, solar wind speed, density (black) and flow pressure (red), temperature, and AE (black) and SYM/H (red) geomagnetic indices. Figure 2F shows the magnitude of north-south component of the geomagnetic field (|BN|) as observed from CMO (black) and predicted by the LSTM models using different time histories as input (dashed lines). Figure 2G shows the time derivative of the north-south geomagnetic field component (d|BN|/dt) observed from CMO (black) and the value derived from the LSTM |BN| models (dashed lines) with gaps in predictions occurring where one or more input variables are missing due to a gap larger than 10 min. Figures 2H–K show other modeling experiments and are discussed later in the paper. In Figure 2F we can see that the performance between the four models for the event are quite close in |BN| predictions. The 1 min time history model shows a drop in |BN| around the 1600 UT mark whereas the other 3 models show a drop occurring an hour later. When looking at the d|BN|/dt values we find the 30 min model shows the most activity, though vastly underestimated, where the other 3 models predict closer to 0.
FIGURE 2. Storm time test predictions of different LSTM configurations for the 07-Sep-2015 storm. Subplots (A–E) show the IMF, solar wind, and geomagnetic indices conditions with the model predictions in the subplots beneath. Subplots (F,G) show the results of predictions utilizing different time histories of input. From these we see the 1-min time history model consistently predicts values less than that of the other 3 models. Subplots (H,I) show the results of predicting for |BN| utilizing the 30-min time history both with and without d|BN|/dt as an input parameter. Likewise, subplots (J,K) utilize the 30-min time history model to predict BN with and without dBN/dt information added as input. From these plots we find that LSTM has an affinity for predicting |BN| with d|BN|/dt information included in the input parameters.
Individual time history model performance is also seen in Figures 3A–D, where the predicted and observed value for the CMO station across the year of 2015 are plotted for correlation. In Figure 3 the dashed red line indicates a perfect correlation between the prediction and observed values and the solid red line shows the line of best fit. The legend in the upper right corner shows the equation for the line of best fit and the Pearson correlation (r) between the predicted and observed data. In these plots we find the Pearson correlation (r) between the 4 models is relatively similar with small increases corresponding with time history supplied to the model. Looking at these values we see an 11% increase in correlation from the 1 min to the 30 min model. As the input data increases we find an 8% increase when doubling the time history from 30 to 60 min and a 2% increase when increasing time history from 60 to 120 min. Looking at the lines of best fit (solid red) we can see a slow shift towards the perfect correlation (dashed red) line as input data increases, from 1 min time history to 120 min. The 30-min time history model was selected as the suitable model due to a reduction in model complexity with nearly similar results as the more complex 120-min model. Further, we find the 30-min time history more suitable for our long term application of GIC predictions due to its performance in event time d|BN|/dt fluctuations as indicated by Figure 2G.
FIGURE 3. Storm time density plots using 2015 as a test set showing the performance of the models shown in Figure 2. The dashed red line indicates a perfect 1:1 correlation between predicted and observed values. The solid red line is a line of best fit correlation between the predicted values and the observed values. Subplots (A–D) show density plots of predicted vs observed values for 2015 CMO when utilizing 1−, 30−, 60−, and 120-minutes of time history. Subplots (E–H) utilize the selected 30-minute time history model and compare predicting |BN| with and without d|BN|/dt information and predicting BN with and without dBN/dt information. We can see in subplots (G) and (H) that LSTM predicting BN provided significantly poor predictions with a correlation of 0 regardless of utilizing previous dBN/dt information as an input to the model. However, when predicting |BN| as seen in subplots (E) and (F) LSTM performs well, with a significant 74% increase in Pearson correlation when including previous d|BN|/dt information as an input parameter to the model.
In Figures 2H,I the 30-min time history model was selected and trained for |BN| with (red line) and without (green line) previous known geomagnetic field time derivative (d|BN|/dt) history as an input parameter. This is also found in Figures 3E,F, which show an increase in Pearson correlation between the predicted and observed CMO geomagnetic field values when d|BN|/dt is included from 0.50 to 0.87. The use of d|BN|/dt as in input variable is spurred by the standard LSTM setup where previous observations variables are passed as a standard input for the upcoming prediction. In the case of our model, passing |BN| at t-1 information creates a feed forward network where the model passes a value with high correlation to |BN| at t-1 as the prediction. The use of d|BN|/dt removes this outcome while providing the model with general information regarding the strength of fluctuations. We acknowledge concern that by using d|BN|/dt as input, the LSTM models may act as a first-order Taylor series expansion (i.e., |BN| (t) = |BN| (t-1) + d|BN|/dt (t-1)*dt. However, our approach differentiates from this expansion by not utilizing |BN| at t-1 as input. Additionally, it uses 30 min of time history of d|BN|/dt and SW/IMF parameters as input. One may find concern that the LSTM model will find high correlation between prediction |BN| (t) and the input feature d|BN|/dt (t-1), however in testing our model does not show such behavior, suggesting the LSTM model learns more complex behaviors than the first order Taylor series expansion.
Lastly, Figures 2J,K show the results of two 30-min time history models trained with (red line) and without (green line) dBN/dt input predicting BN. Both of these models performed poorly, predicting a nearly straight line at 0 nT, with density plots in Figures 3G,H indicating no correlation between the predicted and observed values for CMO in 2015 while predicting BN. From the eight models tested we find that the 30-min model with d|BN|/dt input and |BN| output retains the best performance regarding predictions during storm time while limiting the complexity of the model. The model configuration is then utilized across all 4 stations with the 2000–2015 training set.
Geomagnetic Field Prediction Across the Alaska Chain
The magnetometer chain of FYU, CMO, PKR, and KAV were chosen to create a perpendicular line prediction of the geomagnetic field with respect to the auroral oval. For each of the stations the 30-min time history was chosen for the performance found in Section 3.1. The models employ the same early stopping mechanism described in 3.1 to avoid overfitting of the data. While the models in Section 3.1 were trained from a smaller limited dataset (i.e., 2009–2015), the models trained for each station utilized data from January 2000 to December 2015 for training with 68.75% of the data as the training set, 25% as the validation set, and 6.25% as the test set. Like the models in Section 3.1, the year of 2015 was separated and used as a testing set of the models due to its high prevalence of geomagnetic activity throughout the year. Further, for each station an additional MLR model was made with the same 2000–2015 dataset, utilizing the same input variables and cube root normalization as the LSTM. In some cases the MLR models predicted negative values of |BN|, which would imply an un-physical value of BN. The negative values were removed from the MLR predictions dataset as they are un-physical quantities for |BN|.
In Figure 4 the results of the models are plotted for the 07-Sep-2015 geomagnetic storm and subplots Figures 4A–E follow the same format as Figures 2A–E. In Figures 4F,G we see the predictions for the LSTM model (red) and MLR model (green) compared to the observed values (black) for the FYU station with Figure 4F showing the |BN| predictions and Figure 4G showing d|BN|/dt predictions. Likewise, Figures 4H–M show the |BN| and d|BN|/dt prediction results for PKR, CMO, and KAV, respectively. From these figures we see that the MLR models underestimate |BN| more than the LSTM models. When considering the d|BN|/dt plots we find that the LSTM shows more frequent strong fluctuations, with the caveat they are not always predicted at the correct time. In Figure 5 the LSTM and MLR model predictions are plotted against the observed station values for the year of 2015 in the same format as Figure 3 to test for correlation. The LSTM models in Figures 5A,C,E,G show high Pearson correlation coefficients of 0.80–0.86, while the MLR models in Figures 5B,D,F,H show much lower correlation coefficients of 0.68–0.73. Looking at the average of the Pearson correlation values, LSTM shows a 18.2% increase in performance over the MLR models.
FIGURE 4. Storm time predictions utilizing LSTM (red) and MLR (green) models showing their ability to predict the 07-Sep-2015 event for each of the 4 individual stations across Alaska. Subplots (A–E) show IMF and solar wind properties for the event while subplots (F–M) show the observed and predicted |BN| and d|BN|/dt for the four selected stations. Here we can see that for all stations the LSTM model performs better at matching |BN|.
FIGURE 5. LSTM and MLR density plots using the 2015 test set for all 4 stations. Subplots (A, C, E, G) show prediction results for the LSTM models of each station while subplots (B, D, F, H) show prediction results for the MLR models of each station. Further, we can see that the lines of best fit (solid red) between the LSTM predictions and observed data are much closer to the perfect prediction lines (dashed red) than those of the MLR predictions and the observed data with a 18.2% increase in the average Pearson correlation values.
Polarity and Coupled Model
One of the problems with the LSTM model is its affinity towards predicting |BN| over BN. This means that to obtain the best results we train off the magnitude, thus losing sign information in the process. This lost feature, which we are calling polarity, is a necessary component of the ionospheric current modeling, since ionospheric current directions (i.e., eastward or westward electrojets) can be inferred by the sign of BN. Initially, LSTM seemed to favor predicting when all values in the observation set were positive. To resolve this problem, we applied multiple scaling methods to the dataset when predicting BN. These methods comprised of scaling the data between 0 and 1, −1 and 1, and linearly shifting the BN data above 0 by adding the minimum BN value to each data point. We developed the LSTM models for predicting the scaled BN values and reversed the scaling back to the original scale for the finalized prediction. However, these scaling approaches did not show any significant improvement when compared to the predictions of the original cube root normalized dataset. Therefore, a different approach to polarity utilizing a secondary model was developed. To create this model, the polarity was encoded into a 1 for positive values and 0 for negative values, then this observed polarity was trained for using a LSTM kernel looking at the previous minute of MLT, |BN|, and d|BN|/dt information. With the polarity retained through a secondary model we are able to decode the polarity and multiply it through the geomagnetic field model trained for predicting |BN| and retain BN in a coupled model technique. Figure 1 summarizes the setup of the polarity and coupled models.
Figure 6 shows the results of the coupled model approach applied to all stations for the 09-07-2015 test storm. In the left column we can see the predicted values follow the observed data values, with one false positive prediction around 18 UT for the CMO station. As seen in the right column of Figure 6, the Pearson correlation between predicted and observed geomagnetic field increased for all stations from an average 0.84 to 0.88. This indicates improved performance of a coupled model than the LSTM in Figure 5 while preserving the critical sign information needed for ionospheric current analysis. Further, BN fluctuations observed in Figure 6 show persistence of negative enhancements for the 09-07-2015 storm which the coupled model properly captures with minimal unexpected sign flips. For KAV we can see around 14 UT and 17 UT the coupled model properly flips positive for the short positive durations seen in the observed data during a predominately negative enhancement. This makes the coupled model approach a promising modeling method for our study.
FIGURE 6. Coupled model storm predictions. In the left column the coupled models ability to predict BN for the 07-Sep-2015 storm is plotted against the observed data for each station. In the right column the density plot for each station is shown for the 2015 dataset. The predictions of the coupled models provide better Pearson correlation coefficients than the original LSTM |BN| while following the storm time fluctuations, indicating the Polarity model is providing the necessary information to the |BN| prediction.
Skill Scores and Model Performance
Heidke skill scores (HSS) are a widely accepted method of determining machine learned model performance by testing multiple thresholds to understand model sensitivity and variability at discerning the desired output variable. These scores can be see in Table 1 ranging from 0 (random prediction) to 1 (perfect prediction) and have been split between scores for |BN| and d|BN|/dt sensitivity. The scores are calculated based on whether the predicted and observed values cross the threshold at the same time, with thresholds for d|BN|/dt selected based off of Pulkkinen et al. (Pulkkinen et al., 2013) and thresholds for |BN| selected from the 50, 75, 85, and 99 percentile of the |BN| values over 2000–2014 for the 4 stations. We can see from Figure 5 and Table 1 that the LSTM models show a 18.2% increase in Pearson correlation over the MLR models with the LSTM models achieving HSS above 0.5 for 87.5% of the selected |BN| thresholds for both storm time and 2015 while MLR scored above 0.5 for 18.75% for the same thresholds. When focused on the 2015 HSS of the 99 percentile |BN| threshold we can see the LSTM models outperform the MLR models in all cases, showing a 152.2, 66.6, 95.8, and 215.8% increase for FYU, KAV, PKR, and CMO stations, respectively. Likewise, when focusing the 99 percentile |BN| threshold storm HSS we find a 176.9, 16.1, 510, and 828.6% increase for FYU, KAV, PKR, and CMO, respectively.
TABLE 1. HSS of selected |BN| and d|BN|/dt thresholds for the LSTM and MLR models. Scores are evaluated through direct comparison on a minute by minute basis across the year of 2015 and separately for the 07-Sep-2015 storm. Polarity has been encoded into a value of 0 (−) or 1 (+) and the HSS for this corresponds to accurately assessing a 1.
The HSS results for the d|BN|/dt thresholds show mixed results between the LSTM and MLR models. For the 2015 d|BN|/dt HSS, LSTM models generally perform better than the MLR by showing higher HSS for 10 out of 16 thresholds, while for the storm time d|BN|/dt scores the MLR models perform better by showing higher HSS for 10 out of 16 thresholds. However, our storm time performance testing is limited for only a single geomagnetic storm. For better understanding of storm time model performance, a more comprehensive testing is required with a large number of geomagnetic storms.
The neural networks created are adequate at predicting the overall strength of the field, as indicated by high correlation and good HSS scores for |BN|, but not the variability of the field. This can also be seen in Figures 4G,I,K, and m where the d|BN|/dt of the LSTM models can be seen to exceed that of the observed values at times predating or postdating the observed fluctuations. Additionally, the d|BN|/dt scores for the LSTM models is less than 0.5, regardless of looking at storm time or the full 2015 test year. The poor d|BN|/dt scores of LSTM models are understandable since the models were made to predict |BN| and then derive d|BN|/dt from predicted |BN|. Further, we see a pattern where the scores for d|BN|/dt are decreasing as the threshold increases, which is a pattern found in other recent machine learning studies of geomagnetic fields (Camporeale et al., 2020; Smith et al., 2021), implying it is a common challenge for data-driven models.
Despite the low performance in d|BN|/dt, the performance in predicting |BN| is promising, especially when coupled with the secondary polarity model. The polarity models created scored 0.9 or above at predicting the encoded polarity value, which when multiplied through their respective station the LSTM prediction generally retained or increased in HSS for the |BN| thresholds of the original |BN| model. The HSS within Table 1 show an enhancement within the storm time dBN/dt scores for 14 out of 16 thresholds when using a coupled model approached, while 15 out of 16 2015 threshold scores stay the same or increase. The minimal increase of the coupled models BN HSS from the LSTM models |BN| HSS is understandable because BN is mainly determined by the |BN| models and not by the polarity model. However, the coupled model provides important sign information of |BN|, creating overall larger dBN/dt, and thus boosting the original dBN/dt HSS of the coupled models from the d|BN|/dt HSS of the LSTM models. Additionally, there are more data points for both BN and dBN/dt of the coupled models than the LSTM models, since the dataset for the LSTM models had the data points where the MLR models predicted negative removed. These additional data may play a role in increasing the BN and dBN/dt HSS of coupled models.
There are concerns for the Polarity models to generate false peaks within the d|BN|/dt predictions due to overpredictions in the |BN| models before and after a sign transition or a sign change with incorrect timing. However, the occurrence of this is rare because the current models underpredict |BN| leading to overall lower negative to positive and vice versa jumps within the dataset and thus lead to lower dBN/dt values. A majority of the unexpectedly large d|BN|/dt peaks are due to the original |BN| models predicting a large fluctuation occurring at the wrong time, which persist when multiplied through by Polarity. The Polarity model, while scoring well, fails predominately at times where there is a flip from positive to negative and vice versa. The delayed or preemptive timing of polarity switching may in turn cause unexpected d|BN|/dt patterns, however the vast majority of these patterns will arise during quiet times where the geomagnetic field is fluctuating minutely around 0 nT.
The LSTM models, while promising, have a few different caveats to them in their implementation. The first and foremost caveat is that the models only predict |BN| one time step ahead, which is limiting in advanced GIC prediction. For future work we aim to provide a 10–20 min prediction range thus increasing the model’s practicability. There are two main methods of achieving this, one of which is utilizing the initial t+1 prediction as the starting point of a secondary network that spans the desired prediction range. This setup requires that the secondary network be initialized with the full training dataset to set the internal states of the model. In practice, this is a time consuming approach, which will only increase as more data is used to create the models. This is opposed to our justification for using the 30-min time history model over the 120-min history model, since complexity of the model will play a factor in time to set the states every time a new prediction is made. Another method is to determine the probability of strong d|BN|/dt fluctuations in a set period. This type of approach has been utilized in other studies, such as the one by Maimaiti et al. (2019) to determine the probability of geomagnetic substorm onset and the studies by Smith et al. (2021) and Camporeale et al. (2020) to forecast the probability of specific dB/dt (i.e., surface geomagnetic field time derivative) thresholds. However, this approach for GIC prediction has ambiguity in whether the outcome will occur at the beginning, middle, or end of the window, which may be pertinent information to the end user.
Secondly, our models currently only predict |BN| and in combination with a polarity model, BN, while GICs are influenced by the surface geomagnetic field which is made up of both BN and BE (Ngwira et al., 2008; Bedrosian and Love, 2015; Lotz and Cilliers, 2015). For future work we aim for the models to predict two output variables, |BN| and |BE|. The surface geomagnetic field information, combined with ground conductivities, can be used to determine geomagnetically induced electric fields within the Earth’s surface and GICs in specific electrical systems (Ngwira et al., 2008). The current models are limited in this capacity as local conductivity information is currently unavailable for the Alaska region. Until proper conductivity information is available, the models may still be utilized as an indication system, since GICs oftentime occur with large geomagnetic field perturbations that our models are intended to predict.
Thirdly, our models generally underpredict the geomagnetic field strength and do not properly capture the time variations observed in the data. This underprediction is commonly seen in other Machine Learning models (e.g. Keesee et al. (2020)) and likely due to the choice of training with both quiet and storm time data, where the quiet time data makes a larger portion of the trained set leading to lower overall geomagnetic field predictions. A possible solution around this is to train solely off of storm time data, however this adds an extra layer of complexity to the training process during model training and has been shown to not completely solve the problem on its own (Pinto et al., 2022). However, this approach improves general model performance and is something we plan to implement in the future. Another possible approach would be to create a model that takes the incoming prediction and computes the likely offset for that value to better retain the strength of the observed field data, though such a setup would likely require careful consideration in implementation. Even with the aforementioned approaches, the machine-learning based model may find difficulty in predicting “once-in-a-lifetime” singular events, for example, the Carrington event on 1–2 September 1859 (Green and Boardsen, 2006; Cliver and Dietrich, 2013), because such severe events are very rare. In such a case, a physics-based model is a good alternative for the dB/dt predictions.
Lastly, the internal setup of LSTM requires full time history information of the incoming data, which this study does not have access to due to gaps in data sources. During the training of the models, if data is not split into groups of continuous data, the model will not understand the presence of gaps assuming that the incoming data is continuous. In a real-time situation data outages will occur and the model will continue with these data gaps present in the dataset, unless reset to avoid, which would result in 30 min of time without predictions at a minimum while waiting for a continuous stream of data. With this in mind, we chose to train the model with data gaps included within the dataset, which can be seen in Figures 2, 4 where gaps in data are present, allowing the model to pick up immediately when new information is available rather than requiring a down time to fill a continuous segment. This is done after the linear interpolation by passing the model the dataset as-is and allowing the training algorithms to assume the set is continuous. However, we understand this decision results in performance decreases in datasets with many or expansive gaps in data, such as the PKR dataset which was missing the year of 2010 within the supplied dataset. Despite the known limitations of training from incomplete datasets, we find the models still attain high Pearson correlation coefficients and promising HSS for |BN| and BN prediction, and will benefit from on-going efforts for continuous dataset collection.
While the caveats mentioned above are inherent to the utilization of LSTM and our current efforts in producing a model, there is one limitation within our model that puts a restraint on where these models may be used. The models produced so far rely on past known geomagnetic field information (i.e., d|BN|/dt in the previous 30 min) in the local region, which may not be suitable for all areas looking to perform GIC risk assessment. The Alaska region can accommodate this since the magnetometers chosen also provide real-time data which may be used with the models for real-time predictions. Additionally, the Space Weather UnderGround outreach project initiated at the University of New Hampshire (Smith, 2020) and expanded to the University of Alaska Fairbanks, will build a cost efficient and research-capable array of magnetometers across Alaska and New Hampshire with a 1 nT/s resolution, increasing the spatial resolution of data in these regions. Moving forward our models will utilize the datasets provided by these arrays for forecasting GIC risk.
Potential LSTM Model Use and Future Work
With the inclusion of multi-minute output, the LSTM models will be matured enough to create advanced GIC warnings based on likely dBN/dt thresholds without the need to directly predict GICs. The creation of GICs is a complex problem that is dependent on ground conductivity and the properties of the electrical device that it is being influenced. The land conductivities of Alaska have not been thoroughly studied though conductivity maps exist for the mainland of the United States. Due to the complexity of GIC prediction on every electrical system, it becomes practical to provide GIC warning and geomagnetic field predictions to the end user. In this manner, the end user can apply these predictions coupled with the knowledge of their own system to determine risk.
The increase of model performance via multi time history predictions and closer values to the observed data will allow the models to be utilized in ionospheric current predictions. This possible use case was the reason for pursuing BN predictions instead of |BN| leading to the creation of the polarity model. Previous studies have created modeling techniques to predict the ionospheric current based on local and/or global geomagnetic field patterns and electrodynamics (Lu et al., 1995; Kihn and Ridley, 2005; Vanhamäki and Juusola, 2020). A local model to determine the auroral oval requires, at a minimum, multi-point field values along a line perpendicular to the oval. Our current study is setup for this with the use of PKR, FYU, KAV, and CMO, which roughly lay on a line perpendicular to the auroral oval. With the inclusion of multi time history predictions our LSTM models coupled with ionospheric current modeling have the potential to forecast north-south motion or expansion of the ionospheric currents. Current patterns would be useful in region-based GIC risk assessment and awareness while also providing information to aurora enthusiasts and citizen scientists since strong ionospheric currents have been connected to auroral activity (Akasofu, 1989; Newell et al., 2001).
Summary and Future Work
This study aims to show the progression of LSTM neural networks trained to predict the geomagnetic field at individual stations across Alaska. To achieve this we trained 12 models (4 LSTM, 4 MLR, 4 Polarity) with NASA OMNI IMF and solar wind data coupled with SuperMAG geomagnetic field information from the years 2000–2015 split in a 68.75-25-6.25 training/validation/test set configuration. We produced 8 models to test the configuration of LSTM with our desired inputs and outputs. From these models we chose 30 min of time history utilizing IMF, solar wind, and past d|BN|/dt information as the most effective and applied the configuration to 4 stations across Alaska. We find that the LSTM models generally outperform the MLR models with respect to predicting |BN|, however the results of d|BN|/dt prediction performance are inconsistent between the two modeling methods. Due to the initial model performing best when predicting |BN| a coupled model approach was utilized to retain BN output with performance similar to the original |BN| model it was based on.
The models are limited to single point future predictions and generally underestimate the strength of the geomagnetic field. Future work on these models aims to increase the time history output to 10–20 min of future predictions while also increasing the performance of the models to estimate the geomagnetic field strength and variations in both the North-South and East-West components. With this will come potential integration into ionospheric current models for ionospheric current and auroral activity forecasting. Further, the models will be converted to a regional GIC risk assessment allowing the end user to apply the geomagnetic field predictions to their own electrical systems.
Data Availability Statement
MB was the primary author of the paper, conducted data preparation, model development, and participated in analysis. HC provided guidance on study concept and design, model development and analysis, interpretation of results, and writing. DO contributed to guidance on model development, data preparation, and general guidance. AK contributed to the study design and overall discussion. VP contributed to model development and methodology, overall discussion, general guidance, and writing. MM contributed to model development and methodology, overall discussion, and general guidance. CN contributed to general guidance and discussions regarding GICs. SP contributed to model methodology, data preparation, and overall discussion.
This work was supported by NSF EPSCoR Award OIA-1920965. HC gracefully acknowledges support from the NASA grants, 80NSSC18K1043, 80NSSC20K1670, and 80MSFC20C0019.
Conflict of Interest
CN was employed by the ASTRA LLC.
The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
We thank all members of the MAGICIAN team at UNH and UAF that participated in discussions leading to this article. We also thank USGS, INTERMAGNET, and the Geophysical Institute for supporting and maintaining the magnetometers used in this study as well as promoting high standards of magnetic observatory practice. We thank NASA OMNIweb and SuperMAG for providing organized data access that supported this study.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fspas.2022.846291/full#supplementary-material
Bedrosian, P. A., and Love, J. J. (2015). Mapping Geoelectric Fields during Magnetic Storms: Synthetic Analysis of Empirical united states Impedances. Geophys. Res. Lett. 42, 10,160–10,170. doi:10.1002/2015gl066636
Blake, S. P., Gallagher, P. T., McCauley, J., Jones, A. G., Hogg, C., Campanyà, J., et al. (2016). Geomagnetically Induced Currents in the Irish Power Network during Geomagnetic Storms. Space Weather 14, 1136–1154. doi:10.1002/2016sw001534
Brownlee, J. (2017). Long Short-Term Memory Networks with Python: Develop Sequence Prediction Models with Deep Learning. Jason Brownlee. Available at: https://machinelearningmastery.com/lstms-with-python/.
Camporeale, E., Cash, M. D., Singer, H. J., Balch, C. C., Huang, Z., and Toth, G. (2020). A gray-Box Model for a Probabilistic Estimate of Regional Ground Magnetic Perturbations: Enhancing the Noaa Operational Geospace Model with Machine Learning. J. Geophys. Res. Space Phys. 125, e2019JA027684. doi:10.1029/2019JA027684
Carter, B. A., Yizengaw, E., Pradipta, R., Halford, A. J., Norman, R., and Zhang, K. (2015). Interplanetary Shocks and the Resulting Geomagnetically Induced Currents at the Equator. Geophys. Res. Lett. 42, 6554–6559. doi:10.1002/2015GL065060
Connor, H. K., Zesta, E., Ober, D. M., and Raeder, J. (2014). The Relation between Transpolar Potential and Reconnection Rates during Sudden Enhancement of Solar Wind Dynamic Pressure: Openggcm-Ctim Results. J. Geophys. Res. Space Phys. 119, 3411–3429. doi:10.1002/2013JA019728
Dimmock, A. P., Rosenqvist, L., Hall, J. O., Viljanen, A., Yordanova, E., Honkonen, I., et al. (2019). The GIC and Geomagnetic Response Over Fennoscandia to the 7-8 September 2017 Geomagnetic Storm. Space Weather 17, 989–1010. doi:10.1029/2018SW002132
Dimmock, A. P., Rosenqvist, L., Welling, D. T., Viljanen, A., Honkonen, I., Boynton, R. J., et al. (2020). On the Regional Variability of Db/dt and its Significance to Gic. Space Weather 18, e2020SW002497. doi:10.1029/2020SW002497
Fiori, R. A. D., Boteler, D. H., and Gillies, D. M. (2014). Assessment of Gic Risk Due to Geomagnetic Sudden Commencements and Identification of the Current Systems Responsible. Space Weather 12, 76–91. doi:10.1002/2013sw000967
Keesee, A. M., Pinto, V., Coughlan, M., Lennox, C., Mahmud, M. S., and Connor, H. K. (2020). Comparison of Deep Learning Techniques to Model Connections between Solar Wind and Ground Magnetic Perturbations. Front. Astron. Space Sci. 7, 72. doi:10.3389/fspas.2020.550874
Khanal, K., Adhikari, B., Chapagain, N. P., and Bhattarai, B. (2019). HILDCAA‐Related GIC and Possible Corrosion Hazard in Underground Pipelines: A Comparison Based on Wavelet Transform. Space Weather 17, 238–251. doi:10.1029/2018sw001879
Lu, G., Lyons, L. R., Reiff, P. H., Denig, W. F., de la Beaujardiére, O., Kroehl, H. W., et al. (1995). Characteristics of Ionospheric Convection and Field-Aligned Current in the Dayside Cusp Region. J. Geophys. Res. 100, 11845–11861. doi:10.1029/94JA02665
Maggiolo, R., Hamrin, M., De Keyser, J., Pitkänen, T., Cessateur, G., Gunell, H., et al. (2017). The Delayed Time Response of Geomagnetic Activity to the Solar Wind. J. Geophys. Res. Space Phys. 122, 11,109–11,127. doi:10.1002/2016JA023793
Maimaiti, M., Kunduri, B., Ruohoniemi, J. M., Baker, J. B. H., and House, L. L. (2019). A Deep Learning‐Based Approach to Forecast the Onset of Magnetic Substorms. Space Weather 17, 1534–1552. doi:10.1029/2019SW002251
Ngwira, C. M., Pulkkinen, A. A., Bernabeu, E., Eichner, J., Viljanen, A., and Crowley, G. (2015). Characteristics of Extreme Geoelectric fields and Their Possible Causes: Localized Peak Enhancements. Geophys. Res. Lett. 42, 6916–6921. doi:10.1002/2015GL065061
Ngwira, C. M., Pulkkinen, A., McKinnell, L.-A., and Cilliers, P. J. (2008). Improved Modeling of Geomagnetically Induced Currents in the South African Power Network. Space Weather 6, S11004. doi:10.1029/2008SW000408
Oliveira, D. M., Arel, D., Raeder, J., Zesta, E., Ngwira, C. M., Carter, B. A., et al. (2018). Geomagnetically Induced Currents Caused by Interplanetary Shocks with Different Impact Angles and Speeds. Space Weather 16, 636–647. doi:10.1029/2018SW001880
Pinto, V. A., Keesee, A. M., Coughlan, M., Mukundan, R., Johnson, J. W., Ngwira, C. M., et al. (2022). Revisiting the Ground Magnetic Firld Perturbations challenge: A Machine Learning Perspective. Unpublished.
Pirjol, R. J., Viljanen, A. T., and Pulkkineni, A. A. (2007). “Research of Geomagnetically Induced Currents (GIC) in Finland,” in 2007 7th International Symposium on Electromagnetic Compatibility and Electromagnetic Ecology, St. Petersburg, Russia, 26-29 June 2007. doi:10.1109/EMCECO.2007.4371707
Pirjola, R., Pulkkinen, A., and Viljanen, A. (2003). Studies of Space Weather Effects on the Finnish Natural Gas Pipeline and on the Finnish High-Voltage Power System. Adv. Space Res. 31, 795–805. doi:10.1016/s0273-1177(02)00781-0
Pulkkinen, A., Lindahl, S., Viljanen, A., and Pirjola, R. (2005). Geomagnetic Storm of 29-31 October 2003: Geomagnetically Induced Currents and Their Relation to Problems in the Swedish High-Voltage Power Transmission System. Space Weather 3 (8). doi:10.1029/2004SW000123
Pulkkinen, A., Rastätter, L., Kuznetsova, M., Singer, H., Balch, C., Weimer, D., et al. (2013). Community-Wide Validation of Geospace Model Ground Magnetic Field Perturbation Predictions to Support Model Transition to Operations. Space Weather 11, 369–385. doi:10.1002/swe.20056
Rodger, C. J., Mac Manus, D. H., Dalzell, M., Thomson, A. W. P., Clarke, E., Petersen, T., et al. (2017). Long-Term Geomagnetically Induced Current Observations from New Zealand: Peak Current Estimates for Extreme Geomagnetic Storms. Space Weather 15, 1447–1460. doi:10.1002/2017SW001691
Rogers, N. C., Wild, J. A., Eastoe, E. F., Gjerloev, J. W., and Thomson, A. W. P. (2020). A Global Climatological Model of Extreme Geomagnetic Field Fluctuations. J. Space Weather Space Clim. 10, 5. doi:10.1051/swsc/2020008
Smith, A. W., Forsyth, C., Rae, I. J., Garton, T. M., Bloch, T., Jackman, C. M., et al. (2021). Forecasting the Probability of Large Rates of Change of the Geomagnetic Field in the uk: Timescales, Horizons, and Thresholds. Space Weather 19, e2021SW002788. doi:10.1029/2021SW002788
Wintoft, P., Wik, M., and Viljanen, A. (2015). Solar Wind Driven Empirical Forecast Models of the Time Derivative of the Ground Magnetic Field. J. Space Weather Space Clim. 5, A7. doi:10.1051/swsc/2015008
Keywords: space weather, GIC, geomagnetic storms, ground geomagnetic field, machine learning, neural networks, LSTM
Citation: Blandin M, Connor HK, Öztürk DS, Keesee AM, Pinto V, Mahmud MS, Ngwira C and Priyadarshi S (2022) Multi-Variate LSTM Prediction of Alaska Magnetometer Chain Utilizing a Coupled Model Approach. Front. Astron. Space Sci. 9:846291. doi: 10.3389/fspas.2022.846291
Received: 31 December 2021; Accepted: 24 March 2022;
Published: 02 May 2022.
Edited by:Peter Wintoft, Swedish Institute of Space Physics, Sweden
Reviewed by:Andrew Smith, University College London, United Kingdom
Octav Marghitu, Space Science Institute, Romania
Copyright © 2022 Blandin, Connor, Öztürk, Keesee, Pinto, Mahmud, Ngwira and Priyadarshi. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Matthew Blandin, email@example.com