機械学習･深層学習ミーツ
|天気予報|🌏🌎🌍
グーグル･ディープマインドは、 自社開発した機械学習ベースの気象予測モデル(MLWP)GraphCastが、数値計算ベースの気象予測モデル(NWP)を上回ったと発表した(23年11月14日@同社公式ブログ[*1]。論文[*2]は、サイエンスにて公開されたのが23年11月14日(arXivに第1版が投稿されたのは22年12月24日、第2版は23年8月4日)。比較対象となったNWPは、世界最高レベルと謳われる、欧州中期予報センター(ECMWF)が開発･運用している高解像度予報(High RESolution Forecast：HRES)である。HRESは、緯度/経度 0.1° の解像度で、全球の 10 日間の予報を約1時間で生成する。
GraphCastは、比較対象としてHRES とともに、Pangu-Weatherも意識している。Pangu-Weatherとは、中国(通信機器メーカー大手)ファーウェイ(の子会社ファーウェイ･クラウド)が開発したMLWPである。23年7月20日にnatureで公開された論文[*3](arXivには、22年11月3日投稿)において、HRESよりも優れているとアピールしている。Pangu-Weatherは、FourCastNetを意識して、開発されている。FourCastNetは、NVIDIA、ローレンス･バークレー国立研究所等が開発した、MLWPであり、成果をまとめた論文[*4]は22年2月にarXivに投稿されている。
この一大サーガを、以下Ⅰ-1～Ⅰ-3に整理する。その前に、Ⅰ-0において予め、いくつか整理をしておく。
📖追補📖
新たに、マイクロソフトが開発したモデルAuroraを追加した。Auroraは、NWP及びGraphCastより優秀である！と主張する。
24年7月から、Googleナウキャストが利用可能と報道された(24年6月19日)ので、Met-Net3を追加した。
NeuralGCMを追加した。
GenCastを追加した。
【1】機械学習ベースの気象予測
(0) 注意喚起の為念
｢機械学習ベース(の気象予測)｣という言葉について。英語でもMachine Learning-basedという形容句が付いているので、日本語でも、機械学習ベースとして良いと思われる。しかし、この用語は、誤解を生じさせる可能性があり、深層学習ベースと訳した方が、良いと思われる。というのも、機械学習ベースの気象予測だと、サポートベクターマシンなどの｢機械学習ベース｣の手法を使って、ダウンスケーリングを行った研究が存在するからである。ダウンスケーリングは、元の予報が行われた空間スケールよりも細かい空間スケールにおいて、予報を行うことである。予報時点(時刻)は、同じである。
深層学習ベースの予報は、過去データから学習して、未来の予報を行う。FourCastNet、Pangu-Weather、GraphCastは、どれも深層学習ベースの予報モデルである。実際、[*4]では、畳み込みニューラルネットワークを使った気象予測モデルに対して、DLベース(の気象予測)という表現が使われている。とは言え、平仄を合わせるために、以後(も)、機械学習ベース(MLWP)という言葉を使用する。
(1) 3モデルのアーキテクチャ概要
0⃣ 従前には、深層学習ベース気象予測モデルのアーキテクチャとして、ConvLSTMやオートエンコーダー(自己符号化器)が存在していた[*13]。ConvLSTMは、空間方向における気象データの相関を畳み込みニューラルネットワークで捉え、時間方向の相関をLSTM(Long Short Term Memory；再帰型ニューラルネットワーク)で捉える、という建付けである。トランスフォーマーは、オートエンコーダーの1種である。GraphCastのモデル･アーキテクチャもオートエンコーダーである。
1⃣ FourCastNet は、Fourier ForeCasting Neural Networkの略である(と[*4]に明記されている)。Fourierは、フーリエ･ニューラル演算子を意味している。
2⃣ Pangu-Weather には、ざっくり言うとトランスフォーマーが使われている。少し正確に言うと、ビジョン･トランスフォーマーであり、もう少し詳しく言えば、Swinトランスフォーマーが使われている。Pangu(盤古)とは、中国神話における天地創造の神のことらしい(ファーウェイの大規模言語モデルのPangu、という名称が付けられている)。Pangu-Weatherのパラメータ数は、およそ6,400万である。
3⃣ GraphCastのモデル･アーキテクチャは、名前に仄めかされている通り、｢グラフ･オートエンコーダー｣という種類のグラフ･ニューラルネットワーク(GNN)である。パラメータ数は、およそ3,670万で、比較的少ない。
上記3モデルは、それぞれ深層学習の作法で、｢広域に及ぶ影響を取り込んでいる｣。これは、物理的(気象学的)には、｢非局所的な相互作用を取り込む｣ことに相当している。具体的には、ビジョン･トランスフォーマーの自己注意機構、GNNのメッセージ･パッシングを利用して、｢非局所的な相互作用を取り込んでいる｣。
(2) 気象学的物理量(予報変数)
0⃣ 共通
気象学では、東西方向の風をUと呼ぶ(らしい)。同じく、 南北方向の風をVと呼ぶ(らしい)。 ジオポテンシャルとは、｢海水面からある高度まで単位質量当たりの空気塊を上昇させるのに必要な仕事量｣である。比湿とは、｢湿潤空気中の単位質量当たりの水蒸気の質量｣である。
高層大気の物理量(上空変数)は、特定の等圧面上における値が取り扱われる。特定の等圧面上における値とは、例えば、U500と表記される値である。具体的にU500とは、｢気圧500hPa(ヘクトパスカル)の等圧面上における東西風｣を意味する。また、等圧面は、海面からの高度に相当する。例えば、500hPaは海面から上空5,500mの高度に、850hPaは1,500mの高度に相当する。500と850hPaは特に、代表的な等圧面(高度)のようである。
1⃣ FourCastNet
地表面の物理量(地表変数)として、以下を採用：風速10mのU成分U10、 風速10mのV成分V10。地表から2mの気温T2m(地上気温)、地表気圧sp、平均海面気圧mslp。上空変数として、以下を採用：ジオポテンシャルZ、比湿RH、気温T、U、V。等圧面は、50、500、850、1000hPaの4つ。50hPaではZのみ。1000hPaでは、U、V、Zのみ。500と850hPaでは全て。つまり、{Z、RH、T、U、V}。
さらに、(日平均)鉛直積算水蒸気量TCWV(Total Column Water Vapor)が採用されている。TCWVは、地表面を基点として、その上空を鉛直の柱と考えたとき、その柱に含まれる水蒸気の総量を意味する。つまり降水となり得る、水蒸気、水、氷等全てを含む。
2⃣ Pangu-Weather
地表変数として、以下を採用：平均海面気圧MSLP、地上気温T2、U10、V10。上空変数として、以下を採用：Z、比湿Q、T、U、V。等圧面は、50, 100, 150, 200, 250, 300, 400, 500, 600, 700, 850, 925, 1000hPaの13。
3⃣ GraphCast
地表変数として、以下を採用：平均海面気圧MSL、地上気温2T、風速10mのU成分10U、風速10mのV成分10V、総降水量TP。上空変数として、以下を採用：Z、Q、T、U、V、鉛直風速W。等圧面は、Pangu-Weatherと同じ。
FourCastNet、Pangu-Weather及びGraphCastでは、ERA5を入力(学習データ)として受け取り、ERA5を出力するように学習される。もちろん、入力と出力の間には、時間差がある。この時間差を、(天気予報の文脈では)｢リードタイム｣と呼ぶ。グランド･トルゥースも、ERA5である。
時間間隔は1時間、格子間隔は0.25°で、全球モデルの再解析データの内、最も解像度･時間間隔が細かくなっている。期間は、1979年～となっている。なお、鉛直層は37等圧面(1,2,3,5,7,10,20,30,50,70,100,125,150,175,200,225,250,300,350,400,450,500,550,600,650,700,750,775,800,825,850,875,900,925,950,975,1000hPa)である。
2⃣ 熱帯低気圧の進路予測における｢グランド･トルゥース｣は、 IBTrACS(International Best Track Archive for Climate Stewardship)である[*7]。IBTrACSとは、米国大気海洋庁(NOAA)が公開している世界の熱帯低気圧のデータである。Best Track:ベストトラック(データ)とは、事後解析により最も現実に近いとされるデータのことである。
3⃣ TIGGE(THORPEX Interactive Grand Global Ensemble)データベースは、THORPEXプロジェクトで整備された、中期アンサンブル予報データのデータベースである。THORPEX(THe Observing system Research and Predictability Experiment)は、世界気象機関(WMO)による、1～2週間先の気象予測精度向上を目指したプロジェクトである。世界の数値気象予報機関10機関(気象庁、欧州中期予報センター、米国立環境予測センターなど)のデータが利用できる[*8]。
【2】言葉の整理、予報誤差を表す指標など
気象予測･天気予報の世界では、予報精度等に独特の指標が存在するので、整理した。
⓪ 天気、気象、気候等について、正確な区別は行わない。
① 中期天気予報とは、最大 10 日先までの大気変動の予測を意味する。
② ターゲットという文言は、気象予測･天気予報の文脈では、"評価対象変数"の意味で用いられる。気象予測･天気予報における変数は、やや煩雑である。具体的には、まず、地表変数と上空変数がある。上空変数は、等圧面毎に存在する。さらに、予測期間に応じて、それらの変数は増えていくことになる。例をあげると、地表変数を4つ、上空変数を5つ、等圧面を13とした場合、それだけで、変数は4+5×13＝69個になる。中期予報を考えると、10日間で12時間ごと(つまり1日に2回)の更新なので、変数は20倍になる。つまり、69×20＝1,380変数となる。この1,380個の評価対象変数を｢ターゲット｣と呼ぶ。
1⃣ サイクロン、台風、ハリケーンについて[*9]。
サイクロンとは、低気圧を指す一般的な用語である。台風もそれ以外の低気圧も、すべてサイクロンと呼ぶ。トロピカル･サイクロンとは、熱帯(性)低気圧を指す。ただし、トロピカル･サイクロンの略称として、サイクロンが使われることがある。なお、｢熱帯｣とは、発生する場所を指す用語ではなく、サイクロンの構造を指す用語である。台風やハリケーンなどは、すべて｢強い｣トロピカル・サイクロンである。
台風とは、北西太平洋に位置する｢強い｣トロピカル・サイクロンを指す。ハリケーンとは、北部大西洋、東部北太平洋、中部北太平洋および南東太平洋に位置する｢強い｣トロピカル・サイクロンを指す。以下Ⅰ-1～Ⅰ-3では、｢強い｣トロピカル・サイクロンの総称として、熱帯低気圧という呼称を用いている。
2⃣ 台風の上陸、接近、通過について[*10]。
◇日本における｢上陸｣とは☛ 台風の中心が、北海道・本州・四国・九州の海岸に達した場合を言う。
◇日本における｢接近｣とは☛ ㊀ある地点への台風の接近：台風の中心が、その地点を中心とする半径300km以内の域内に入ること。㊁ある広がりをもった地域(地方予報区など)への台風の接近：台風が、その地域の地理的な境界線(海岸線、県境など)から半径300km以内の域内に入ること。
◇日本における｢通過｣とは☛ 台風の中心が、小さい島や小さい半島を横切って、短時間で再び海上に出る場合を言う。
3⃣ アンサンブル予報について[*11]
複数の予測の集合をアンサンブル、個々の予測をメンバーと呼ぶ。摂動を加えているメンバーを｢摂動ラン(あるいは摂動予測)｣、摂動を加えていないメンバーを｢コントロールラン(あるいは制御予測)｣と呼ぶ。
4⃣ 欧州中期予報センター(ECMWF)。そして、IFSとHRES
⓪ ECMWF及び、その数値予報モデル概要
ECMWFは、欧州35カ国が支援する政府間組織で、2025年に設立50周年を迎えた。数値気象予報の世界的リーダーとして広く認められている。ECMWFの数値予報モデルは数百万行に及ぶ巨大な計算プログラムであり、｢大気と、大気と海洋、陸地、海氷との相互作用を支配する物理法則｣に関するあらゆる知識を網羅している。この数値モデルによって生成される予報は世界中の、｢気象観測所、海洋ブイ、衛星からの観測データ｣に基づいて定期的に調整される。全観測データの95%以上は衛星データで占められており、欧州宇宙機関(ESA)との緊密なパートナーシップを維持している。出典：https://www.nature.com/articles/s41467-025-65837-2
ECMWFのアンサンブル予報は、50の予測計算を並列に実行している。出典：ibid
① IFS
ECMWFが開発し、現業予報で使用している全球気象予測モデルがIFS(the Integrated Forecasting System)である。IFSの力学コア🌏は、静水圧近似(静力学平衡とも呼ばれる)🌏、予報時間の間隔が2時間、半陰解法、セミ･ラグランジュ法🌏である。物理過程🌏は格子点法🌏で行い(つまりパラメタリゼーション)、力学過程🌏はスペクトル法🌏を用いる。垂直方向では、モデルは有限要素法を用いて離散化される。水平方向では、適合ガウス格子🌏が使用される。
出所：https://www.ecmwf.int/en/research/modelling-and-prediction/atmospheric-dynamics
② HRES
高分解能のIFSがHRES(Highest-RESolution)。水平解像度0.1°、鉛直解像度137層。スケジュールはやや複雑だが、例えば0～90時間であれば、｢00UTC、06UTC、12UTC、18UTC｣全てで1時間ごとに実行される。出所：https://www.ecmwf.int/en/forecasts/datasets/set-i
🌏 こちらの各項を参照。もしくは、Ctrl＋Fで検索して頂くか･･･
③ 為参考：Destination Earth(DestinE)
DestinEは、欧州の環境プログラムをデジタルで実現するためのイニシアチブ。地球の高精度なデジタル･レプリカ(デジタルツイン)を開発し、自然現象、災害、そしてそれらに影響を及ぼす関連する人間活動をモデル化およびシミュレートすることを目的としている。ECMWFは、2つの優先度の高いデジタル･ツイン(気候変動適応🐾と極端気象に関するデジタル･ツイン)と、デジタルツインの実現に必要なソフトウェアおよびデータ環境であるデジタル･ツイン･エンジンの実装を担っている。
🐾 気候変動適応とは、気候変動を所与とした上で、その影響を軽減しつつ、柔軟に対応する取り組みと説明される。加えて、気候変動を逆手にとり、機会と捉えること(蜜柑の代わりにアボカドを栽培する等)も気候変動適応とされる。
(1) 予報誤差を表す基本的な指標[*11]
予報誤差を表す基本的な指標として平均誤差(Mean Error：ME、バイアスとも呼ぶ)並びに、自乗平均平方根誤差(Root Mean Square Error：RMSE)がある。ME は予測値の実況値からの偏りの平均であり、0 に近いほど実況からのずれが小さいことを示す。RMSE は最小値の 0 に近いほど予測が実況に近いことを示す。実況とは、特定時刻の観測データを意味する。通常は、最新時刻のデータを意味し、最新時刻を持ってリアルデータと見做している。
さらに、アノマリー相関係数(ACC)がある。これは、予測値の基準値からの偏差(気象学的には大胆にも、これをアノマリーと呼んでいる)と、実況値の基準値からの偏差との相関係数である。基準値としては、気候値を用いる場合が多い。気候値(平年値とも呼ばれる)は長期平均値のことであり、通常は、30年間の平均値である。
FourCastNet、Pangu-Weather及び、GraphCastでも、NWPとの比較にRMSEと私たちの数値モデルは、数百万行に及ぶコードを持つ巨大な計算プログラムであり、大気と、大気と海洋、陸地、海氷との相互作用を支配する物理法則に関するあらゆる知識を網羅しています。このモデルは、スーパーコンピュータ（現在、1秒あたり30兆兆回の計算能力を持つ）を用いて、1日に数回計算を行っています。が使われている。
(2) スキルスコア[*11]
スキルスコアは、気候学的確率などによる予測の難易を取り除いて、予測の技術力を評価する指数である。代表的なスキルスコアは 、Heidke のスキルスコア(HSS)で、以下の式で求められる適中率である。適中率とは、(文字通り)予測が適中した割合である。
HSS＝(FO＋XXーS)/(NーS)
S＝Pc(FO + FX) + Px(XO + XX)
Px=X/N＝(FX＋XX)/(FO＋XX＋FX＋XO)
(”適中”と表現される)FOは、｢あり(例えば、雨が降る)｣と予測して、実際に｢あり(例えば、雨が降った)｣だった事例の数。(同じく、”適中”と表現される)XXは、｢なし｣と予測して、実際に｢なかった｣事例の数。(”空振り”と表現される)FXは、｢あり｣と予測して、実際は｢なかった｣事例の数。(”見逃し”と表現される)XOは、｢なし｣と予測して、実際には｢あった｣事例の数。
Pcは、｢気候学的出現率｣と呼ばれ、M/N＝(FO＋XO)/(FO＋XX＋FX＋XO)で表される。「現象あり」の平均的な出現確率であり、この量は実況のみから決まり、予測の精度にはよらない。Pxは、「現象なし」の平均的な出現確率であり、１ーPcである。
スキルスコアには、他にも、ギルバート･スキルスコア(エクイタブル･スレットスコア)やブライア･スキルスコア、フラクション･スキルスコア、ROC 面積スキルスコアなどがある。
【3】感想など・・・
(0) 数値天気予報で、日本は世界でも進んでいたという印象があった。しかしMLWPモデルでは、プレゼンスがないように見える。
➡ 新聞報道(25年5月14日付け日経)によると、気象庁は、機械学習(深層学習)ベースの気象予測の導入を検討しているらしい。モデル･アーキテクチャは、グラフ･トランスフォーマーか？
❚追記1❚
MetNet-3ベースの雨量予測モデルを24年7月から日本で運用している株式会社ウェザーニューズは、｢2kmメッシュ/15分毎の細かさで、欧州全域の気温、気圧、湿度、降水量、降雪量、風向、風速、日射量などの予測データを提供｣することを発表した(24年7月17日)[*35]。
❚追記3❚
資料↓を読む限り、MLWPモデルでプレゼンスを発揮するのは、難しい気がする。
(資料) 日本気象予報士会東京支部、機械学習GraphCastと気象予報の将来、2024年8月18日(第75回東京支部例会)、https://www.yoho.jp/shibu/tokyo/2408tange&yamazaki.pdf
❚追記4❚
ECMWFもMLWPモデル｢AIFS｣を開発している(論文は[*47]で、24年8月7日付け)。アーキテクチャは、エンコーダー・プロセッサ・デコーダー。プロセッサは、GNNベースのトランスフォーマー。
👉 AIFSは、｢熱帯低気圧の進路を含む多くの指標において、最先端の物理ベース･モデル(数値計算ベースの気象予測モデル)を最大20%上回る性能を発揮した｣と、ECMWFが報告している。予報1回あたりのエネルギー消費量は、約1/1,000に削減した、とする。
(1) 予報変数の予測で、MLWP ＞NWPという議論は、あまり意味がない。MLWPの入力データとグランド･トルゥースが同じであり、MLWPの精度が上がって出力とグランド･トルゥースの差が小さくなった結果、MLWP ＞NWPになっただけである。ここでのポイントは、NWPの入力とMLWPの入力は異なるということである。
(2) 現時点の評価は、｢中期予報の範囲において予測できる程度に、気象現象の力学過程及び物理過程を、MLWPがやっと表現できるようになった｣であろう。MLWP ＞NWPの議論が意味を持つのは、同じ入力(観測データ)を用いた上で、MLWP ＞NWPとなった場合であろう。
(3) 極端気象現象の高精度予測
それとは別に、MLWPには、大きな期待がかけられている。それは、｢アンサンブル予報による極端気象現象の高精度予測が、迅速に行える｣期待である。ここで言う極端気象現象とは、熱帯低気圧(≃台風)や、(線状降水帯発生に伴うゲリラ的)集中豪雨等を指している。
熱帯低気圧の進路予測や、集中豪雨の予測に、アンサンブル予報が有効であることが知られているが、それはアンサンブル･メンバー数が多い事が必須条件である。ただし、アンサンブル･メンバー数を増やすと、それだけ計算時間を要する。(事前学習済の)MLWPは、NWPに比べて桁違いに高速なので、一刻を争う極端気象現象の予測におけるMLWPの期待は、大きいと考えられる。
熱帯低気圧(台風)の数値予報で言えば、進路予測の精度は確実に向上していると、広く認識されている。一方で、強度予測は、精度の向上が足踏みしていると目されている。その理由の一つとして、雲モデル(積雲対流スキーム)があげられていた。つまり、雲のモデル化が不十分であるため、と考えられていた。
¦為念－FYI¦ 24年8月発生した台風10号の進路予測において"結果的に"、気象庁モデルの精度は、ECMWFやJTWC🌀の精度に劣ったとされる。実際の進路は、気象庁が予測した進路より西側(大陸側)に寄っていた。もっとも、台風経路は、｢太平洋10年規模変動(気候システムの自然変動)と同期しながら、今世紀に入って大陸側にシフトしている｣との指摘がある[*48]。従って、24年台風10号の進路予測が外れた理由は、構造的な要因に由来するかもしれない。
🌀 JTWC(Joint Typhoon Warning Center：米軍合同台風警報センター)は、米国海軍と米国空軍が、ハワイ州真珠湾海軍基地に合同で設置した、米国防総省の機関。JTWCは、あくまで軍事作戦の一環として、台風警報を発令するため、非公開となることもある。ちなみに、｢JTWCは軌道予測にどのような数値モデルを使用していますか？｣という質問に対する回答は、削除されている(https://www.metoc.navy.mil/jtwc/jtwc.html?faq#models)
〖参考〗JAMSTEC(海洋研究開発機構)は、新たな雲モデルを開発し、｢熱帯低気圧の(強度を含む)再現性を向上させた｣と発表(24年3月15日)[*28]。新しい雲モデルは、｢統計的に扱うことで、計算格子内に異なる種類の複数の雲が存在｣しているようにモデル化する。再現性が向上した理由は、｢熱帯低気圧中心の下層付近における、上昇を伴う雲による加熱率を、より正確に再現できた｣ため。
❚追記2❚
東京大学大気海洋研究所は、｢機械学習モデルと数値予報モデルによるハイブリッドモデルが、熱帯低気圧の予測において、従来手法を上回った｣と発表した(24年7月19日。論文は、[*36])。予測の対象は、熱帯低気圧の進路と強度。グランドトルゥースは、IBTrACS-WMO ver.4🐾0。比較対象モデルは、㊀各国該当機関🐾1の全球モデル、㊁Pangu-weather🐾2、㊂ERA5を側面境界値とする数値予報モデル(WRF🐾3モデル)🐾4、及びPanguの出力を側面境界値とする数値予報モデル🐾5。機械学習モデルの出力を数値予報モデルの入力として使用し、モデル全体の精度をあげた→進路と強度の誤差は、対㊀比では59％と32％減少。対㊁比では2％と59％減少。対㊂比では32％と23％減少。
ハイブリッドモデル🐾6における、機械学習モデルは、中国が開発した｢Pangu-Weather｣。数値予報モデル(大気モデル)はWRFモデルが使用された。なお、機械学習モデル＋数値予報モデルのハイブリットという座組は、グーグルのNueralGCMと同じである。
🐾0 米海洋大気庁(NOAA)に所属する国立環境情報センター(NCEI)が管理している熱帯低気圧ベスト･トラック･データ(の第4版)。
🐾1 中国気象局CMA、カナダ気象センターCMC、欧州中期予報センターECMWF、米環境予測センターNCEP。
🐾2 ERA5で学習したバージョンと、GFSのデータで学習したバージョンがある。GFSはNOAAによる気象予報モデル。
🐾3 Weather Research and Forecasting
🐾4 海洋モデルを結合させたバージョンが含まれる。
🐾5 (Pangu-weatherの出力を側面境界値として用いる場合は)spectrul nudgingという技法を使って側面境界値の不一致を解消したバージョンを含む。
🐾6 GDS2024世界デジタルサミットで、ヒューレット･パッカード･エンタープライズのエン･リム･ゴー氏は｢AIによる(気象)予測では、物理方程式を使わない｣と発言している。これまでは、それでtrueであるが、今後はfalseであろう。つまり、物理方程式も使って、さらに精度を高める方向に向かうだろう。
(4) (2)での議論も被ってくるが、MLWPの学習データは再解析データなので、観測データを入力データとすることは、現時点では出来ない(👇)。しかし、その解決を、力技で行うことも可能ではないか、と思われる。生成モデルを使うことで、再解析データを生成することも、できるのではないだろうか(言うは易し、ではあるが)。あるいは、ニューラル演算子を使ったベイズ的能動学習というアプローチもあるだろうか。
👉 長谷川直之(2025)†によれば、｢観測データをそのまま使って学習する純粋なAI予報モデルの研究も進められています｣(232頁)。長谷川氏は、第27代気象庁長官。
† 長谷川直之、天気予報はなぜあたるようになったのか、集英社インターナショナル、2025
(5) 生成モデルという側面でも、再帰型ニューラルネットワークという側面でも、(制約付き)ボルツマン･マシンが気象予測に用いられても良いと思えるが、適用例はほとんどないようである。
(6) また、MLWPをPhysics-Informedで補強するというアプローチは、ありだと思えるが、未だ具体的には行われていないようである。
👉 フィンランド･アールト大学の研究者によるCLIMODE: CLIMATE AND WEATHER FORECASTING WITH PHYSICS-INFORMED NEURAL ODESという論文[*46]が出た(24年4月15日@arXiv)。※知らないだけで、他にもあるかもしれない。
(7) Physics-Informedニューラルネットワーク(PINN)のソルバーであったSimNetを前身とするNVIDIA Modulusは、 GraphCast アーキテクチャも含んでいる[*14]。GraphCastはグーグルが開発したモデルで、グーグルのTPU(v4)を使って学習されている。NVIDIAは、それをも対応しているという事実は、NVIDIAがAIプラットフォーマーの”絶対王者”であることを示しているように思える。
(8) 国立研究開発法人海洋研究開発機構が、地域気候サービスのための生成AI基盤モデルを開発するらしい。具体的には、｢TCFDレポート生成モデル｣と｢自治体の災害対策立案用モデル｣を開発するらしい[*49]。
GPUベンダー(であり、AIプラットフォーマー)米NVIDIA、米ローレンス･バークレー国立研究所等[*12]は、高解像度のMLWPモデル｢ForeCasNet｣(arXivにて論文[*4](以下、本論文))を発表した(23年2月22日)。
モデル開発の背景には、従来のMLWPでは、降雨予測や地表面風速の予測が行えないという問題意識がある。降雨予測は、災害対応。風速予測は、 風力エネルギーの資源計画改善という、具体的な問題設定が、それぞれある。技術的には、従来のMLWPでは、解像度が低すぎて、対応できなかった。従って、高解像度(具体的には、0.25°＝赤道付近では約30km×30kmの空間解像度に相当する)気象予測が可能な、MLWPモデルを開発することが目標である。
物理現象として気象予測を考えたとき、｢非局所的な相互作用｣を捉えることが、本質的に重要である。ニューラルネットワーク･アーキテクチャで、それを可能にするアプローチは、いくつも考えられる。本論文では、ビジョン･トランスフォーマー(ViT)を採用している。先行研究では、畳み込みニューラルネットワーク(CNN)が採用されていた。画像認識タスクのみならず、広範な応用分野で、ViT ＞ CNNが観測されたのだから、発想としては自然であろう。
【1】本論文の主張
本論文は、次のように主張している。
1⃣ 降水量と地表面風速の予測が可能な、最初の高解像度MLWPモデルを開発した。
2⃣ NWPと同程度の精度かつ、NWPより遥かに高速なモデルを開発した。
3⃣ 数千のアンサンブルメンバーによるアンサンブル予測が可能なモデルを開発した。
なお、本論文(及び本稿)で言及されるNWPは、ECMWFの全球気象予報モデル(IFS)である。IFSは、米国立環境予測センターの全球気象予報モデルよりも、優れているとされている。
【2】事前整理
(1) ビジョン･トランスフォーマー
ビジョン･トランスフォーマー(ViT)は、長距離の依存関係を適切にモデル化できるために、優れた性能を発揮すると理解されている。そのカラクリは、トークン混合である。ViTの文脈では、｢特徴量を成分として持つベクトル｣をトークンと呼ぶ。つまり、トークン＝(特徴量1、特徴量2、･･･、特徴量n)Tである(上付き添え字Tは、転置の意味)。この例では、ベクトルの成分はn個である。この数nをチャネルの次元と呼ぶ。もっとも、特徴量1や特徴量2,･･･は、ニュアンスとしては、特徴量候補1、特徴量候補2,･･･、と呼ぶ方が正しいだろう。
具体的に述べると、長距離の依存関係を適切にモデル化するためにViTは、トークンの相互作用を計算する。算術的には、ベクトルの内積を計算して、相互作用の強さを定量化する。トークン間の相互作用を計算する仕組みは、自己注意機構(Self Attention Mechanism)と呼ばれる。
なおトランスフォーマーやViTは、豊富なデータで学習しなければ、優れた性能を発揮できないと
認識されている。もっとも、気象予測に関して、その心配はないだろう。
(2) トークン混合
自己注意機構では、トークンにグラフ構造を課し、トークン間の”グラフの類似性”を学習する、(効果的な)トークン混合が行われる。トークン混合(token mixing)は、トークンを空間方向に混ぜる｢空間混合(space mixing)｣と、次元方向に混ぜる｢チャネル混合(channel mixing)｣で構成される。ただし、空間混合をトークン混合と呼ぶこともある。
空間混合は、異なる空間位置におけるトークンの混合である(つまり、トークンの場所を変える)。 チャネル混合は、トークンの中の特徴量の混合である(つまり、ベクトルの成分の場所を変える)。混合することで、遠くにあるトークンの影響も取り込むことが、”可能な”仕組みとなっている。可能ではあるが、それが効率的であるかは、別問題である。
実際、トークン混合は、トークン数に応じて、2次関数的にスケールするため、高解像度入力では実行困難となる。そこで、効率的なトークン混合を如何に実行するか、が実用上の課題となる。
(3) 適応的フーリエ･ニューラル演算子
本論文では、ニューラル演算子を使った｢効率的な｣トークン混合を提案･実施している。正確には、適応的(adaptive)フーリエ･ニューラル演算子(AFNO)[*15]を使った空間混合が行われる。チャネル混合は、自己注意機構と同じ仕組みを採用する。なお、｢適応的｣は、正確な和訳ではない。
具体的には、トークンを関数空間の連続要素として扱い、空間混合を畳み込み積分としてモデル化する。そのために、ニューラル演算子を抜擢した。その数学的な背景として[*15]では、｢自己注意機構は、カーネル積分として記述することできる｣ことを上げている。
ニューラル演算子は、無限次元空間内の連続関数間の写像を学習する、ニューラルネットワークである。本論文では、ニューラル演算子の中でも、フーリエ･ニューラル演算子(FNO)を採用している。
FNOでは、入力データに並進不変性を課すことで、カーネル積分を(非局所的な)畳み込み積分に置き換え(制限)する。そうすることで、高速フーリエ変換の適用を可能とし、それが推論の高速化を可能とする。FNOは推論が高速であることに、最大の特色がある。ニューラル演算子及びフーリエ･ニューラル演算子の詳細は、こちらを参照。
AFNO)は、元々は、画像を取り扱えるように拡張したFNOである。FourCastNetでは画像を扱っているわけではないが、数値データ(今の場合、ERA5の再解析データ)をピクセル値と考えれば、画像アナロジーで捉えることが可能である。FNOは入出力が連続関数であるが、画像には(離散的な入力データにも)不連続性があるため、不連続性に適応する必要がある。このため、①チャネル混合重みにブロック対角構造を課す、②トークン間で重みを共有するなどの、いくつかの改修が施されている。また、AFNOでは(本末転倒にも思えるが、必然的に)、離散フーリエ変換を使用している。
(ただ･･･)結果として、AFNOの計算複雑性は、O(NlogN)に低減される。自己注意機構の計算複雑性O(N2)に比べて、ほぼ2次加速である。ここで、Nはトークンの数である。
【3】FourCastNetモデルの詳細
(0) 改めて概要
FourCastNetのモデル･アーキテクチャは、あくまでビジョントランスフォーマー(ViT)である。フーリエ･ニューラル演算子(FNO)ではない。その意味で、FourCastNetという名称は、ややミスリードである。FNOはあくまで、トークン混合を効率的に行う施策として、採用されている。効率的とは、計算複雑性を低減させる、という意味である。位置埋め込みや残差接続は、普通に行われる。
FourCastNetはViTとFNOベースのトークン混合を組み合わせたことで、｢長距離の依存関係を適切に取り込む｣こと、並びに｢(計算複雑性が低いため)高解像度に対応すること｣ができるモデルである。
(1) データセット
学習データセットは、1979 年～2015 年(37年分)のERA5データセットを使用。検証データセットは、2016 年と 2017 年の2年分。テスト データセットは、2018 年以降のデータで構成される。各変数は、721 × 1440ピクセルの 2次元場として表される。ERA5の時間分解能は1時間であるが、FourCastNetではサブサンプリングを行い、間隔を6時間としている。
ちなみに、エンドツーエンドの学習には、64 個の Nvidia A100 GPU のクラスターで約 16 時間を要した。また、FourCastNet の学習に必要なエネルギーは、NWPにおいて50個のアンサンブル･メンバーで 10 日間の予測を生成するのに必要なエネルギーとほぼ同じ、であった。
(2) ハイパーパラメータ等
① バッチサイズ → 64
② 学習率(事前学習/再学習/総降水量モデル) → 5×10−4 /1×10−4 /2.5×10−4(総降水量モデルは、別モデルである)。
③ 学習率スケジューリング → コサイン･スケジューリング(cosine decay)
④ パッチサイズ → 8×8
⑤ ドロップアウト率 → 0
⑥ フーリエ･ニューラル演算子のブロック数 → 8
⑦ フーリエニューラル演算子の埋め込み次元 → 768
⑧ 活性化関数 → GELU
(3) 総降水量モデル
総降水量(TP)の予測は、困難であることが知られている。TPの確率分布はゼロで強くピークに達し、正の値に向かって長い裾が広がる。つまり、TP は他の予報変数よりも、疎らな空間特徴を示す。さらに、TP は、大気の動的時間発展を導く変数(風、圧力、温度など)に大きな影響を与えず、NWP で正確に捕捉するには、相変化などのプロセスの複雑なパラメーター化が必要になる。
これらの理由から、本論文では、ERA5の再解析データを使って学習したモデル(以下、便宜上、基本モデルと呼ぶ)の出力を使用して、TP用のモデルを別途学習する。このTP用モデルは、基本モデルと同様にAFNOアーキテクチャを採用する。さらに、2次元畳み込み層(周期的なパディング付き)と、ReLU活性化層が追加されている。基本モデルが 6 時間単位で予測を行うため、TP用モデルも、6 時間ごとの累積総降水量を予測するように学習する。さらにTPは、対数変換される: TP’ = log(1 +TP/ε)、ε = 1 × 10−5。こうすることで、降水量ゼロの予測が妨げられ、値の分布の偏りが少なくなる。
【4】NWPとの比較結果
(1) 予報変数のNWPとの比較
比較がフェアではないという意味で、比較に意味がないので、個別評価は割愛する。代わりに、㈠アンサンブルしたFourCastNetの結果と、㈡アンサンブルしていない単独のNWP結果との比較を評価する。アンサンブルは、白色雑音を加えて初期条件を摂動することで、実行された。雑音をスケールする係数は0.3を使用している。対象とする予報変数は、Z500とU10である。経時に伴い、RMSEが増加しない、ACCが低下しない、と精度が高いと定義している。
全てにおいて、すなわちZ500のRMSEとACC、U10のRMSEとACCにおいて、ほぼ完全に㈡が㈠を凌駕している。残念ながら、FourCastNetの精度は、その程度である。
(2) 熱帯低気圧の進路予測
0⃣ 前振り
まず｢以前のモデルは、熱帯低気圧の(進路)予測を考慮するのに十分な、長い予測リードタイムを備えた風速やその他の重要な予報変数の、正確な予測を生成できなかった｣とした上で、｢FourCastNetは、かなり優れた解像度を備えており、熱帯低気圧の発生と進路を予測できる変数の、正確な中範囲予測を生成した｣と主張している。そしてケーススタディとして、2018 年に発生したハリケーン･マイケルを検討している(Pangu-Weatherでも取り上げられている)。
1⃣ 予測ロジック、評価指標等
❶熱帯低気圧の進路予測と、❷中心気圧の低下(熱帯低気圧の発達)を予測している。
❶ 平均海面気圧の極小値を”台風の目”として、”台風の目”を追跡することで、熱帯低気圧の進路を予測した。アンサンブル･メンバー数が100とやや少ない、FourCastNetのアンサンブル予報が出力した平均海面気圧極小値の平均位置と、グランド･トルゥースを比較する(グランド･トルゥースは、ERA5の再解析データである)。精度は、やや物足りないように思えるが、グランド･トルゥースは緯度経度に関する90%ileの幅には収まっている。
❷ 本論文でも、｢(発生から)36 時間後に始まり48 時間後の時点に及ぶ、気圧の急激な低下を完全には予測できていない｣と述べられているように、熱帯低気圧の強さの予測は精度が低い。この理由として、｢対流および放射プロセスを考慮していない｣ことを推測している。
(3) 大気河川の予測
FourCastNetによって推論される(日平均)鉛直積算水蒸気量と、グランド･トルゥース(ERA5の再解析データ)を比較している。本論文では、｢8日を超えるとアノマリー相関係数ACC＞0.6となるので、予測精度が非常に優れている｣と主張している。ACCは1に近づくにつれて、実況値と相関が高いことを意味する。0.6で優れているのか、という感じではある。
(4) 極端気象現象の捕捉
極端気象現象の捕捉能力を、NWPと競っている。例によって、グランド･トルゥースは、ERA5の再解析データである。極端気象現象として、極端な降水量と強風を取り上げている。予測値の裾が、どの程度の極端値に相当するか、で極端気象現象の捕捉能力を評価している。具体的内容は、以下の各項で述べる。
1⃣ 降水量
例えば予測値の99%ile値とグランド･トルゥースの99%ile値を比較して、どの程度の極端な降水量(正確には、降水強度と思われる)を捕捉できているかを評価する。グランド･トルゥースの99.99%ile値はおよそ50mmである(最近の例と照らすと、やや物足りない)。これに対して、NWPは、およそ40mm程度。FourCastNetは、およそ30mm程度であり、どちらのモデルも、極端な降水を過小評価しているが、NWP＞FourCastNetである。モデルの予測に系統的な偏りがあるかどうかを確認する｢相対分位誤差｣と呼ばれる指標で比較しても、同じ結果である。
2⃣ 風速
U10に関するグランド･トルゥースの99.99%ile値は、およそ20m/sである。NWPはわずかに過大評価しており、FourCastNetはわずかに過小評価している。風速に関しては、NWPと互角と考えても良いだろう。ただし、 相対分位誤差で測ると、ややNWPに軍配があがる。
【5】考察
(1) FourCastNetは、NWPと同程度の解像度を有するMLWPモデルが(より優れてはいないものの)、NWPと同程度の精度を叩き出すことを目標としている。故に、その目標は達成された、と評価できるかもしれない。
(2) FourCastNetは地表風速(並びに降水量)の予測に拘っているが、MLWPの目指すべきは、個別予報変数の予測ではないだろう。メンバー数を増やしたアンサンブル予報で、台風の進路予測や集中豪雨予測などの、極端気象現象の精度を向上させることが本命であろう。アンサンブルメンバー数1000の熱帯低気圧進路予測が示されても、良かったように思う。
(3) ViTをモデル･アーキテクチャとして選択するのは、自然と思われるが、効率的なトークン･ミキサーのアーキテクチャとして、フーリエ･ニューラル演算子を選ぶというのは、一見あまりピンと来ない。仕掛けが大袈裟･大仰過ぎるように感じるし、矛盾も感じる。
矛盾というのは、次のような意味である：元々、(フーリエ･)ニューラル演算子は、対象が無限次元＝連続系であるにも関わらず、有限次元しか扱えないニューラルネットワークで表現させることには無理がある。そこで、無限次元を扱えるニューラル演算子が導入された。しかし、本論文のトークン･ミキサーでは、不連続な画像を扱えるように改良した、AFNOという、(フーリエ･)ニューラル演算子の不連続バージョンを導入している。これでは、わざわざ、ニューラル演算子を導入した意味がないように感じられる。
ただ結果として、効率的な(計算複雑性を削減した)トークン･ミキサーが構築できたことは、間違いない(ので、結果オーライということで良いのだろう)。
(4) 本論文で指摘されている通り、FourCastNetは、物理ベースの NWP モデルと組み合わせることも可能である。つまり、Physics-Informed FourCastNetも可能である。
中国のファーウェイ(の子会社ファーウェイ･クラウド)が開発した Pangu-Weatherは、先行MLWPモデルであるFourCastNetの精度は、不十分(例えば、5 日間のZ500予報のRMSEは484.5で、IFSの333.7よりも”はるかに悪い”)と主張して、FourCastNetの改善を目指している。以下、本稿では[*3]を、本論文と呼ぶ。
ビジョン･トランスフォーマー(ViT)をモデル･アーキテクチャとして選んだFourCastNetに対して、Pangu-WeatherはSwinトランスフォーマーを選んでいる。Swinトランスフォーマーは、中国にあるマイクロソフトの研究所で開発されたトランスフォーマーである。もちろん、メイドイン中国であるから採用したわけではないだろう。
【1】本論文の主張
(1) Pangu-Weatherは、FourCastNetと比較して、大幅な精度向上が得られた。
(2) Pangu-Weatherは、極端気象現象の予測に優れている。
【2】事前整理
(1) Swinトランスフォーマー
ビジョン･トランスフォーマーから、Swinトランスフォーマーへの改良点は、以下の2点であると理解されている[*16]。
1⃣ ウィンドウ(窓)、ずらしウィンドウ(shifted window)
ビジョン･トランスフォーマーは、画像を複数のトークンに分割し、トークン間の類似性を内積計算によって定量化することで、トークン毎の特徴量を抽出する。トークンの数が増えると、計算コストが2次関数的に増大してしまうことが欠点であった。
Swinトランスフォーマーでは、トークンを｢ウィンドウ(窓)｣という複数のグループに分け、トークン間の類似性を、ウィンドウ内に限定して計算することで、計算コストを抑える。
このとき、(運悪く)トークンがウィンドウによって分断され、本来高いはずの類似性が、考慮されない可能性がある。そのため、ウィンドウをずらした、ずらしウィンドウ(shifted windowに対する正式な和訳ではない)に対しても内積計算を行う。当然、そうすることで、計算コストの抑制効果は減ずることになる。
2⃣ 階層型構造
Swinトランスフォーマーでは、隣り合う複数のトークンの特徴量をまとめ、より大きなトークンの特徴量に変換(集約)することで、サイズが異なるトークンの特徴量が考慮される階層型構造を採用している。この仕組みにより、Swinトランスフォーマーでは、全体の特徴量から細部の特徴量まで、サイズが異なる特徴量を考慮することができる。
この性質は、流体解析に向いていると考えられる。ViTでは、非局所的な相互作用を考慮することはできても、サイズの異なるスケール間の相互作用を考慮することは、難しいと思われるからである。大小様々なスケール間の相互作用を考慮することは、流体解析には、本質的に重要である。
👉 マイクロソフトの気象予測基盤モデルAuroraに関する論文[*29]には、Pangu-WeatherにおけるSwinトランスフォーマーに由来するこの階層構造が、｢グラフ･ニューラルネットワークのメッセージ･パッシングに相当する｣旨の記述がある。Auroraは"プロセッサー"にSwinトランスフォーマーを採用している。
【3】Pangu-Weatherの詳細
(1) モデルのオリジナリティ
1⃣ 3次元化
Pangu-Weatherでは、深層ニューラル ネットワークの入力と出力を 3 次元で概念化できるように、高さ情報を新しい次元に埋め込んでいる。3次元モデルは、異なる気圧レベルの大気状態間の関係を捉える能力があり、｢2 次元モデルと比較して、大幅な精度向上が得られる｣としている。それ自体、納得性は高い。
2⃣ 階層型時間集計アルゴリズム
予測リードタイムを増加させて、一連のモデルを学習する階層型時間集計アルゴリズムを適用した。例えば、リードタイムが56時間の場合、24時間予測モデルを2回、6時間予測モデルを1回、1時間予測モデルを2回実行する。こうすることで、｢反復回数が削減され、累積予測誤差が軽減された｣としている。
3⃣ 地球固有の位置バイアス
地球固有の事前分布を深いネットワークに注入するために、Swin の元の相対位置バイアスを置き換える地球固有の位置バイアスを設計した。ただし、この変更により、バイアス･パラメータの数が 527 倍に増加し、約 6,400 万個のパラメータが含まれることとなった。
(2) アーキテクチャの詳細
0⃣ エンコーダー層、デコーダー層
8つ(2+6)のエンコーダー層と 8つ(6+2)のデコーダー層を持つ、標準のエンコーダー･デコーダー･アーキテクチャを採用している。時間とメモリの両方の複雑さを軽減するため、標準の Swin トランスフォーマーよりも大幅に少なくなっている。より大きなメモリを備えた、より強力なハードウェアが準備された場合、ネットワークの深さを増やせるため、さらに精度が向上する(と期待される)。
1⃣ パッチ埋め込み
次元削減には、パッチ埋め込みを使用する。 ViTの文脈でパッチ埋め込みを説明すると、以下のようになる：2次元である画像データをパッチごとに1次元のシークエンス･データに変換し、線形射影して得られた出力を、パッチ埋め込みという。
標準のViTに従って、パッチ埋め込みには、GELU(Gaussian Error Linear Unit)を活性化関数とする線形層を使用した。
2⃣ ダウンサンプリングおよびアップサンプリング
Swin トランスフォーマー の実装に従い、ダウンサンプリングおよびアップサンプリング操作を使用して、異なる解像度の隣接層を接続する。ダウンサンプリングでは、4 つのトークンが1 つに埋め込まれ、線形層を実行して次元を削減した。アップサンプリングの場合は、逆の操作が実行される。
3⃣ 3次元地球専用のアーキテクチャ
各エンコーダ層とデコーダ層は、地球の形状に合わせて特別に設計されている。自己注意機構は、ViTの標準的な機構を使用している。計算コストを削減するために、特徴マップをウィンドウに分割するウィンドウ注意機構を(Swinトランスフォーマーから)継承している。
ずらし窓機構は、以下のように調整されて適用されている。層ごとに格子分割が、前の層とウィンドウ サイズの半分だけ異なる。経度方向の座標は周期的であるため、左右端の半分のウィンドウが 1 つの完全なウィンドウに埋め込まれる。緯度方向の座標は周期的ではないため、緯度方向に沿っては、埋め込みは実行されない。
4⃣ 地球固有の位置バイアス
Swin トランスフォーマーは、相対位置バイアスを使用して注意の並進不変コンポーネントを表現した。バイアスは各ウィンドウの相対座標に基づいて計算された。しかし、各トークンは地球の座標系上の絶対位置に対応する。一部の気象状態(予報変数)は、絶対位置と密接に関係している。
ジオポテンシャル、風速、温度の特性を捉えるために、地球固有の位置バイアスを導入した。これは、(相対座標ではなく)絶対座標に基づいて、各トークンに位置バイアスを追加することで機能する。
(3) データセット
”FourCastNetとの比較を公平にするため”、学習データは1979 年～2017 年(39 年分)のERA5データセット。検証データは2019 年のデータで、テストデータは2018 年のデータである。公平の意味は、よくわからない。FourCastNetの学習データは、1979 年～2015年(37年分)。検証データは16年と17年。テストデータは、18年以降、となっている。
等圧面は、13(50、100、150、200、250、300、400、500、600、700、850、925および 1,000 hPa)が選択されている。各変数は、721 × 1440ピクセルの 2次元場として表される(これはFourCastNetと同じ)。Pangu-Weatherの予報間隔(予報時間の最小単位)は、ERA5と同じ1時間(FourCastNetは6時間)である。
ちなみに、エンドツーエンドの学習には、192 個の NVIDIA Tesla-V100 GPU のクラスターで約 16 日かかる。
(4) ハイパーパラメータ等
① オプティマイザー → Adam
② エポック数 → 100
③ コスト関数 → 平均絶対誤差損失
④ バッチ サイズ → 192
⑤ 学習率 → 0.0005から始まって、0まで減少させる
⑥ 学習率スケジューリング → コサイン･スケジューリング(cosine decay)
⑦ 重み減衰 → 3 × 10−6
⑧ ドロップアウト率 → 0.2 (ScheduledDropPath)
⑨ 過学習の軽減策(⑦及び⑧に加えて) → 学習データのすべての開始時点を、各エポックでランダムに並べ替えた
⑩ 正規化 → 平均値を引いてから、標準偏差で割ることで行った。各変数の平均と標準偏差は、1979 年から 2017 年の気象データに基づいて計算された
⑪ 各変数の重み → 平均損失値に反比例し、これらの変数による寄与の等価性を容易にするように設計された
【4】比較結果
(1) 予報変数の比較
1⃣ FourCastNetとの比較
比較指標は、RSMEである。単一メンバーの予報(つまりアンサンブル予報ではない)で、Pangu-Weatherの 5 日間 Z500 予報のRMSE は 296.7である。対して、FourCastNetのRMSEは 462.5であり、大幅に改善されていると主張している。
2⃣ NWPとの比較
NWPは、ECMWFのIFSである。単一メンバーの予報(つまりアンサンブル予報ではない)では、Pangu-Weather の 5 日間 Z500 予報の RMSEは、先述の通り296.7である。NWPのRMSEは333.7であり、Pangu-Weatherの方が優れている、と主張している。
3⃣ 予測時間ゲイン
本論文では、Pangu-Weatherの優位性を実証するために、｢予測時間ゲイン｣と呼ぶ概念を導入した。これは、Pangu-Weather と他モデルが同じ精度を報告した場合のリードタイム間の平均差に相当する。
Pangu-Weather は通常、NWPよりも 10～15時間の予測時間の増加を示し、比湿度などの一部の変数では、その増加は 24時間を超える。この理由を、｢従来の NWP 手法では特定の変数を予測するのが難しい一方で、AI ベースの手法では、豊富な学習データから効果的なパターンを学習することでメリットが得られることを意味している｣と解釈している。FourCastNet との比較では、Pangu-Weather の予測時間の増加は 40時間と大きい。
4⃣ NWPとPangu-Weatherの予測において、質的な違いが存在するとの主張
Z500とT850、T2m、10 メートルの風速(U10 or V10)の結果を、NWP(IFS)および グラウンド･トゥルース(ERA5)と比較した。比較した結果、本論文は、｢Pangu-Weatherの出力とNWPの出力は、グラウンド･トゥルースに十分に近いが、Pangu-WeatherとNWPには明らかな違いがある｣と主張する。具体的には、｢Pangu-Weatherでは、より滑らかな等高線が生成され、NWPはそれほど滑らかではない。これは、Pangu-Weatherが隣接する地域に対して、同様の値を予測する傾向があることを示唆している｣とする。そして、その理由を、｢初期条件を使用して、偏微分方程式系を解くことによって、各グリッド セルで単一の推定値を計算する一方で、天候の混沌とした性質と、初期条件と格子点スケール未満の力学過程･物理過程に関する、必然的に不正確な知識が存在するため｣と推測している。
(2) 熱帯低気圧の進路予測
1⃣ 概要
本論文では、掲題に関するPangu-Weatherの主な利点として、｢初期段階で進路予測できること｣を上げている。地表変数の 1 つである平均海面気圧(MSLP)の極小値を見つけることにより、熱帯低気圧の進路予測において、高い精度を達成している。
2⃣ グランド･トルゥース、比較対象モデル、比較した指標
❶ グランド･トルゥース は、IBTrACSのデータである(IBTrACSは、Ⅰ-0【1】(3)2⃣を参照)。
❷ 比較対象モデル はECMWF-HRESである。
❸ 比較した指標は、発生から3 日後と 5 日後の｢台風の目｣の平均直接位置誤差である。
3⃣ 詳細アルゴリズム
開始時点とそれに対応する”台風の目”の初期位置を考慮して、6 時間予測アルゴリズムを繰り返し呼び出し、次の条件を満たす MSLP の極小値を探した(この条件は、ECMWFの熱帯低気圧進路予測プロトコルである。[*19]を参照)。
㈠ 北半球では半径 278 km 以内に 5×10−5より大きい最大850hPa 相対渦度があり、南半球では-5×10−5より小さい最小相対渦度があります。
㈡ 温帯低気圧の場合、半径 278 km 以内で最大の層厚†は、850hPa から 200hPa になる。
㈢ 熱帯低気圧が陸上にある場合、半径 278 km 以内では10mの最大風速が8 m/sより大きくなる。
㈣ ”台風の目”の目の位置が特定されると、進路予測アルゴリズムは引き続き 445 km 付近で次の位置を見つた。㈠～㈢の条件を満たすMSLP の極小値が見つからない場合、進路予測アルゴリズムは終了する。
† 気象学における層厚は、 二つの等圧面の間のジオポテンシャルの差である。
4⃣ 結果詳細
比較対象熱帯低気圧として、IBTrACS と ECMWF-HRES の両方に出現する、2018 年の名前付きの熱帯低気圧 88 個を選択した。なお、リードタイムは、6時間の倍数に設定された。
Pangu-Weather の 3 日目と 5 日目の平均直接位置誤差は、120.29 kmと195.65 kmである。ECMWF-HRESは、162.28 kmと272.10 kmであり、どちらもPangu-Weather の方が小さい(つまり、予測精度において優れている)。
さらに(2018年に発生した)個別の熱帯低気圧：台風コンレイ(2018年の台風25号)、台風ユトゥ、台風マーオン、ハリケーン･マイケルについて、HRESと比較し、以下のように総括している。コンレイについて、HRESは中国に上陸すると予測したが、Pangu-Weatherは上陸しないと予測。実際は上陸しなかった(台風の｢上陸｣、｢接近｣、｢通過｣について、Ⅰ-0【2】(0)2⃣を参照)。ユトゥについては、Pangu-WeatherはHRESよりも48時間前に正しい進路を予測した。ハリケーン･マイケルについては、実際に上陸した時間と予測した時間の差が、Pangu-Weatherは3時間で、HRESは18時間だった。さらに上陸地点が、HRESは東にずれていた。台風マーオンは、Pangu-Weatherは上陸地点を正しく予測したのに対して、HRESは誤っていた。
さらに、Pangu-Weatherは、アルゴリズム内部で、層厚と渦度を、ジオポテンシャルと風速から導いていた。本論文では、これをもって、｢Pangu-Weatherの出力が、説明可能であることを示している｣と主張している。
5⃣ 注意事項
本論文では、華々しい結果を訴えるとともに、ECMWF-HRESとの比較は、公正ではないことも併記している。つまり、HRESと Pangu-Weatherでは、入力データが異なる(Pangu-Weather は再解析データを使用している)ので、フェアではないと主張している。これは、その通りであろう。
(3) 極端気象現象の捕捉
FourCastNetで行われている(すなわち[*4]で行われている)ので、本論文でも掲題の作業が、軽く行われている。グランド･トルゥースはERA5、NWPはECMWFのIFS。対象となっている予報変数は、以下の通り：Z(500hPaと850hPa)、T(同じ)、U500、V500、Q500、U10、V10、T2M。
NWPと区別がつかないのは、Z500、Z850、T500、Q500。ほぼ区別がつかない(99.9%ile以降でのみ差が生じている)のは、T850、U500、V500、T2M。NWPのほうが若干精度が良い。U10とV10は、99%ile以降で差が生じている(最も分が悪い)。基本的に、Pangu-Weatherの予測は、全て、過小評価となっている。
(4) 為念：推論速度
Pangu-Weather の推論速度は、FourCastNet の推論速度に匹敵する。システムレベルの比較では、FourCastNet は Tesla-A100 GPU(312テラFLOPS)で、24 時間の予報を推測するのに 0.28 秒を必要とする。Pangu-Weather は Tesla-V100 GPU (120テラFLOPS)で、1.4 秒を必要とする。GPU のパフォーマンスを考慮すると、Pangu-Weather は FourCastNet よりも約 50% 遅い。
一方で、数百のノードを備えたスパコンで数時間を必要とするNWP(IFS)よりも 10,000 倍以上高速である。
【5】アンサンブル予報の研究
アンサンブル予報を研究するために、FourCastNetを調査するとともに、メンバー数100のアンサンブル予報を生成した。具体的には、99 個のランダムな摂動を生成し、それらを摂動のない初期状態に追加した。これらの予測結果を単純に平均することによって、100 メンバーのアンサンブル予報を生成した。なお、摂動には、パーリン･ノイズが含まれている(らしい)。パーリン･ノイズは、単純な乱数に比べて、滑らかに変化する自然なノイズを生成できる。複数周波数のノイズを統合することで、自然界に見られる自己相似性を再現できる、といった特徴がある[*17]。
結果は、短期(たとえば 1 日)予報では、単一メンバー法＞アンサンブル予報であり、リードタイムが5 ～ 7日の場合は、アンサンブル予報＞単一メンバー法である。これは FourCastNetとも一致している。
[結果1] アンサンブル予報は、短期予測では、予期せぬノイズを導入するリスクがある。
[結果2] アンサンブル予測は、Q500や U10などの滑らかでない変数に対して、より多くの利点をもたらす。
[結果3] Pangu-Weather の｢スプレッド･スキルの関係｣は 1 より小さい。アンサンブル予報が上手くいくためには、アンサンブル･スプレッドが、アンサンブル平均予報値の予報誤差(RSME)と同程度であることが望ましい。つまり、スプレッド･スキルの関係が1になることが望ましい。スプレッド･スキルの関係が1より小さいということは、スプレッドが小さい、つまり、アンサンブル･メンバーの分散が小さいことを意味している[*18]。
【6】考察
(0) ネイチャー･インデックス(Nature Index)によると、中国(のテクノロジー)企業の社員が、ネイチャー誌の論文において、単独著者となったのは初めてのことらしい。ネイチャー･インデックスは、高品質な科学ジャーナル82誌に掲載された原著論文を対象に、研究成果を国･機関別にプロファイリングしたデータベースである。
(1) FourCastNetとの比較で、Pangu-Weatherの方が精度が高い理由は、3次元化が大きいのではないかと推測する。
(2) 効率的なトークン･ミキサーについては、FourCastNetでは、ViTに外付けする形でフーリエ･ニューラル演算子が適用された。Pangu-Weatherでは、Swinトランスフォーマーが元々、有している性能で効率的な空間混合を実現している。結果として、FourCastNetの方が速い(より効率的に空間混合を実現している)が、Pangu-Weatherも、実用的に遜色ないスピードを実現していると考えられる。
(3) 熱帯低気圧の進路予測において、Pangu-Weatherは、｢層厚と渦度を、ジオポテンシャルと風速から導いている｣と、本論文で述べられている。これは、Pangu-Weatherの内部で正しく、力学過程がモデル化されている傍証であろう。Swinトランスフォーマーと流体解析の相性の良さを
(4) Swinトランスフォーマーと流体解析の相性は良いと考えられる。前項(3)は、その一つの実例と考えられる。相性が良いのであれば、Swinトランスフォーマーを使った流体解析例が、存在するはずである。しかし、適用例は、ほとんどないようである(知らないだけ？)。もし、相性は良いものの、適用例がほとんどないのであれば、それは、なぜだろうか？
グーグル･ディープマインドの研究者は、従来の数値予報(NWP)モデル及び、機械学習ベースの気象予測(MLWP)モデルを上回るモデルGraphCastを構築した、と主張する論文[*2](以下、本論文)を発表した(23年11月14日@Science)。従来のMLWPモデルとは、中国のファーウェイが開発したPangu-Weatherを指している。
先行モデルとの対比で言うと、FourCastNet及びPangu-Weatherのモデル･アーキテクチャは、ビジョン･トランスフォーマーである。一方、 GraphCastは、グラフ･ニューラルネットワーク(GNN)である。
🌏 [*2]の筆頭著者Rémi Lamは、2024年の｢nature's 10｣の一人(https://www.nature.com/immersive/d41586-024-03890-5/index.html)。nature's 10とは、当該年の｢科学を形作った｣とnatureが選んだ人。Rémi Lamは、GenCastの著者グループにも名を連ねる。
❚為参考❚
スタンフォード大学人間中心AI研究所(HAI)は、2017年から毎年｢AI Index Report｣という調査報告書を公開している(2020年はコロナ禍のため、例外)。2024年版第5章科学と医療[*55]には、5.1｢注目すべき科学的マイルストーン｣というセクションがある。6個のモデルが紹介されているが、その内2個が、気象関連である。2個の中にGraphCastが含まれている。
GraphCastは、以下のように紹介されている：GraphCastは、高精度の10日間天気予報を1分以内に提供する新しい天気予報システムである。グラフ･ニューラル･ネットワークと機械学習を活用して、GraphCastは膨大なデータセットを処理し、気温、風速、大気の状態などを予測する。GraphCastは、気象パターンを解読し、異常気象への備えを強化し、地球規模の気候研究に貢献する貴重なツールとなり得る。
【1】本論論文の主張
(1) 中期予報(最大10日先の天気予報)において、NWPよりも、概ね高精度である。
(2) 予報変数の予測において、Pangu-Weatherよりも、概ね高精度である。
(3) 熱帯低気圧の進路予測において、NWPよりも高精度である。
【2】事前整理
(1) ECMWFの熱帯低気圧進路予測プロトコル(トラッカー)
熱帯低気圧の進路予測は、かなり複雑なロジックのもとで、実行されている。132ページに及ぶ本論文のSupplementary Materials[*19]によれば、ECMWFのプロトコルは以下のようなものである。Pangu-Weatherは進路予測に、このプロトコルを、そのまま採用している。一方、GraphCastのプロトコルは異なる。
複数の時間ステップにわたる圧力レベル 200、500、700、850および 1000 hPaでの変数 10U、10V、MSL および U、V、Z のモデルの予測が与えられると、ECMWFのトラッカーは各時間ステップを順次処理して、 進路全体にわたる熱帯低気圧の位置を繰り返し予測する。トラッカーの各 6 時間予測には 2 つの主要なステップがある。最初のステップでは、熱帯低気圧の現在位置に基づいて、トラッカーは 6 時間先の位置推定値を計算する。2 番目のステップでは、新しい推定値の近くで、熱帯低気圧の中心に特徴的ないくつかの条件を満たす場所を探す。
新しいの熱帯低気圧の位置推定値が計算されると、トラッカーはこの推定値から 445 km 以内の平均海面気圧(MSL)のすべての極小値を調べる。それから、以下3 つの条件を満たす、現在の推定値に最も近い極小値を検索する。
⓵ 渦度チェック・・・極小値から 278 km 以内の 850 hPa での最大渦度は、北半球では 5×10ｰ5/sより大きく、南半球では5×10ｰ5/sより小さい。 渦度は水平風 (U および V) から求める。
⓶ 風速チェック・・・熱帯低気圧が上陸している場合、(中心から)半径278 km以内のU10(地表から上空10mにおける東西風)の最大風速が8m/sを超える。
⓷ 層厚のチェック・・・低気圧が温帯性の場合、半径 278 km 以内に、850hPaから200hPaの間に、最大層厚がある(層厚は、Z850 ～ Z200 として定義される)。
⓸ これらすべての条件を満たす最小値がない場合、トラッカーは、熱帯低気圧が存在しないとみなす。
このトラッカーを使用すると、熱帯低気圧が特殊な条件下で一時的に消滅し、その後再び出現することができる。一方、GraphCastが採用したプロトコルでは、熱帯低気圧はいったん消えると、復活はしない。
【3】GraphCastの詳細
以下、ことさら断らない場合でも、[*19]に含まれる情報を適宜、追加している。
(0) 改めて概要
1⃣ MLWPモデルの概要
GraphCastの入力は、空間解像度が緯度/経度0.25°である。これは、FourCastNetやPangu-Weatherと同じ。予測間隔(リードタイム)は、6時間である。これはFourCastNetと同じであるが、Pangu-Weatherは1時間である。詳しく言うと、最新の2時刻(現在時刻と 6 時間前)における気象データを入力として受け取り、6 時間先の気象状態を予測する。
単体の｢Google Cloud TPU v4 デバイス｣で、1 分以内に正確な10日間の天気予報を生成する。ECMWFのHRES(空間分解能0.1°)は、11,664 コアのクラスター上で実行され、約1時間を要する。
2⃣ アーキテクチャ概要
(先述あるいは既述の通り)GraphCastは、GNNベースのMLWPである。気象データはグラフとしてモデル化され、ノードは、地球のグリッド･セルを表す。GraphCastは、このグラフベースの表現により、｢データ内の『局所的(短距離)並びに、非局所的(長距離)』相互作用を、どちらも捉えることができる｣ところに特色がある。もちろん、トランスフォーマー(ビジョン･トランスフォーマー)も自己注意機構により、非局所的相互作用を取り込むことが可能である。しかし、[*19]によると、｢トランスフォーマーでは、計算の複雑さを軽減するために、全ての相互作用は考慮されない｣。
GNNを選択した理由として、｢これまでの GNN ベースの学習シミュレータは、偏微分方程式によってモデル化された流体およびその他のシステムの複雑な力学を学習するのに非常に効果的であり、気象力学のモデル化への適合性をサポートしている｣が上げられている。
(1) アーキテクチャの詳細
0⃣ 前説
GraphCastのモデル･アーキテクチャは、グラフ･オートエンコーダである。本論文によれば、｢エンコーダ＋プロセッサ＋デコーダ｣構成のGNN に基づいたニューラル ネットワーク･アーキテクチャとして実装されている。部外者(?)が書いた[*20]には、[*2]では全く示されていない、エンベッダーという名称のコンポーネントが追加されている。[*19]にも、その文言はないが、そのような解釈をしたほうが、見通しはよくなるだろう(し、一般化して捉えることができるだろう)。[*20]によると、その機能は｢入力された特徴量を潜在空間に埋め込む｣ことである。
1⃣ エンコーダ(符号化器)
符号化器 の目的は、マルチ･メッシュ上で排他的に実行されるプロセッサ用に、潜在空間に埋め込まれた特徴量を準備することである。まず、グリッド･ノード、メッシュ･ノード、メッシュ･エッジ、｢グリッドからメッシュ･エッジ｣、および｢メッシュからグリッド･エッジ｣の5つの特徴量を、5つの多層パーセプトロンを使って、固定サイズの潜在空間に埋め込む。
潜在空間への埋め込みは、マルチ･メッシュを扱うという意味で重要な処方と思われる。それは、マルチ･メッシュで扱うことが物理的に相応しい系を、深層学習の枠組みで実装する、という意味において、である(あくまで、物理の文脈ではなく、深層学習の文脈で重要という意味)。潜在空間へ埋め込むことで、大きさの異なるメッシュ(スケールの異なる物理現象)間の相互作用を、小さいスケール間の相互作用と同じ枠組みで処理することを、可能にしていると考えられる。
2⃣ マルチ･メッシュ
GraphCastでは、地球はマルチ･メッシュで表現される。実は地味に、マルチ･メッシュは大きな働きをしている、と理解している([*19]でマルチ･メッシュのアブレーション分析が行われていることから判断して、その理解は正しいだろう)。マルチ･メッシュのおかげで、❶メッセージパッシングの数を減らすことができ、❷大小さまざまなスケール(具体的には、6つのスケール)、つまりマルチ･スケールの相互作用を考慮することができる。これは同じ事象(メリット)を、❶機械学習の文脈、❷流体力学の文脈で述べている。なお、❶については、下記(2)過剰平滑化対策、を参照。
❷は、気象現象を流体現象と捉えた場合、本質的に重要である。Pangu-Weatherでは、Swinトランスフォーマーのアーキテクチャ特性によって、❷を実現していると考えられる。
マルチメッシュは、解像度の高い正20面体🐾メッシュの集合であり、正20面体(ノード 12 個、面 20 個、エッジ 30 個)を6 回繰り返し微細化することによって定義される。各回の微細化で、三角形が 4 つの小さな三角形に分割され、面とエッジの数が4 倍になる。最終的には、ノード数40,962と 面の数81,920を持つ、正二十面体メッシュノードが構築される。ちなみに46かける20は81,920となるが、46かける12は、49,152で40,962にはならない。
🐾 全球モデルのメッシュは、正20面体。
3⃣ プロセッサ
(GraphCastの)プロセッサは、少ないメッセージパッシング･ステップで、各マルチメッシュ･ノードを更新する。具体的には、各層の MLP に非共有ニューラル ネットワークの重みを付けた、16の非共有 GNN 層を使用して、局所的および非局所的な情報を、効率的に伝播する。GraphCast内の全MLP のニューラル パラメーターは、層間では共有されない。ただし、特定の層内では、常に空間的位置間で共有される。
GNN 層は、最初に隣接ノードの情報を使用して、各メッシュ･エッジを更新する。次に、各メッシュ･ノードを更新し、そのメッシュ･ノードに到着する、全てのエッジからの情報を集約する。そして両方を更新した後、特徴表現は、残差接続を使用して更新される。
4⃣ デコーダ(復号器)
復号器は、更新処理されたマルチ･メッシュ上の特徴量を、グリッド表現にマッピングする。
(2)過剰平滑化(Over-Smoothing)対策
GNNは、メッセージ･パッシング(MP)を用いて、情報を取り込む。MP自体は、近隣ノードの情報(つまり局所的情報)を取り込むだけなので、非局所的な情報を取り込むには、MPを何度も繰り返す(反復する)必要がある。しかし、MPを多数回行うと、｢過剰平滑化｣と呼ばれる現象が発生し、非局所的な情報の取り込みがうまく行われない。
GraphCastは、サイズが異なるマルチ･メッシュを用いている。大きなサイズのメッシュを使うことで、長距離＝非局所的な情報を、比較的少ない回数で取り込むことを可能としている、と考えられる。比較的という意味は、均一の小さいサイズのメッシュを使って、MPを多数回行うことに比較して、という意味である。
(3) データセット
ERA5の再解析データセットを使用する。学習データは、1979～2015年(37年分)のデータを使用。検証データは2016～2017年、テストデータは、2018～2021年を使用。格子点(グリッド･ポイント)は、721 × 1440である。
(4) ハイパーパラメータなど諸々
1⃣ 最適化関係
勾配消失を防ぐために、L2ノルムを用いた、勾配クリッピングを適用している。しきい値は、32である。
① コスト関数 → 平均自乗誤差(MSE)
② オプティマイザー → AdamW(パラメータβ1＝0.9、β2＝0.95)
③ 重み減衰率 → 0.1
④ 活性化関数 → Swish
2⃣ カリキュラム学習
カリキュラム学習を採用している。具体的には、学習率と自己回帰ステップ数を変化させた、3段階のカリキュラム学習を実施している。第1段階は1000回の勾配降下更新で、1回の自己回帰ステップ。学習率は0から1.0×10-3まで直線的に増加する。
第2段階は、299,000回の勾配降下更新と、1回の自己回帰ステップ。ハーフコサイン減衰関数で0に戻る、学習率スケジュールからなる。第3段階は11,000回の勾配降下更新、自己回帰ステップ数は2から12まで、1000回の勾配降下更新ごとに1ずつ増加する。学習率は3.0×10-7に固定された。
3⃣ 正規化
❶ 入力の正規化・・・すべての入力を正規化した。まず各変数について、1979 年から2015 年までの等圧面ごとの平均と標準偏差を計算した。それを使用して、平均0と分散1に正規化しました。相対的なエッジの距離と長さについては、特徴量を、最長のエッジの長さに正規化した。
❷ 出力の正規化・・・モデルの出力は、Ytである。ここでYtは、時刻t+1における変数Xt+1と、時刻tにおける変数Xtの差である。つまり、Yt= Xt+1ーXtである。
各等圧面ごとに、Ytの標準偏差を計算することで、モデル出力を正規化した。具体的には、モデル出力に、この標準偏差を乗算して Y^t+1を求める。そして、XtにY^t+1を加えて、X^t+1を求める。
4⃣ ニューラルネットワークのパラメータ化
GraphCast 内のニューラル ネットワークは、すべて多層パーセプトロン(MLP)であり、1 つの隠れ層を持つ。隠れ層と出力層のサイズは512である。ただしデコーダの MLPの最終層のサイズは227 である。これは、予測変数の数と一致させるためである(地表変数5＋上空変数6×等圧面37＝227)。
【4】比較結果
(0) HRESのグランド･トルゥース[*19]
1⃣ 概要
GraphCastはERA5のデータを予測するように学習されており、ERA5のデータを入力として使用する。しかし、HRES の予測は HRES の解析に基づいて初期化されている。その理由は、｢一般に、モデルを自身の解析結果と照らし合わせて検証することで、最良の予測精度が得られる｣からである。そこで、ERA5のグランドトゥルースに対してHRES 予測を評価するのではなく、HRES 予測の第0ステップでさえ誤差がゼロでないことを意味するHRES-fc0データセットを構築し、将来の初期化におけるHRES 予測の初期時間ステップを格納した。HRES-fc0 を HRES 予測を評価するためのグランドトゥルースとして使用する。
2⃣ 初期化の違い
公平な比較を行うために、GraphCast の ERA5 初期条件が、HRES で使用される条件よりもさらに先の将来を見据えた同化ウィンドウから導出されたことを確認する必要がある。HRES初期化は、常に 3 時間後の観測値を同化するが、ERA5初期化は 9 時間後および 3 時間先の観測値を同化する。最大 3.75日までのパフォーマンスを比較する場合、HRES に対して同じ初期化を使用する。さらに、HRES アーカイブされた予測は、｢0時、12時｣初期化からのみ利用可能である。
(1) 比較モデルと比較指標
比較するベースラインは、HRESであり、比較指標はRSMEとACC。(0)で既述している通り、グランド･トルゥースは、GraphCast(MLWP)とHRES(ベースラインであるNWP)で異なる。
要するに、GraphCastとNWPの比較はシンプルではなく、ややこしいが、比較はかなり厳密だと思われる。つまり、NWPに比べて高精度は、確からしいと思われる。
(2) 比較結果
1⃣ HRESとの比較
比較対象となった変数は、地表変数が4つ(MSL、2T、10U、10V)で上空変数が5つ(Z、Q、T、U、V)で等圧面が13。従って、4+5×13＝69。12時間ごとに、(中期予報なので)10日間の期間にわたり比較している。つまり合計で69×2×10＝1,380個の変数(ターゲット)を比較していることになる。結果は、GraphCast は 1,380 ターゲットの 90.3% で HRES を上回り、ターゲットの89.9%で HRES を大幅に上回った。[*19]では、p値及びt値で有意性検定も行っている。
なお、HRESがGraphCast よりも優れたパフォーマンスを示したケースは、成層圏に局在していることを指摘している。そして、50 hPa レベルを除くと96.9%、50および 100 hPaを除外すると99.7%で HRESよりも高精度である、と述べている。さらに、この結果は、地域によらず、世界中で一般的に当てはまることも指摘している。
2⃣ Pangu-Weatherとの比較
比較した指標は、Z500、T500、T850、Q500、U500、V500、2T(Pangu-Weatherでは、T2M)、10U(Pangu-Weatherでは、U10)の8つ。6時間ごとの予測値を10日間にわたり比較した。ただし、[*19]に示された図S9は7日間まで。RMSEを指標として、ほとんど全てで、GraphCastが上回っている(RSMEが小さい)。Pangu-Weather が GraphCast を上回ったケースは、リードタイム6 時間と12時間の Z500だけ(つまり２つだけ)である。割合として、250/252＝99.2%という数字を出している。
(3) 熱帯低気圧の進路予測
極端気象現象の予測として、熱帯低気圧の進路予測のみを取り上げた(本論文では、大気河川と極端な高温･低温も取り上げられているが、割愛)。
[*19]を含めても、細かい記述が抜けているように感じる。さすがに熱帯低気圧の位置は、”台風の目の位置”で捕捉している、という理解で良いと思う(㊀明示はされていない)。グランド･トルゥースは、 IBTrACS(Pangu-Weatherと同じ)。ベースラインは、TIGGEアーカイブに保存されているECMWFのプロトコルを採用した、HRESの 0.1° 予報から取得した進路予測と記載されている。㊁再解析データが用いられているか否かは明示されていない(が、他所の記述を見る限り、再解析データは使われていないだろう)。㊂アンサンブル予報に関する情報の記載は、ない。
予測した熱帯低気圧の位置と、グランドトルゥースとの差の｢中央値｣を、評価指標としている。Pangu-Weatherは平均値であった。結果として、GraphCastは、予測期間5日にわたる全期間において、HRESよりも精度が高かった(Pangu-Weatherは、3日目と5日目をピンポイントで比較していたが、それは重要な差ではない)。なお、GraphCastとECMWFでは、進路予測のプロトコルが同一ではない。
【5】アブレーション分析[*19]
マルチ･メッシュが、GraphCast のパフォーマンスにどのような影響を与えるかをよりよく理解するために、GraphCast のパフォーマンスが、マルチ･メッシュなしで学習されたモデルと比較されている。マルチ･メッシュなしモデルのアーキテクチャは 、GraphCast と同じ(同じエンコーダとデコーダ、同じ数のノードを含む)である。マルチ･メッシュなしモデルは、短距離の情報のみを伝播することができる(GraphCastは長距離の情報も伝播できる)。
GraphCast は、50hPaで 5 日を超えるリードタイムを除いて、すべての予報変数に対してマルチ･メッシュ構造の恩恵を受けている。この改善は、リードタイムが 5 日未満の場合、すべての等圧面にわたるジオポテンシャルと平均海面圧力で特に顕著である。 また、5 日未満のリードタイムで GraphCast が、ジオポテンシャルでHRES を上回るパフォーマンスを発揮するには、マルチメッシュが不可欠であることを示された。
【6】考察
(0) トランスフォーマーを作ったグーグルは、MLWPにトランスフォーマーを使わない。日本なら、他社がトランスフォーマーを使って高精度なモデルを作ると、家元が黙っていて良いのか･･･という議論になりそう。
(1) Pangu-Weatherに勝っている理由は、モデル･アーキテクチャに由来するわけではないだろう。つまり、ViTよりもGNNが勝っている、ということではないだろう。精度は、ハイパーパラメータのチューニングといった工夫で変わるのだろう。
(2) MLWP＞NWPという結論になっている。しかし、それはMLWPがERA5の(大量の)再解析データを学習データとして、与えられているからである。つまり、MLWPがグランド･トルゥースを高精度に近似できるようになった、ということであって、NWPとの比較で優位になった、というのは少し論点が異なるだろう。
(3) 予報変数の予測、熱帯低気圧の進路予測や集中豪雨の発生予測において、今後は、再解析データを入力とするのではなく、観測データを入力として、NWPと優劣を比較するフェーズに移行するのだろう。その場合、再解析データを使って事前学習されたMLWPに、観測データから別のニューラルネットワークが作った再解析データを入力する、という枠組みになると思われる。生成モデルを使う、あるいはニューラル演算子を使ったベイズ的能動学習を使う、という手があるだろうか。
(4) アブレーション分析やスペクトル分析がある一方で、アンサンブル予報に関する情報がないのは寂しい。
【0】はじめに
マイクロソフト他†1の研究者は、｢基盤モデルが、最先端のNWP(数値計算ベースの気象予測モデル)を凌駕することを示した｣と主張する論文[*29](以下、本論文)を発表した(24年5月28日@arXiv)。その基盤モデルが、多様な気象･気候データを使って、100万時間以上の事前学習を実施したAuroraである。なお、グーグルのGraphCastを最高性能MLWPとして取り上げ、タイマン勝負の結果、Auroraが優位であるとも主張している。
†1 墺ヨハネス･ケプラー大学、英ケンブリッジ大学、蘭アムステルダム大学、Polyコーポレーション。
❚為参考❚
スタンフォード大学人間中心AI研究所(HAI)は、2017年から毎年｢AI Index Report｣という調査報告書を公開している(2020年はコロナ禍のため、例外)。最新の2025年版第5章科学と医療[*54]に、｢注目すべきモデル｣という囲み記事がある。9個のモデルが紹介されているが、その内5個が、気象関連である。5個の中にはAuroraも含まれている。
AI Index Report2025では、次のように簡潔に紹介されている(ドメインは、earth science)：Auroraは、100万時間を超える地球システムデータに基づいて学習された大規模な基盤モデルであり、大気質、海洋波浪、サイクロンの進路、高解像度気象に関する最先端の予測を提供する。Auroraは、従来のシステムを凌駕する性能を持ちながら、わずかな計算コストで動作し、最小限のリソースでドメイン間の再学習が可能である。
【1】本論文の主張
本論文は、以下を主張する：
『基盤モデルの原則が、気象予測にも当てはまるという証拠を、初めて示した』。
基盤モデルの原則とは、｢データとモデルがスケーリングされると、パフォーマンスが大幅かつ予測可能な形で向上する｣ことを指している。データのスケーリングとは、㊀ERA5以外のデータを考慮すること、㊁データの解像度を0.25°から0.1°に上げること、を指している。モデルのスケーリングとは、｢パラメーター数が 2倍になるごとに性能が約 5% 向上する｣ことを指している。
【2】本論文のセットアップ1
(0) Auroraの概要
Auroraは基盤モデルであり、事前学習＋再学習＋追加学習を実行している。なお、再学習は、❶短いリードタイム再学習と、❷ロールアウト再学習、の2種類を実行している。事前学習の検証には、0.25°解像度のHRES予測の 1 年分(2020 年)を使用する。テスト年は、データセットに応じて2022年と 2023 年である。
(1) 学習データセット
1⃣ 事前学習
Auroraは、ERA5を含めて以下10個のデータセットを用いて、事前学習を行った。
① ERA5は、ECMWF(欧州中期予報センター) の気象と気候に関するグローバル再解析データセット。
② HRES予測は、ECMWF が実行する運用NWP予測モデルの高解像度バージョンを指す。これは最も正確な NWP 予測モデルと考えられており、0.1°の解像度で実行される。本論文では、0.1°の解像度データと、解像度0.25°に再グリッド化されたデータの両方を使用する。
③ HRES解析は、ECWMF の公式解析データを表し、HRES-T0解析の上に追加の同化ステップが含まれている。解像度は、0.45°。HRES-T0解析は、HRES 予測を初期化するために使用される初期条件を提供し、予測の質を評価するためのグラウンド･トゥルースと見なされることがよくある。
④ IFS ENSは、IFS のアンサンブル予報データである。IFS ENSでは、50のアンサンブル･メンバーが18 km†2の粗い解像度で実行される。確率モデルで表される物理法則を適用し、初期条件に摂動を与えて、アンサンブル･メンバーは生成される。WeatherBench2 リポジトリのデータセットが使用されるが、このリポジトリには、500、700、850hPa の等圧面しかない。
⑤ IFS ENS 平均には、IFS ENS に基づく各変数の平均予測データが含まれる。WeatherBench2によって提供され、500、700、850hPa の等圧面しかない。
⑥ NOAAによるGFS予測は、基本解像度 18 km の運用予測データを提供する。本論文では、0.25°に再グリッド化されたデータを使用する。これらの予測のゼロ時間(初期化)は、IFS-HRES ゼロ時間と同様に導出されたリアルタイム運用解析である。
⑦ NOAAによるGFS-T0 解析は、GFS予測のゼロ時間(初期化)から得られる、リアルタイムの運用解析を指す。
⑧ NOAAによるGEFS 再予測は、限られたカバレッジ、機器または観測システムのバイアス、およびモデル自体の制限など、入力データに内在する不確実性に対処するために、21 のアンサンブル･メンバーに基づいている。実際には、このような大量のデータは、過去数か月分しかアーカイブされていない。そのため、2000 年から 2019 年までの再予測データを使用し、毎日 00 UTC(協定世界時)に再解析初期条件で初期化する。この設定では、アンサンブル･メンバーは 5 つだけであり、すべてが含まれている。GEFS には 6 つの圧力レベルがあり、WeatherBench2 の 850、925、1000 と一致する 3 つの等圧面を使用する。
⑨ CMIP6†3のデータセットは、CMCC-CM2-VHR4と ECMWF-IFS-HR の2 つである。それぞれ 7 つの等圧面(50、250、500、600、700、850、925hPa)がある。
⑩ MERRA-2は、米航空宇宙局(NASA)GMAO†4による大気再解析データセットで、宇宙からのエアロゾル観測を組み込んでいる。解像度は、0.625°×0.5°である。WeatherBench2 の13 レベルに対応する等圧面を使用する。
†2 2023 年6月27日以前。27日以後は、解像度が9 kmになった。
†3 CMIP6(Coupled Model Inter-comparison Project phase 6:第六期結合モデル相互比較プロジェクト) は気候モデル相互比較プロジェクトであり、陸、海、大気、エアロゾル変数を含むさまざまな気候モデリング実験を組み合わせている。
†4 Global Modeling and Assimilation Office
2⃣ 追加学習
㊀ CAMS(コペルニクス大気監視サービス)は、大気質、気候パターン、紫外線レベルを監視するヨーロッパを代表するサービスである†5。具体的には、世界および地域の大気質、インベントリ･ベースの排出量、観測ベースの温室効果ガスおよびバイオマス燃焼、太陽エネルギー、オゾン、紫外線による地表フラックスに関する情報を1日2回提供している。
CAMSのモデルは頻繁に更新される(とされている)。モデルのデータ解像度は 0.4°であり、地球全体の大気組成の予測、解析、再解析データを作成する。WeatherBench2の13等圧面を使用する。
㊁ CAMSRAは、CAMSからの第 4 世代 ECMWF大気構成グローバル再解析(EAC4)を指す。データ解像度は 0.75°で、CAMS と同様に、気象変数のほか、大気の構成を表す変数が含まれている。WeatherBench2 の13 等圧面を使用する。
†5 ちなみに、｢コペルニクス気候変動サービス｣は、欧州および世界全体の気候の過去、現在、そして未来に関する品質管理された情報を提供している。｢コペルニクス緊急管理サービス｣は、自然災害の早期警報を提供している。
【3】本論文のセットアップ2
(1) モデル･アーキテクチャ
0⃣ 概要
基本アーキテクチャは、オートエンコーダ(エンコーダ･デコーダ)＋プロセッサーという構成で、これは、既存モデル(FourCastNet、Pangu-Weather及びGraphCast)と同じである。プロセッサーは、(Swin)トランスフォーマー。ただしAuroraは、多種多様なデータセットを用いて事前学習するため、Auroraのトランスフォーマーは、マルチ･モダリティに対応するトランスフォーマーPerceiver(及びPerceiver IO)†6である。グーグル･ディープマインドが開発(2021年発表)したPerceiverは、｢クロス･アテンションとlatentトランスフォーマー｣という基本モジュールの繰り返し、で構成される。
†6 以降では、PerceiverとPerceiver IOを区別しない。なお、Perceiver IOでは、入力データと同じモダリティの出力を得ることができる(ので、正確にはAuroraはPerceiver IOであろう)。
1⃣ エンコーダー(符号化器)
⓪ 概要
Auroraは、多種多様なデータで事前学習しているところに、大きな特徴がある。多様なデータセットは、【2】(1)で示している(ERA5を含めた10個のデータセット)。多種多様なデータセットを(事前)学習することは、学習モデルにとってチャレンジである。Auroraでは、エンコーダーの設計を工夫することで対応した。具体的には、多種多様なデータセットを、標準化された3次元テンソルにマッピングする、柔軟なエンコーダーを設計した。標準化された3次元テンソルが、モデル入力である。
① 入力
エンコーダーは、すべての変数を、通常の緯度経度グリッド上の H(高さ)×W(幅)画像として扱う。各変数について、現時刻 t の状態と時刻 t−1 の状態を含める。これにより、T†7×H×W テンソルが生成される。等圧面をC、大気変数をVAで表すと、大気変数のデータセットは、 VA×C×T×H×W テンソルと表される。同様に、地表変数をVSで表すと、地表変数のデータセットは、VS×T×H×Wテンソル(等圧面は1000hPaであり、C=1)として表される†8。
†7 Tは、時間次元である。データは12時間毎更新なので、T＝2となる。
†8 実際には、すべての計算がバッチ処理されるため、これらのテンソルの前に追加のバッチ次元が追加される。
② 等圧面における埋め込み(embedding)
Auroraは、標準のヴィジョン･トランスフォーマーと同様に、H×W の各画像を P×P パッチに分割する。各等圧面のパッチは、線形層によって ℝDのベクトル†9にマッピングされる(埋め込まれる)。つまり、VA×C×T×P×P → C×D およびVS×T×P×P → 1×D。異なる変数を持つデータセットに対応するために、この線形変換は、各変数 v に対して、その変数に固有の重み Wvのセットを使用して、動的に構築される。
†9 言うまでもなく、ℝDのベクトルは、D個の実数値を成分として持つベクトルである：(x1,x2,･･･xD)T。さらに言うまでもなく、ここで現れた(･･･)TのTは転置(Transpose)の意味であり、横を縦にする、という意味である。
③ 集約(aggregation)
次のステップは、データセット間で数が変化する可能性がある大気中の｢物理的な圧力面｣を、固定された｢潜在的(latent)な圧力面｣にマッピングする(結果として、減らす)ことである。物理的な圧力面は、上記②の等圧面である。
既述(0⃣)の通り、Auroraのモデル･アーキテクチャは、Perceiverである。Perceiver への入力は(Perceiverもトランスフォーマーなので)、潜在的クエリ･ベクトルQと、潜在的キー･ベクトルKと潜在的バリュー･ベクトルV†10である。出力は、大気変数の潜在状態を符号化する 3×Dテンソルである†11。地表変数の埋め込みは、残差 MLP(多層パーセプトロン)に渡されるだけである。この地表変数の潜在状態は、垂直方向の大気変数の潜在状態と結合され、各パッチの位置における気象状態の(3 + 1)× D 潜在表現を生成する。
†10 クエリ(query)ベクトル、キー(key)ベクトル、バリュー(value)ベクトル という文言は、トランスフォーマーで用いられている文言であり、正式な訳語は(おそらく)ない。
†11 言うまでもなく、この場合のDは、ℝDに埋め込まれた等圧面におけるデータの個数である。つまり、物理的な空間次元を意味していない。また３×Dの3は、入力が(Q,K,V)であることから来ている。
④ 情報を追加してデータセットを完成させる
各パッチの潜在表現を集約すると、3 × H/P × W/Pテンソルという3次元テンソルが得られる(H×W の各画像を P×P パッチに分割したのであった)。このデータ(トランスフォーマー的に言えば、トークン)に、緯度と経度の座標、および物理的サイズ†12に関する情報を追加する。最後に、各パッチの絶対時間†13情報を含める。
†12 パッチあたりの面積(単位km2)。ただし、緯度･経度に依存する。
†13 1970 年1月1日を起点とする時間を指している(ようである)。
2⃣ プロセッサー
Swinトランスフォーマーを用いる。Pangu-WeatherもSwinトランスフォーマーを使用している。ただしPangu-Weatherで行われた、固定空間解像度の入力を必要とする、地球固有の位置バイアスは、Auroraでは用いられていない。代わりに、エンコーダーで位置符号化とスケール符号化を選択する。こうすることで、入力に(解像度等の)制約が課せられない。
Auroraでは学習全体の安定性を高めるために、Resポスト･ノルム･レイヤー正規化❚補足1❚を使用するが、(Swinバージョン2の)コサイン･アテンション❚補足2❚ではなく、(トランスフォーマーの標準である)ドット積アテンションを選択している。
❚補足1❚ Resポスト･ノルム･レイヤー正規化
オリジナルのトランスフォーマーにおける｢レイヤー正規化(LN)｣は、ポスト･ノルムと呼ばれ、マルチ･ヘッド･アテンション層の後に(LNを実装した層が)配置された。その後、LNをアテンション層と残差MLP(多層パーセプトロン)の前に、それぞれ配置するプレ･ノルムが採用され、ポスト･ノルムよりも優れた性能を示した。Swinバージョン1もプレ･ノルムである。Swinバージョン2では、アテンション層と残差MLPの後ろに、それぞれLNを置いた。これを、Resポスト･ノルムと呼ぶ。Resポスト･ノルムは、ポスト･ノルムのようにメイン分岐で勾配を止めず、プレ･ノルムよりも勾配爆発を避けるために、アテンション層と残差MLPの出力を、より良く制御する(とされる)。
レイヤー正規化は、可変長系列データ用バッチ正規化、という理解で良い(だろう)。言うまでもなく、バッチ正規化は、学習の｢安定化と高速化｣を目的としている[*30]。
❚補足2❚ コサイン･アテンション
オリジナルの自己注意(セルフ･アテンション)計算では、画素対の類似項はクエリー･ベクトルとキー･ベクトルのドット積(内積)として計算される。このアプローチを大規模な視覚モデルで用いると、特にResポスト･ノルム構成において、少数の画素ペアに支配されることが多い。この問題を緩和するために、Swinバージョン2では、画素ペアの注目度ロジットをスケーリングされた余弦関数(いわゆるcos)で計算するアプローチを提案した。これを、コサイン･アテンションと呼ぶ[*30(再)]。
3⃣ デコーダー(復号器)
デコーダーは、潜在的な圧力面における出力を、通常の緯度経度グリッド上の画像にマッピングする。エンコーダーのパッチ埋め込み層と同様に、出力パッチを構築する線形層は、各変数に関連付けられた重みを選択することにより動的に構築される。この全体的なアーキテクチャにより、デコーダーは任意の変数セットに対して、任意の等圧面で予測を出力できる。
4⃣ 正規化
Aurora は、エンコーダーで処理する前にすべての変数を正規化し、デコーダーの出力を非正規化して最終予測を生成する。全ての地表変数及び、全ての等圧面における全ての大気変数は、スケール(一定の空間的係数)とセンター(中心値)によって個別に正規化される。センターは、ERA5学習データ全体にわたって計算された、経験的期待値によって推定される。スケールは、ERA5学習データ全体にわたって計算された、経験的標準偏差によって推定される。スケールとセンターは、すべてのデータセットに使用される。最終予測は、デコーダーの出力を非正規化する(元に戻す)ことによって生成される。
(2) 学習方法、ハイパーパラメータ等
損失関数は、予測値とグランドトルゥースとの間の、平均絶対誤差(MAE)で与える。なお、損失関数における重みは、変数毎に細かく設定されている。例えば、事前学習における、平均海面気圧(MSLP)の重み＝1.5、風速10mのU成分U10の重み＝0.77、風速10mのV成分V10の重み＝0.66、地表から2mの気温Tの重み＝3.0等である。再学習だと1.5→1.6、0.77→0.77(変わっていない)、0.66→0.66(変わっていない)、3.0→3.5等となっている。
1⃣ 事前学習
32 個のGPUを使って、ステップ数150,000で、事前学習が行われた。バッチサイズは、GPU毎に 1。学習率は、1,000 ステップに対してゼロから線形ウォームアップで増大させた後、ハーフコサイン減衰関数で減衰させる。基本学習率は 5×10－4。使用するオプティマイザーは AdamW。AdamW の重み減衰は、 5×10－6。正則化の手法としては、ドロップパス †14を採用している。ドロップ確率は 0.2。
モデルをメモリに収めるために、バックボーン層にアクティベーション･チェックポイント†15を使用し、すべてのモデル勾配をGPU 間でシャーディング†16する。
†14 ドロップアウトは、中間層のノード出力を一定の割合でランダムに0とすることで、結合を欠落させる正則化手法であった。これは、横方向(幅方向)にネットワークを小さくしていると考えられる。これに対し、ドロップパス(stochastic depthとも呼ばれる)は、層数を変更することで縦方向(深さ方向)にネットワークを小さくする正則化手法である。
†15 アクティベーション･チェックポイント(勾配チェックポイントとも呼ばれる)はメモリ使用量を削減する手法であり、具体的には、次のように学習が実行される：フォワード処理では、一部の層の演算結果だけを保持する。バックワード処理では、演算結果が保持されていない層に関してはフォワード処理を再実行する。
†16 シャーディング(データ並列処理)は、メモリを節約する分散学習法である。具体的には、モデルの状態(モデルパラメータ、勾配、およびオプティマイザの状態)をデータ並列グループ内のGPU間で分割することで、メモリを節約する。
2⃣ 短いリードタイム用の再学習
㈠ HRES 0.25°解析
8個のGPUを使って、ステップ数8,000で、モデル全体の重みを再学習する。バッチ サイズは GPU毎に1。各反復で、2 つのロールアウト･ステップを実行し、これらの両方のステップを逆伝播する。モデルは、両方のロールアウト･ステップで平均されたMAE損失を、最小化するように最適化される。
この再学習における学習率は、まず1,000ステップに対して線形ウォームアップを使用し、その後は、5×10−5の一定学習率を使用する。重み減衰は、事前学習と同じ値(つまり、5×10－6)を使用し、ドロップ･パスを無効にする(正則化手法を使用しない)。モデルが 2 つのロールアウト･ステップで、メモリに収まるように、事前学習と同様に、エンコーダーとデコーダーのアクティベーション･チェックポイントと勾配シャーディングを使用する。
㈡ HRES 0.1°解析
8個のGPUを使って、ステップ数12,500で、モデル全体の重みを再学習する。バッチ サイズは GPU毎に1。学習率は、1,000ステップに対して線形ウォームアップを使用し、その後 2×10−4の一定学習率を使用する。重み減衰は0 に設定し、ドロップ･パスを無効にする。アクティベーション･チェックポイントを使用し、｢重みと勾配｣にシャーディングを使用する。
㈢ CAMS 0.4°解析
まず、『学習』について。
シングル･ステップ予測(この場合は 12 時間)で『学習』し、バッチ サイズは GPU毎に1に固定する。100ステップに対して線形ウォームアップを使用するが、その後は学習率スケジュールを使用しない。重み減衰も使用しない。ドロップ･パスは無効にする。アクティベーション･チェックポイントを使用して、｢重みと勾配｣にシャーディングを使用する。
次に、『追加学習』について。
この追加学習は 2 つのステップで進行する。最初のステップでは、16 個のGPU を使用して、高学習率†17を使ってステップ数22,000 、次に低学習率†18を使ってステップ数14,500 で 、CAMS 再解析データを追加学習する。CAMS 再解析データから CAMS解析データへの転送を最大限にするために、CAMS 再解析データは CAMS 解析データの解像度0.4°に再グリッド化される。
2 番目のステップでは、8 個の GPU を使用して、高学習率を使ってステップ数7,500、低学習率を使ってステップ数5,500で、CAMS解析データを追加学習する。最終モデルは合計ステップ数49,500で追加学習される。
†17 高学習率とは、新しい汚染変数のみのエンコーダー･パッチ埋め込みでは、1×10−3、ネットワークの残りの部分では 1×10−4を意味する(らしい)。
†18 低学習率とは、新しい汚染変数のみのエンコーダー･パッチ埋め込みでは、1×10−4、ネットワークの残りの部分では 1×10−4を意味する(らしい)。
3⃣ ロールアウト(運用開始)のための再学習
長期にわたるマルチ･ステップ･ダイナミクスを確保するために、AIモデルは通常、ロールアウト専用にモデルを再学習する(らしい)。ロールアウトの再学習には、バックボーンの自己注意(セルフ･アテンション)操作に関係するすべての線形層に、低ランクAdaption(LoRA)層†19を使用する。
さらに、シングル･ステップの再学習と比較してメモリの増加を回避するために、プッシュフォワード･トリック†20を(”表面的に”)使用する。このトリックでは、勾配は最後のロールアウト･ステップを通じてのみ伝播される。ただし、各学習ステップで長いロールアウトを生成することによる遅延を回避するために、深層強化学習で使用される方法と同様に、メモリ内のリプレイ･バッファーを使用してこれを大規模に実行する。各学習ステップで、モデルはリプレイ･バッファーから初期条件をサンプリングし、次のタイムステップでの予測を計算し、この予測をリプレイバッファーに追加する。定期的に、新しい初期条件がデータセットから取得され、リプレイバッファーに追加される。
この手順により、モデルは余分なメモリや速度のペナルティなしで、すべてのロールアウト･ステップで再学習できる。
㈠ HRES 0.25°解析
0.25°20 個の GPU を使用して、ステップ数13,000のLoRA層を再学習する。GPU毎のバッファー･サイズは200。これにより、バッファーの合計サイズは 4,000 になる。データセットのサンプリング期間は 10。モデルが、初期のステップを上手く予測できるようにするため、最初の5,000ステップでは、バッファーに最大 4 日先の予測のみを保持するスケジュールを使用する。4～10日のリード タイムでは、5,000ステップ以降にのみ、バッファーに保持される。学習率は、5×10−5で一定。
㈡ HRES 0.1° 解析
0.1°解析の再学習は 、32個の GPU を使用して行われる。0.25°データの 6.25 倍大きいため、GPU毎のバッファー サイズは20。データセットのサンプリング周期は、0.25°解析と同じ10。学習率は、5×10−5で一定。ステップ数は6,250で、LoRAの重みを使って再学習する。
㈢ CAMS 0.4°解析
この追加学習は、16 個のGPU を使用し行われる。GPU毎のバッファー サイズは 200。データセットのサンプリング周期は、他と同じく10。学習率は、5×10−5で一定。ステップ数は6,500で、LoRAの重みを使って再学習する。
†19 LoRAは、大規模モデルを高速かつ低消費メモリで再学習する手法と位置づけられる。LoRAでは、重み行列の差分が学習対象となる。正確には、この差分行列を”低ランク行列”に分解して、(再)学習が実行される。
†20 プッシュフォワード･トリック自体は、分布シフト問題の解決策(の一つ)である。分布シフト問題は、(少なくとも狭義には)強化学習において現れる問題で、｢学習時に用いるデータの分布と、推論時に用いるデータの分布における差異が、学習を阻害する(かもしれない)問題｣を意味する。本質的には、一般的な教師あり学習においても、学習データと推論時に用いる(本番)データとの分布における差が、推論性能を悪化させる問題、という認識で良いだろう。
ただし、ここでは、あくまで｢メモリ増加回避策｣として、実施されている(という理解で良いだろう)。
【4】比較結果
(0) 予備的なIFSとの比較ー計算コスト
リードタイムが最大 15 日の場合、統合予報システムIFSは、ゴールド スタンダードで最先端の数値予報システムである。ただし、このシステムは相当な計算コストがかかる。10 日間の予報を作成するには、それぞれ 36 コアの 352 個のハイエンド CPU ノードで約 65 分かかり、リードタイムは 1 時間あたり約 5,720ノード秒に相当する。
Aurora は、｢単体の A100 GPU でリードタイム 1 時間あたり約 1.1 秒で予測できるため、IFS に比べて約 5,000 倍の高速化が実現する｣とする。なおGraphCastは、｢単体のGoogle Cloud TPU v4 デバイスで、1 分以内に(正確な)10日間の天気予報を生成する｣と主張している。リードタイム10日(240時間)とすると、Aurora≒264秒、GraphCast≒60秒となり、GraphCastが速い。
(1) 大気汚染ーAuroraとCAMSの比較
0⃣ まとめ
Aurora を CAMS解析データに合わせて追加学習する†21ことで、CAMS 予測と同等か、それを上回る運用予測を、桁違いに低い計算コストで作成できることを実証した(と主張)。
†21 CAMS解析データによるAI基盤モデルの追加学習は、以下の理由から(一般的には、)困難と主張している：㈠CAMS は頻繁に更新され、データ分布に大きな影響を与える。㈡CAMS解析データは 2015 年までしか遡れず、少ない。㈢大気汚染変数は、気象変数とは異なり、大きなダイナミック･レンジを持つ濃度値である。これらの変数は不均一で、疎らに偏っていることがよくある。㈣大気汚染変数は、工場の排出などの人為的要因に大きく依存する。
1⃣ 変数
世界保健機関(WHO)によると、大気汚染の主な原因は、6つの大気汚染物質：一酸化炭素CO、窒素酸化物NO、二酸化窒素NO2、二酸化硫黄SO2、オゾンO3、および粒子状物質(PM)†22である。Aurora は、この6つの大気汚染物質を変数として、モデル化する。
†22 PMは、大きさに応じてPM1、PM2.5、PM10 等と呼ばれる。大気質の警告は、通常、PM2.5 および/または PM10 のしきい値に基づいている。
2⃣ データのセットアップ
2017年10月から 2022 年 5 月までの CAMS解析データでAuroraを追加学習し、2022 年 5 月から 2022 年 11 月までの CAMS 解析データでテストする。CAMS解析データは、時間的範囲が非常に限られているため、汚染物質の動態をよりよく理解するために、2003年1 月から2021年12 月までの低解像度データと、CAMS再解析データ†23も追加学習プロセスに組み込む。
†23 本論文では、｢CAMS 再解析データは、解像度が低く、IFS のかなり古いバージョンを使用しているため、品質が低いと考えられる｣と書かれている。
3⃣ 結果
すべてのターゲット†24の95%で CAMS と競合し(RMSE†2520% 以内)、すべてのターゲットの 74%で RMSEに関して CAMS 予測と同等か、それを上回る運用予測を達成した。また、Aurora は、最上層大気のオゾンと下層大気のすべての種の、12 時間予測において CAMS より優れている。
†24 気象予測･天気予報の文脈では、｢評価対象変数｣を意味する。Ⅰ-0【2】(0)0⃣②を参照。
†25 自乗平均平方根誤差。平均誤差(Mean Error：ME、バイアスとも呼ぶ)と並んで、予報誤差を表す基本的な指標。
(2) 中期気象予測ーAuroraとIFS-HRESとの比較
0⃣ 概要
Auroraは、最先端の数値中期天気予報†26システムである統合予報システムIFS-HRESと同じ空間解像度0.1°†27で中期天気予報を実行し、IFS-HRES の予測精度を初めて上回った(と主張)。これまで、最先端の機械学習ベース天気予報モデル†28の解像度は、0.25°†29であった。解像度を、0.25°から0.1°に引き上げるための最大障害は、高解像度の学習データが不足していることである(と、本論文は主張している)。
†26 中期天気予報は、10日先までの天気予報を指す。
†27 赤道付近において、約11km×11kmの空間解像度に相当する。データ量は、データポイントあたり 1.78 GBに相当する。
†28 例えば、ForeCasNet、Pangu-Weather、GraphCastを指している。
†29 赤道付近において、約30km×30kmの空間解像度に相当する。
1⃣ 変数
変数は{U,V,T,Q,Z}及び10U及び10V、MSL、2Tである。気象学では、東西方向の風をU、南北方向の風をVと呼称する。Tは気温、Qは比湿†30、Zはジオポテンシャル†31である。風速10mのU成分は10U、風速10mのV成分は10Vと呼称される。MSLは平均海面気圧で、MSLPとも表記される。2Tは、地上2mの気温を指し、地上気温と呼ばれる。表記は、T2mやT2などいくつか存在する。
等圧面は、500, 600, 700, 850, 925hPa†32が採用されている†33。リードタイムは、(中期予報なので)12時間～10日である。
†30 湿潤空気中の単位質量当たりの水蒸気の質量。
†31 海水面から、ある高度まで単位質量当たりの空気塊を上昇させるのに必要な仕事量。
†32 500hPaは、海面から上空5,500mの高度に、850hPaは1,500mの高度に相当する。500hPaと850hPaは特に、代表的な等圧面(高度)である。
†33 {U,Q,Z}は500hPaでの比較結果が、詳しく示されている。{T}は、500hPaと850hPaでの比較結果が、詳しく示されている。
2⃣ データのセットアップ
0.25°データは、1950 年まで遡り、解析データ及び再解析データ†34が存在する。0.1°解像度の解析データは、2016 年以降のみ利用可能であり、再解析データは存在していない。 Aurora は、0.25°データで事前学習した後、0.1°データで追加学習することで、高解像度の予測機能を実現できる。2016～2022 年の IFS-HRES解析データで事前学習済Aurora モデルを追加学習することで『0.25°事前学習＋0.1°追加学習→0.1°予測』の実現を実証する†35。
†34 再解析データとは、最新の数値予報システムと過去の観測データを活用して、過去の大気の状況を｢空間3次元+時間の4次元データ｣として再現したものである。
†35 解像度の向上に対処するために、事前学習後に、モデルにいくつかの小さな変更を加えている。評価では、AuroraをIFS解析で初期化し、運用設定を模倣した IFS解析に対して、予測を評価する。
3⃣ 結果
Aurora は、大多数のターゲットでIFS-HRESよりも RMSE が低い。パフォーマンスの向上は、12 時間を超えるリードタイムで顕著であり、長いリードタイムでは、RMSE が最大 60% 減少する。温度や風速成分などで、改善(RMSEの減少)が顕著である。最小の改善は、ジオポテンシャルで見られ、IFS-HRES は高大気レベルで、より正確な予測を出す傾向がある。
(3) 極端イベント
気象学的極端イベントとして、Ciarán(キアラン)を取り上げている。キアラン は、2023年後半に北西ヨーロッパで発生した嵐(暴風雨)であり、英国気象庁によれば、英国で11 月に記録された史上最低気圧であった。機械学習ベースの気象予測モデル(MLWPモデル)は共通して、｢2023 年11月2日 00 UTC に発生する最大 10 m 風速の急上昇を捉えることができない｣との研究成果を受けて、Auroraに対して検証を行った。
結果は、「Auroraは、FourCastNet、Pangu-Weather並びにGraphCastとは異なり、最大 10 m の風速の急激な上昇を正確に予測｣でき、｢IFSの結果とほぼ一致した｣。
(4) グーグルGraphCastとのタイマン勝負
0⃣ 概要
現在 0.25°で最も優れたMLWPとされるグーグル･ディープマインドが開発したGraphCastと、Auroraを比較する。 GraphCastは、0.25°解像度のERA5のみで学習されている。([*2]を読む限り、10日ではないかと思われるが･･･)本論文では、GraphCastのリードタイムは最大5日と述べている。
1⃣ 実験セットアップ
グラウンド･トゥルースは、HRES-T0解析データである(再解析データ†34(再)ではない(はず))。Auroraは、GraphCastの学習データと同じ｢2016～2021年のHRES-T0解析データ(解像度0.25°)｣で追加学習される。GraphCastとAuroraはともに、WeatherBench 2†36によって提供されるHRES-T0データセットの00Z/12Z初期化を使用して評価される。
比較指標には、RMSE及びアノマリー相関係数(ACC)を使用する†37。
†36 さまざまな天気予報モデルを評価および比較するためのフレームワーク。WeatherBench 2 は、オープンソースの評価コード、グラウンドトゥルースおよびベースライン データセットで構成されており、ERA5 データセットの包括的なコピーが含まれている。
†37 先行MLWPである｢FourCastNet、Pangu-Weather及び、GraphCast｣でも、NWPとの比較にRMSEとACCが使われている。
2⃣ RMSEによる比較結果
Aurora は、94% のターゲット変数で、GraphCast と同等かそれ以上の性能を発揮している。AuroraとGraphCastは、GraphCastがロールアウト時に再学習された2 ～ 3 日のリードタイムかつ下層大気†38で、互いに最も似ている。Aurora は上層大気†39において、GraphCast比最大40%RMSEが向上している。さらに、短いリード タイムと長いリード タイムで最大 10 ～ 15% の大幅な改善が見られる。
一方GraphCast は、最大 5 日間、比湿(Q) のほとんどの等圧面†40でわずかに優れたパフォーマンスを示している。
†38 本論文の図5に対する目視(定量的な数値での比較ではなく、図から読み取ったという意味)で、850～1000hPa。850hPa＝1500ｍ上空。1000hPa≒地表。
†39 上層大気では、GraphCastの性能が低いことが、知られているらしい。
†40 本論文の図5に対する目視(同)で、150hPa～1000hPa。
3⃣ 極端な値の予測能力に関する比較
風速や気温などの、地表変数(地表面の物理量)の”極端な値”を正確に予測することは、生命への影響を軽減するための将来計画において、非常に重要である†41。結論を簡潔に述べると、以下の通り：㊀MLWPの予測性能はNWPを上回る、㊁リードタイムが長くなるにつれて、その傾向は顕著になる、㊂風速では、AuroraがGraphCastを上回る。
補足すると、次のようになる。2022 年IFS-HRES 0.25°データセットの 06Z/18Z 初期化を使用して、地表変数分布の極値(裾)を予測するAurora の性能を、GraphCast および IFS-HRES と比較した。
地表面での風速予測では、Aurora が GraphCast よりも優れており、リードタイムが長くなるにつれて、地表変数分布全体で IFS-HRES に対するMLWP モデルの性能が向上する。ただし、2Tは分布の暖かい部分と冷たい部分で動作が異なる。この理由は、次のように説明されている：
冬季の極端現象に対する IFS-HRES の性能には偏りがあり、冬季の予測精度が低下することが報告されている。MLWP モデルには、物理的制約内での単なる前向きシミュレーションではなく、このような偏りのある予測データが含まれているため、IFS などの NWP モデルと比較すると、これらの偏りが現れる傾向がある。
†41 晴原柑九朗(performed by 山下智久;ブルーモーメント)の自論とも符号する。
【5】改善の余地
(1) 確率論的予測
現時点で、Auroraは、決定論的な予測しか生成できない。確率論的予測は、降水量など、その挙動が本質的に確率的である変数の場合に特に重要である。将来的には、モデルを確率バージョンに再学習するか、異なるデータソースで学習された可能性のある決定論的 Aurora モデルのアンサンブルを使用することで、これに対処できる。
(2) さらなるデータの多様性
Auroraは、データの多様性を押し広げたが、それでもグローバル･データセットでのみ学習されている。HRRR†42やCONUS404†43などの多くのローカル高解像度データセットが利用可能であり、そのようなデータを活用することは、将来の研究にとって有望な手段である。
†42 HRRR(the High-Resolution Rapid Refresh)はNOAA(米国海洋大気庁∊商務省)のリアルタイム3km分解能、1時間毎に更新される雲分解、対流を考慮した大気モデルで、3kmレーダー同化を伴う3kmグリッドで初期化される。レーダーデータは15分毎に、1時間の時間間隔にわたるデータがHRRRに同化され、13kmレーダー強化ラピッド･リフレッシュ(※)による1時間毎のデータ同化によって提供される詳細な情報を追加する。出所：https://rapidrefresh.noaa.gov/hrrr/
※ ラピッド･リフレッシュ(RAP)は、NCEP(米国立環境予測センター∊国立気象局∊NOAA)で運用されている大陸規模のNOAA毎時更新同化･モデル化システムである。RAPは北米をカバーし、主に数値予報モデルとそのモデルを初期化する解析/同化システムで構成されている。出所：https://rapidrefresh.noaa.gov/
†43 CONUS404は、米国本土の水文モデルや気象解析に適した、高解像度の水文気象データセット。40年以上にわたって4kmの解像度で、米国大陸(CONUS:contiguous United States or continental United States)をカバーしていることからこの名前が付けられた(つまりCONUSは、アラスカ州とハワイ州を除いた米国になる)。USGS(米国地質調査所∊内務省) Water Mission Areaとの共同研究の一環としてNCAR(米国大気研究センター∊大気研究大学連合←米国立科学財団)が実行したWeather Research and Forecasting(WRF)モデルシミュレーションによって作成された。出所：https://rda.ucar.edu/datasets/ds559-0/
【6】考察
(0) ザックリ言うと、基盤モデルは、やはりスゴい･･･ということになるだろうか。
(1) 本論文で述べられている通り、｢気象予測(予報)でも、基盤モデルが、性能向上に資する｣ことを実証したところが、最大の訴求点である。もう少し正確な記述をするならば、気象予報の基盤モデルAuroraを、適正に再学習した結果、最先端数値予報の精度を凌駕したとの主張が展開されている。数理的に｢Auroraが数値予報を凌駕した理由｣を考察すれば、Auroraの予測が、結果的にアンサンブル予測になっていると考えるしかないだろう。結果的にという文言は、｢多種多様な(事前学習では10種類の)データを使った事前学習＋再学習｣という枠組みを使うと、結果的にアンサンブル予測になるのだろう、という意味で使っている。
👉 本論文(Supplementary MaterialsのG、特にG.8)には、ロールアウト用再学習(LoRA層を使用した再学習)を実施したAuroraの出力(予測)は、アンサンブル予測と見做せる旨の記述がある。
(2) なお、0.1°の再解析データを作成すれば、数値予報の精度がAuroraを逆転するだろうが、それは本質的ではない。再解析データを作成する手間をかけずに、高度な予測･予報ができるという意味が大きい。つまり、適宜、局所的にデータを細かくするなどして、高精度かつ迅速な予測が可能になるはずである。しかも、リアルタイム(に近い)データを、反映させられる可能性もあるので、さらに有望であろう。
(3) データのスケーリング則も示されている。パラメータ2倍で性能5%アップが、どの程度間尺に合うかは、よく分からないものの、スケールで勝ち負けが決まるゲームは、米国好みではあるだろう。
(4) Auroraは多種多様なデータを学習データとするため、必然的に、エンコーダーに工夫を凝らしている。また、事前学習･再学習･追加学習においても、様々なテクニックを駆使している。全体的に、かなり大変な作業であっただろうと思われる。
【0】はじめに
グーグル日本法人は、2024年6月19日、該社公式ブログにおいて、以下のような控えめな発表を行った[*31]：日本において、株式会社ウェザーニューズとのパートナーシップを通して、MetNet-3 のモデルを基に構築された新しい雨量予測モデルを開発した。このモデルにより、最大12時間先まで、5 分毎の降水量予測であるナウキャストが、(検索画面に)表示されるようになる。 MetNet-3の予測においては、継続的で”タイムリー”なデータ供給がキモであるから、ウェザーニューズの貢献は(実は)大きい。
MetNet-3は、MetNet🐾1モデルの第三世代である。グーグル･ディープマインドとグーグル･リサーチの研究者によるMetNet-3を記述した論文[*32](以下、本論文)自体は、2023年7月6日に公開されている(@arXiv)。MetNet-3もMLWPであり、その意味ではⅠ-1～Ⅰ-4で記述したFourCastNet、Pangu-Weather、GraphCast、Aurora(以下、FPGAとする)と同じである。また、主張も、｢数値計算ベースの予測モデル(NWP)を精度で上回る｣なので、同じである。しかし、以下の点で異なる：
㈠ FPGAは中期予報である。MetNet-3は最大24時間先までの”超短期”予報🐾2である。
㈡ FPGAは全球モデルである。MetNet-3は、対象範囲が局所的🐾3なモデルである。
㈢ FPGAのモデル出力は、大気変数を含む。MetNet-3のモデル出力は、降水量と地表変数のみである。
🐾1 meteorology(気象学)＋network、であると合理的に推測できる。
🐾2 気象予測･予報において、中期は最大10日先までを意味する。短期は、最大2日(48時間)先までを意味する。
🐾3 MetNet-3は、米国本土(CONUS：アラスカ州とハワイ州を除く、地続きの米国)を対象としている。
【1】本論文の主張
本論文は、以下のように主張する：
(1) MetNet-3は、ほぼリアルタイムなデータをモデルに反映させることで、降水強度🐾4を｢時間解像度2分、空間解像度1km｣で、地表変数🐾5を｢時間解像度5分、空間解像度4km｣で予測する。
(2) 最大24時間先の予報において、既存NWP(☛【2】(0))より高精度である。
(3) モデルの実行時間は、約1秒。
🐾4 降雨という文言では、雨のみが対象となる。降水には、雨以外の雪、霰、雹なども全て含む。降水強度とは、1時間あたりの降水量(単位:mm/時間)である。もちろん、降雨強度という文言もある(というより、降雨強度の方がメジャー)。ちなみに空港では、層状性の雲から降る雨(地雨)と、対流性の雲から降る雨(しゅう雨)を区別する。
🐾5 ☛【3】(1)2⃣。
【2】事前整理
(0) 既存NWP・・・天気予報モデル
1⃣ HRRR(the High-Resolution Rapid Refresh)
HRRRはNOAA(米国海洋大気庁)の単一メンバー予報モデル(アンサンブル予報ではない、という意味)。空間解像度は、3km。時間解像度は１時間。予報は1時間間隔。1時間毎に更新される雲分解、対流が考慮されている。3kmレーダー同化を伴う、3kmグリッドで初期化される。レーダーデータは、15分毎に1時間の時間間隔にわたるデータがHRRRに同化され、13kmレーダー強化ラピッド･リフレッシュ🐾6による1時間毎のデータ同化によって提供される詳細な情報を追加する。
🐾6 ラピッド･リフレッシュ(RAP)は、NCEP(米国立環境予測センター)で運用されている大陸規模のNOAA毎時更新同化･モデル化システム。RAPは北米をカバーし、主に数値予報モデルとそのモデルを初期化する解析/同化システムで構成されている。出所：https://rapidrefresh.noaa.gov/
2⃣ HREF(High Resolution Ensemble Forecast system)
NOAAの高解像度マルチ･メンバー予報(つまり、アンサンブル予報)。アンサンブル･メンバー数は、10。空間解像度は、3km。時間解像度は１時間。予報は6時間間隔。
3⃣ HRES(High RESolution)
欧州中期予報センター(ECMWF)が提供する、単一メンバー中期気象予測モデル(つまり、アンサンブル予報ではない)。最も正確なNWPと言われている。空間解像度は、0.1°である。0.1°は、7～11kmの空間解像度に相当する(緯度によって変わる)。時間解像度は、1時間。予報は6時間間隔。
4⃣ ENS(Ensemble)
欧州中期予報センター(ECMWF)のアンサンブル予報。アンサンブル･メンバー数は、50。空間解像度は、0.2°である。0.2°は、14～22kmの空間解像度に相当する(緯度によって変わる)。時間解像度は、1時間。予報は6時間間隔。
(1) U-Net
U-Net は元々、バイオメディカル画像セグメンテーション用畳み込みニューラルネットワークとして提案されたモデルである。モデル･アーキテクチャは、オート･エンコーダー(エンコーダー･デコーダー)🐾7。2 次元の畳み込みと逆畳み込みを含む4 つのブロックで構成される、画像間回帰タスクの一般的な選択肢と評される。気象予測においてU-Netを使う価値は以下の通り：U-Netは｢気象場の地域特性を学習できる｣[*33]。
🐾7 もともとオート･エンコーダー(AE、自己符号化器)は、勾配消失と過学習を避けるために、開発された。
(2) MaxVit(Multi-axis Vision transformer)
MaxVitのモティベーションは、長距離相互作用をキャプチャ(捕捉)することである。アーキテクチャ的には、 MaxVitブロック(MB)の連なりと表現できる。MBは、MB畳み込み(MBConv)＋ブロック･アテンション＋グリッド･アテンションで構成される。ブロック･アテンションが短距離相互作用をキャプチャする。グリッド･アテンションで、距離相互作用をキャプチャ(捕捉)する。この｢スケールの異なる相互作用をキャプチャする｣ことが、流体力学が支配する物理現象の解析には必須である。
ブロック･アテンション ≃(オリジナルの)Swinトランスフォーマー、グリッド･アテンション≃ずらし窓の Swinトランスフォーマー。つまり、プロセッサーのコンポーネントは、Pangu-Weatherと同じと考えられる。思想にまで範囲を拡げると、GraphCastにおけるGNNのメッセージパッシングと同じである。
(3) スクイーズ＆エキサイテーション･ブロック
スクイーズ＆エキサイテーション･ブロックは、大域的な情報を処理するために導入される(ので、MaxVitには不可欠)。本論文の文脈で言うと、大域的な情報＝長距離相互作用である。スクイーズ･ブロックでグローバル平均プーリングをかけることで、大域的な情報を処理する(らしい)。
ボトルネック比については、原論文[*34]において、注意機構とゲート機構という概念の中で説明されている。注意機構は、利用可能な計算資源を、信号の最も有益な構成要素に偏らせる手段と解釈できる。これに対して、ゲート機構は、計算効率の高い方法でチャネルごとの関係をモデル化することで、ネットワークの表現力を高めることに重点を置いている。ボトルネック比は、ゲート機構におけるハイパーパラメータで、注意機構とゲート機構のバランスをとる役目を果たしていると考えられる。
※ なお[*34]では、削減比として、2,4,8,16,32が上げられているが、本論文のボトルネック比0.25は、削減比4に該当(1/4=0.25)という理解で良いのだろうか？
(4) 重みのイニシャライザ(初期化子)
深層学習モデルの学習プロセスにおいて重みの分布に偏りが発生すると、学習が安定しない。また、収束も遅くなる。そこで、重みをランダムに初期化する手法が採用される。広く知られている初期化手法は、以下の３つである。
Xavierイニシャライザ(もしくは、GlorotNormalイニシャライザと呼ばれる)は、｢√(2/(n+N))を標準偏差とする正規分布で、ノードの重みを初期化する｣。ここで、nは前層(入力)のノード数、Nは後層(出力)のノード数である。活性化関数がシグモイド関数の場合に使用される、イニシャライザの第一世代。HeNormalイニシャライザは、｢√(2/n)標準偏差とする正規分布で、ノードの重みを初期化する｣。Xavierイニシャライザ の改良版という位置づけで、活性化関数がReLUの場合に使用される。LeCunNormalイニシャライザは、｢√(1/n)標準偏差とする正規分布で、ノードの重みを初期化する｣。
【3】MetNetｰ3の詳細
(1) 前説
1⃣ 概要
MetNet-3は、地形埋め込み、U-Net、MaxVitの 3 つで構成されている。ネットワーク全体では、2 億 2,700 万の学習可能なパラメーターがある。ここで言う｢地形埋め込み｣とは、人手に頼らず、地形情報を(ネットワークが)自動的に検出し、埋め込み処理を自動的に行うことを指している。
2⃣ モデルの変数
変数は、瞬間降水量・降水強度と地表変数である。瞬間降水量は、マルチレーダー/マルチシステム(MRMS)から取得され、2分間隔で推定値が取得される(値は、1時間あたり降水量に換算される)。降水強度は、 レーダーと地上雨量計の両方から取得される(時間間隔は、60分＝1時間)。地表変数は、地表上10mの風速及び風向、地表上2mの温度、地表上2mの露点、である。地表変数はCONUSに広がる、942か所の気象観測所(OMOステーション)から得られる。
なお、モデルの出力(予測する変数)には、同化された気象状態も含まれる。
3⃣ データセット
利用可能なデータは、2017 年 7 月から 2022 年 9 月までの期間にわたる。データには、過去 90 分間のレーダーデータ(降水量)、過去6 時間の OMOステーションからのデータ、GOES衛星🐾8からの画像、同化された気象状態、緯度と経度の情報、高度情報、現在の時刻が含まれる。学習、検証、テストの各データ セットは、重複することなく、順番に期間から生成される。
19 日間の学習データ、1 日間のブラックアウト、2 日間の検証データ、1 日間のブラックアウト、2.5 日間のテスト データ、1 日間のブラックアウトの連続期間を使用して、ブラックアウト期間にはサンプリングを行わずに、それぞれ学習、検証、テスト データをサンプリングする。
学習サンプルの数を増やすために、正確なリードタイムの観測が利用できない場合はいつでも、線形補間を使用してトレイン スプリット内のターゲットを一時的に補間する。空間的には、ターゲット パッチは、経度と緯度で 0.5 度間隔で CONUS 領域上のグリッド上の交差点からランダムにサンプリングされる。
地表変数はOMO ステーションの測定値を取得し、ステーションがある 4 km x 4 km ピクセルにマッピングする。特定の地域に複数のステーションがある場合は、それらの測定値の平均を取得する。 OMOステーションは、757 の学習ステーションと 185 のテスト･ステーションの 2 つのグループに分けられる。テスト･ステーションのデータは学習中には一切使用されず、テスト ステーションの結果のみが報告される。
🐾8 NASAとNOAAの共同プロジェクトである、静止軌道からの実用気象観測衛星。GOES衛星は、1975年以降センサーの改良を重ねながら、地球の気象現象の観測･予測、大気の測定、太陽活動の監視等に関する画像とデータを提供し続けている。出所：https://www.restec.or.jp/satellite/goes-series-all.html
(2) ハイパーパラメータ等
オプティマイザ・・・AdamW(β1🐾9＝0.9🐾10、β2🐾11＝0.999🐾12)
学習率・・・8.0×10－5
重み減衰率・・・0.1
活性化関数・・・MBConv内＝GELU、その他＝ReLU
Polyak減衰率(忘却係数)🐾13・・・0.9999
バッチサイズ・・・64
学習ステップ数・・・26万
再学習ステップ数・・・8万
🐾9 勾配の移動平均の減衰率。
🐾10 β1の値としては、いわゆる｢規定値(デフォルト)｣である。ほとんどのタスクにおいて、既定値で良い結果が得られる、とされる。
🐾11 勾配の二乗の移動平均の減衰率。
🐾12 β2の値としては、いわゆる｢規定値(デフォルト)｣の一つである。0.9 or 0.99 or 0.999が使われる。それぞれ、平均化の長さ＝10 回、100 回、1000 回のパラメーター更新に対応している。
🐾13 一般的な理解として、Polyak平均と忘却係数は、サンプリングのばらつきを抑制するために、用いられる。忘却係数をαとすると、t+1番目のサンプル＝(1－α)×t番目のサンプル＋α×Polyak平均である。
(3) MetNet-3における”修正された”U-NetとMaxVit
1⃣ U-Net
U-Netのダウンサンプリング･パスでは、2つの畳み込みResNetブロックを適用し、4 kmと8 kmの解像度レベルで最大プーリングを使用して2倍にダウンサンプリングする。アップサンプリング･パスでは、16 km レベルと 8 km レベルの両方でカーネル (2, 2) とストライド (2, 2) を使用した、逆畳み込みを使用してアップサンプリングし、2 つの畳み込み ResNet ブロックを適用する。4 km 解像度から 1 km 解像度へのアップサンプリングは、4 x 4 ピクセルの正方形で各アクティベーションを繰り返し、2 つの ResNet ブロックを再度適用することによって実行される。
2⃣ MaxVit
12個の修正されたMaxVitブロック(MB)を使用する。元のアーキテクチャと比較して、次の変更を導入した。元のMaxVitアーキテクチャに存在していたMLPサブ･ブロックを削除し、正規化されたキーとクエリを使用し、各MBの出力からMaxVitの出力への｢スキップ接続｣を導入した。すべてのアテンションウィンドウのサイズは8 x 8で、32個のアテンションヘッドを使用する。MBConvは、拡張率🐾144と、ボトルネック比(削減率)0.25のスクイーズ＆エキサイテーションを使用する。
🐾14 拡張率とは、MBConvの最初の畳み込みでチャネル数を何倍にするか、を表す係数のことらしい。また、4は、拡張率として通常使われる範囲の値(らしい)。出所https://qiita.com/omiita/items/1d96eae2b15e49235110
(4) リードタイム･コンディショニング
MetNet-3では、リードタイムは2分間隔で、0～24時間を表す0～721までのインデックスを持つワンホット埋め込みとして符号化し、連続的な32次元表現にマッピングされる。ただし、気象予測タスクは、リードタイム🐾15が長くなるにつれて著しく難しくなり、学習に悪影響を与える可能性がある。これに対処するため、MetNet-3では、t = 24時をt = 0時の10倍少ない頻度でサンプリングする。この偏ったサンプリング･スキームにより、サンプリング頻度が低い、長いリードタイムを含む、すべてのリードタイムで結果が改善される。
偏ったサンプリング･スキーム実現のために、MetNet-3では、リードタイムに加法乗法コンディショニング(FiLM)が適用されている(と理解している)。今の文脈で説明するとFiLMは、リードタイムに応じたアフィン変換を適用することで、アクティベーション分布🐾16をリードタイムごとに異なるにように調整する。より具体的には、リードタイムを入力として｢スケーリング係数｣および｢シフト係数｣を出力する関数を、それぞれ学習していると推量される。線形結合の文言を使うと、スケーリング係数は｢傾き｣で、シフト係数は｢切片｣である。
🐾15 気象予測の文脈におけるリードタイムとは、予測間隔を意味し、モデルに入力されるデータが取得された時刻と、モデルの出力時刻(例：明日9時の雨量であれば、”明日の9時”が、出力時刻)との差がリードタイムである。
🐾16 アクティベーション分布＝隠れ層各ノードにおける出力の分布。
(5) 損失関数及び、損失関数の再スケーリング
損失関数には、交差エントロピーを用いる。ただし、HRRR 同化状態の予測に対してのみ、平均二乗誤差(MSE)を損失関数とする。これらの損失関数は、大きさが大きく異なる可能性があるため、大きさが同程度になるように再スケーリングする必要がある。
再スケーリング手法は、各変数の平均がほぼ 0、標準偏差が 1 になるように、MSE損失のすべてのターゲットを再スケーリングする方法以外に、各入出力サンプルの勾配を動的に再スケーリングする新しい手法も導入する。より正確には、各サンプルのモデル出力に対する MSE 損失の勾配を計算した後、サンプルの勾配の全体的な大きさを変更せずに、各出力チャネルで同じ L1 ノルムを持つように再スケーリングする。このスケーリングにより、各出力チャネルの影響が制限され、したがって、ターゲット チャネルのごく一部が破損したとしても、モデルへの影響が制限されることが保証される。
(6) 正規化・正則化・初期化
1⃣ 正規化
ネットワーク全体で、プレ･レイヤー正規化を使用する(つまり、サブブロックの入力にレイヤー正規化を適用する)。また、指定されたサブブロックの最後の畳み込みではない、各畳み込みの後にレイヤー正規化を適用する。
2⃣ 正則化
各サブモジュールの出力を残差ブランチに追加する前に、各ResNetブロックの最初の畳み込みの後に、ドロップアウトを適用する。ドロップアウト率は0.1。またMaxVitでは、ドロップパス(あるいはstochastic depthと呼ばれる)も使用する。ドロップパス率は、ネットワーク全体で、0から0.2まで線形に増加する。
ドロップアウトは、中間層のノード出力を一定の割合でランダムに0とすることで、結合を欠落させる正則化手法。これは、横方向(幅方向)にネットワークを小さくしていると考えられる。これに対し、ドロップパスは、層数を変更することで縦方向(深さ方向)にネットワークを小さくする正則化手法である。
3⃣ 初期化
LeCunNormalイニシャライザ(初期化子)を使用する。さらに、MaxVit 内の各サブ･ブロックの最後の線形層の初期化を 1/√N で再スケールする。ここで、N は指定された残差接続で、出力が追加されるサブ･ブロックの数。
(7) 高密度化(デンシフィケーション)
入力データが取得できない地点でも、予測値を出力できるようにすることを、MetNet-3では、高密度化(デンシフィケーション)と呼んでいる。これを、データ同化と呼んでも間違いではない🐾17。地表変数のグランドトルゥースは、CONUSに疎らに配置された気象観測所(OMOステーション)でのみ観測される。地表変数をOMOステーションが存在しない場所でも予測できるようにするため、学習中に各 OMOステーションを｢25% の確率で、ランダムにマスクする｣🐾18。これにより、特定の場所に入力 OMO 変数がない場合でも、モデルが OMO 変数を予測するように学習される。
この高密度化は、マスク言語モデリングあるいは、画像識別タスクにおけるcutout等のアナロジーとして捉えることが可能であろうか。
🐾17 データ同化の主目的の一つとして、｢データがない場所や時間の物理量を求める｣がある。
🐾18 この処置は、ホールド･アウトとは別のプロセスであることに注意。
(8) 再学習
降水量と地表変数の予測の品質には、トレードオフがある。このトレードオフを利用した。最初に降水量に使用するモデルを学習し、その後、降水量モデルと比較して地表変数における損失の重みを 100 倍に増やしてモデルを再学習する。さらに、この地表変数モデルでは地形埋め込みを無効にする。これは、地形埋め込みは、異なる場所間の転送を妨げる可能性があるため。
【4】比較結果
(0) 前説
1⃣ 比較対象モデル
HREF及びENSと比較する。HRRR、HRESは参考用という位置づけ。4つをまとめてベースラインと呼んでいる。HREFとHRRRはCONUSを対象とするローカルモデルである。ENSとHRESはグローバルモデル。HREFとENSはアンサンブル予報であり、HRRRとHRESは単一メンバー予報である。
2⃣ 比較指標
連続ランク確率スコア(CRPS)、重要成功指数(CSI。スレット･スコアとも呼ばれる)、および平均絶対誤差 (MAE) 。 CRPS は、確率的予測が関与する場合に、最も広く使用される指標の1つ(であり、アンサンブル･モデルに対して、より優れた評価指標)である。故に、｢HREF及びENSとの比較に特に適している｣(と、本論文に書かれている)。
CRPS → 予測値が確率的である場合に、MAEを拡張した指標
CSI＝TP/(TP＋FP＋FN) → CSIが高い＝稀な気象を適中させる能力が高い
(1) 結果
1⃣ 降水量
① 瞬間降水量
まず、CRPSで評価する：MetNet-3のCRPSがENSのCRPSよりも、24 時間の全リードタイム範囲にわたって、小さい。t=1時でCRPSの差は最大である。t=24時に近づくにつれて、差は縮まっていく。CRPSの経時変化を曲線と見做すと、その形状は、t=1～4時を除いて、MetNet-3とENSでほぼ同じ。
次に、CSIで評価する：MetNet-3 は、最初の15 時間のリードタイムで小雨(1 mm/h)に対して ENS よりも優れている。t=1時でMetNet-3 のCSIは0.5を越えている(CSI＝1で完全な予報)。t=15時でCSI＝0.275程度。
中雨(4mm/h)に対しては、最初の20時間のリードタイムでENSより優れている。ちなみに、t＝1時でMetNet-3 のCSIは0.3程度(当然、小雨より低い)。大雨(8mm/h)に対しては、24 時間の全リードタイム範囲でENS よりも優れている。ちなみに、t＝1時でMetNet-3 のCSIは0.2程度。
② 降水強度
まずCRPSで評価する：リードタイムの最初の19 時間まで、ENSとHREFの両方を上回っている。t=1時でCRPSの差は最大である。t=19時に近づくにつれて、差は縮まっていく。ちなみに、曲線形状はt=5時以降は、ほぼ同じ。
次にCSIで評価する：小雨(1mm/h)最初の14時間のリードタイムでENSより優れている。t=1時でMetNet-3 のCSIは0.65を越えている。t=14時のCSIは0.3程度。中雨(4 mm/h)は、最初の18 時間のリード タイムまでベースラインを上回っている。t=18時のCSIは、0.15程度。大雨(8mm/h)は、最初の19時間のリードタイムまでベースラインを上回っている。t=19時のCSIは、0.1を下回っている。
2⃣ 地表変数
MetNet-3 は、24 時間の全リードタイム範囲で、すべての地表変数について、ENS よりも優れた CRPS と MAEを達成している。ただ、t=0→24時に対するCRPSとMAEの”悪化度”は、MetNet-3の方が急である。長いリードタイムに対して、MLWPの予測が、本質的に困難であることを示唆しているのだろう。
さらに本論文では、異なる指標を使って、MetNet-3がENSよりも優れていることを主張する。具体的には、MetNet-3 とENSによる気温と風速の予測分布を、グランドトルゥースと比較している。グランドトルゥース＝OMOステーションの観測値で、リードタイム範囲は24時間である。観測値のほとんどは、MetNet-3 の80% 信頼区間🐾19に収まるが、ENS の80% 信頼区間には収まらない。つまり、ENSは、アンサンブル予報として分散が不十分で、MetNet-3に劣るという主張である。
🐾19 予測の10番目と90番目の四分位数の範囲を意味する。
【5】考察
(0) 24年7月から、利用可能なので、どういうパフォーマンスを示してくれるか楽しみ。
(1) タイムリーなデータを高速で処理することによって、MLWPがNWPを越えたという内容。それ自体にサプライズはない。ただ、この枠組みでは、線状降水帯によるゲリラ豪雨の｢降水量予測｣はできないだろう。線状降水帯によるゲリラ豪雨のデータがない(だろう)し、おそらく降水データ以外の物理データを併せないと、高精度な予測は難しいだろう。
(2) MetNet-3における工夫としては、リードタイム･コンディショニングが、肝要であると思われる。
【0】はじめに
グーグル🐾1・米MIT・欧州中期予報センター(ECMWF)の研究者は、｢1～10日間の予報で機械学習モデルに匹敵し、1～15日間の予報で欧州中期予報センターのアンサンブル予測に匹敵し、数十年にわたる気候指標を正確に追跡できる｣モデルを開発した、と主張する論文[*37](以下、本論文)を発表した(24年7月22日@nature)。当該モデルは、NeuralGCMと呼ばれる。
気象予測の分野で、これまで世間を賑わせてきた機械学習ベースのモデル(FourCastNet、Pangu-Weather、GraphCast)🐾2は、あくまで天気予報モデルであった。つまり、粗っぽく言うと、中期(つまり最大10日先)の風速や気温を予測しようとするモデルであった。これに対して、NeuralGCMは、大気(循環)モデルである。大気モデル(General Circulation Model：GCM)は、地球全体を対象として、地球大気のダイナミクスを表現するモデルである。解析期間が長期にわたる地球温暖化傾向等をシミュレートするような場合には、GCMが必要になる。
NeuralGCMは、数値解析モデルとニューラルネットワーク･モデルを連成させたハイブリッド･モデルである。ハイブリッド･モデルは、無理ゲーであるパラメタリゼーション(☛【2】(1)2⃣)を回避できるが、別の障害が存在した。NeuralGCMは、その障害を首尾良く回避できた(☛【2】(1)2⃣及び【3】(6)1⃣)ので、良い成果に繋がった、と考えられる。
🐾1 正確には、グーグル･リサーチとグーグル･ディープマインド。
🐾2 AuroraとMetNet-3は、毛色が違う。ざっくり言うと、Auroraは｢基盤モデル＋リードタイムに合わせた再学習｣という枠組みを採用したモデル。MetNet-3は、降水量の超短期予測モデルである。
❚為参考❚
スタンフォード大学人間中心AI研究所(HAI)は、2017年から毎年｢AI Index Report｣という調査報告書を公開している(2020年はコロナ禍のため、例外)。最新の2025年版第5章科学と医療[*54]に、｢注目すべきモデル｣という囲み記事がある。9個のモデルが紹介されているが、その内5個が、気象関連である。5個の中にはNeuralGCMも含まれている。
AI Index Report2025では、次のように簡潔に紹介されている(ドメインは、weather forecasting)：NeuralGCMは、微分可能な物理ベースのソルバーと、気象と気候の両方をシミュレートする機械学習コンポーネントを組み合わせたハイブリッドモデルである。NeuralGCMは、短期･中期予測において、主要な機械学習ベース･モデルや物理ベース(数値計算ベース)･モデルと同等かそれ以上であり、数十年にわたって気候指標を正確に追跡し、熱帯低気圧のような複雑な現象を捉える。
【1】本論文の主張
本論文は、以下を主張する：NeuralGCMは、
(1) 最先端の物理ベース･モデル(具体的には、ECMWF-ENS)よりも優れたCRPS🐾3を備えた、正確なアンサンブル天気予報を行う最初の機械学習ベースのモデルである。
(2) 全球雲解像大気モデル(具体的には、X-SHiELD)に匹敵する空間バイアスを実現し、現実的な熱帯低気圧の進路をシミュレートし、現実的な過去の気温傾向を使用してAMIP🐾4のようなシミュレーションを実行できる最初のハイブリッド モデルである。
(3) 既存の機械学習ベースモデル(具体的には、Pangu-Weather、GraphCast)よりも、｢地衡風と非地衡風の鉛直構造及び、その比率｣を正確に描写する。
(4) 従来の全球モデル(GCM)より8～40 倍粗い水平解像度で実行しても、同等程度の精度を達成する。このため、計算リソースを 3 ～ 5 桁節約できる。
🐾3 CRPS＝連続ランク確率スコア。詳細は、【3】(9)を参照。
🐾4 AMIP＝Atmospheric Model Intercomparison Project：大気大循環モデル相互比較実験。MIP(モデル相互比較プロジェクト)とは、モデル誤差を低減するために、シミュレーション結果と観測データとの比較検証、並びにシミュレーション結果間の比較検証を行うプロジェクトであった。1991年に、大気循環モデル(Atmospheric Model)で始まったので、AMIP(Atmospheric Model Intercomparison Project)と呼ばれた。
【2】事前整理
(1) 全球モデル(Global Spectral Model：GSM)
1⃣ 概要
予報モデルには、用途によって、GSM(大気循環モデルとも呼ばれる)、全球アンサンブル予報モデル、領域モデル🐾5等🐾6が存在する。これは、日本のみならず世界共通である(為念)。GSMとは、地球全体を予報領域とした数値予報モデルである。数日より先の予報には、地球全体をカバーするGSMを考慮する必要がある。なぜなら、日本から遠く離れた地域における大気の状態であっても、数日後には日本に影響を与えるからである。静力学平衡🐾7の仮定をした静力学方程式系を基礎方程式として、スペクトル･モデル(☛下記(2)を参照)を採用して、日本では1988年に実用化された。日本のGSMは水平格子間隔が約13km(2023年3月～)、鉛直層は128層(2021年3月～)に分割し、予報期間は11日間である。GSMの予測精度が良いと言われているのは欧州中期予報センター(ECMWF)、米国、英国及び、日本らしい[*38]。ちなみに、米国のGSMは13km間隔・127層・予報期間16日間、欧州中期予報センター(ECMWF)は9km・137層・10日間である。
🐾5 少なくとも日本の領域モデルには、局地モデル(Local Forecast Model：LFM)、メソモデル、メソ･アンサンブル予報システム(Meso-scale Ensemble Prediction System：MEPS)が含まれる。
🐾6 日本だと他には、季節アンサンブル予報システム(JMA/MRI-CPS2)がある。
🐾7 静力学平衡 (静水圧近似とも呼ばれる)とは、鉛直の気圧傾度力＋重力を0とおく仮定(近似)である。発達した積乱雲等でなければ、かなりよい精度で成り立つ(出所：https://www.jma.go.jp/jma/kishou/books/nwptext/45/1_chapter4.pdf)。静力学平衡の下では、鉛直速度の時間変化率
を予報する必要がないため、計算量が少なくなる。気圧傾度力とは、気圧の高いところから低いところに(空気塊を、)動かす力を言う。
2⃣ パラメタリゼーション[*39]
GCM内部では、｢雲の形成、放射輸送、降水｣等の未解明物理過程(☛下記(2)を参照)は、パラメタリゼーションによって表される。グリッド･サイズ未満の捕捉できない物理量をa、グリッド上の物理量をbとする。言い換えると、aは未解明物理過程における物理量である。bを使って、aによってbに起きる変化を定式化することを、パラメタリゼーションと呼ぶ。つまりパラメタリゼーションとは、bで捕捉できない物理量をbで定式化するという”無理ゲー”である。そのため、GCMにおける誤差要因の一つである。
本論文のNeuralGCMは、データ(ERA5)を使い、未解明物理過程をニューラルネットワークで表現することによって、この無料ゲーを避けている。そう解釈することができる。
(2) 力学過程と物理過程[*40]
気象予報の文脈における力学過程(力学フレームとも呼ばれる)とは、ザックリ言うと、流体力学に基づく方程式を解く部分である。噛み砕くと、｢数値予報モデルの基礎方程式に含まれる移流や気圧傾度力の時間変化率(＝時間微分計算)を求める部分｣と｢実際に時間積分を行うところをあわせた部分｣を指す。
力学過程以外の｢外力、非断熱加熱、相変化に伴う加湿の効果を計算する部分｣と｢それらの計算に必要な大気以外とのやりとりや内部的な変化を考慮する部分｣等は物理過程と呼ばれる。[*41]では、｢放射過程、乱流過程、雲過程｣という分解が成されている。
(3) スペクトルモデル(Spectral Model)[*42],[*43]
GSMはGlobal Spectral Modelの略語であった。Spectral Model(スペクトルモデル)＝スペクトル法を採用したモデルであり、スペクトル法＝水平方向の離散化手法の一つ、である。数値解析における普通の感覚だと、水平方向の離散化とは、解析対象を碁盤目状に切って格子を作り、連続的な物理量(気象学的物理量あるいは、予報変数)の適当な平均値を、格子点の値としてしまう。これは、格子点法と呼ばれる。
GCMは、解析対象である地球を平面ではなく、球面と(正しく)みなす。球面とみなした上で、予報変数を球面調和函数で展開する🐾8。直感的な表現を使うと、波の重ね合わせとして、予報変数を表現する。スペクトル法では時間微分計算を解析的に実行できるため、力学過程(☛上記(2)参照)の計算精度は、スペクトル法の方が高い。
🐾8 経度方向は、フーリエ級数で表現される。緯度方向は、ルジャンドル陪関数で表現される。
(4) スペクトルモデル のT○○、TL○○という表記[*42],[*43]
スペクトル法(スペクトルモデル)では、波の重ね合わせとして、予報変数を表現した。従って、スペクトルモデルでは、重ね合わせる波の数を増やせば、近似精度が上がる。これは格子間隔を狭くすれば近似精度が上がることと同じ理屈である。球面調和関数の最大合計波数を切断波数と呼ぶ。重ね合わせる波の中で、波長が最短の波の数が、切断波数となる。スペクトルモデルでは、例えば、T213やTL319のように、この切断波数を明示することで、モデル解像度を表現する。213や319が切断波数である。
Tは、東西波数及び、全波数(＝東西波数＋南北波数)の空間で、波数の切断を三角形(Triangular)型で行う、ことを意味している。TLは、結果的に、セミ･ラグランジュ法を採用していることを指す(Lは線形項のみを考慮したこと🐾9を意味する)。セミ･ラグランジュ法は、クーラン条件(CFL条件)による積分時間間隔の上限を回避するために開発された。セミ・ラグランジュ法では、時間Δt後に、｢風上側の｣どの情報が運ばれてくるかを考えて、出発点を決める。つまり、出発点の値が、Δt後の格子点の値となる。
🐾9 セミ･ラグランジュ法では、非線形項である移流項が時間微分の中に隠れる。そのため、非線形項を考慮することなく、結果的に線形項のみを考慮する。
(5) ガウス格子(ガウス･グリッド)[*42],[*43]
スペクトルモデルは、｢予報変数を球面調和函数で展開する｣が、これは力学過程の話である。スペクトルモデルでも、物理過程は格子点法を用いて計算する。それは、｢陸面・海面の効果や雲の生成｣等は、｢波の重ねあわせと考えるよりは、局所的な効果や変化と考えたほうが都合が良い｣からと説明される。つまりスペクトルモデルでは、物理過程に現れる物理量は、格子上で計算された後、計算結果を、波数毎の物理量に変換するという処理が行われる。この、スペクトルモデルにおいて物理過程を計算するために用いる格子を、ガウス格子と呼ぶ。
気象庁のGCMでは、適合ガウス格子(reduced Gauss grid)が用いられている。ECMWFのGCMでは、三次八面体逓減ガウス格子(cubic octahedral reduced Gaussian grid)なるガウス格子を使用している。適合ガウス格子では、精度に影響がない範囲で、中高緯度の格子が標準ガウス格子(regular Gauss grid)よりも少なくなっている。このため、物理過程の計算量が減少するというメリットがある。
【3】NeuralGCMの詳細
(1) 概要
1⃣ モデル･アーキテクチャ
NeuralGCM では、力学コアで考慮されていない物理過程の影響と、力学コア内の計算エラーがニューラル ネットワークによって近似される。数式的に説明すれば、ニューラルネットワークで表現された物理過程(と計算エラー)が、プリミティブ方程式(☛(3)2⃣参照)の右辺に付加される。これを便宜上、修正されたプリミティブ方程式と呼ぼう(本論文では、そういった文言はない)。プリミティブ方程式の左辺は予報変数の時間微分であるから、予報変数の時間発展を駆動する、”考慮されていなかった(ある意味、未知の)”外力あるいはポテンシャルみたいなものが、正しく付け加えられた、と解釈される。修正されたプリミティブ方程式の結果は、ニューラルネットワークに還流するから、力学コアとニューラルネットワークは、連成(coupling)していることになる。
改めて、構成的・構造的な説明をすれば、NeuralGCMは、㊀力学コアと、㊁物理モジュールという、2つのコンポーネントで構成される。㊁が、上記でニューラルネットワークと言い表したものである。言わずもがな、㊀が力学過程を担当し、㊁が物理過程を担当する。そして、重複するが、｢㊀と㊁が連成している｣ところにNueralGCMの特徴がある(☛2⃣を参照)。
㊀ 力学コアは、離散化された支配的な動的方程式を解くコンポーネントであり、重力とコリオリ力の影響下で大規模な流体運動と熱力学をシミュレートする。㊁物理モジュールは、ニューラル･ネットワーク(NN)を使用して、｢雲の形成、放射による熱輸送、降水、サブグリッド･スケールのダイナミクス🐾10｣等の未解明物理過程が、シミュレーション結果に与える影響を推測する。
力学コアの詳細は、下記(3)を参照。物理モジュールの詳細は、下記(4)を参照。
🐾10 グリッドのサイズよりも細かい現象は、捉えることができない。サブグリッド･スケールのダイナミクスは、その意味での、｢未解明｣物理過程である。
2⃣ 特徴ーなぜ、NeuralGCMはoutperformeするのか？
一言で述べれば、GCMの大きな誤差要因である｢未解明物理過程｣を、ニューラルネットワークで精度よく近似できたから、であろう。もちろん、高品質かつ豊富な再解析データERA5の存在が大きい。
機械学習モデルとGCMを組み合わせたハイブリッドモデルは、これまで、未解明の物理過程を力学過程とは無関係に学習していた(オフライン学習)。その後、機械学習モデルが、既存の GCMに挿入される。しかし、学習中に機械学習モデルと支配方程式が適当に連成していないと、数値的不安定性(イニシャル･ショック)や気候ドリフト🐾11などの深刻な問題が発生する可能性がある。
NeuralGCMはハイブリッドモデルにおいて障害となりうる数値的不安定性(イニシャル･ショック)や気候ドリフトを、うまく回避したハイブリッドモデル、と評価することができるのだろう。回避策の一つとしては、損失関数の工夫があげられるだろう(☛(7)及び(8)を参照)。イニシャル･ショックの回避策は、エンコーダーにおいて工夫がなされている(☛(6)1⃣を参照)。また、【5】気候シミュレーションに関する評価結果、も参照。
🐾11 合成の誤謬という言葉は、それほど適当ではないが、ザックリと言えば、そういう意味である。今の文脈では、異なるモデルを結合させた場合に、非現実的な結果が出力されることを指す。例えば、大気モデルと海洋モデルを結合させたモデル(いわゆる大気･海洋結合モデル)が、現実的な気候を再現できないような場合を指す。
3⃣ モデルのバリアント
決定論的NeuralGCMと確率論的NeuralGCMがある。確率論的NeuralGCMとは、つまり、アンサンブルモデル(NeuralGCM-ENS)である。決定論的モデルには、水平解像度(2.8°、1.4°、0.7°)に応じて、3種類のバリアントが存在する。アンサンブル･モデルの水平解像度は、1.4°である。
水平解像度2.8°、1.4°、0.7°は、それぞれTL63、TL127、TL255が対応する(☛TL○○の意味は【2】(4)を参照)。
(2) ニューラルネットワーク(物理モジュール)の学習
1⃣ 入出力
ニューラルネットワークへの入力には、発散(の鉛直構造)、渦度、風ベクトル、温度、比湿、雲氷混合比🐾12、雲水混合比🐾13、地表気圧(の対数値)が含まれる。出力は、関連する大気状態の変化率である。
🐾12 湿潤空気中の単位質量当たりの雲氷の質量。湿潤空気＝乾燥空気＋水蒸気＋雲水＋雲氷＋雨＋雪。雲に含まれる水分は、水and/or氷である。
🐾13 湿潤空気中の単位質量当たりの雲水の質量。
2⃣ 決定論的モデルの学習
NeuralGCMは、ERA5からサンプリングされた最大5日間の予報変数を予測するように学習されている。NeuralGCMモデルを 2.8°、1.4°、0.7° の解像度で学習するために、ERA5 データを対応するガウス格子に再グリッド化🐾14する。最終的な 2.8° および 1.4°モデルは、1979 年から 2017 年までのデータで学習され、2018 年のデータで評価された。 0.7°モデルは、1979 年から 2019 年までのデータを使用して学習された。
学習中は、学習の反復の関数としてリードタイムの長さも変更し、学習の過程でモデルが改善されるにつれて予測範囲を広げる。
🐾14 相対的な領域の重なりによって重み付けされた寄与を線形に集計する、保守的な再グリッド化スキームを使用する(らしい)。
3⃣ 確率論的モデル(アンサンブル･モデル)の学習
学習した空間的および時間的相関を持つガウスランダム場からサンプリングすることにより、学習済エンコーダーと学習済物理モジュールに、ノイズを注入する。学習では、予測ごとに2つのアンサンブルメンバーを生成する。
確率論的モデルは、1979 年から 2019 年までのデータを使用して学習された。
(3) 力学コア
1⃣ 離散化
力学コアは、計算領域を離散化するためにガウス･グリッド(【2】(5)参照)とσ座標🐾15を使用する。ガウス･グリッドを使用すると、グリッド空間表現と球面調和関数基底間の高速で正確な変換が可能になる。その結果、ガウス積分法によって定義される等角経線と不等間隔の緯度が得られる。
本論文では、さまざまな水平解像度(2.8°、1.4°、0.7°)で予測を行う一連のモデルを学習した。これらは、TL63、TL127、TL255 に対応する(☛TL○○は【2】(4)を参照)が、モデル方程式を解く際は、T62、T125、および T253 を使用する(☛T○○は【2】(4)を参照)。これらは、エイリアシング🐾16エラーを回避し、Google TPU でコストがかかる 128 の倍数を超える配列次元の増加の必要性を最小限に抑える。また、すべてのモデルは、鉛直離散化に 32 の等距離σ座標を使用する。
🐾15 気圧座標p(x,y,z,t)を地表面気圧ps(x,y,t)で規格化した変数σ＝p/psは、z座標の1価関数であるから、鉛直座標として用いることができる。これがσ座標である。
🐾16 エイリアシングとは、格子数の不足により高波数が低波数に投影される現象。
2⃣ プリミティブ方程式
NeuralGCM の力学コアは、｢㊀運動量方程式、㊁熱力学の第 2 法則、㊂熱力学状態方程式(理想気体)、㊃連続の式、および ㊄静力学平衡🐾7(再)↑｣の組み合わせを表す、プリミティブ方程式を解く。方程式を解くために、水平方向の風に対する発散渦度表現を使用し、次の 7 つの予報変数の方程式を導く: ❶発散 δ(∇と速度ベクトルの内積)、❷渦度 ζ(∇と速度ベクトルの外積)、❸温度 T、❹地表気圧psの対数値、および 3 つの水分種(❺比湿q、❻雲氷混合比qci、❼雲水混合比qcl)。
パラメータは、鉛直軸に平行な上向きの単位ベクトル k、水平速度ベクトル u🐾17、コリオリのパラメータ f、理想気体定数 R、定圧比熱容量 Cp、κ = R/Cp、σ座標での診断された🐾18鉛直速度 σ˙、診断された流体粒子の圧力変化 ω ≡ dp/dt、診断されたジオポテンシャル Φ、診断された仮温度 Tνである。
🐾17 ❶発散 δ、❷渦度 ζ及び、鉛直軸に平行な上向きの単位ベクトル k(並びに∇とラプラシアン∆)を使って表される。
🐾18 基礎方程式に基づき、時間積分されて未来の値が求められる量を｢予報変数｣と呼ぶ。気象予報の文脈では、予報変数から時間積分を経ずに求められる量を｢診断(された)量｣と呼ぶ(ことがある)。
(4) 物理モジュール
物理モジュールでは、関連する大気状態の変化率が、2つのステップで計算される：㈠1⃣特徴量の抽出と2⃣正規化、及び㈡3⃣ニューラルネットワークでの演算処理と4⃣再スケーリング、である。
1⃣ 特徴量の抽出
NeuralGCM では、メイン･ネットワークによって入力機能として使用される有用な特徴量を抽出することを目的とした 2 つの埋め込みモジュール：鉛直埋め込みネットワークと表面埋め込みネットワーク、を使用する。鉛直埋め込みネットワークは、大気の鉛直構造全体にわたって共通の特徴量を抽出することを目的とする。表面埋め込みネットワークは、大気と表面(陸地、海面、海氷面)との境界の状態を推定することを目的とする。これを実現するために、海面温度(SST)と海氷密接度(sea ice concentration)🐾19など、追加の表面関連入力を受け取る。
従来の気象予報プロトコルでは、提供されるSSTと海氷密接度は、予報開始の最初の時刻の前日から取得された値を使用して、各予報を通じて一定のままであった。このため、モデルは予報の初期化時に利用可能なデータのみに依存するようになる。NeuralGCMでは、ERA5データからSSTと海氷密度を指定し、2.8°解像度モデルの場合は 6 時間間隔、1.4°解像度モデルの場合は 12 時間間隔で更新する。
🐾19 氷に覆われている海面の占める割合。単位は、%。0%の場合は海氷が全く存在せず、100%の場合は全体が海氷に覆われている状態を示す。
2⃣ 特徴量の正規化
学習ダイナミクスを改善するために、すべての入力特徴量を平均 0、分散 1 で近似的に分布するように正規化する。これは、初期化時にすべての特徴量の平均と標準偏差を推定し、ニューラルネットワークに特徴量を入力する前に適切なシフト変換と再スケール変換を追加することで行う。気圧ごとに正規化される比湿q を除き、すべての特徴量は全ての大気圧で均一に正規化された。
3⃣ ニューラルネットワークでの演算処理ーネットワーク･アーキテクチャ
NeuralGCMは、エンコード-プロセス-デコード・アーキテクチャを採用する(これは、既存モデル＝FourCastNet、Pangu-Weather、GraphCast及びAuroraと同じである)。
プロセス・コンポーネントには、5 層の全結合MLP(多層パーセプトロン)ブロックを使用する。すべての入力機能は、入力機能をサイズ 384 の潜在ベクトルにマッピングする線形層である｢エンコード｣層に渡される前に連結される。すべてのプロセス・ブロックは、384 個の隠れユニットを持つ 3 層 MLP ブロックを使用して、潜在ベクトルの更新を計算する。最後に、線形｢デコード｣層は、サイズ 384 のベクトルを出力ベクトルにマッピングし、その後、さまざまな変数の気圧ごとの値に分割する。
4⃣ 再スケーリング
学習した予報変数の変化率を予測する場合、ニューラルネットワークの出力は、スケーリングされる。具体的には、有限差分法を使用して1 時間間隔で、ERA5 データから推定される各変数の変化率の0.01標準偏差で、スケーリングされる。標準偏差は、地球全体で平均化された 10 個のスナップショットを使用して推定される。
(5) ハイパーパラメータ等(オプティマイザー周り)
0⃣ 学習率の最大値
水平解像度によって異なる。2.8°と1.4°の場合、2×10-3。0.7°の場合は、1×10-3。 アンサンブル･モデルの場合も、1×10-3。
1⃣ デコーダーの再学習以外
オプティマイザー・・・Adam：β1 = 0.9、β2 = 0.95、ε = 10-6 )
学習率のウォームアップ・・・最初の2,000ステップで最大値になるように直線的に増加。
学習率のスケジューリング・・・2,000ステップ～15,000ステップまで一定。その後、減衰時間を 10,000 ステップに設定して、減衰率0.5で指数関数的に減衰する。
2⃣ デコーダーの再学習
学習率のウォームアップ・・・最初の1,000ステップで最大値になるように直線的に増加。
学習率のスケジューリング・・・1,000ステップから2,000ステップまで一定。2,000ステップ以降、1,000ステップごとに、減衰率0.5で指数関数的に減衰する。
(6) エンコーダー、デコーダー
大気のσ座標表現を使用する NeuralGCMを ERA5データとインターフェースするために、エンコーダー･モジュールとデコーダー･モジュールを使用する。各モジュールは、学習された補正と組み合わせたσ座標と気圧座標間の再グリッドに基づいている。
1⃣ エンコーダー
エンコーダーを使用して、気圧座標のERA5データからσ座標のモデル状態を 3 つの手順で取得する。まず、ジオポテンシャルが地表気圧と一致する圧力を計算して、各経度緯度の地表気圧を計算する。次に、関連するすべての大気変数を NeuralGCM のσ座標に線形補間する。最後に、学習された物理モジュールと同じ構造のニューラルネットワークを使用して計算された補正を補間結果に追加する。最後のステップは、イニシャル･ショック🐾20を大幅に軽減する。これは、数値計算的解釈では数値振動、大気物理学的解釈では大気重力波を除去することに相当する。
線形補間と組み合わせる前に、ネットワーク出力は、対応する変数の 0.02 標準偏差でスケーリングされる。
🐾20 力学的バランスが崩れたデータ(物理的に不整合なデータ)を数値計算の入力値とすることで、数値振動が発生して起こる数値不安定性を、イニシャル･ショックと呼ぶ(らしい)。数値気象予報の文脈では、しばしばデータ同化において、観察される。致命的でないイニシャル･ショックは、モデルの散逸過程によって時間とともに散逸する[*39]。
2⃣ デコーダー
デコーダーを使用して、σ座標上のモデル状態を気圧座標上のERA5データに 3 つのステップでマッピングする。
㈠まず、温度と湿度からジオポテンシャルを診断する(☛”診断”という文言については、🐾18を参照)。㈡次に、結果を ERA5の気圧座標に補間する。地表気圧については、線形補間を使用する。σ座標から気圧座標にデータを外挿するには、ジオポテンシャルと温度に、線形外挿を使用する。他のすべての変数については、外挿に境界値を使用する。これらの処理は、ECMWF が圧力外挿に使用するより複雑な手順を近似することを目的としている。㈢最後に、学習済物理モジュールと同じ構造のニューラル ネットワークを使用して計算された補正を補間結果に追加する。
補間の結果と組み合わせる前に、ネットワーク出力は対応する変数の 0.02 標準偏差でスケーリングされる。
(7) 決定論的モデルの損失関数ー概要
0⃣ 概要
決定論的モデルは、損失関数の目的と最適化されるパラメータが異なる2 つの段階で学習される：1⃣プライマリ学習と2⃣デコーダーの再学習、である。1⃣と2⃣の損失関数は、どちらも、距離の定義が異なる、再スケール変数🐾21の平均二乗誤差(MSE)である。
🐾21 損失関数を計算する前に、大気変数の自然スケールの違い及び、リードタイムの関数としてのグラウンドトゥルースからの偏差の増加に対処するために、再スケーリングを行う。
1⃣ プライマリ学習
プライマリ学習では、NeuralGCM予測と ERA5 データ間のさまざまなタイプの不一致を考慮するために、｢精度(accuracy)、シャープネス、バイアス｣の 3 タイプの組み合わせを使用した。3タイプの損失関数は、さらに、距離の定義が異なる(☛下記(7)2⃣～4⃣参照)。バイアス損失関数は、ペナルティ項として機能する。
また、バイアス損失関数を除くすべての損失関数は、気圧座標で表示される｢データdata｣とσ座標で表示される｢モデルmodel｣の両方の表現に適用される。これにより、総損失関数は、㊀＋㊁＋㊂＋㊃＋㊄となる(☛下記(7)1⃣参照)。 なお、㊀～㊄は全て、λ×Mという形式で表される。λは係数で、Mは距離の二乗の平均値で表される。
2⃣ デコーダーの再学習
この学習フェーズでは、デコーダー･モジュールのパラメーターのみを更新する。デコーダーの再学習では、球面調和関数から格子点空間への変換の切り捨てエラーから生じる、スペクトル･アーティファクトを削除する。このアーティ･ファクトに対処するため、エンコーダー・モジュール及び学習済物理モジュールを、すべてフリーズする。そして、6,12,18,24 時間🐾22の予測に対して、グリッド空間表現で計算されたMSE(平均二乗誤差)を使用して、デコーダー・パラメーターを最適化する。
4,000 の勾配降下ステップを使用するが、通常、損失と評価指標は、最初の1,000ステップ後に横ばいになる。
🐾22 NeuralGCM-1.4°の場合は、12 時間および 24 時間。
(8) 決定論的モデルの損失関数ー詳細
0⃣ 距離について
距離は、予測値とグランドトルゥースとの間で計算する。3タイプの損失関数(下記2⃣～4⃣)は、距離の定義が異なる。なお、グランドトルゥースの時刻はt＋τとする。ここでτはリードタイム🐾23。予測値は時刻tで初期化された🐾24、時刻t＋τの計算値である。
🐾23 リードタイムτ(単位は、hour)は{6, 12, 18, 24, 36, 48, 60, 72, 96, 120}から選ばれる。
🐾24 時刻tにおけるグランドトルゥースを初期値として与えられた、という意味。
1⃣ 総損失関数の概要
㊀＝λdata×MData、㊁＝λspec×MDataSpec、㊂＝λmodel×MModel、㊃＝λspec×MModelSpec、㊄＝λbias×MMSBiasである。
係数は、λdata = 20、λspec = 0.1、λmodel = 1、λbias = 2は経験的に選択された。
㊀と㊂が｢2⃣精度(accuracy)｣に該当する。ただし、㊀がデータ(気圧座標表示)で、㊂がモデル(σ座標表示)である。 ㊁と㊃が｢3⃣シャープネス｣に該当する。ただし、㊁がデータ(気圧座標表示)で、㊃がモデル(σ座標表示)である。 ㊄が｢4⃣バイアス｣に該当する。｢3⃣シャープネス｣および｢4⃣バイアス｣に対応する損失関数は、｢総損失関数と比較して小さいものの、NeuralGCM予測の鮮明さにプラスの効果があることがわかった｣とのコメントが付されている。
2⃣ 精度(accuracy)損失関数
距離2＝∑(予測値ーグランドトルゥース)2。∑は、全波数と東西波数に対してとる。ただし予測値とグランドトルゥースには、フィルタリングを適用する。フィルタリングの目的は、カオス ダイナミクスによって生じる｢予測可能性の範囲を超える｣成分を削除することである。具体的には、まず、精度損失関数に含めるべき全波数の最大値(しきい値)を、各リードタイムに対して推定する。次に、しきい値よりも高い総波数を除去する。
3⃣ シャープネス損失関数
距離2＝∑(S(予測値)ーS(グランドトルゥース))2。ここでスペクトル･エネルギーS(＊)は、＊に対して東西波数で総和をとった量である。∑は、全波数に対してとるが、上限値までとする。この上限値は、水平解像度2.8°、1.4°、0.7°に対して、 42、80、120が適用された。
4⃣ バイアス損失関数
距離の定義は、2⃣と同じであるが、距離計算する前に、バッチとリードタイムの次元で平均化する点が異なる(またフィルタリングは行わない)。バイアス損失関数は、ペナルティ項として機能する。
(9) 確率論的モデルの損失関数
アンサンブルモデル(NeuralGCM-ENS)の損失関数である。連続ランク確率スコア(CRPS)は、確率的予測が関与する場合に、最も広く使用される指標の1つであり、アンサンブル･モデルに対して、より優れた評価指標であった。NueralGCM-ENSの損失関数は、多変量CRPSであり、2種類のCRPSの和(足し算)で構成される。一つは、精度を高めるために用いられる。もう一つは、アンサンブルの広がりを高めるために用いられる。
CRPSは、6時間から5日間のロールアウト長(リードタイム？)で計算される。また、各変数について、最大カットオフ波数未満のグリッド空間のCRPSと球面調和基底のCRPSの合計を使用する。
【4】中期予報に関するNueralGCMの評価
(0) 前説
0⃣ベースラインと評価指標
ECMWFのアンサンブル･モデル(ECMWF-ENS)をベースラインとして使用する。評価指標には、二乗平均平方根誤差(RMSE)、二乗平均平方根バイアス(RMSB)、連続ランク確率スコア(CRPS)及び、スプレッド･スキルの関係(spread-skill ratio)🐾25、を使用する。
🐾25 アンサンブル･スプレッドと、アンサンブル平均予報値の予報誤差(RSME)の比。アンサンブル予報が上手くいくためには、スプレッド･スキルの関係が1になることが望ましい。[*18]を参照。
1⃣ 総評
㈠ 決定論的モデル×短リードタイム
与えられた初期条件に対して単一の天気予報を生成する決定論的モデルは、短いリードタイムで RMSE スキルを使用して効果的に比較できる。最初の 1 ～ 3 日間は、大気変数に応じて、天気パターンの時間発展を正確に追跡する予報によって RMSE が最小化される。この時間スケールでは、NeuralGCM-0.7° と GraphCast(どちらもmade by グーグル！)が最良の結果を達成し、さまざまな変数間でわずかなばらつきがあることがわかった。
㈡ 決定論的モデル×長リードタイム
リードタイムが長くなると、カオス的発散によりRMSE が急速に増加するため、決定論的モデルにとってRMSEの有用性が低下する。RMSB は、時間の経過に伴う永続的なエラーを計算し、モデルがはるかに長いリードタイムでどのように機能するかを示す。ここでも、NeuralGCM モデルは以前のアプローチと比較して有利であり、熱帯地方の特定の湿度に対するバイアスが大幅に少なくなっている。
㈢ アンサンブル･モデル
アンサンブル･モデルは、特にリードタイムが長い場合に、天気予報固有の不確実性を捉えるために不可欠である。約7 日間を超えると、ECMWF-ENS と NeuralGCM-ENS のアンサンブル平均は、決定論的モデルよりも RMSE がかなり低くなっており、これらのモデルが起こり得る天候の平均をよりよく捉えていることを示している。
NeuralGCM-ENS は、アンサンブル平均 RMSE、RSMB、CRPS について、ほぼすべての変数、リード タイム、鉛直レベルで ECMWF-ENS と比較して誤差が低く、スキルの空間パターンも同様である。ECMWF-ENS と同様に、NeuralGCM-ENS の｢スプレッド･スキルの関係｣は約 1 である。
(1) 熱帯低気圧、大気河川、熱帯収束帯
1⃣ 概要
本論文では、熱帯低気圧、大気河川、熱帯収束帯🐾26に対するNeuralGCM のパフォーマンスが示されている。NeuralCGMの比較対象となっている機械学習ベースの予報(MLWP)モデルは、Pangu-WeatherとGraphCastである。この二つは、ERA5及びECMWF-HRES予測よりも、大幅に｢ぼやけた予測｣を行っている。対して、NeuralCGMは解像度が粗い🐾27にもかかわらず、優れている。
NeuralGCM と ECMWF の両方からのアンサンブル平均予測(つまり、NueralGCM-ENSとECMWF-ENSの予測)は、平均的な意味で ERA5 に近いため、長いリード タイムでも本質的に滑らかである。
🐾26 南北両半球からの貿易風が合流する帯状の境界。赤道前線は、同じ意味。出所：https://www.jma.go.jp/jma/kishou/know/yougo_hp/haichi3.html
🐾27 GraphCast および Panguの水平解像度は 0.25°、NeuralCGMは最も細かい0.7°を使用。ちなみに、ERA5の水平解像度は0.25°で、ECMWF-HRESの水平解像度は0.1°。
2⃣ ぼやけた予測について
ぼやけた予測は、物理的に一貫性のない大気条件に対応し、極端な気象を誤って表現する。さまざまな予測モデルのぼやけ具合は、パワースペクトルで定量化できる。NeuralCGM-0.7° のパワースペクトルは、他のMLWPモデル(Pangu-WeatherとGraphCast)よりも一貫して ERA5 に近いことを示しているが、ECMWF-HRESと比べると、ぼやけている。
NeuralGCM のスペクトルは解像度が上がるにつれて精度が上がる。これは、より高い解像度で学習された NeuralGCM モデルを、さらに改善できる可能性を示唆している。
(2) 水収支
NeuralGCM では、移流は力学コアによって処理され、機械学習のパラメータ化は大気の鉛直柱内の局所的過程をモデル化する。したがって、純粋な機械学習手法とは異なり、水平輸送やその他の解決されたダイナミクスにより、局所的なソース(湧き出し：流入)とシンク(吸い込み：流出)を分離できる。これにより、結果の解釈が容易になり、水収支の診断が容易になる。
具体的には、機械学習ベースのアプローチのように直接予測するのではなく、降水量から蒸発量を引いたものを診断する(☛”診断”という文言については、🐾18を参照)。短期天気予報では、降水量から蒸発量を引いた平均は、ERA5 データに非常に近い現実的な空間分布を示す。NeuralGCM-0.7°の降水量から蒸発量を引いた分布は、熱帯地方の極端な事象を過小評価しているが、温帯では ERA5 の分布とほぼ一致している。NeuralGCM の現バージョンでは、大気柱における大気変数の変化率を直接予測するため、降水量と蒸発量を区別できない。
(3) 地衡風バランス
本論文では、NeuralGCM、GraphCast、ECMWF-HRES が地衡風バランス、つまり中緯度における大規模なダイナミクスを駆動する支配的な力の間のほぼ平衡をどの程度捉えているか、が扱われている。ECMWFの数値予報モデルの研究者が、MLWPモデルの評価を行った論文[*44]では、Pangu が地衡風と非地衡風の鉛直構造を誤って表現していることが強調され、リードタイムが長くなると劣化が見られることが指摘されている。同様に、GraphCast の誤差はリードタイムとともに悪化することが観察されている。
対照的に、NeuralGCM は、ERA5データと比較した場合、さまざまなロールアウト🐾28で GraphCast と比較して、地衡風と非地衡風の鉛直構造とその比率をより正確に描写している。ただし、ECMWF-HRES は、NeuralGCM よりも ERA5 データにわずかに近くなっている。 NeuralGCM では、地衡風の垂直構造の表現は最初の数日間はわずかに劣化するだけで、その後は特に 5 日目以降は目立った変化は見られない。
🐾28 ここではリードタイムと同じ意味(だと思われる)。具体的には、1日、5日、10日が選択されている。
【5】気候シミュレーションに関する評価結果
(0) 前説
NeuralGCM モデルは、最大3日先の天気を予測するように学習されているが、一般的には中期の時間スケールをはるかに超えて、大気をシミュレートできる。NeuralGCM を使用した気候シミュレーションでは、2.8° および 1.4° の決定論的モデルを使用する。これらのモデルは学習に比較的費用がかからず、より広いパラメーター空間を探索して安定したモデルを見つけることができる。
以前の研究では、数値的不安定性と気候ドリフトのため、ハイブリッド モデルを使用した(時間スケールを伸ばす)拡張シミュレーションの実行は困難であることがわかった。本論文では、選択したモデルの安定性を定量化するために、複数の初期条件を実行し、そのうちのいくつが不安定にならずに終了したか、が報告されている。
(1) 季節サイクル
0⃣ 概略説明
NeuralGCM が季節サイクルのさまざまな側面をシミュレートする能力を評価した。モデルは、NeuralGCM-1.4°を採用した。具体的には、10 日間隔の37 種類の初期条件で、2019年から2 年間のシミュレーションを実行した。これらの 37 種類の初期条件のうち、35 種類は不安定性なく 2 年間を正常に完了した。
1⃣ 世界平均気温
NeuralGCMの35のシミュレーションで捕捉された2020 年までの世界平均気温の時間的変動を、ERA5 再解析およびstandard climatology benchmark🐾29と比較した。NeuralGCM の世界平均気温の季節性と変動性は、ERA5 で観測されたものと定量的に類似している。 NeuralGCMのアンサンブル平均気温RMSEは、ERA5をベンチマークした場合0.16℃となり、standard climatology benchmarkのRMSE 0.45℃を大幅に上回っている。
🐾29 よくわからない。
2⃣ 降水量の年間バイアス
NeuralGCMとX-SHiELDを比較した。X-SHiELDは、個々の対流セルを明示的にシミュレートするために、十分に高い全球解像度で実行される、実験的全球雲解像モデルである[*45]。X-SHiELDなどの全球雲解像モデルは、その解像度が深層対流を解像できるため、特に水循環のシミュレーションでは最先端のモデルと考えられている。比較期間は、X-SHiELDデータが利用可能な2020年1月19日から2021年1月17日である。
NeuralGCMの降水量の年間バイアス(RMSE 1.09 mm)は、X-SHiELD(RMSE 1.74 mm)と標準気候学ベンチマーク(RMSE 1.36 mm)よりも小さい。さらに、NeuralGCM は X-SHiELD よりも上部および下部対流圏の温度バイアスが低い。また、X-SHiELD の降水量バイアスを、 NeuralGCM-1.4° の｢降水量－蒸発量｣バイアス(為念：ーはマイナスの演算記号。つまり、蒸気量を控除した降水量バイアス)と間接的に比較すると、NeuralGCM の方がわずかに大きいバイアスとグリッド スケールのアーティファクトが示されている。
3⃣ 熱帯低気圧
NeuralGCMが熱帯低気圧を再現する能力を評価するために、熱帯低気圧トラッカー TempestExtremes🐾30を使用した。NeuralGCMは、1.4° の粗い解像度でも熱帯低気圧の現実的な軌跡と数を再現する(対応する期間の ERA5＝86個に対して、83個)。これに対して、X-SHiELDを1.4°の水平解像度に再グリッド化すると熱帯低気圧の数が大幅に過小評価される(40個)。
🐾30 TempestExtremes(TE)は、長方形格子または非構造格子/ネイティブ格子上の地域的または全球的な、地球システムデータセットの特徴検出、追跡、科学的分析のための多面的なフレームワーク。TEフレームワーク･バージョン2.1では、熱帯･亜熱帯低気圧、モンスーン性低気圧･低気圧、大気河川、大気ブロック、降水クラスター、熱波など、を調べるための広範なサポートを提供する(らしい)。出典：https://gmd.copernicus.org/articles/14/5023/2021/
5⃣ その他
NeuralGCMは、全球可降水量の年間サイクルや全球総運動エネルギーなどの指標からもわかるように、季節サイクルを正確にシミュレートしている。さらに、NeuralGCMは、ハドレー循環🐾31や帯状平均🐾32帯状風🐾33などの、重要な大気ダイナミクス、さまざまな季節における渦運動エネルギーの空間パターン、モンスーン循環の独特の季節的挙動を捉えている。
🐾31 大気(大)循環の根源は、赤道付近の温かい空気と極付近の冷たい空気との間で起こる熱対流である。ただし、地球自転の影響を受け、緯度に応じて3つの循環が生まれる。3つとは、赤道付近の低緯度の｢ハドレー循環｣、極付近の｢極循環｣及び、その中間緯度の｢フェレル循環｣である。ハドレー循環は、地表面で｢貿易風｣を発生させる。ハドレーという名称は、英国の弁護士＆気象学者のジョージ･ハドレーに由来する。
🐾32 気象学の文脈で、帯状平均とは、｢東西方向にわたる平均｣を意味する。平均をとる範囲は地球1周分であり、緯度を特定して平均をとる。
🐾33 対流圏上部または圏界面(対流圏と成層圏の境界面)付近の狭い領域に集中して吹いている、帯状の"非常に強い"風。気象学の文脈で、帯状とは、東西方向という意味である。"非常に強い"を定量的に表現すると、秒速100m(出典：https://www.nagare.or.jp/download/noauth.html?d=30-5tokushu3.pdf&dir=135)。ちなみに観測史上最大の台風(2013年台風30号)の最大瞬間風速が、105m。
(2) 10年スパンのシミュレーション
1⃣ 気温変動
NeuralGCM が過去の気温傾向をシミュレートする能力を評価するために、NeuralGCM-2.8° を使用して 40 年間の AMIP🐾4(再)のようなシミュレーションを実施した。1980 年に 10 日間隔で初期条件を設定して実行した 37 回の異なる実行のうち、22 回のシミュレーションは 40 年間を通じて安定していた。その結果を、CMIP6🐾34の22 回のシミュレーションと比較する。
⓵22回のNeuralGCMシミュレーション全てと、⓶22 回のCMIPシミュレーションの平均は、ERA5 データで観測された地球温暖化の傾向を正確に捉えている。年毎の気温の傾向とERA5 データには強い相関関係があり、NeuralGCM が SST(海面温度)強制の気候への影響を、効果的に捉えていることを示唆している。1981～2014 年の平均空間バイアスを比較すると、⓵が ⓶よりも、バイアスが小さいことがわかった。
🐾34 CMIPとは、世界気候研究計画(WCRP)によって開始された、相互比較プロジェクトプロジェクトである。CMIP6(Coupled Model Inter-comparison Project phase 6:第6期結合モデル相互比較プロジェクト)は、気候モデル相互比較プロジェクトであり、｢陸、海、大気、エアロゾル変数｣を含むさまざまな気候モデリング実験を組み合わせている。
2⃣ 熱帯温暖化傾向の鉛直構造
気候モデルが上層圏で過大評価する傾向がある、熱帯温暖化傾向の鉛直構造を調査した。線形回帰で計算された NeuralGCMの傾向は、CMIP6よりも ERA5 に近い。特に、上層圏のバイアスは減少している。ただし、NeuralGCM は、規定の SSTによって温度が通常より制約される地表近くでも、CMIP6よりも予測のばらつきが大きくなっている。
3⃣ ゼロショット予測
SST(海面温度)の上昇を伴うAMIPシミュレーションを実施することで、NeuralGCM がこれまでに見たことのない温暖な気候に汎化できるかどうかを評価した。NeuralGCM は、中程度の SST上昇(+1℃および+2℃) に対しては、地球温暖化の特徴をいくつか示している。ただし、より大幅な上昇 (+4℃)に対しては、地球温暖化の特徴を示せていない。さらに、SST の上昇を伴うAMIP シミュレーションでは気候ドリフトが見られ、このコンテキストでの NeuralGCM の限界が強調された。
【6】まとめ・考察
(0) MLWPモデルにとっては"敵方総大将"にあたる、欧州中期予報センター(ECMWF)の研究者が、著者として名を連ねているところに信頼感がある(という印象を与える)。
余談1：グーグルは、MLWPモデルにおいて(高精度の)、超短期･中期･超長期のモデルを開発したことになる・・・スゴい。
余談2：力学過程を担う数値計算モジュールに、ニューラルネットワークの出力を付加するという枠組みは、マテリアルズ･インフォマティクスにおける｢機械学習ポテンシャルを使った第一原理計算｣をイメージさせる。
(1) 改めてまとめると、パラメタリゼーションという名を持つ無理ゲーの回避策は、物理過程をニューラルネットワークで表現したことである。そして重要なことは、数値計算ベースで解析する力学過程とニューラルネットワークを上手く連成させた、ことである。上手く連成させないと、数値振動や気候ドリフトが発生し、安定したシミュレーションは実行できない。NeuralGCMは(決定論的モデル・確率論的モデルともに)、損失関数を工夫することで、気候ドリフトを回避している(と理解している)。その為、損失関数は、やや複雑になっている。数値振動(イニシャル･ショック)は、エンコーダーを工夫することで、回避している。ただし、初期条件によっては、安定したシミュレーションを実現できない。
ちなみに本論文では、｢雲のフィードバックなど、天気の時間スケールに微妙な影響を与える気候の重要なプロセスを学習するには、おそらく別の学習戦略が必要になるでしょう｣との記述がある。
(2) 複数の主張がなされているが、最も重要な主張は、NeuralGCMが、｢ECMWF-ENSよりも優れたCRPSを備えた、正確なアンサンブル予報を行う最初の機械学習ベースのモデル｣というものである。つまり、アンサンブル･モデル(確率論的モデル)が重要である🐾35。アンサンブル･モデルと決定論的モデルとの違いは、本質的には、損失関数の違いのみである(と理解している)。アンサンブル･モデルの損失関数は、アンサンブル･モデルに対して｢より優れた評価指標｣とされる CRPSで構成したところがポイントであろう。しかも、2つのCRPSを使っている。これは、高い精度を求めつつ、アンサンブルの広がりを高めるためである。
🐾35 例えば、中期予報は、オマケみたいなものであろう。中期予報について、Pangu-WeatherやGraphCastよりも高精度と主張しているわけであるが、Auroraには勝てないだろう(アプローチが異なるので、ある意味、当然)。
(3) 具体的な計算リソースの削減について、本論文では、次の記述がある： NeuralGCM-1.4°では、1つのTPU(テンソル処理ユニット)を用いて、24時間で70,000日のシミュレーションを行った。これに対して、X-SHiELDでは13,824のCPUコアで19日のシミュレーションを行った。👉これは、公平な比較になっているのだろうか？
複雑系における将来予測🐾1である気象予報は、本来、確率でしか与えることはできない。確率で予報を示すことは、当たらなかった時の”保険”をかけているわけではない。そうすることが、正解である。気象予報のコミュニティでは、確率で示した予報を｢アンサンブル予報｣と呼ぶ。
グーグル･ディープマインドの研究者は、｢アンサンブル予報で、NWP🐾2を上回るMLWP🐾3を開発した｣と主張する論文[*50](以下、本論文)を発表した(24年12月4日@nature)。これまで、決定論的な予報(別の表現を使うと、単一メンバーによる予報)でNWPを超えた、と主張するMLWPはいくつも存在する(本頁他項で掲示した通り)。
アンサンブル予報については、｢ForeCasNet(NVIDIA)、Pangu-Weather(中国ファーウェイ)、MetNet3(グーグル)及びNeuralGCM(グーグル)｣がチャレンジしているが、｢GraphCast(グーグル)及びAurora(マイクロソフト)｣は、チャレンジしていない。そして、MetNet3及びNeuralGCMは、部分的にであるが、既にECMWF-ENSを超えたと主張している。部分的の意味は、MetNet3は超短期🐾4であり、NeuralGCMはECMWF-ENSより解像度が粗い🐾5。今般グーグルが開発したMLWP『GenCast』は、ガチンコ勝負で、ECMWF-ENSを超えたと主張する。もちろん、ガチンコ勝負とは、土俵が同じという意味である。
🐾1 企業の事業価値評価算出も、複雑系における将来予測に該当する。従って、”幅”でしか、事業価値を与えることはできない。しかし、これは、価値を価格に置き換えても成立する。つまり、(新聞記者は勘違いしているようであるが、)価格も幅でしか、与えることはできない。
🐾2 為念：数値計算ベースの気象予測モデル
🐾3 為念：機械学習(深層学習)ベースの気象予測モデル
🐾4 MetNet3は、24時間先までの予報しかできない。
🐾5 決定論的NeuralGCMの水平解像度は、2.8°、1.4°、0.7°であるが、確率論的NeuralGCM(つまり、アンサンブルモデル)の解像度は1.4°である。対するECMWF-ENSの解像度は、0.25°。
❚為参考❚
スタンフォード大学人間中心AI研究所(HAI)は、2017年から毎年｢AI Index Report｣という調査報告書を公開している(2020年はコロナ禍のため、例外)。最新の2025年版第5章科学と医療[*54]に、｢注目すべきモデル｣という囲み記事がある。9個のモデルが紹介されているが、その内5個が、気象関連である。5個の中にはGenCastも含まれている。
AI Index Report2025では、次のように簡潔に紹介されている(ドメインは、weather prediction)：GenCastはAIを搭載した気象モデルで、拡散ベースのアプローチを用いて高精度の15日予報を提供し、ほぼすべての指標でENSのような従来のシステムを凌駕している。数時間単位ではなく数分単位で予測を生成し、災害対応、再生可能エネルギー、農業など幅広い分野で応用されている。
【1】本論文の主張
本論文は、以下を主張する：
｢0.25°の水平解像度で、15 日間のアンサンブル予報を生成するGenCastは、世界最高の運用中期天気予報である欧州中期予報センター(ECMWF)のアンサンブル予報、ECMWF-ENSを超えた｣。少し詳しくに述べると、次の通り(主語は、GenCast)。
(1) アンサンブル･スキル(☞詳細は【4】(1)参照)
｢1,320個あるターゲット🐾6の97.2% に対してECMWF-ENS よりも優れたCRPS🐾7を示した｣
｢2つの地表変数･3つの上空変数に対して、スプレッド･スキルの関係🐾7は、同等｣
｢2つの地表変数に対して、ECMWF-ENS よりも優れたREV🐾7を示した｣。
(2) 熱帯低気圧の進路予測
｢アンサンブル平均･熱帯低気圧軌道の位置誤差で、ECMWF-ENS よりも優れている」
｢進路確率🐾8を使ったREVで、ECMWF-ENS よりも優れている｣
(3) その他(☞【4】(3)を参照)
｢気象状態に固有の空間依存性構造(の一部)を、ECMWF-ENS よりも上手く捉えている｣
🐾6 ターゲットという文言は、気象予測･天気予報の文脈では、"評価対象変数"の意味で用いられる。
🐾7 【2】(3)を参照
🐾8 【4】(2)1⃣を参照
【2】事前整理
(1) 拡散確率生成モデル
0⃣ 時間のインデックス
広く知られているように、拡散確率生成モデルには、いわゆる順過程と逆過程が存在する。そのため、時間ステップのインデックスを導入すると説明が容易になる。順過程における時間tのインデックスは、0→T、逆過程はT→0とする。
1⃣ 推定の枠組み
拡散確率生成モデルは、一般的には、拡散モデルと呼ばれている。｢推定の枠組み｣という観点で見ると、拡散モデルは、ベイズ推定を使っていると考えて良いだろう。つまり、事前確率分布から事後確率分布を逐次的に推定していく枠組み、と解釈して良いだろう。
2⃣ 事後確率分布からのサンプリング
あくまで、GenCastの最終アウトプットは、確率分布である。確率分布を求めるには、確率変数を扱う必要がある。確率変数を作るために、拡散モデルでは、決定論的な観測データ🐾9にガウス雑音を加え続け、最終的には、ほぼ正規分布に従うデータ＝ノイズ(ほぼ白色雑音)にしてしまう。この正規分布に従うデータを、ここでは便宜上、最終母集団(t=T)と呼ぼう。なぜ、正規分布するかというと、正規分布が扱いやすい確率分布であるからに他ならない。この正規分布が、”最初の”事前確率分布(t=T)となる。
これから推定するという意味で"未知の"事後確率分布からサンプリングする方法には、MCMC(マルコフ連鎖モンテカルロ法)がある。性能が低いプレーンなMCMCを使うかは別として、”未知の”事後確率分布からサンプリングした値に、最終母集団からノイズを除去した母集団(t=T－1)との差が最小になる条件を課すことで、事後確率分布(t=T－1)が推定できる。ただ、本論文は、その手法は用いていない(☞(2)2⃣参照)。
🐾9 本稿では、アプリケーションをあからさまに、気象予報に全振りしているので、”観測”データとしている。もちろん、一般的には、”観測”データに限定されない。
3⃣ 条件付き確率の積
事後確率分布(t=T－1)は、最終母集団(t=T)を所与とした場合に、母集団(t=T－1)が従う条件付き確率p(xT－1|xT)である。xT－1 ∊ 母集団(t=T－1)、xT ∊ 母集団(t=T)である。
上記の諸手続きを繰り返すことで、元の決定論的な観測データが、｢実は、確率変数であった｣とみなした場合の、事後確率分布(言うなれば、始祖の事後確率分布)を求めることができる。もちろん、これは、順過程･逆過程がマルコフ過程であると仮定した場合に成立する。白色雑音を順次加えていく順過程は何の問題もないが、逆過程はマルコフ過程になるように、ノイズを除去していくことが求められる。もっとも、逆過程を(拡散確率)微分方程式で記述する枠組みを採用していれば、そもそも懸念は生じない(はず)。
4⃣ ノイズの除去(デノイジング･ネットワーク)ー頭出しのみ
拡散モデルは、ノイズを除去するノイズ除去器が、ニューラルネットワークを使って実装できるため実現(実装)可能となった。拡散モデルを使ったMLWPはGenCastだけではない(が、ほぼない)。例えば、米フロリダ国際大学他[*51]によるCoDiCastは、クロス注意機構ブロックとU-Net🐾10で、 デノイジング･ネットワークを構築している(ので、新味はない)。GenCastが具体的に、どのように実装しているか(どのようなアーキテクチャを採用しているか)については、【3】(3)で述べる。
🐾10 U-Net は元々、バイオメディカル画像セグメンテーション用畳み込みニューラルネットワークとして提案されたモデルである。モデル･アーキテクチャは、オート･エンコーダー(エンコーダー･デコーダー)。気象予測においてU-Netを使う価値は、U-Netは｢気象場の地域特性を学習できる｣[*33]ところにある。気象場の地域特性とは、例えば、｢風が吹いた後の気温の下がり方｣や｢地面から空気への熱の伝わり方｣が地上変数・上空変数に与える影響を指している(と思われる)。
5⃣ 冷静に考えてみれば・・・
冷静に考えてみれば、未知関数の勾配を求めるタスクで、｢行って・来い｣自体は、不思議でも何でもない。機械学習･深層学習の文脈であれば、誤差逆伝播法が該当する。気象予報の文脈では、高精度な初期値の作成のために使われる、4次元変分法におけるアジョイント法が、該当する。t＝Tから逆向きに解くことによって、(拡散モデルでは、スコア関数と呼ばれる)未知関数の勾配を求めることができる。｢行って・来い｣は、数学的観点から、むしろ標準的と言って良いのかもしれない。拡散モデルでは、順過程・逆過程ともに、(各過程を離散化されたものと見做して、連続化することによって)確率微分方程式で定式化することができる。この確率微分方程式は、驚くことに常微分方程式に変換できる。
拡散モデルの特徴として、元の大きい問題(観測データが従う確率分布の推定)を細かい工程に分解して実行することで、安定的な学習が可能、が上げられる。ただし、それらの代償として、推定に時間がかかる。
(2) 未知の事後確率分布を求める
1⃣ スコア関数
拡散モデルにおける"未知の"事後確率分布は、スコア関数と呼ばれる。正確には、スコア関数＝∇log(未知の事後確率分布)である。未知の事後確率分布を使う代わりに、スコア関数を使う理由は、正規化定数を扱う必要がないからである。ここで言う正規化定数とは、統計力学で言えば分配関数である。分配関数を求めるには、煩雑な🐾11積分計算が必要であり、計算コストが高いので忌避すべき作業とされる。
なお、スコア関数による推定と、変分法による推定は一致することが(数学的には、ギルサノフの定理によって)証明されている[*52]。わざわざ標語的に述べれば、拡散モデルの枠組みでは、スコア関数を正確に推定することがキモである。
2⃣ 微分方程式を解く・・・
本論文では、前述の通り、MCMCあるいはランジュバン･モンテカルロ法🐾12といったsampling手法を使って"未知の"事後確率関数を推定するというアプローチを採用しない。ここでは、下記との混同を避けるため、samplingという表記にした。sampling手法は遅いのが、不採用の理由であろう。
本論文の枠組みでは、拡散モデルにおける順過程･逆過程を微分方程式で記述している。GenCastは、逆過程の微分方程式を解くことで、新しいサンプル(地表変数及び上空変数の予測値)を生成する。つまり、微分方程式の答えが、地表変数及び上空変数の予測値である。言い方を変えれば、｢微分方程式を解く行為が、未知の確率分布から新しいサンプルをサンプリングすること｣に等しい。なお、微分方程式を解いてサンプリングするとアイデアは本論文のオリジナルというわけではない。
xが従う事後確率p(x)を推定、p(x)からサンプリングーーーーーーーー>{x1,x2,･･･,}
xに関する(p(x)を"パラメータ"として含む)微分方程式を解いた結果 ー>{x1,x2,･･･,}
例えば、拡散方程式の拡散係数を変化させれば、拡散方程式の解(解の軌跡)は変わる。常微分方程式のパラメータに、ランダムネスの種を仕込んでおけば、解は変化する。ランダムネスの種が、適当な確率分布に従っていれば、微分方程式の解は、その確率分布の影響を受ける。要するに、微分方程式の解が適当な確率分布からのサンプリングに相当するように、微分方程式を仕込めば良いということになる。拡散モデルの場合、それがスコア関数ということになる。正確には、スコア関数を少し変形させる(というか掛け算を行う)。スコア関数はノイズ除去器として表現され、ノイズ除去器はネットワークに学習器として実装される(実装できる)。つまり、ノイズ除去器＝スコア関数は、データによって学習される。未知の事後確率分布は、あくまで陰的(implicit)に求められている。
解く微分方程式は、拡散確率微分方程式(拡散SDE)もしくは拡散常微分方程式(拡散ODE)である。SDE と比較すると、ODE はランダム性がないため、より大きなステップ サイズで解くことができるため、ODEのほうが計算量が少ない(ODEは2次加速する!)[*52]ので、後者が選択される。具体的に本論文では、2次精度のODEソルバー｢DPM-Solver++2S ソルバー｣を使って拡散ODEを解いている(＋様々な工夫を凝らしている)。
🐾11 英語ではintractableという文言が用いられる。｢手に負えない｣という日本語があてられる。
🐾12 スコア関数を使うことで、高性能なサンプリングが行うことができる。
(3) 評価指標
1⃣ CRPS(Continuous Ranked Probability Score：連続ランク確率スコア)
正しく、数式で書くとやや煩雑であるが、本質的には、次の㊀と㊁の差が、CRPSである：㊀アンサンブル予測値とグランドトルゥースとの差(平均値)、㊁異なるアンサンブルメンバーにおける予測値の差(平均値)。平均は、アンサンブル･メンバー全てに渡ってとる。
意味合いとしては、｢予測値が確率的である場合に、平均絶対誤差(MAE)を拡張した指標｣である(から、そういう理解で良いだろう)。アンサンブル予測(確率的な予測を)することで生じる｢ばらつき｣を控除することで、アンサンブル予測のデメリットを補償し、決定論的な予測における評価指標と環境を近づけている、と解釈できるだろう。
2⃣ スプレッド･スキルの関係(spread/skill ratio)
アンサンブル予報が上手くいくためには、スプレッド･スキルの関係が1になることが望ましい。スプレッド･スキルの関係は、次の式で定義される：
スプレッド･スキルの関係＝アンサンブル･スプレッド/アンサンブルRMSE
ここで、 アンサンブル･スプレッド＝√アンサンブル分散、アンサンブル分散＝アンサンブル予測値とアンサンブル平均値とで計算した分散である。アンサンブルRMSE(二乗平均平方根誤差)は、アンサンブル平均値とグランドトルゥースとでとったRMSEである。
3⃣ REV(Relative Economic Value：正式な訳語があるかは分からないが、相対的経済価値)
REVは、｢概念上、損得勘定を加味した｣気象予報に対する評価指標である。概念上という意味は、REVが、実際に金銭的価値を表現しているわけではない、という意味である。REVは、次のように表現できる。
REV＝(min(コスト、真の期待損失)ー(予測期待コスト＋実際の損失))/(min(コスト、真の期待損失)ー真の期待コスト)
ここで、損失とは、気象に起因する災害による損失である。コストとは、損失を避けるために要する費用といった意味である。REVは、本質的には、｢(予測期待コスト＋実際の損失)と、真の期待コスト｣を比較している。真の期待値と予測期待値の差が、実際の損失とバランスしているかを測っていると考えても良い。min(コスト、真の期待損失)は、基準として設定してある(との解釈で良いだろう)。
『真の期待値、予測期待値』という言葉の意味は、二値分類における混同行列に現れる、{TP、FN、FP、TN}を使って表現した方が分かり易い。
真の期待値＝(TP＋FN)×(裸の値)
予測期待値＝(TP＋FP)×(裸の値)
実際の損失とは、FN×(裸の損失)である。重ねての注意喚起になるが、REVは金銭的価値を表現しているわけではない。つまり、名は体を表していない。あくまで、比率である。そして、比較のために使う場合、損失とかコストは適当な値で構わない(もちろん、コストが損失を上回ることはない。保険の掛け金が、損失額を上回らないのと同じ)。
【3】GenCastの説明
(0) 概要的なこと
0⃣ アンサンブル予測(予報)について
アンサンブル予測において、複数の予測の集合をアンサンブル、個々の予測をメンバーと呼ぶ。またアンサンブル予測では、元のデータ(具体的に名をあげれば、ERA5の生データということになる)に摂動を加えたデータを使って行われる。ちなみに、摂動を加えているメンバーを｢摂動ラン(あるいは摂動予測)｣、摂動を加えていないメンバーを｢コントロールラン(あるいは制御予測)｣と呼ぶ。
冒頭【0】でも上げたように、ForeCasNet、Pangu-Weather及びNeuralGCMはアンサンブル予測を行っているが、ForeCasNetとPangu-Weatherは摂動ランを使っている。ECMWF-ENSも同様である。ちなみに、ForeCasNetはERA5に白色雑音を加えで摂動データを作っている。Pangu-Weatherは、パーリン･ノイズという雑音を加えている。対して、NeuralGCMはCRPSを損失関数とすることで、アンサンブル予測を可能にしている。
1⃣ 改めてGenCastとは
GenCastは、拡散確率生成モデルを基にしたMLWPである。生成モデルは、確率分布を生成する学習モデルである。GenCastは、大気状態(＝地表変数及び上空変数)の確率分布をモデル化し、新しい｢地表変数及び上空変数｣を生成できる(条件付き)拡散モデルとして実装されている。生成するモデルは、変分自己符号化器(変分オートエンコーダ：VAE)や敵対的生成ネットワーク(GAN)を始め、数々ある。2015年に提案された拡散モデルは、当初、GANに太刀打ちできなかった。しかし、(極めて雑に表現すると)種々の革新の結果、GANを超えるまでになった。ということで、GenCastは拡散モデルを採用したのだろう。
拡散モデルは、反復的な改良のプロセスを通じて機能する。GenCastにおける将来の大気状態Xt+1は、純粋なノイズとして初期化された残差Zt+10を、時間ステップにおける2つ前の大気状態(Xt、Xt－1)に基づいて反復的に改良することで生成される(反復プロセスについては、下記(1)を参照)。
なお、ECMWF-ENSとの比較を公正に行うためであろうが、摂動ランを使ったアンサンブル･モデルも別途用意されている。それが、GenCast-Perturbedである。ただし、本稿ではGenCast-Perturbedの結果は割愛する。
2⃣⃣ 学習スケジュール、学習環境、推論時間など
GenCastの出力は、2 段階(ステージ1、ステージ2)のプロセスで学習されたモデルによって生成された。ステージ 1 は事前学習ステージで、200 万学習ステップを要した。このステージでは、グラウンドトゥルース･データセットが、水平解像度が0.25° から 1° にダウンサンプリングされた。この学習ステージは、32 個の TPUv5を使用して 3.5 日強を要する。ステージ1が完了した後、ステージ 2 が実行され、モデルは0.25° で再学習される。再学習には6万4,000 学習ステップを要し、32 個の TPUv5で1.5 日弱を要する。
Cloud TPUv5 デバイスで 15 日間の GenCast 予報を 1 つ生成するには、約8分かかる。
3⃣⃣ アンサンブル･メンバー数
ECMWF-ENSには 50 個の摂動アンサンブル メンバーが含まれているため、GenCastも50 メンバーの アンサンブルを使用した。
(1) GenCastの基本的な枠組み
1⃣ [*53]と同じところ
基本の枠組みは、本論文における参考文献21[*53]を踏襲している。ODEソルバーの各離散求解ステップrθが、時間ステップのインデックスtを付けた残差Z＝ZtをNステップかけて、改良していく。
Zti+1＝rθ(Zti；Xt－1、Xt－2、σi+1、σi)
ZtN＝Zt＝XtーXt－1である。σiは、ノイズレベル･スケジュールと呼ばれるもので、これも[*53]を踏襲している。
σi＝(σmax1/ρ＋i(σmin1/ρーσmax1/ρ)/(N－1))ρ
σmax＝80(サンプリング)、88(学習)
σmin＝0.03(サンプリング)、0.02(学習)
ρ＝7
i＝{0,･･･,N－1}
N＝20
である。
2⃣ [*53]と違うところ
① ノイズ除去器の学習における期待値計算を行う場合の分散が異なる。☞(3)2⃣参照。
② ノイズ除去器の学習における期待値計算を行う場合の重みが異なる。変数レベルごとの損失重みと、そのグリッド セルの正規化された面積で重み付けしている。☞(3)2⃣参照。
(2) オプティマイザー・ハイパーパラメータ
オプティマイザー・・・AdamW
学習率スケジューリング関数・・・コサイン
1⃣ ステージ1
バッチサイズ・・・32
ウォームアップ･ステップ数・・・1000
総学習ステップ数・・・200万
ピーク学習率・・・1.0×10-3
重み減衰率・・・0.1
2⃣ ステージ2
バッチサイズ・・・32
ウォームアップ･ステップ数・・・5000
総学習ステップ数・・・6万4000
ピーク学習率・・・1.0×10-4
重み減衰率・・・0.1
(3) デノイジング･ネットワーク
拡散モデルにおいて、ノイズの除去は重要なプロセスである。ノイズ除去器は、ニュラルネットワークとして実装することができる(つまり、学習することができる)。 基本的なアーキテクチャは、本論文における参考文献21[*53]で提示されているアーキテクチャを踏襲している。
1⃣ ノイズ除去器
ノイズ除去器Dθのアーキテクチャは、自己符号化器である。エンコーダ＋プロセッサ＋デコーダ。プロセッサ部分が、グラフ･トランスフォーマーである。Dθは、[*53]では次のように、実装されている。
Dθ(残差)＝係数1×残差＋係数2×fθ(係数3×残差、係数4)
ここで、係数1＝1/(1＋σ2)、係数2＝σ/√(1＋σ2)、係数3＝1/√(1＋σ2)、係数4＝ln(σ)/4である。σは、ln(σ)～正規分布(平均=-1.2、分散＝1.22)で定義される。時間ステップのインデックスtを付けた残差Zは、Zt＝XtーXt－1である。fθは、ニューラルネットワーク関数であり、本論文では、グラフ･トランスフォーマーにて実装される(という理解で良いだろう)。
2⃣ ノイズ除去器の学習
∑tE(･)を最小化するように学習を行う。E(･)は、次の値の期待値である。
λ(σ)∑i∑jwjai(ZtーYt)2
ここで、
λ(σ)＝1＋(1/σ)2
ln(σ)～正規分布(平均=-1.2、分散＝1.22)
wj：変数レベルごとの損失重みで、GraphCastと同様に設定される。
ai：緯度経度グリッドセルの面積。ただし、グリッド全体で正規化されている。
Yt＝Dθ(Zt＋ε)
Zt：残差
ε ～ Pnoise(·|σ)
i：グリッド内の位置(緯度と経度の座標)をインデックスする。
j：変数をインデックスし、上空変数の場合は等圧面をインデックスする。
t：時間ステップをインデックスする。
である。なお、Pnoiseは通常、独立同一分布(i.i.d.)に従うガウス分布になるように選択される。ただし、地球規模の気象変数の球面形状をより適切に考慮したノイズ分布を使用する方が有益であることがわかった(らしい)。そこで、本論文では、緯度経度グリッドで i.i.d. 白色雑音をサンプリングするのではなく、球面上で等方的(isotropic)ガウス白色雑音をサンプリングし、それをグリッドに射影する。
期待値計算は、σ train、ε ～Pnoise(・;σtrain) で行う。ここで、
σ train＝(σmax1/ρ＋u(σmin1/ρーσmax1/ρ)ρ
u～一様分布[0,1]
である。
【4】比較結果
(0) データセット
1979 年から 2019 年の期間をカバーしている、ECMWFの ERA5 アーカイブから構築されたデータセットで、GenCast は学習された。このデータセットには、利用可能な変数の部分セット、13 の等圧面、0.25°グリッドに対する ERA5 の最良推定解析が含まれている。時間解像度は、1 時間から 6 時間まで、サブサンプリングした。このデータセットから、12 時間の時間解像度でシーケンス を抽出して GenCast を学習した。
GenCast の開発フェーズでは、1979 年から 2017 年までの日付を使用して学習を行い、2018 年に結果を検証した。テスト･フェーズを開始する前に、すべてのモデルと学習の選択を固定し、1979 年から 2018 年までのデータでモデルを再学習し、2019 年に結果を評価した。
(1) アンサンブル･スキルの評価
気象予報の文脈では、予報能力のことをスキルと呼ぶ。スキルを評価する指標は、スキル･スコアと呼ばれる。
1⃣ 評価指標
GenCastとECMWF-ENSについて、㈠CRPS(連続ランク確率スコア)、㈡｢スプレッド･スキルの関係(spread-skill ratio)｣並びに、㈢REVを使って比較評価した。㈠と㈡は、アンサンブル･モデルの評価指標として標準的である。為念：㈠は、小さい方が優れている。㈡は1に近いと優れている。㈢は０以上1以下で、大きい方(1に近い方)が優れている。
2⃣ ターゲット(評価対象変数)
㈠ CRPSの対象
4つの地表変数及び、8つの等圧面における5つの上空変数が対象である。リードタイムは、1～15日(つまり全範囲)。GenCastは、15日間にわたり12時間間隔で予測するので、30セットの予報が出力される。従って、(4＋5×8)×30＝44×30＝1,320個のターゲットが出力される。
地表変数は、①地表から上空2mにおける気温、②平均海面気圧、地表から上空10mにおける風速の③東西方向成分(気象学では、Uと呼称･表記される)及び、④南北方向成分(同Vと呼称･表記される)。8つの等圧面は、200,250,300,500,700,850,925,1000hPaである。上空変数は、⑤風速U、⑥風速V、⑦ジオポテンシャル、⑧気温、⑨比湿、である。
㈡ ｢スプレッド･スキルの関係｣の対象
500hPaのジオポテンシャル、700hPaの比湿、850hPaの気温、地表2mの気温、地表10mの風速の東西成分(つまりU)が対象である。リードタイムは、1～15日(つまり全範囲)。
㈢ REVの対象
99.99 パーセンタイルを超える｢地表2mの気温、地表10mの風速｣が対象である。リードタイムは、1日、5日、7日。
3⃣ 結果
㈠ CRPS
1,320個のターゲットの97.2%(単純に計算すると1,283個)で、GenCast＞ECMWF-ENSと主張している。つまり、GenCastのCRPSが、より小さい。短いリードタイム(12時間及び24時間)×｢⑤U、⑥V及び⑧気温｣×等圧面｢ 200,250,300 ｣でのみ、GenCast＜ECMWF-ENSである。
㈡ スプレッド･スキルの関係
両者とも、ほぼ1であり、本論文では”引き分け”というニュアンスである。ただ細かくみると、①地表から上空2mにおける気温×｢リードタイム～5日まで｣、におけるECMWF-ENSの｢スプレッド･スキルの関係｣は、1より十分小さい(ように思える)。
㈢ REV
｢地表2mの気温、地表10mの風速｣×｢リードタイム1日、5日、7日｣の全てにおいて、GenCast＞ECMWF-ENSである(つまりREVが大きい)。
(2) 熱帯低気圧の進路予測
1⃣ 評価指標
㈠アンサンブル平均･熱帯低気圧軌道の位置誤差及び、㈡進路確率を使ったREVで評価する。㈠のリードタイムは、1～5日、㈡のリードタイムは、{1,3,5}日。位置誤差におけるグランドトルゥースは、TempestExtremes 熱帯低気圧トラッカー🐾13である。ちなみに、ECMWF-ENS予報を 6 時間から 12 時間の解像度にダウンサンプリングして、GenCast との公平な比較を行っている。
進路確率は、特定の時間に特定の 1° グリッド ボックスを、熱帯低気圧の中心が通過すると予測するアンサンブル メンバーの割合として計算される。1° を選択したのは、赤道で 111 km に相当するためである。これは、熱帯低気圧進路確率を定義するために使用される一般的な半径 120 km に近い。
🐾13 TempestExtremes(TE)は、長方形格子または非構造格子/ネイティブ格子上の地域的または全球的な、地球システムデータセットの特徴検出、追跡、科学的分析のための多面的なフレームワーク。TEフレームワーク･バージョン2.1では、熱帯･亜熱帯低気圧、モンスーン性低気圧･低気圧、大気河川、大気ブロック、降水クラスター、熱波など、を調べるための広範なサポートを提供する(らしい)。出典：https://gmd.copernicus.org/articles/14/5023/2021/
2⃣ 結果
㈠ 位置誤差
全てのリードタイム(1～5日)において、GenCastが優れている。図から読み取るしかないので、目安に過ぎないが、10km程度は、誤差が小さい。
㈡ 進路確率
全てのリードタイム(1、３，5日)において、GenCastが優れている。
(3) その他
1⃣ 本論文では、｢予測モデルが、気象状態に固有の空間依存性構造(の一部)を、どの程度うまく捉えているか｣も評価している。CRPSを使ってGenCastとECMWF-ENS予報を比較し、GenCastが優れていた。詳細は割愛。
2⃣ 本論文では、地域風力発電予測という形でも比較を行っている。これは、予測値＝地表から上空10mの風速を使って、最終的に風力発電量をCRPSで比較している。結果は、GenCast＞ECMWF-ENS予報である。これは、風速の予測のみならず、上記1⃣の空間依存性構造の把握を含めて、GenCastが優れている、ことを主張している。
【5】考察
(0) コアな部分は、NVIDIAが22年10月@arXivに公開した論文[*52]に依存している。実用的な機械学習･深層学習モデルの開発力に関するNVIDIAの地力は、相当高いことが伺える。
(1) とは言え、気象予報に(おそらく)初めて拡散モデルを適用し、ECMWF-ENSを超えているのであるから、グーグルは褒められるべきであろう。
(2) GenCastの性能は、白色雑音を基にした摂動ランによるアンサンブル･モデルGenCast-Perturbedの性能を上回っている。つまり、観測データの確率分布を考慮したアンサンブル予報は、白色雑音ベースのアンサンブル予報よりも性能が高い。当然の帰結であろう。そして、これが、GenCastがECMWF-ENS予報に勝った理由ということになるだろう。
海洋研究開発機構(JAMSTEC)は、｢大気変数を、高精度に予測できる学習モデルを構築することに成功した｣と主張する論文(以下、本論文[*21])を発表した(24年3月20日)。本論文の学習モデルは、シグネチャを”説明変数”として使っている。
シグネチャを使った学習モデルは、時系列解析以外では、あまり用いられていない(シグネチャが何であるかは後述する)。数少ない(?)例外は、同じJAMSTECによる先行研究[*22]である。本論文は、[*22]の枠組みを、大気に適用した研究と位置づけられる。
❚為念1❚ 米国の量子コンピュータ･ベンダーRigetti Computingは、該社公式ブログに、｢シグネチャを使った景気後退予測モデルの量子版は、古典版より精度が高い｣という主張を投稿している(23年4月)。その投稿では、シグネチャを使った予測モデル･古典版がプロビット･モデルと比較されている(シグネチャを使ったモデルが、やや性能が高い。こちらを参照)。
❚為念2❚ シグネチャという文言は、多様な分野で使用されている(が、もちろん意味は、全く異なる)。例1：サイバーセキュリティの分野で、シグネチャと言えば｢マルウェアの帰属を特定･判別するために用いられるデータ｣を指す。シグネチャ･ベースのセキュリティ、といった表現が多く見られる。例2：生物学･医学の分野では、｢健常者群、がん患者群などの2群間を分類する際に使用される、特徴的な遺伝子｣をシグネチャと呼ぶようである。
【1】本論文の主張
(1) シグネチャを使用した学習モデルは、温度と水蒸気混合比†1の絶対平均差が最小限で、高レベルの精度を示した†2。
(2) シグネチャを使用した学習モデルは、激しい降雨時であっても、水蒸気と温度の急激な変動を含む、鉛直構造と大気の不安定性をうまく捉えることができた。
(3) シグネチャを使用した学習モデルは、生の説明変数を使用した学習モデルよりも、精度が高かった†3。
†1 水蒸気混合比[単位：g/Kg]＝水蒸気質量[g]/乾燥空気質量[Kg]
†2 グランドトルゥースは、気象庁の運用メソスケール･モデルの解析値である。
†3 精度の計測指標は、差と二乗平均平方根誤差(Root Mean Squared Error：RSME)である。
【2】事前整理
(0) 前振り：ラフパス理論とシグネチャ
ラフパス理論は、確率論における比較的新しい概念と紹介される。数学的に言えば、｢確率微分方程式の求解法が、微分方程式の求解法と同様に扱える枠組みを提供する｣。曰く、ラフパスの空間に、適当な位相を導入すると、線積分が測度とは関係なく定義できる[*23]。
以下では数学ではなく、あくまで機械学習モデルを議論の対象とする。さらに、機械学習モデルは線形回帰モデルを考える。誤解を恐れずに言うと、線形回帰モデルの範疇であれば、ラフパス理論におけるパスとは、”説明変数と目的変数のセット”に過ぎない。さらに線形回帰モデルの範疇では、シグネチャは、説明変数から新たに作られる”｢合成｣説明変数のセット”に過ぎない。
ここで、シグネチャがセットであることについて少しだけ述べる。シグネチャは、次のように表される(詳細は、(2)1⃣を参照)：{1、S(1)(X)、S(2)(X)、S(1,1)(X)、S(1,2)(X)、･･･}。ここで、Xはパスであり、S(0)(X)＝1である。もっとも、S(0)(X)という表記はされない。S(1)(X)などは、シグネチャの｢各項｣といった呼ばれ方をする[*24]。
なお、パスの和訳は経路(あるいは、道、路)であるが、経路と訳してしまうと、遥かに有名な｢ファインマンの経路積分｣と混同しかねないので、日本語を使わずにパスとする。
(1) シグネチャによる線形回帰モデルとは、何ぞや？
｢生の｣説明変数による線形回帰モデルは、
目的変数＝∑係数†1×｢生の｣説明変数＋定数†2
という形をとる。シグネチャを使うと、
目的変数＝∑係数×｢合成｣説明変数＋定数
となるだけである。もちろん、シグネチャを顕に書けば
目的変数＝∑係数×シグネチャ＋定数
となる(だけである)。では、なぜ、こんなことをするかと言うと、シグネチャという合成説明変数が、｢生の｣説明変数より、説明能力が高いと考えられる(場合がある)からである。機械学習･深層学習寄りの言葉で言うと、特徴量エンジニアリングをしている、とイメージしても良いだろう。
違った例でイメージを表すと、正規分布における、平均と分散を考えると良いかもしれない。正規分布は、1次モーメント＝平均と2次モーメント＝分散を指定すれば、完全に表現できる(3次モーメント＝歪度、4次モーメント＝尖度、･･･、の全ての高次モーメントは、1次モーメントと2次モーメントから計算できる)。つまり、生のデータから平均と分散という合成量を生成することで、データの分布(正規分布)を完全に表現できる(👉下記(2)2⃣も参照)。この例では、平均と分散を使った回帰モデルでデータの分布(正規分布)を表現するわけではないから、かえって混乱するかもしれないものの、説明能力が完璧な合成量という意味では、正規分布における平均と分散は、良い例だと思う†3。
†1 この係数は、線形回帰の文脈では回帰係数と呼ばれる。深層学習モデル(ニューラルネットワーク)だと、｢重み｣と呼ばれる。
†2 機械学習･深層学習の文脈では、バイアス項と呼ばれる。
†3 混乱ついでに、情報が集約された”合成量”という意味では、テンソルネットワーク(今の場合であれば、もっともプリミティブな行列積状態を考えれば良い)を上げることもできるだろう。
(2) 数学的な補足
1⃣ シグネチャと積分
シグネチャの各項(以下、本稿では、｢シグネチャ各項｣あるいは｢シグネチャ項｣という文言と区別しない)は、パスの反復積分として定義される。この積分は、リーマン･スティルチェス積分の一種(特殊な形)である†1。先の表記S(1)(X)等、を使うと、段数†21の場合であれば、
S(1)(X) = ∫dX(1)
と1重積分として表現される。幾何学的に言うと、これは、線分である。積分区間は明示していないが、[a,b]のように適当にとる。実際は、S(1)a,b(X)のように下付き添え字で表記する。積分記号にも、∫a＜s＜bとついて、その場合には、積分測度dXにもdXsと付く(ので、見た目が煩わしい)。
段数2だと(そして、ここでは、添え字をフルに付けると)
S(1,2)a,b(X) ＝∫a＜s＜bS(1)a,s(X)dX(2)s
＝∫a＜r＜s＜bdX(1)rdX(2)s
＝∫ab∫asdX(1)rdX(2)s
と二重積分として表現される。段数2のシグネチャ項の幾何学的解釈は、直接(ストレートには)、面積とはならないが、段数2のシグネチャ項を使って面積を表現することは容易である(なお、段数がより高いシグネチャ項の幾何学的解釈は、残念ながら、全く直感的ではない)。S(1,2)(X)の計算に、S(1)(X)が使われているところが、『反復』積分という名前の所以である。
†1 より一般的には、ヤング積分に拡張できる。リーマン･スティルチェス積分の範囲では、有界変動を持つパスに対する反復積分のみが定義可能である。ヤング積分にまで拡張すると、”もっと不規則な”パスに対する反復積分まで定義可能となる。なお、リーマン･スティルチェス積分を簡潔に述べると、∫f(x)dg(x)と表させる積分のことである。
†2 文言は、[*24]に従った。
2⃣ シグネチャとモーメント
(1)で、シグネチャ項と、｢平均、分散｣(統計的モーメント)とのイメージ的なアナロジーについて触れた。実は、シグネチャ項を使って統計的モーメントを簡単に計算することができる。しかも、それは2次までのモーメントに留まらず、高次のモーメントにまで当てはまる。詳しくは、[*25]を参照。
3⃣ シグネチャと基底関数
(0)では、線形回帰の議論の範疇であれば、シグネチャ各項は、｢合成｣説明変数に過ぎないと書いた。
一応、次のような問題を考えよう：説明変数と目的変数との間には確かに何らかの関係があって、回帰タスクにより、その関係は顕示可能とする。その前提下であっても、勝手に作った｢合成｣説明変数と目的変数との間に、何らかの関係が確かに存在すると期待しても良いのか。
シグネチャ項を｢合成｣説明変数とする場合の答えは、イエスである。シグネチャ項は、パスの関数が作る空間における基底(関数)と見做すことができる。つまり、パスXの連続関数f(X)†1は、シグネチャ項の線形結合で近似†2できる。もっと平たく言えば、説明変数X、目的変数 ＝Xの適当な関数値f(X)であるから、目的変数は、シグネチャ項の線形結合により、近似できる。連続関数f(X)は、もちろん非線形で構わない。多変量回帰も、非線形関数を線形結合で表現する枠組みであるが、シグネチャ項を使った表現は、多変量回帰よりも優れている場合があると期待される(できる)。該当する”場合”については、下記(3)を参照。
|蛇足|・・・再度、混同しかねない正規分布の平均と分散で例えると、平均と分散という合成量(統計量)は説明能力が完璧であったが、正規分布に対してのみ有効である。一方、シグネチャ(各項)は、任意の(連続)関数に対して有効である。故に、基底ということになる。
†1 正確に言うと、有界変動を持つパスの集合のコンパクト部分集合上の連続関数。有界変動≒区分的に微分可能、である。先述の通り、有界変動の制限は、パスの反復積分＝シグネチャが、リーマン･スティルチェス積分として存在することを保証するための制限であった。
†2 任意の精度で近似可能。
(3) シグネチャを使う意味
データの持ち方で学習モデルの性能が変わるという意味では、グラフ(構造)を使ったニューラルネットワーク｢グラフ･ニューラルネットワーク(GNN)｣をイメージするといいかもしれない。GNNは、構造化データを適切に処理できると考えられている。例えば、分子シミュレーションの代理モデルであれば、分子をグラフ構造で捉えて、NNで学習することにより、他よりも高性能な学習モデルとなる(と言われている)。つまり、分子(が示す様々な物性値等)をNNの入力(説明変数)として考えた場合、グラフ構造で表現することがより望ましい(だろう)という”ドメイン知識”を使って、データを洗練させると、学習モデルの性能はあがるという例である。
すると、シグネチャは、どういう”ドメイン知識”が適用できるかという話になる。シグネチャはパスの幾何学的形状の影響を受けやすい[*25]。この性質を、分かり易く例えれば、｢”カクカクしたデータ”を扱うケースでは、シグネチャを説明変数として使った学習モデルは、性能が高いと期待できる｣と表現できるだろう。具体的には、しばしば引き合いに出される時系列データあるいは、もっと直截的に株価データを扱うケースで、シグネチャを使った学習モデルは高性能であることが期待できる。
先行研究[*22]では、水深でスライスした観測プロファイル(圧力、塩分濃度、水温)†が、”カクカク”しているという”ドメイン知識”をもとに、シグネチャを使った回帰モデルを、予測モデルとして採用している(と思われる)。
† ちなみに、観測範囲は、水深約2,000mから海表面まで。
(4) ペナルティ項
パスは、そのシグネチャによって完全に決定されるのか?･･･という問いに対する答えは、実は、ノーである(が、大きな問題はないとされている)。また、シグネチャ項には、共線性の問題が発生する。これは、回帰タスクを実行する際に障害となるので、解決しなければならない。標準的な解決策の一つは、ペナルティ項の導入(正則化)である。具体的には、ラッソ、リッジ、エラスティックネット正則化が行われる。先行研究[*22]では、ラッソ正則化が採用されている。本研究では、リッジ正則化が採用されている。
【3】枠組み
(0) 先行研究の枠組み
先行研究[*22]のモチベーションは、観測プロファイルの品質管理を自動化したい、というものであった。品質管理は、合格･不合格の二値管理である。シグネチャ項の線形和で、スコアを算出し、設定したしきい値を越える･越えないで、合格･不合格とする。
シグネチャ項を使った学習モデルは、主成分分析法などの従来手法より、精度が高いと結論されている。本研究は、この先行研究に触発されて、大気プロファイルに対して、 シグネチャ項を使った回帰モデルを適用している。
(1) データのセットアップ
大気データは、水平解像度10 km を採用する気象庁の運用メソスケール･モデルの解析値である。地表面及び16の等圧面†1における、3時間毎の｢圧力、温度、相対湿度†2｣で構成される。地点は、福岡県(北緯 33.6 度、東経 130.4 度)。期間は、2019年から2020年で、合計5,848サンプルから構成されている(2年×365日×24時間/3時間＝5,840)。データは、無次元値化されている†3。
また、偏りのない分析にするため、データは時間軸に沿ってランダムにシャッフルされた。データセットは、80%と20%の割合で分割された。前者(80%)は学習データとして、後者(20%)は検証データセットとして使用された。
†1 1000、975、950、925、900、850、800、700、600、500、400、300、250、200、150、100hPa
†2 ただし、250hPaを超える等圧面における相対湿度は、ゼロに設定されている。
†3 1000hPa、100hPa、1hPaの値を使用して無次元値化された、と書かれている。
(2) モデルのセットアップ
ニューラルネットワークは、全結合の隠れ層(2つ。最初の層は32 ノード、2 番目の層は128ノード) と出力層(364 ノード)で構成。シグネチャは、Pythonライブラリesig†を使用して計算され、反復積分の次数は、5に設定された。最適化計算に使用したソフトウェア･ライブラリは、SciPy。オプティマイザは、BFGSを使用した。
† シグネチャを計算するPythonのライブラリには、esigの他に、iisignature[*26]、Signatory[*27]がある。
(3) 本研究の枠組み
2年間(2019年及び2020年)のデータで学習済のモデルに対して、2021年のデータを入力し、2021年の出力値を予測値として、グランドトルゥースと比較する。2021年のデータも3時間ごとなので2,920サンプルで構成される(1年×365日×24時間/3時間＝2,920)。
(4) 為念：｢合成｣説明変数としてのシグネチャ
【2】(0)でシグネチャ項は、｢合成｣説明変数に過ぎない、と書いた。本研究における説明変数は、気象学的に意味のある物理量であるから、｢合成｣説明変数も、気象学的に意味のある物理量でなければ平仄が合わない。幸い本論文には、シグネチャ項(反復積分)の持つ物理的な意味が、いくつか述べられている。
⇒ 圧力と水蒸気(相対湿度か?)からの2次反復積分(シグネチャ項)は、｢降水可能水蒸気｣を示す。そして、圧力と温度からの反復積分(シグネチャ項)は、大気中の総熱量を示す。
【4】シグネチャを使用した学習モデルの検証結果
(0) グランドトルゥース及び比較指標
グランドトルゥースは(既述通り)、気象庁の運用メソスケール･モデルの解析値である。
比較指標は(こちらも、既述通り)、差と二乗平均平方根誤差(Root Mean Squared Error：RSME)である。
(1) 比較対象モデル
❶シグネチャを説明変数とした深層学習モデルと、❷ベースライン：｢生の(raw)｣大気プロファイル値を説明変数とした学習モデルが、比較された。生の大気プロファイル値は、｢圧力、温度、相対湿度｣である。ベースラインもモデル･アーキテクチャは、同じである。すなわち、全結合の｢隠れ層×2と出力層｣で構成されている。隠れ層のノード数は、16と32。出力層は、48ノードである。
目的変数は、気温と水蒸気混合比である。2021年8月における❶と❷それぞれの予測値が比較された。
(2) シグネチャ項を使った学習モデルの精度
2021年における、㊀年間平均値、㊁夏の期間に渡る平均値、㊂冬の期間にわたる平均値を、気温と水蒸気混合比において算出している。その上で、差とRSMEで比較している。概ね、学習モデルの予測精度は、高いと言えるであろう。
差とRSMEで比較するとRSMEの方が(グランドトルゥースとの)乖離は大きいが、水蒸気混合比は気温に比べて、その乖離幅は小さい。乖離幅に関して、気温は㊀～㊂で大きな違いはないが、水蒸気混合比は㊂が明らかに小さい。気温は、高高度になるほど乖離が大きくなる(250hPaを超える等圧面における相対湿度は、ゼロに設定されていることから、水蒸気混合比は高高度ではゼロに収束している)。
気温の乖離は、高度の上昇に伴って、長らくマイナス方向に振れた後、短くプラス方向に振れる。その後、再びマイナスに振れる。この傾向は㊀～㊂で同じである。ただし、実際にプラスの乖離を示すのは、㊂のみである。また、㊀と㊂は800hPa等圧面まで、乖離がほぼ0。㊁は600hPa等圧面まで、乖離がほぼ0。
(3) ❶シグネチャ項を使った学習モデルと、❷ベースラインとの比較
結論として、全体的に見て、❶は❷より精度が高い、と言える。折角なので意地悪な視点で、重箱の隅を突いてみよう。気温で、❶と❷を比べると、高高度かつRSMEでは、むしろ❷の方が精度が高いと言える。しかし、❷は地表面付近では、(グランドトルゥースとの)乖離が顕著に大きい。水蒸気混合比は、全体的に❶の精度が高い。地表面付近で❷の乖離が顕著に大きいという傾向は、同じである。本論文では、｢中層(975hPa～800hPa)では持続的な正のバイアスが観測された｣と指摘されている。そして、その原因について、｢入力ベクトルと、表層,大気境界層,自由大気層の値との間に、❷モデル内で相関があることに起因していると考えられる｣と推論している。
【5】考察
(0) 日本の気象予報は、数値予報に関しては、世界でも高水準にあると思われる。ただ、機械学習･深層学習(を使用･適用して、付加価値を付けるあるいは性能向上を図る)に関しては、遅れていたように思う。本研究は、嬉しい例外の一つとして上げられるだろう。
(1) シグネチャを使った学習モデルが、金融や経済分野における時系列データ解析以外で使われることは珍しい。気象分野で使われたのは、本研究が初めてらしい。見事にハマった、という感じであろうか。
(2) 【0】でも取り上げた通り、シグネチャを使った学習モデルは、既に量子版を検討するフェーズにある。量子化することで、性能が上がると期待できる合理的な理由があるからである。気象分野にあてはめたこの事例も、量子化を期待したい。
【尾注】
*1 https://deepmind.google/discover/blog/graphcast-ai-model-for-faster-and-more-accurate-global-weather-forecasting/
*2 Remi Lam et al.、Learning skillful medium-range global weather forecasting、https://www.science.org/doi/reader/10.1126/science.adi2336
*3 Kaifeng Bi et al.、Accurate medium-range global weather forecasting with 3D neural networks、https://www.nature.com/articles/s41586-023-06185-3
*4 Jaideep Pathak et al.、FOURCASTNET: A GLOBAL DATA-DRIVEN HIGH-RESOLUTION WEATHER MODEL USING ADAPTIVE FOURIER NEURAL OPERATORS、https://arxiv.org/pdf/2202.11214.pdf
*5 https://www.jstage.jst.go.jp/article/jwea/45/2/45_261/_pdf
*6 https://climcore.rcast.u-tokyo.ac.jp/reanalysis/
*7 https://www.ncei.noaa.gov/products/international-best-track-archive
*8 松枝未遠･中澤哲夫、TIGGE データを利用した顕著現象発生予測プロダクトの開発とその評価、https://www.dpac.dpri.kyoto-u.ac.jp/workshop/2014/proceedings/07-matsueda.pdf
*9 http://agora.ex.nii.ac.jp/digital-typhoon/help/world.html.ja
*10 http://agora.ex.nii.ac.jp/digital-typhoon/help/landfall.html.ja
*11 以下を参考にした：気象庁･数値予報解説資料集(令和4年度) 4.7 表記と統計的検証に用いる代表的な指標)：https://www.jma.go.jp/jma/kishou/books/nwpkaisetu/R4/4_7.pdf
*12 他には、ミシガン大学、ライス大学、カリフォリニア工科大学、パーデュー大学。全て米国の機関である。
*13 https://www.jma.go.jp/jma/kishou/books/nwpreport/64/chapter5.pdf
*14 https://developer.nvidia.com/ja-jp/blog/develop-physics-informed-machine-learning-models-with-graph-neural-networks/
*15 John Guibas et al.、ADAPTIVE FOURIER NEURAL OPERATORS: EFFICIENT TOKEN MIXERS FOR TRANSFORMERS、https://arxiv.org/pdf/2111.13587.pdf
*16 https://www.mizuho-rt.co.jp/publication/column/2022/infocomm0317.html
*17 関口智大、アーティスト制御可能なオーロラシミュレーションに関する研究、https://core.ac.uk/download/pdf/188016534.pdf
*18 池田翔他、気象庁全球週間アンサンブル予報のダウンスケールデータを用いた相対湿度および葉面濡れの確率予報実験と検証、日本気象学会機関誌｢天気｣69巻(2022年)3号、pp3-18、https://www.jstage.jst.go.jp/article/tenki/69/3/69_133/_pdf/-char/ja
*19 https://www.science.org/doi/suppl/10.1126/science.adi2336/suppl_file/science.adi2336_sm.pdf
*20 https://developer.nvidia.com/ja-jp/blog/develop-physics-informed-machine-learning-models-with-graph-neural-networks/
*21 Mikiko Fujita et al.、Prediction of Atmospheric Profiles With Machine Learning Using the Signature Method、https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023GL106403
*22 Nozomi Sugiura、Machine Learning Technique Using the Signature Method for Automated Quality Control of Argo Profiles、https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2019EA001019
*23 稲浜譲、論説|ラフパス理論と確率解析、数学67巻(2015) 第3号、pp. 291-313、https://www.jstage.jst.go.jp/article/sugaku/67/3/67_0673291/_pdf/-char/ja
*24 杉浦望実、特集「データ同化の方法」［研究詳解］シグネチャ法入門、統計数理(2022) 第70巻 第2号、pp.251–267、https://www.ism.ac.jp/editsec/toukei/pdf/70-2-251.pdf
*25 Ilya Chevyreva ＆ Andrey Kormilitzin、A Primer on the Signature Method in Machine Learning、https://arxiv.org/pdf/1603.03788.pdf
*26 Jeremy Reizenstein、The iisignature library: efficient calculation of iterated-integral signatures and log signatures、https://arxiv.org/pdf/1802.08252.pdf
*27 Patrick Kidger and Terry Lyons、SIGNATORY: DIFFERENTIABLE COMPUTATIONS OF THE SIGNATURE AND LOGSIGNATURE TRANSFORMS, ON BOTH CPU AND GPU、https://openreview.net/pdf?id=lqU2cs3Zca
*28 https://www.jamstec.go.jp/j/about/press_release/20240315/
*29 Cristian Bodnaret al.、AURORA: A FOUNDATION MODEL OF THE ATMOSPHERE、https://arxiv.org/pdf/2405.13063
*30 Ze Liu et al.、Swin Transformer V2: Scaling Up Capacity and Resolution、https://openaccess.thecvf.com/content/CVPR2022/papers/Liu_Swin_Transformer_V2_Scaling_Up_Capacity_and_Resolution_CVPR_2022_paper.pdf
*31 https://blog.google/intl/ja-jp/products/explore-get-answers/search-japan-nowcast-hashtag/
*32 Marcin Andrychowicz et al.、Deep Learning for Day Forecasts from Sparse Observations、https://arxiv.org/pdf/2306.06079
*33 関山剛、「深層学習を使った気象シミュレーション代理モデルの可能性」、https://cesd.aori.u-tokyo.ac.jp/fugaku/pdf/20220311_sekiyama_pub.pdf
*34 Jie Hu et al.、Squeeze-and-Excitation Networks、https://arxiv.org/pdf/1709.01507
*35 https://jp.weathernews.com/news/47907/
*36 Hao-Yan Liu et al.、A Hybrid Machine Learning/Physics-Based Modeling Framework for 2-Week Extended Prediction of Tropical Cyclones、https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2024JH000207
*37 Dmitrii Kochkov et al.、Neural general circulation models for weather and climate、https://www.nature.com/articles/s41586-024-07744-y
*38 気象庁、令和5年度数値予報解説資料集｜第1章基礎編1.7.2全球モデル、https://www.jma.go.jp/jma/kishou/books/nwpkaisetu/latest/1_7_2.pdf
*39 三好建正、解説｜ビッグデータ同化と気象予測、応用物理 第90巻 第8号(2021)、pp.470-475、https://www.jstage.jst.go.jp/article/oubutsu/90/8/90_470/_pdf
*40 https://www.jma.go.jp/jma/kishou/books/nwptext/51/2_chapter4.pdf
*41 富田浩文、全球非静力学大気モデルの計算方法－最近の動向と今後の方向－、https://www2.ccs.tsukuba.ac.jp/people/kuramasi/algorithm10/slides/tomita_110208.pdf
*42 https://www.jma.go.jp/jma/kishou/books/nwptext/45/1_chapter4.pdf
*43 榎本剛、新用語解説｜三次八面体逓減ガウス格子、日本気象学会機関誌｢天気｣64巻(2017年)、No.3、pp.189-191(☛https://www.metsoc.jp/tenki/result.php。通し番号だとそうなのだろうか。PDFだとpp.57-59)、https://www.metsoc.jp/tenki/pdf/2017/2017_03_0057.pdf
*44 Massimo Bonavita、On Some Limitations of Current Machine Learning Weather Prediction Models、https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023GL107377
*45 https://www.gfdl.noaa.gov/shield/#X-SHiELD
*46 Yogesh Verma et al.、CLIMODE: CLIMATE AND WEATHER FORECASTING WITH PHYSICS-INFORMED NEURAL ODES、https://arxiv.org/pdf/2404.10024
*47 Simon Lang et al.、AIFS - ECMWF’S DATA-DRIVEN FORECASTING SYSTEM、https://arxiv.org/pdf/2406.01465
*48 https://www.kyushu-u.ac.jp/ja/researches/view/1137
*49 https://image.itmedia.co.jp/l/im/aiplus/articles/2410/10/l_dy_n_13.jpg
*50 Ilan Price et al.、Probabilistic weather forecasting with machine learning、https://www.nature.com/articles/s41586-024-08252-9.pdf
*51 Jimeng Shi et al.、CoDiCast: Conditional Diffusion Model for Weather Prediction with Uncertainty Quantification、https://arxiv.org/pdf/2409.05975
*52 鈴木大慈、学術変革領域研究(A)｢学習物理学の創成｣領域会議 拡散モデルチュートリアル(+α)、2023年9月25日、https://speakerdeck.com/taiji_suzuki/kuo-san-moderutiyutoriaru?slide=55
*53 Tero Karras et al.、Elucidating the Design Space of Diffusion-Based Generative Models、https://arxiv.org/pdf/2206.00364
*54 全体レポートは巨大(456頁)なので、各章ごとに分割してアクセスできる。第5章は、こちら：https://hai-production.s3.amazonaws.com/files/hai_ai-index-report-2025_chapter5_final.pdf
ちなみに、全体はこちら：https://hai-production.s3.amazonaws.com/files/hai_ai_index_report_2025.pdf
*55 https://hai-production.s3.amazonaws.com/files/hai_ai-index-report-2024_chapter5.pdf
