(1) 革新的データ同化の創出を目指して

科学研究を進める上において，物理・化学法則等に基づく数値モデルと，観測・実験に基づくデータの比較が重要であることは論を待たない．しかしながら，近年の巨大スパコンの登場や大規模地球観測網・実験設備等の整備に伴い，大規模数値モデルと大容量観測データを突き合わせることすら容易ではなくなってきた．数値モデルと観測データをベイズ統計学の枠組みで統融合するための計算技術であるデータ同化は，時々刻々と入力する観測データに基づいて各時刻における状態の逐次推定を行う「逐次データ同化」と，予め決められた時間窓において観測データと最も整合する状態を探索する「非逐次データ同化」とに大別される．逐次データ同化で行う状態推定においては，アンサンブルカルマンフィルタや粒子フィルタに代表される逐次ベイズフィルタが用いられるが，ナイーブな実装をしてしまうと，世界最大のスパコンを以てしても計算機資源が全く不足するという事態が簡単に起こる．大規模数値モデルへデータ同化を実装する際には，4次元変分法を始めとする非逐次データ同化を用いるのが常套であり，例えば気象予報は主に4次元変分法に基づいて行われている．実装の手間は逐次ベイズフィルタよりも大きくなるが，必要な計算コストは格段に小さくて済む．
従来の4次元変分法は，事後分布の局所最大を与える状態を推定するのみであり，その不確実性を推定することが原理的に不可能であるという大きな欠点があった．例えば，台風の進路予測においては，中心位置の不確実性を「予報円」によって表現するが，これは従来の4次元変分法は算出することができず，アンサンブル計算によって求められている．すなわち，現在の天気予報では，非逐次／逐次データ同化をアドホックに組み合わせて実施されているのが実情である．我々は，2nd-order adjoint法を採り入れることにより，不確実性評価が可能な4次元変分法を開発することにより，これを解決した（後述(2)を参照）．このようにして得られた不確実性は，観測デザイン最適化のためのフィードバックともなる極めて重要な情報である．
それでもなお，現行の4次元変分法には，初期解に強く依存した局所最適解しか得られないという問題がまだ残されており，固体地球科学の実問題に適用するためには，これを乗り越えるようなアルゴリズム開発が求められる．本年度は，後述(3)の研究でも用いたレプリカ交換モンテカルロ法をプラグインすることにより，「革新的4次元変分法」とも言うべき，初期解に依存しない大域的最適解を求めることが可能な4次元変分法のアルゴリズム開発に着手した．

(2) 4次元変分法による大規模データ同化に基づく予測不確実性評価法の開発

非逐次型データ同化法である4次元変分法は，フォワード計算を実施するための数値モデルから「アジョイントモデル」と呼ばれる方程式系を構成する必要があるため，アンサンブルカルマンフィルタや粒子フィルタのような逐次ベイズフィルタを用いたデータ同化手法と比較して，実装は容易ではないものの，必要なメモリ容量と計算時間を数値モデルと同等に抑える事ができるため，連続体計算のような自由度の大きい数値モデルを用いるデータ同化を実施する際には，極めて有効である．しかしながら，従来の4次元変分法に基づくデータ同化では，勾配法による最適化を行うため，逐次型データ同化法では自然に得ることができる推定値の不確実性を評価することが不可能であった．昨年度本センターでは，その4次元変分法の弱点を解決する新アルゴリズムを開発し，大自由度系の数値モデルに対する推定値の不確実性の高速かつ高精度な評価を可能した．しかしながら、この手法は変数の初期値に対する不確実性評価法であり，その得られた不確実性が未来の状態の予測へどの程度影響するかは未知であった．この未来の状態予測への影響を評価することにより，予測の高精度化・観測へのフィードバック等へ利用できるため，その評価法を確立することは喫緊の課題である．そこで，本年度は昨年度の手法をさらに高度化させ，得られた推定値の不確実性が系の将来予測に及ぼす影響（予測不確実性）を調べるアルゴリズムを開発した．一般に，不確実性の伝搬を評価するためには，事後分布の時間発展を計算する必要があるが，それを直接評価すると自由度の指数関数オーダーに比例する計算コストを必要とするため現実的ではない．そこで本センターでは，事後分布を直接評価しない代わりに，事後分布を最大にする解の時間発展を考えることで，計算コストを劇的に抑える計算手法を開発した．具体的には本提案手法では，まず，事後分布の変数をモデル変数のベクトルZ(t)とパラメータベクトルmに分解する．そしてパラメータベクトルmの推定・不確実性評価を昨年度開発した手法により行ない，最適推定値m’とその不確実性?mを得たのち，それらを用いて3つの解，Z(t|m’)，Z(t|m’+m)，Z(t|m’-m)の時間発展を計算する．事後分布の時間発展を直接評価することは困難であるが，これらの解の時間発展は元の数値モデルと同程度の計算コストで計算できる．これらの時間発展を追うことで，パラメータの不確実性がモデル変数の時間発展に及ぼす影響を定量的に評価することが可能となる．この手法は系のパラメータの自由度がモデル変数の自由度に比べて少ない場合には，事後分布の直接評価に比べて，圧倒的に計算コストを軽減できる．本提案手法を多相系の大規模フェーズフィールドモデルに適用したところ，モデルパラメータの不確実性に依存した系の時間発展の予測不確実性を高速かつ高精度に評価できることを確認した．本提案手法は，一般の自励系モデルに対して適用できるため，浅水方程式などの津波モデルや，断層運動を取り扱う弾性体モデルなど，固体地球科学で用いられる様々な大規模シミュレーションモデルへの応用が可能である．

(3) データ駆動型モデリングに基づく首都圏を中心とした地震動イメージング技術の開発

本センターは，防災科学技術研究所データプラットフォーム拠点形成事業（防災分野）「首都圏を中心としたレジリエンス総合プロジェクト（首都圏レジリエンスプロジェクト）」に参画し，稠密な観測網に基づく都市の地震被害評価技術の開発を実施している．
巨大地震発生時に，首都圏を中心とした都市部における構造物の地震応答を即時的に評価することは，構造物の被害推定だけでなく，地震後の迅速な救助活動や二次災害の軽減に繋がる．この地震応答の計算の際には，初期条件として構造物直下における入力地震動が必要となるが，すべての構造物において地震動を直接観測することは現実的でない．そこで，本センターではこれまで，観測された地震データから，地震観測網よりも圧倒的に密に分布した構造物直下における入力地震動のイメージング手法の開発を行ってきた．
昨年度までに，観測データから入力地震動を推定するスパースモデリングに基づくデータ駆動型モデリング手法に加え，波動方程式を物理的な拘束条件として課すことで観測データをより定量的に説明する入力地震動をイメージングする手法を開発した．具体的には，波動方程式に含まれる地震波速度や層厚といった地下構造に関するパラメータや，震源に関するパラメータを未知変数とし，観測データと数値モデルの定量的な適合度を示す事後分布を定義する．一般に，この事後分布は極めて複雑な形をしており，解析的に最適化を行うことは容易ではない．そこで，マルコフ連鎖モンテカルロ（MCMC）法に基づく事後分布からのサンプリングを行い，未知パラメータの最適化を図りながら，同時に入力地震動のイメージングを実行するという方針をとった．MCMC法としては，多峰性を持つ分布からでも効率的にサンプリングすることが可能なレプリカ交換モンテカルロ（REMC）法を採用した．レプリカ交換モンテカルロ法は，サンプリングを実行すべき分布に加え，それを滑らかにした形状を持つ複数の分布から同時にサンプリングを行う．このとき，適当なタイミングでサンプルの交換を分布間で行うことにより，広範囲かつ効率的な探索を行うことが可能である．
本年度はREMC法を用いて，関東地方の地震動イメージングを行い，速度応答スペクトルを用いて構造物の揺れの簡易的な評価を行った．関東地方では，首都圏における地震像の解明を目的として，2007年度以降，都心部を中心に数kmの観測点間隔でおよそ300点の地震計が設置されている（首都圏地震観測網，MeSO-net）．本年度は，1923年大正関東地震時に被害が大きかった地域の一つである東京都北東部の10km四方の領域を対象に，MeSO-netで得られた観測波形を用いて地震動イメージングを行った．その結果，高層建築物で卓越する周期5-10秒の長周期地震動に対して，観測波形の大部分を説明可能な地震動のイメージングに成功した．このような長周期の帯域では，観測波形が再現されている上，推定された応答スペクトルも観測から得られる応答スペクトルと良い一致を示した．一方，中小規模の建物を含む一般的な構造物は0.2-0.5 秒程度の周期帯が卓越する．しかしながら，被害を受けた構造物の卓越周期は長くなることから，一般的な構造物の大規模被害のみを想定する場合は周期1秒以上の地震動を評価すれば十分であると考えられる．そこで，様々な規模の構造物の揺れの評価に向けて，周期1-10秒の地震動イメージングを行ったところ，振幅の大きな成分の直達波の地震動がある程度再構築でき，また応答スペクトルを再現することに成功した．この結果は，構造物の応答評価という観点において，地震動イメージング手法が1秒程度の短周期帯まで適用可能であることを示している．また，得られた観測データから，REMC法による地震動イメージング，また速度応答スペクトルまでの一連の計算を可能とする地震動イメージングモジュールの開発を行っている．今後，地震動イメージング手法の更なる高度化や高速化により，将来的に地震発生時の即時的な被害推定や二次災害の軽減に貢献することが期待される．