\section{Introduction}

 Many ground systems are influenced or disrupted by extreme space weather, such as power infrastructure, satellite communication, television and radio links, as well as GPS navigation systems \citep{allen1989effects, kan1979energy, ding2007large, board2009severe}. Forecasts of space weather events would provide early warning to these systems, allowing for precautions to be taken to mitigate the effects of otherwise damaging activity. One area of particular interest is the modeling of variations in Earth's magnetic field, or geomagnetic activity. There exist many different classes of geomagnetic activity, such as magnetic storms, magnetospheric substorms, quiet day variation, etc \citep{chapman1962geomagnetism}. The physical processes behind many of these phenomena begin in the interior of the sun, where pressure builds to a release of a steady stream of radiation and energetic particles known as the solar wind. Geomagnetic storms at Earth are often the result of periods of drastically intensified solar wind known as coronal mass ejections. In recent years, the problem of forecasting such storms has gained much attention. Geomagnetic indices were developed to measure the general intensity of space weather from the ground \cite{kivelson1995introduction}. Recent works focus on modeling a particular class of geomagnetic activity by studying the relationship between solar wind measurements obtained from space weather datasets or derived from tertiary datasets, and a proxy for magnetic activity such as magnetic indices which summarize specific geomagnetic phenomena.
 
  In this work, we seek to predict multiple magnetic indices such as the auroral electrojet indices (AE, AU, AL), the disturbance time index (Dst), the Solar Radio Flux at 10.7cm (F10.7) index, and the Planetary K (Kp) index several hours in advance. Where the majority of works have focused on predicting a single magnetic index, we demonstrate an approach to simultaneously forecast several indices multiple hours in advance. In our work, we frame the problem as a multivariate sequence-to-sequence learning task \citep{sutskever2014sequence} where the objective is to predict the four magnetic indices using solar wind measurements. We utilize long short-term memory networks which have been applied to myriad applications involving sequential data and are known for addressing vanishing gradient issues encountered with vanilla recurrent neural networks \citep{hochreiter1997long}. 

\vspace{-.1in}
\section{Related Work}

\vspace{-.1in}
\subsection{Magnetic Index Forecasting}

Many works focus on modeling a particular class of geomagnetic activity by examining a single magnetic index. One example area of research is that of studying solar wind coupling where the objective is to compute some coupling function and solve a least-squares regression problem to predict a particular magnetic index using the time-lagged coupling function. \citet{mcpherron2015optimum} introduced a coupling function which partially linearizes the relationship between a set of commonly used solar wind variables and solved a least-squares regression problem to predict the AL index one hour in advance. They reported an average coefficient of determination ($R^2$) of $.67$ over two-month periods spanning 1963-2014 which improves upon other solar wind coupling functions \citep{newell2007nearly, rostoker1972response, clauer1986relationship}.

Other works which utilize neural networks have emerged in recent years. \citet{revallo2014neural} investigated forecasting the Dst index one hour in advance using neural networks, and \citet{boberg2000real} sought to predict the Kp index in real time using solar wind measurements. In \citet{maimaiti2019deep}, the authors derived a substorm classification dataset by applying the substorm identification criteria \citep{newell2011substorm} to the SuperMAG SML index \citep{gjerloev2012supermag}. The forecasting approach involves training variations of (previously) state-of-the-art convolutional networks, such as ResNet \citep{he2016identity} to predict the onset of a magnetospheric substorm within the next hour.

\vspace{-.1in}
\subsection{Solar Wind Forecasting}
Another area of research is the prediction of Solar Wind measurements at L1. NASA's solar dynamics observatory curates a database of ultraviolet images of the solar corona \citep{galvez2019machine} which has been used by many works to forecast measurements at L1. Works typically model complex coupling processes which occur between the Sun and the Earth's magnetic field in order to predict solar wind speed (V) or interplanetary magnetic field strength (B). \citet{upendran2020solar} used a CNN-LSTM architecture leveraging a GoogLeNet \citep{szegedy2015going} as a feature extractor to forecast daily averages of the solar wind multiple days in advance. \citet{chandorkar2019dynamic} proposed a new regression framework to model the varying time lag implicit to the coupling processes between the solar magnetic field and solar wind propagation around the heliosphere.

\vspace{-.1in}
\section{Data and Preprocessing}

\vspace{-.1in}
\subsection{Solar wind and IMF Measurements}
 We use measurements of the solar wind speed $V$, proton density $n$, and the components of the interplanetary magnetic field (IMF) $B_x, B_y, B_z$.  These measurements are taken by the ACE and WIND satellites located at L1. We also utilize historical values of the magnetic indices AE, AL, AU, Dst, F10.7, and Kp. This data is curated by NASA's Space Physics Data Facility and is available through the OMNIWeb database\footnote{OMNIWeb data available at https://omniweb.gsfc.nasa.gov/}. As additional features, we include calendar year, day of year, and hour of day signals. These signals are converted to sinusoids with periods 11 years (roughly equivalent to the average length of a solar cycle), 365 days, and 24 hours in order to account for cyclical variations in the data. All data available on the OMNIWeb database has been preprocessed and resampled to several different time resolutions. For our experiments, we utilize hourly resolution data. 

\vspace{-.1in}
\subsection{SuperDARN Measurements}

We also utilize measurements derived from the SuperDARN network of radars \citep{greenwald1995darn}. Specifically, we used the spherical harmonic fit technique developed by \citet{ruohoniemi1989mapping} to produce a database of convection patterns covering the years from 2013 to 2017 (inclusive) on a five-minute cadence. These measurements were subsequently averaged and resampled to a 1 hour resolution. From these patterns we extracted two parameters, the Cross-Polar-Cap Potential (CCP), which parameterizes the degree to which solar wind energy is coupled into the Earth’s magnetosphere, (hence is one of the best parameters available for describing the level of global-scale magnetospheric activity), and the Polar Cap Radius (PCR), which is the average colatitude of the Convection Reversal Boundary. These measurements provide the best proxy for magnetic conditions at the poles aside from the auroral electrojet indices themselves. For a more detailed discussion, the reader is referred to \citep{bristow2007superposed}.

\vspace{-.1in}
\section{Model Training}

\vspace{-.1in}
\subsection{Dataset Assembly}
We treat this forecasting problem as a multivariate sequence-to-sequence learning problem. At some time $t$, the inputs to our model comprise the samples of the time series of all input features in the interval $[t - T_h, t]$ where $T_h \geq 0$ is the number of samples in the past considered. The target is the samples of the time series of the four magnetic indices in the interval $(t, t + T_p]$ where $T_p \geq 1$ is the number of samples in the future predicted. Figure 1 illustrates the structure of input/output pairs on a subset of features. To assemble our datasets, we simply slide this structure through the entire time series one sample at a time. For all of our experiments, we discard all datapoints with missing values. To avoid violating causality, we split the data into training, validation, and testing sets sequentially with sizes being roughly 60, 10, and 30 percent, respectively. We center and scale the input features and training labels by the mean and standard deviation of each indiviual feature in the training set. We perform no additional preprocessing or filtering by the magnitude of the magnetic indices.

\begin{figure}[t]
  \centering
  \includegraphics[width=0.9 \textwidth]{pairStructure2.png}
  \caption{Structure of a single datum in our dataset. The blue shaded region between $t-T_h$ and $t$ comprises all samples in the input time series and the green shaded region between $t$ and $t + T_p$ comprises all samples in the target time series. In this example, $T_p = 6, T_h = 5$.}
\end{figure}

\vspace{-.1in}
\subsection{Model Training and Hyperparameter Tuning}
In this work, we consider a simple neural network consisting of a LSTM which embeds the input time series into a low dimensional space followed by a linear layer which performs a linear transformation on the final output hidden state of the LSTM to the dimensionality of the output time series. We use Mean Squared Error as our loss function, and optimize the weights of the network using the Adam optimizer \citep{kingma2014adam} with weight decay. We conduct all experiments using the PyTorch deep learning framework \citep{NEURIPS2019_9015}. We treat the optimizer learning rate, weight decay coefficient, batch size, hidden dimension, and the number of layers of the LSTM as training hyperparameters. For each experiment, we conduct a random search \citep{bergstra2012random} over 100 total models and train the models for 1350 epochs using the median early stopping rule. We select the best model ranked by the performance on the validation set for all experiments.

\vspace{-.1in}
\section{Results}
\textit{20 Year OMNIWeb Data} In this experiment, we use 20 contiguous years of OMNIWeb Data (without including SuperDARN data) ranging from the years 2000-2019. We use a time history and a lead time of $T_h = T_p = 6$ hours. Our training, validation, and testing datasets correspond roughly to the years 2000-2011, 2012-2013, and 2014-2019, respectively. Our baseline for these exeriments is a persistence forecast in which the prediction for the magnetic indices for all $t$ is equivalent to the most recent available observation in the historical indices. In this experiment, we report performance in terms of pearson corellation coefficient $\rho$ computed over the testing set. Table 1 illustrates the performance drop off over time for our network (mNet) and a persistence forecast (pers). The purpose of this experiment is to compare to a baseline persistent forecast, but it should be noted that we also compare to \citep{mcpherron2015optimum} by computing the coefficient of determination (or explained variance) for a 1 hour AL forecast using predictions from the same network. We find that we achieve an $R^2$ of $.744$ over our test set ranging from 2014-2019. For this experiment, the best performing model uses a one layer LSTM, a batch size of 502, hidden dimension of 58, an initial learning rate of .00368, and a weight decay parameter of 9.7E-8.

\begin{table}[h]
\small
\caption{Pearson correlation coefficent for mNet and baseline (pers) predictions with test labels.}
\label{sample-table}
\begin{tabular}{ccccccccccccc}
\hline
 & \multicolumn{2}{c}{\textbf{AE}} & \multicolumn{2}{c}{\textbf{AU}} & \multicolumn{2}{c}{\textbf{AL}} & \multicolumn{2}{c}{\textbf{Dst}} & \multicolumn{2}{c}{\textbf{f10.7}} & \multicolumn{2}{c}{\textbf{Kp}} \\ \hline
\textbf{Hrs} & \textbf{mNet} & \textbf{pers} & \textbf{mNet} & \textbf{pers} & \textbf{mNet} & \textbf{pers} & \textbf{mNet} & \textbf{pers} & \textbf{mNet} & \textbf{pers} & \textbf{mNet} & \textbf{pers} \\
\hline
1 & .901 & .478 & .885 & .491 & .874 & .422 & .976 & .790 & .995 & .993 & .937 & .658 \\
2 & .782 & .452 & .789 & .462 & .743 & .397 & .949 & .759 & .994 & .991 & .882 & .632 \\
3 & .696 & .429 & .715 & .437 & .656 & .375 & .918 & .730 & .993 & .990 & .825 & .605 \\
4 & .650 & .406 & .665 & .414 & .616 & .352 & .892 & .703 & .992 & .989 & .790 & .578 \\
5 & .622 & .382 & .627 & .394 & .594 & .328 & .870 & .678 & .991 & .988 & .757 & .557 \\
6 & .598 & .364 & .598 & .379 & .574 & .309 & .850 & .655 & .990 & .987 & .726 & .535 \\
\hline
\end{tabular}
\end{table}

\textit{SuperDARN measurements} We compare the performance of our model on three different datasets to examine the added predictive capability of SuperDARN measurements. We start with the original OMNIWeb dataset as described in the previous experiment (base), create an additional dataset by replacing the historical values of the magnetic indices with the SuperDARN measurements (sdrn), and finally create a baseline dataset where we only use solar wind measurements (sw). The hyperparameters for the selected models are given in Table 2. For a fair comparison, we discard all datum corresponding to hours which are not represented in all three datasets. Similar to the first experiment, we use a time history and a lead time of $T_h = T_p = 6$ hours. Figure 2 illustrates the performance of the three models in terms of the pearson correlation coefficient of the predictions with test labels. 

\begin{figure}[b]
  \centering
  \includegraphics[scale=.28]{decorrFig.png}
  \caption{Pearson correlation for the predictions of the base model (blue/solid), sdrn model (orange/dotted), and sw model(red/dashed) with test labels.}
\end{figure}

We find that the SuperDARN measurements alone do not provide a substantial improvement over the base solar wind model, and thus we cannot conclude on their efficacy as a sufficient proxy for the magnetic indices. We also remark that the inconclusive results may be due to the limited data availability of SuperDARN measurements. To this end, we hope to increase the amount of SuperDARN measurements in future work. 


\begin{table}[t]
\small
\centering
\caption{Hyperparameters for the models in the \textit{SuperDARN measurements} experiment.}
\begin{tabular}{cccccc}
\hline
\textbf{Model} & \textbf{LSTM Layers} & \textbf{Weight Decay} & \textbf{Hidden Dimension} & \textbf{Batch Size} & \textbf{Learning Rate} \\
\hline
\textbf{base} & 1 & 4.21E-8 & 59 & 1731 & 4.73E-2 \\
\textbf{sw} & 1 & 2.31E-7 & 44 & 1032 & 4.79E-2 \\
\textbf{sdrn} & 1 & 4.21E-2 & 48 & 1060 & 2.87E-6 \\ \hline
\end{tabular}
\end{table}





\vspace{-.1in}
\section{Conclusion}
\vspace{-.1in}

We have presented a sequence-to-sequence learning approach to forecasting several geomagnetic indices simultaneously. These results demonstrate the efficacy of such an approach to forecasting these magnetic indices, an improvement over the current state-of-the-art for magnetospheric substorm prediction as well as an ability to outperform a persistence forecast over a longer time horizon. We find that the addition of the currently available SuperDARN data does not improve the predictive power of our models substantially. We seek to include additional results after more SuperDARN data has been processed in the future. Much additional future work is planned. For example, a larger model class may result in improved forecasts. Many recent works have emerged for similar time series forecasting problems including the use of attention, one-dimensional convolutional networks, and other recurrent neural networks. We seek to investigate these by performing a model architecture search \citep{elsken2018neural} or by using pretrained models similar to the approach used in \citep{maimaiti2019deep}. Additionally, prior works have used the ultraviolet coronal image dataset curated by NASA to forecast the solar wind and magnetic conditions at L1 orbit point at a coarse granularity. One potential future direction is to utilize a similar approach to forecast the magnetic indices directly and model the Earth-sun system in an end-to-end fashion. Such an approach may benefit from a pretrained convolutional network. Lastly, we seek to investigate the application of physics-informed neural networks to this problem by exploiting domain knowledge as a regularizer, such as the relationship between the AE, AL, and AU indices.

\vspace{-.1in}
\section*{Broader Impact}
\vspace{-.1in}
Space weather is a broad term used to describe many potentially catastrophic phenomena for which we currently have limited forecasting capability. We know such events will happen, we know the consequences can be disastrous, but we do not have the sophistication to predict when or how a particular event will evolve. One example is a famous 1859 `Carrington Event' that occurred during which northern lights were visible from Central America. Communication infrastructure was not nearly as advanced as it is today, but early telegraph stations caught on fire from induced electrical activity. The same thing could happen today to our electric power grid. The National Academies of Sciences estimates that an event of similar magnitude occurring today could cause several trillion dollars, knocking out power for many countries for many months. Currently, alerts to power utility companies about potential space weather disruptions are based exclusively on these geomagnetic indices, thus improving forecasts of these indices would immediately provide a tangible benefit to society. Space weather forecasting remains in its infancy, in part due to the cacophony of physical processes that cascade from the Sun to the Earth when a solar storm occurs. For decades, physicists have investigated this problem and have made substantial advances in modeling individual components of these complex systems, but have struggled to model the system with analytical techniques in an end-to-end fashion. As such, the broad impacts of effectively applying machine learning to the space weather forecasting problem could be huge, as ML techniques have a unique ability to describe systems that cannot be modeled analytically.
    
\vspace{-.1in}
\begin{ack}
This work benefited directly from conversations with Jacob Bortnik and Ryan McGranaghan. We would like to thank them for their contributions.

This work was supported by the Defense Advanced Research Projects Agency (DARPA) under cooperative agreement D19AC00009 from the Department of the Interior to Georgia Tech. The first author is supported by the NDSEG fellowship.
\end{ack}

\bibliographystyle{unsrtnat}
