天気予報とデータサイエンス

サマリー

岩波データサイエンスシリーズ「時系列解析−状態空間モデル・因果解析・ビジネス応用」より。前回は機械学習のシミュレーションとの融合技術として、データ同化とエミュレーションについて述べた。今回は天気予報とデータサイエンスについて述べる。

天気予報とシミュレーション

気象庁などの機関が天気予報を作る際には、いくつかの方法がある。そのうち主要なものの一つは現在の雨の分布から外装によって将来の雨の分布を予測する運動学的な方法となる。現在では気象レーダーは日本全域をカバーしており、気象庁のサイトで降雨の状況をほぼリアルタイムで把握できる。そこで、時間方向の差分から雨雲の動きを判断し、この雨雲が今後どこに移動していくかを予測できる。

気象庁が短時間予測としてウェブサイトで公開している予測はこのようにして作られている。一方、運動学的な予報が有効なのは、数時間以内のごく短時間の予報に限られている。というのも、雨雲を流す風自身がさまざまな物理的プロセスによって時々刻々と変化していくからである。そのため、数時間以上先の予報は、大気の流れを記述する流体力学の基礎方程式などの偏微分方程式系を離散化した物理モデルを構築し、数値計算によって近似解を得ることで実施されている。ここではこれを運動学的な方法と区別して物理学的な手法と呼ぶ。

物理学的手法では、物理法則を模した数値シミュレーションを行う。ただし、数値シミュレーションを開始するにあたって、現在の大気の状態を計算の初期値としてシステムに代入する必要がある。この初期値を作成するために、まず、観測データを収集することから始める。世界各国の機関は、地上観測・レーダー観測・船舶による観測・航空観測・衛星観測といった様々な観測データはGTS(世界希少通信回線)と呼ばれる回線を通じて迅速に共有される枠組みが出来上がっていて、気象庁などの機関は、これらの観測の中から怪しいものを品質管理によって取り除く作業をしている。

さらに、観測データを初期値とするだけではもったいないので、データ同化というプロセスを噛ませる。データどうかは、観測値やシミュレーション結果に見られる物理量の誤差の分散や共分散を利用し、観測点ではない場所にも観測データの情報を波及させるシステムとなる。データ同化プロセスは人工知能の一種とみなすことができ、機動的観測の基礎理論にも関わり、このプロセスにより初期値が生成され、数値シミュレーションを始める条件が整う。

数値シミュレーションに用いる物理モデルは、ニュートンの運動法則・質量保存・期待の状態方程式・エネルギー保存則・相変化などの物理法則を表現するものとなる。これらの物理法則は高い精度で成り立つため、一般的には統計的に構築されたモデルよりも正確な結果をもたらすことが多いとされている。

ただし、物理モデルは紙と鉛筆で解ける方程式系ではないため、偏微分を空間・時間方向への差分で四則演算に置き換え、近似的に数値計算を行う必要がある。すなわち、記述する方程式を空間・時間方向に離散化した巨大な連立方程式を作るものとなる。これを概念的に以下のように書く。

\[\mathbf{x}_t=M(\mathbf{x}_0)\]

ここで下付きの添字は時刻を表しており、Mは物理法則に基づく非線形を含んだ時間発展演算子となる。xは離散化された地点ごとの温度や風速がすべて収められており、状態変数ベクトルと呼ばれる非常に長いベクトルになっている。ここでは、x₀が初期値、x_tが物理モデルによる予測結果となる。

下図に物理モデルを格子に分割した際のイメージを示す。気象庁では、地球全体の大気の状態を予測する全球モデルGSMや日本付近の状態を細かく分割する領域モデルMSMを用いて日々の天気予報をおこなっている。

待機は水平方向につながっているので、全球モデルの場合は、地球上を約20km相当の水平格子間隔、鉛直方向には乳状から上空約80kmまでを100層に分割し計算をおこなっている。格子点の下図は約1億3000万で、各格子点に4つずつの物理量が配置されているため、予報すべき変数、すなわち、状態変数ベクトルxのナカザは約5億になる。

物理法則は高い精度で成り立つが、粗い空間離散化を適用しているので、完璧なものにはならない。というのも、格子点間隔より細かいスケールの現象を陽に再現できないからである。例えば、雲は大きさが数百mなので、GSMやMSMでは雲の一つ一つを改造することはできない。しかも、天気予報を出すためには、タイムリミットがあるため、安易に格子点を増やすということもできない。

さらに現実的な層変化のモデル化は大変難しいため(たとえば、現実には気温が0℃以下になっても水はすぐに氷になるわけではない)、大胆な矜持が用いられる。このような不完全を補うためには、高度な計算機を導入して格子点間隔を小さくするとともに、物理法則の理解を深めていくことが重要となるが、まだまだ現実世界を表現するには道半ばといったところとなる。

数値シミュレーションによって得られる予報値は、先述した解像度不足や物理モデルの不完全性のため、現実世界とは少し異なったものとなる。そのため、予報値に対し、過去の実績に基づいてニューラルネットワークやカルマンフィルタなどを用いた統計的補正をかける。この作業をガイダンスと呼ぶ。最終的には、ガイダンスの結果を予報官が見て最終判断をし、説明の文をつけて一般に向けて発表することとなる。

まとめると天気予報は(a)観測データを入手する、(b)データ同化によって観測データを最大限有効活用して初期値を生成する、(c)初期値を物理モデルに与えて予測シミュレーションを行う、(d)(c)で得られた結果をガイダンスによって補正し、予報官が確認する、という手順を経て届けられるということになる。

データ同化

ここで、データ同化のプロセスについて詳しく述べる。物理モデルを用いて良い予報を出すには、得られた観測データを利用して、もっともらしい現在の状態を推定し、これを予報の初期値として与える必要がある。

これは天気予報を偏微分方程式の初期問題と捉えることに由来するものとなる。初期値を生成するには観測データを用いるが、観測データの情報だけを使うのではなく、データ同化と呼ばれる仕組みで物理モデルの結果と組み合わせて初期値を生成する。それは、観測データの数が不足しており、観測データ自身も誤差を持っているからとなる。

観測データがすべての状態変数を復元するのに十分なほど得られているのであれば、最小二乗法などによって初期値を決めることも可能だが、全球モデルにおいては、約5億という膨大な状態変数の下図に対し、手に入る観測データの下図は数百万個のオーダーとなる。そのため、観測データだけから、直接的に状態変数を推定しようとしても、粗っぽい近似になってしまう。

特に海上や上空、顕著現象周辺域のデータは衛星などのリモートセンシングに限られることがほとんどになるので、観測を密にすることは不可能となる。数値計算を開始するにあたっては、全格子点のすべての部律令に対して初期値を与えなければならない。そこで、過去に出した予報値をまずは第一推定値(事前情報)とし、観測点に近いところでは観測値に近づけるとともに、観測点から少し離れたところでは物理量どうしの共分散を利用して適切な修正をかけることを考える。さらに、観測データが不確実性の高いものであれば、観測値よりも過去に出した予報値に近づける必要がある。このようなことを統合的に扱うプラットフォームがデータ同化システムとなる。

データ同化では観測データと物理モデルから得られる共分散の情報を組み合わせて使うが、その例として、台風の中心気圧だけが観測値として得られたと考えてみる。

もし、観測された中心気圧が過去に出した予報値よりも低ければ、データ同化システムは中心気圧を下げるように初期値を修正する。しかし、容易に想像がつくように、中心気圧が低い時には周囲の風速も強いはずである。なので、中心気圧を下げる補正をかける場合には周囲の風速を強めるようにする補正も同時にかける。

このことは数学的に言えば、中心気圧と風速の誤差の共分散は負であるので、中心気圧の観測値から風速も補正できた、ということになる。台風の中で中心気圧と風速の間に負の相関があることは多くの人が知っているが、実際にはさまざまな物理量同士の関係性を使う。これをコンピューターに計算させるというのが、データ同化システムの肝となる。

ここからデータ同化を数学的に解説する。仮に、時刻tにおいて観測データが収集できたとして、時刻tを初期時刻とする天気予報を行うための初期値xtの最適推定値を得ることを目指す。まずは、さきほどの離散化された時間発展方程式を用いて、誤差の時間発展を記述する。

\[\Delta\mathbf{x}_i=M(\mathbf{x}_0+\Delta\mathbf{x}_0)-M(\mathbf{x}_0)=M\Delta\mathbf{x}_0\]

ここで、Δxは真値からの　差として定義される(不可知な)誤差、M=∂M/∂xは物理モデルMのヤコビ行列を表す。また、最後の部分では、Δxの2次以上の項がΔxの1次の項に比べて十分に小さいとしてテイラー展開の1次の近似を用いている。また誤差の共分散を表す行列とその時間発展は以下のように書くことができる。

\begin{eqnarray}\mathbf{B}_0&=&<\Delta\mathbf{x}_0\Delta\mathbf{x}_0^T>\\\mathbf{B}_t&=&<\Delta\mathbf{x}_t\Delta\mathbf{x}_t^T>=<\mathbf{M}\Delta\mathbf{x}_0\Delta\mathbf{x}_0^T\mathbf{M}^T>=\mathbf{MB}_0\mathbf{m}^T\end{eqnarray}

< >は期待値、上付きのTは転置を表し、B₀は状態変数の誤差共分散であることを表す。これらの式を導入するのは、Δxは確率的に生成されるものであり、我々はせいぜいその共分散B₀しか推定できないという想定を暗にふくんでいる。また、この式は状態変数の誤差共分散が物理モデルの特徴を反映しながら時々刻々と変化していくことを表している。

4次元変分法

予報誤差共分散の時間変化の計算は、さきほどの中心気圧と風速の関係で見たように、観測データがまばらな海上や上空の状態をより正確に復元し、正確な予報を行うために重要となる。しかし、ここで一つ問題が出てくる。要素変数が約5億個あるとすると、予報誤差共分散B₀やB_tの要素数は約25京個になるというものである。いかに計算性能の良いスーパーコンピューターといえども、このように大規模な系の共分散を陽に予報することは不可能である。そこで、気象庁では4次元変分法と呼ばれるアルゴリズムをデータ同化手法として実装し、計算の高速化を図っている。

簡単のために、誤差の時間発展が線形的であり、観測データが時刻tでのみ得られたと仮定する。このとき、4次元変分法では、以下のように状態変数の初期値ベクトルを独立変数とする評価関数を定義し、最適化問題を解く。

\begin{eqnarray}J(\mathbf{x}_0)&=&\frac{1}{2}(\mathbf{x}_0-\mathbf{x}_{0,fg})^T\mathbf{B}_0^{-1}(\mathbf{x}_0-\mathbf{x}_{0,fg})\\& &+\frac{1}{2}(\mathbf{H}_tM(\mathbf{x}_0)-\mathbf{y}_t)^T\mathbf{R}_t^{-1}(\mathbf{H}_tM(\mathbf{x}_0)-\mathbf{y}_t)\end{eqnarray}

ここで、x_0,fgは過去に出した予報に基づく状態変数で第一推定値と呼ばれている。そして、時刻tにおける観測値をまとめてベクトルにしたものをy_t、観測値の誤差共分散行列をR_t、物理モデルの格子点から観測値が得られた地点への内装などを表す行列をH_tとしている。

この式では、第二項が主要項であり観測値と数値モデルで用いる状態変数の初期値との差異を定量化している。第一項は過去に出した予測値から離れすぎてはいけないという拘束を表している。これは、ベイズ統計学的に見ると事前分布を与えることに相当している。気象学的な観点でいうと、いくら観測データと状態変数がフィットするからといって、観測がない地点で地上気温が100℃になったらダメでしょ、という拘束をかけていることになる。

この評価関数の勾配はM=∂M/∂x₀であることを使って以下の式のように計算できる。

\[\frac{\partial J(\mathbf{x}_0)}{\partial \mathbf{x}_0}=\mathbf{B}_0^{-1}(\mathbf{x}_0-\mathbf{x}_{0,tg})+\mathbf{M}^T\mathbf{H}_t^T\mathbf{R}_t^{-1}(\mathbf{H}_t\mathbf{M}(\mathbf{x}_t)-\mathbf{y}_t)\]

つまり転置行列M^Tを用いて、この最適化問題を反復計算で解くことにより、予報計算の数十倍程度の計算時間で状態変数の初期値が決定される。ここには、巨大な行列の逆行列を求めるという演算が含まれるが、実質上は\(\mathbf{x}_0=\mathbf{B}_0^{1/2}\mathbf{v}_0\)という変数変換によって定義されるv₀を独立変数とすることで逆行列の計算を回避したり、特異値と特異ベクトルを求めて一般化逆行列に置き換えたりする。

ここで∂J(x₀)/∂x₀=0となるx0をx_0,optと置き、誤差が正規分布をしているとして、最適化された時刻tにおける状態変数をx_t,optと書くこととすると、以下の式のようになる。

\[\mathbf{x}_{t,opt}=\mathbf{Mx}_{0,opt}mathbf{x}_{t,fg}+\mathbf{MB}_0\mathbf{m}^T\mathbf{H}_t^T(\mathbf{R}+\mathbf{HMB}_0\mathbf{M}^T\mathbf{H}^T)^{-1}(\mathbf{y}_t-\mathbf{H}_t\mathbf{x}_{t,fg})\]

この式は、4次元変分法がx_t,fgをアップデートするために、観測値と物理モデルから作られた第一推定値の残差y_t–H_tx_t,fgとB₀の時間発展MB₀M^Tを計算し、それを観測誤差と比較する(RとMB₀M^Tの大きさに応じて修正幅が決まる)という作業をしていることを表す、B₀の時間発展であるB_t=MB₀M^Tは中心気圧と風速といった物理量どうしの共分散を表すので、これが観測データから得られる情報を波及させることに役立つ。

4次元変分法の原理は、統計的推定に深く根ざしている。データどうかでは物理法則の矜持式を拘束条件として、各格子点における気温・風速・湿度といった状態を日々最適化する問題を解いている。なので、日々の天気予報では誰もが知らず内にデータサイエンスの恩恵にあずかっていることになる。

気象学の分野でも、多層ニューラルネットワークが現在のシステムに取って代わると考える人が増えてはいるが、ガイダンスの高度化が進むことはあったとしても、多層ニューラルネットワークが既存のデータ同化システムを超えることは難しい。

というのも、現状において、すでにニュートンの運動法則や質量保存といった精度の良い物理の基礎法則を拘束条件として、観測データとの差異を最小するような問題を解いているからである。逆にいうと、気象学で現代的なデータサイエンスのアプローチによる利得を最大限に得たいならば、物理モデルが苦手としているところ、たとえば、雲物理であるとか、乱流であるとか、数値モデルが苦手としている現象を狙うべきである。

次回はタンパク質の3次元形状(ミスフォールディング)のシミュレーションと機械学習による解析について述べる。