(1) シミュレーション／データ両駆動型データ同化の創出を目指して

データ同化は，数値シミュレーションモデル（以下，数値モデル）と観測データをベイズ統計学の枠組みで統融合するための計算基盤技術であり，高精度かつ高詳細な予報が可能となった現代の気象予報は，データ同化の賜物と言っても過言ではない．データ同化は，本センターが取り扱っているような大規模な数値モデルを対象とする場合において，特にその威力を発揮する．しかしながら，そのような大規模数値モデルに対しても，データ同化の本来の目的である確率密度関数（以下，分布）の推定をきちんと行うためには，従来のデータ同化手法では不十分であった．本センターでは，(2)で述べる4次元変分法に基づくデータ同化手法の革新的アルゴリズムの新規開発に成功した．
また，これまでのデータ同化の枠組みは，所与の数値モデルに大きく依存する「シミュレーション駆動型モデリング」であったが，これが効果的に機能するためには，その数値モデルが起こり得る現象を精緻に記述したものであることが要求される．しかしながら，これが可能な研究分野は極めて少数派であり，例えば地震関連の観測データには，地下構造に関する情報が不十分であること等に起因して，数値モデルでは再現しきれない現象が非常に多く含まれていると思われる．そこで本センターでは，他の研究部門や研究センターとも連携して「データ駆動型シミュレーション プロジェクト室」を2014年度に発足し，従来のシミュレーション駆動型データ同化の枠組みに，スパースモデリングに代表されるデータ駆動型モデリングの手法をプラグインすることにより，シミュレーション／データ両駆動型のデータ同化法の創出を目指している．地震データに適用可能なデータ駆動型モデリング手法の第一歩として，(3)で述べるように，限られた地震観測データから空間的に密な波動場を推定する「地震動イメージング法」の開発を実施した．

(2) 大規模データ同化のための革新的4次元変分法の開発

非逐次型データ同化法である4次元変分法は，フォワード計算を実施するための数値モデルから「アジョイントモデル」と呼ばれる方程式系を構成する必要があるため，アンサンブルカルマンフィルタや粒子フィルタのような逐次ベイズフィルタを用いたデータ同化手法と比較して，実装は容易ではない．しかしながら，必要なメモリ容量と計算時間を数値モデルと同等に抑える事ができるため，連続体計算のような自由度の大きい数値モデルを用いるデータ同化を実施する際には，極めて有効である．ただし，従来の4次元変分法に基づくデータ同化では，勾配法による最適化を行うため，逐次型データ同化法では自然に得ることができる推定値の不確実性を評価することが不可能であった．例えば，気象予報においては主に4次元変分法が用いられているが，台風の中心位置予測に関する不確実性を示す「予報円」を描く際には，4次元変分法と逐次型データ同化手法をアドホックに組み合わせて求めているのが実情である．
本センターでは，大自由度系の数値モデルに対する推定値の不確実性の評価を可能にする革新的4次元変分法を開発した．4次元変分法で最適化される評価関数は，一般には最適値の近傍ではモデルパラメータを含む状態変数の2次形式で記述可能であるため，状態変数が従う分布は多変量正規分布によって近似できる．この場合，着目する状態変数の不確実性を示す標準偏差は，その分布の周辺化によって得られるが，それは評価関数の2階導関数からなる行列（ヘッセ行列）の逆行列から導かれる．しかし，ヘッセ行列の逆行列の計算は，自由度の2乗に比例するメモリ容量と，自由度の3乗に比例する演算を必要とするため，通常の数値微分を用いてそれを評価することは現実的ではない．そこで，我々はsecond-orderアジョイント法を導入することにより，ヘッセ行列の逆行列を評価する新しい手法を提案した．second-orderアジョイント法は，任意のベクトルとヘッセ行列の積を計算する手法であり，所与の数値モデルから接線形モデルを，アジョイントモデルからsecond-orderアジョイントモデルを構成する必要があるが，それぞれ数値モデルと同等なメモリ容量と計算時間で実行可能であるため，たとえ自由度が大きい場合でも，十分にヘッセ行列を評価することができる．second-orderアジョイント法と反復法を組み合わせることにより，求めたい推定値の不確実性を示すヘッセ行列の逆行列の対角成分を直接計算することが可能になった．
数値モデルとして相界面移動を記述するフェーズフィールドモデルを題材に，本提案手法を擬似データに基づく「双子実験」により，その性能を評価した．本センターは，科学技術振興機構 戦略的イノベーション創造プログラム 革新的構造材料「マテリアルズインテグレーションシステムの開発」に参画し，構造材料分野で広く用いられているフェーズフィールド法に資するデータ同化技術の開発を行っている．フェーズフィールドモデルは，相の存在確率を連続場の変数として扱うモデルであるが，相界面の分解能をある程度保証する必要があるため，一般に空間格子数は極めて多く，自由度も非常に大きくなる傾向にある．擬似データをノイズが付加された相の存在確率場の時系列とし，提案手法を適用することによって，相界面の移動速度を特徴付けるモデルパラメータの真値が再現されるかどうかを確認した．その際，擬似データを取得する時間の長さ，取得する時間の間隔，あるいは擬似データに含まれるノイズの大きさを変化させた場合のそれぞれについて，パラメータの標準偏差の推定値の変化を定量的に調べた．その結果，データの量やノイズの大きさに依存した推定値の標準偏差が得られ，本提案手法は推定器として非常に良い性質を持つことが明らかになった．また, 場の初期状態推定についても，最初に設定した真の初期場を復元できることがわかった．さらに，本モデルの持つ非線形性に起因する臨界的性質より，推定される初期場に不確定性が励起されることが明らかになった．この不確定性は，ノイズが大きすぎる場合や，推定したい初期時刻から最初の観測時刻までの時間が長すぎる場合に顕著となるが，本モデルの臨界的性質を前もって知っておくことができれば, 除去することが可能である．

(3) 構造物の即時被害予測に資するデータ駆動型モデリングに基づく地震動イメージング技術の開発

本センターは，文部科学省研究開発局の科学技術振興費「都市の脆弱性が引き起こす激甚災害の軽減化プロジェクト」に参画し，観測に基づく都市の地震被害評価技術の開発を実施している．大地震発生時に，都市の構造物における地震応答を数値モデルで即時に評価することで，救助活動の円滑化や二次災害の軽減に貢献することが可能になる．この地震応答の計算の際には，初期条件として構造物直下における入力地震動が必要となる．本センターでは，観測された地震データから，地震観測網よりも圧倒的に密に分布した構造物直下における入力地震動のイメージング手法の開発を行っている．
本年度は，昨年度までに開発した観測データから入力地震動を推定するスパースモデリングに基づくデータ駆動型モデリング手法に加え，波動方程式を物理的な拘束条件として課すことで観測データをより定量的に説明する入力地震動をイメージングする手法を開発した．具体的には，波動方程式に含まれる地震波速度や層厚といった地下構造に関するパラメータや，震源に関するパラメータを未知変数とし，観測データと数値モデルの定量的な適合度を示す事後分布を定義する．一般に，この事後分布は極めて複雑な形をしており，解析的に最適化を行うことは容易ではない．そこで，マルコフ連鎖モンテカルロ（MCMC）法に基づく事後分布からのサンプリングを行い，未知パラメータの最適化を図りながら，同時に入力地震動のイメージングを実行するという方針をとった．MCMC法としては，多峰性を持つ分布からでも効率的にサンプリングすることが可能なレプリカ交換モンテカルロ法を採用した．レプリカ交換モンテカルロ法は，サンプリングを実行すべき分布に加え、それを滑らかにした形状を持つ複数の分布から同時にサンプリングを行う。このとき，適当なタイミングでサンプルの交換を分布間で行うことにより，広範囲かつ効率的な探索を行うことが可能である．地下構造を仮定して生成した擬似データに対して適用する双子実験により，提案手法の性能評価を行った．はじめに，半無限構造を仮定し双子実験を行ったところ，仮定した地下構造を反映するモデルパラメータ値周辺で峰を持つ分布が推定された．一方で，事後分布のみからサンプリングを行うメトロポリス法と呼ばれる通常のMCMC法では，局所的な峰からのサンプルのみが生成される結果となり，本手法が事後分布から大域的かつ効率的にサンプルを得るのに有効であることが示された．次に3層の水平成層構造を仮定した地下構造で同様の実験を行ったところ，事後分布の局所的な峰に捕捉され，真値に近い値がサンプリングされない場合が存在した．これは半無限構造の場合よりも推定すべきパラメータ数が増加したため，推定が困難となっていることに起因しているものと思われる．しかしながら，このようにすべてのパラメータが必ずしも正確に決まらない場合においても，仮定した地震動に近いイメージングを得ることが可能であるという結果が得られた．提案手法の性能評価のために，観測データのみを用いた古典的な空間補間法であるクリギング法との比較を行ったところ，クリギング法の場合は0.30Hz以下の低周波域においても仮定した波動場をきちんと再現できず，物理的な拘束条件の導入は波動場のイメージングにおける必要不可欠な要素であることが示された．我々が開発した地震動イメージング法を首都圏地震観測網MeSO-netの実データに適用したところ，1.0Hz以下の周波数帯域においては，観測データの速度スペクトルと良く一致する推定結果を得ることに成功した．
また即時的な被害予測のためには，地震動イメージング手法の高速化が必須である．本センターでは，地震動イメージングに要する計算時間およびメモリ容量の削減を目的として，各構造物の入力地震動推定に真に有効な観測点を選択する手法の開発を進めている．具体例として，まずはクリギング法と非凸最適化問題を組み合わせることから始めた．クリギング法で用いられる目的関数に，正則化項としてL1−L2ノルムやL1−largest-Kノルムを加味することによって，DC計画問題と呼ばれる最適化問題に帰着させ，恣意性が少ない観測点選択を行うアルゴリズムを開発した．数値実験を用いて本アルゴリズムを検証したところ，上記のどの正則化項を用いても有効な観測点が選択されることが確認された．