導入 ~宇宙天気とDst指数~

　さて、みなさんは宇宙天気という言葉を耳にしたことがあるでしょうか。宇宙天気とは、一般的に太陽から地球までの宇宙空間環境の状態を表す言葉です。この宇宙空間は常に一定というわけではなく、太陽活動に大きく影響を受けます。とくに、太陽で巨大な爆発が発生し、それによって生じた太陽風が地球磁気圏と呼ばれる領域に到達した場合は、地球周辺の宇宙環境に大きな影響を与えます。日本では、とくに情報通信研究機構で盛んに研究が行われております。詳しくはこちらのリンクを御覧ください。

　宇宙天気の状態を表す指数というものが存在し、その中の１つにDst指数というものがあります。Dst指数は、グローバルな地磁気の擾乱具合を表すパラメータとして知られています。以下の時系列データのように、Dst指数が大幅に減少した場合(一般的には -100nT 以下まで達した場合)に磁気嵐と呼ばれます。非常に大規模な磁気嵐が発生すると、人工衛星の故障、通信環境の混乱、電力網の破壊などが発生する可能性があります(参考:1989年3月の磁気嵐)。ですので、このようなDst指数の大幅な減少を予測するということは非常に重要となります。

　この磁気嵐の規模を示すDst指数は、太陽風内の物理量(磁場・速度・プラズマ密度等)と大きく関連しています。これは、太陽風のエネルギーの地球磁気圏への流入量がその物理量で決まり、その蓄積されたエネルギーが解放される過程が磁気嵐と考えることができるためです。それでは、太陽風とDst指数の時系列データを用いて、数時間後のDst指数の予測を試みていきたいと思います。なお、今回の予測の大部分は Gruet et al. (2018)の論文を踏襲しています。興味のある方は(専門用語が難しいですが)論文を参照してください。

データについて

データ取得元

　使用するデータは、Dst指数の時系列データと太陽風物理量(南向き磁場・磁場強度・太陽風速度・プラズマ密度)です。以下にその時系列データの一例を示します。

　Dst指数のデータは、京都大学大学院理学研究科附属地磁気世界資料解析センターのものを使用しています。なお、Dst指数は1時間値となっています。また、太陽風データはOMNIと呼ばれる、地球磁気圏付近に到達したときの太陽風データの5分値を用いています。

データセット(学習データ・テストデータ)

　データセットは、2000年1月1日~2018年12月31日のうち、Dst指数が -100nT を下回った周辺の時間帯を集めて作成しました。そのようなイベントは67イベント存在し、それを学習データ49イベントとテストデータ18イベントに分けました。学習データ全体で約3900時間ほど、テストデータ全体で約1700時間ほどのデータとなっています。

予測モデルの構築

　それでは、いよいよ予測モデルの構築を行いたいと思います。予測モデルは、LSTMを用いて以下の図のようにしました。
　
　現在時刻 $T$ から、その時刻を含んで過去5時間までのデータを利用します。太陽風データについては、5分値となっているので、1時間分の平均値を使います。6時間分の各物理量のデータを順にLSTM Cellに流し込み、最後の出力を全結合層の入力とします。その出力をもう１度全結合層の入力として、その出力を $p$ 時間後のDstの予測値とします。今回は、以下のようなパラメータで実験しました。

• $p = 4$ [hours]
• Output dim: 125(LSTM), 25(Full Connected Layer 1)
• Dropout ratio: 0.5
• Batch size: 32
• Epochs: 100

上記の条件で学習を行ったモデルによる予測結果の例を以下の図に示します。

図の最上段の4時間先のDst指数の予測値(赤線)を見ると、それなりにうまく予測できているように見えます。次からは、さらに有用な予測とするためにもう一歩踏み込んでいきたいと思います。

ガウス過程による確信区間を含んだ予測

　上記の結果でも一定程度予測できていると言えなくもないですが、予測には不確実性がつきものです。そのような不確実性を考慮に入れることを目標として、ガウス過程を用いて予測値とその確信区間を同時に算出していきます。ガウス過程は様々なところで説明されていますが、例えばこのスライドが詳しいです。
　ガウス過程による予測値分布は、以下の式で表されます。

\begin{align}
 p(y^* | {\boldsymbol x}^*, D) &= N(m({\boldsymbol x}^*) +{\boldsymbol k}_*^TK^{-1}{\boldsymbol y},\, k_{**}-{\boldsymbol k}_*^TK^{-1}{\boldsymbol k}_*) \\
\end{align}

ただし、

\begin{align}
K &=
\begin{pmatrix}
k({\boldsymbol x}_1, {\boldsymbol x}_1) & \cdots & k({\boldsymbol x}_1, {\boldsymbol x}_N)\\
\vdots & \ddots & \vdots \\
k({\boldsymbol x}_N, {\boldsymbol x}_1) & \cdots & k({\boldsymbol x}_N, {\boldsymbol x}_N)
\end{pmatrix} 
\\
{\boldsymbol k}_* &= (k({\boldsymbol x}^*, {\boldsymbol x}_1),\; k({\boldsymbol x}^*, {\boldsymbol x}_2),\; \ldots,\; k({\boldsymbol x}^*, {\boldsymbol x}_N))^T \\
k_{**} &= k({\boldsymbol x}^*, {\boldsymbol x}^*) \\
k({\boldsymbol x}, {\boldsymbol x}') &= \theta_1 \exp\left(-\frac{|{\boldsymbol x}-{\boldsymbol x}'|^2}{\theta_2}\right)
\end{align}

です。$N$は学習データの数で、最後の $k({\boldsymbol x}, {\boldsymbol x}')$ の定義はRBFカーネルと呼ばれるものです。また、事前分布の平均値 $m({\boldsymbol x}^*)$ を何らかの形で定義しなければなりませんが、ここでは、これまでに構築してきたニューラルネットモデルの出力を用います。
　
　以上の条件でガウス過程回帰を行い、Dst指数の予測値をその確信区間とともに出力したものが以下の図になります。

確信区間(薄青色で塗りつぶされた区間)が、実際の観測値(赤線)をほぼをカバーしていることがわかると思います。これにより、観測値はおよそこの幅の間に収まるだろうと考えられるわけです。

　なお、今回はガウス過程回帰のためにGpyTorchを使用しました。詳しくは紹介しませんが、PyTorchに慣れている方は、非常に使いやすいと思います。

考察

　今回は、$p=4$ 時間後のDst指数の予測を行いました。上記の例を見ると、結果として、それなりに予測がうまくいっていると言えると思いますが、ではこれを$p=5, 6, ...$と伸ばしていくとうまくいくのかというと、そうはなりません。なぜなら、太陽風の影響がDst指数として現れるまでの時間スケールがある程度決まっているため、その時間の範囲内のものしか予測できないと考えられるからです。実際に、Gruet et al. (2018) の論文でも、$p=6$のときに予測が大きくずれるという事象が起こっていました。
　
　また、今回は予測がうまくいった例を挙げていますが、あまりうまくいかなかった例ももちろんありました。モデル選択やハイパーパラメータの調整などはあまり真面目にやっていないため、このあたりは改良の余地は十分にあるかと思います。例えば、Transformer などを使うと精度が向上するのではないかという期待をもっています。

まとめと今後について

　今回は、宇宙天気予報の一環として、Dst指数の予測に挑戦してみました。非常にマニアックな分野ではありますが、宇宙天気という分野に興味を持つ方々が増えると非常に嬉しく思います。なお、今回の記事の開発中のコードはGitHubに上げておりますので、興味のある方はご覧ください(ただし、コードの解説などはとくにないことをご了承ください)。
　
　最終的には、いくつかの指数を予測できるようにして、そのような指数を予測するシステムとしてウェブ上に公開できるようにしたいと考えています。

クレジット

The geomagnetic field data used in this article was provided by the WDC for Geomagnetism, Kyoto (http://wdc.kugi.kyoto-u.ac.jp/wdc/Sec3.html).