津波シミュレーションモデルと潮位データの同化におけるモデリングと解析

津波シミュレーションモデルと潮位データの同化におけるモデリ

ングと解析

Modeling and analysis on assimilation of tsunami simulation model

and tide gauge data

中村和幸

1,2∗

樋口知之

1,2

広瀬直毅

3

Kazuyuki Nakamura

1,2

, Tomoyuki Higuchi

1,2

, and Naoki Hirose

3

1

_{統計数理研究所}

1

_{The Institute of Statistical Mathematics}

2

_{科学技術振興機構 戦略的創造研究推進事業}

2

_{Japan Science and Technology Agency CREST}

3

_{九州大学応用力学研究所}

3

_{Research Institute for Applied Mechanics, Kyushu University}

Abstract: We designed a new data assimilation approach for a tsunami simulation model and

tide gauge records, and obtained corrected bottom topography of the Japan Sea. In this data assimilation problem, the number of the variables in the simulation model is very large, whereas observed data set is limited, which results in an ill-posed inversion problem. Therefore, some techniques are required in order to analyze the system, such as introducing appropriate prior informations to the system model. In this paper, we discuss the new approach for constructing the state space model of the tsunami simulation model and the tide gauge records. It is also discussed in view of the graphical model. We show the result of the numerical experiment for validation of the method. We also show the result of the analysis using the tide gauge records of Okushiri Tsunami.

1

序論

津波については，2004 年に発生したインド洋津波以 来，幅広く研究が広がっている分野である．津波は，急 激な海水面高低差が発生することにより引き起こされ る．地震に伴う津波は，地震発生の際に海底地形が急 激に変動し，その上部の海水面に高低差が生じること が原因であり，それが伝播して沿岸に到達することで 被害を生ずることになる． 津波に関する研究については，フィールドワークに よる調査，たとえば最大波高の調査や，潮汐のための 検潮データの解析によるもの，過去の津波の痕跡をも とにした研究などが行われているが，発生頻度の問題 から測定のみをもとにした研究は容易ではない．その ため，フィールドワークや計測データの解析とは別に， 物理過程から導かれる数値シミュレーションコードに より実際の津波を模擬し，その結果と実測との比較に ∗_{連絡先：統計数理研究所}       〒106-8569 東京都港区南麻布 4-6-7       E-mail: [email protected] よってダイナミクスを議論する形式での研究が，海岸 工学や地球物理学の分野においてよく行われている． これらのシミュレーションにおいて，境界条件の一 つである海底地形は，既知のものとして扱われている． しかし，海底地形データセットには誤差があることが わかっている [1]．日本海周辺についても，米海軍海洋 学局，米地球物理学データセンター，日本水路協会，韓 国成均館大学などのデータを比較すると，大きな差異 があることが知られている．これら異なるデータをも とに津波のシミュレーションを行なうと，結果も無視で きないほどに異なる．海底地形の探索については，音波 のみが利用可能であるなどの物理的制約が影響し，観 測量を増やして精緻にするという方向では限界がある ため，より正確な津波シミュレーションを行うには，海 底地形を修正し誤差の評価をする何らかの方法が必要 となる．そもそも一般的に，シミュレーションモデル は，モデル化されない部分，境界条件・初期条件，未 知あるいは人工的に導入したパラメータ等の不確実性 を含むものである．この不確実性を減少させるために 観測情報を導入するようになったのが，データ同化で

人工知能学会研究会資料

SIG-DMSM-A701-02 (7/25)

ある． データ同化は，主に気象・海洋学の分野で発展して きたものである．通常の統計的推定の問題等と大きく 異なる点として，モデルの規模が挙げられる．地球物 理学のデータ同化においては，シミュレーションモデ ルの変数の数が最大で数百万のオーダーに及ぶものを 対象とする．また，地球物理学・シミュレーション科 学両面における蓄積の結果，シミュレーションモデル はある程度精度の高いものができている．しかしなが ら，シミュレーションモデル単体では現実の情報を反 映していないため，実際に起こっていることの再現に は限界がある．その一方で，観測はさまざまな物理的・ 社会的制約により，得られる情報が制限される．デー タ同化は，これらの欠点を相互に補う方法であり，予 報精度向上のためにシミュレーションモデルの初期値 を整備することからはじまり，推定・同定による科学 的発見，さらにはシミュレーションモデルの向上のた めの解析などを目的として発展してきている． 同化の手法としては，逐次型と非逐次型の手法があ る．本発表では，逐次型の手法について扱う．非逐次 型の例として，4 次元変分法 [2] が挙げられる．逐次型 の手法としては，カルマンフィルタ，拡張カルマンフィ ルタ，アンサンブルカルマンフィルタ [3, 4]，粒子フィ ルタ [5, 6] などが用いられている．本研究では同化計 算コスト，特にフィルタリングステップでのメモリと 計算のコストの理由により，粒子フィルタを用いた． 津波シミュレーションの分野において従来から用い られている，観測情報とシミュレーションモデルの統 合の手法として，グリーン関数法による，初期分布の インバージョン解析 [7, 8] が挙げられる．この手法で はまず，津波波源周辺での単位波源を用意してシミュ レーションモデルにより伝播させ，各観測点での観測 系列への応答を求める．このような単位波源を場所を 変えて複数用意し，その最小二乗フィッティングを行 うことにより，津波波源の空間的な形状を得る手法で ある．この手法も強力ではあるものの，線形フィッティ ングの方法であることから，非線形応答に対して，必 ずしも十分な手法とはいえない．システム理論の観点 から見ると，初期分布構成による方法は，応答関数を 固定した状態での入力の再現を行っていることになる． この場合，応答関数そのものに誤差（境界条件やモデ ル化されない力学などに相当）が含まれる場合，その 誤差は初期分布に含まれる誤差として表現されること となる．その結果，ある程度の予測性能向上が見込ま れる一方で，過渡応答全体に対するフィッティングには 限界があることから，応答系列内で非線形性が顕著に なり始める前での適切なタイミングでの打ち切りが必 要となる．一方，今回使用する海底地形修正による方 法の場合，応答関数を条件に合わせて変更することに 対応する．この場合，内部パラメータにより応答が変 わるようなシステムのモデル同定，あるいはエラーを 最小にするようなパラメータ推定の問題となる．その 結果，非線形モデルに対してより柔軟に適応可能であ ると想定される．以上のことから，初期分布インバー ジョンと海底地形補正は並立する方法であり，互いの 欠点を補うことが可能であるが，両者の間には不定と なる誤差項が存在し得ることから，同時に行うのは容 易ではない．以下では後述の理由により，海底地形補 正に絞って議論することにする． 本研究においては，津波のシミュレーションモデル として，海洋学において用いられる浅水波モデルを用 いる．このシミュレーションモデルに含まれる変数は， 海水面高，海水の流速ベクトル，水深などである．今 回，実データの解析に用いたシミュレーションモデル の変数の数は，9× 104程度である．一方観測として用 いるデータは，設置点 18 点の海水面高を各時点におい て計測した沿岸潮位計データである．以下では，これ らのシステムに対して，非線形非ガウス型状態空間モ デルでの定式化を行い，粒子フィルタを適用すること で，逐次型データ同化による解析を行なう．同化の際 に注目する変数は，各地点における水深である．理由 は，前述したように水深データの計測には誤差があり， 同化することによりシミュレーションモデルおよび観 測データに適合した水深データを得られるためである． このことにより，海底地形に対する新たな知見の獲得， 海洋シミュレーションの精度向上のための適切な境界 条件設定，津波の正確な予測などが期待できる．以上 の理由から，本研究においては，海底地形補正を行う ことを目指す．

2

逐次データ同化の枠組みとモデリ

ング

以下では，シミュレーションモデルと観測データに ついて，非線形非ガウス型状態空間モデルの構成法と， その推定について記述する．

2.1

非線形非ガウス型状態空間モデルの構成

2.1.1 システムモデルの構成 はじめに，一般のシミュレーションモデルからどの ように非線形状態空間モデルのシステムモデルとなる かを議論する．以下では，シミュレーションモデル内 で時間発展をするシステムを考える．まず，状態変数 xtとしてシミュレーションモデル内のタイムステップ t での全変数を設定する．すると，数値シミュレーショ ンモデルは xt= ˜_{f (xt−1}) (1)

というように形式的に書くことができる．すなわち，シ ミュレーションモデル内のあるタイムステップでの全 変数を決定すれば，次のタイムステップにおけるシミュ レーションモデル内の変数は一意に決まることになる． この時点では，確率変数に対応させる誤差項は存在せ ず，単なる有限差分方程式系である．ここで，モデル 化されない誤差や不確実性をもつ境界条件等を示す項 として，vtを導入し，上式と組み合わせて実際のシス テムの不確実性を近似することを考える．すると，こ れに従うシステムは xt= f (xt−1, vt) (2) と修正される．ここで f は不確実性の近似の形式によっ て決まり，もし vtが加法的であるとしたならば f (xt−1, vt) = ˜_{f (xt−1) + vt} (3) となる． 今回用いるシミュレーションモデルは，浅水波方程 式 [9] を時間的に離散化，空間的に格子化したものであ る．ここでは構造格子を用いており，以下ではその格 子点数を M とする．津波シミュレーションモデルに含 まれる主要変数は，各格子点 m 毎に海水面高 ηm，流 速ベクトルの緯度方向／経度方向の各成分 umならび に vm，各点の水深 dmとなる． 海底地形補正を行う際は，以下のようにモデルを設 定する．ここで決める必要があるのは，システムノイ ズ vt と初期状態 x0である．初期状態については，以 下の通りとした．時刻 0 での流速ベクトルについては， 全成分 0 で固定とし，海水面高 η·は，地震波形のイン バージョンから決まる諸パラメータをもとに，海底地 形の変動の形と同様の形が地上に形成されるという仮 定 [10] に従って，海水面高を設定する．一方，注目し ている海底地形については，津波伝播中の海底地形は 不変とみなすのが自然であることより，それを表現す る水深成分について，x0すなわち時刻 0 での状態ベク トルの水深成分についてのみ誤差分布を仮定し，残り の時刻においては，変動しないものとした．一方モデ ル化されない誤差項に相当する vtであるが，ここでは 観測が少ないことを考慮に入れて，モデル設定の段階 では常に 0 であるとした．ただし，後述する粒子フィ ルタの問題により，後に水深に対応する成分のみにシ ステムノイズを微小に付加する． 2.1.2 観測モデルの構成 今回用いる観測データは，沿岸に設置されている潮 位計データである．データの時間間隔は，30 秒から数 分程度である．潮位計は，潮汐変動を記録することを 目的として津波を計測することを意図していないもの 1 |1 − − t t x x_{t| −t 1} t y 䉲䊚䊠䊧䊷䉲䊢䊮 t t x| xt+|1t 䉲䊚䊠䊧䊷䉲䊢䊮 1 − t y 1 + t y ᷹ⷰ䊂䊷䉺 ᷹ⷰ䊂䊷䉺 䉲䊚䊠䊧䊷䉲䊢䊮䊝䊂䊦 ਄䈱䉺䉟䊛䉴䊁䉾䊒 t _t t+1 ୃᱜ ୃᱜ 䋨䊐䉞䊦䉺䊥䊮䉫䋩 䋨䊐䉞䊦䉺䊥䊮䉫䋩 1  図 1: 逐次データ同化の手続き の，周期が数分以上になる津波に対しては十分な時間 分解能があり，津波による海水面変動をとらえること が可能である．従来は他に有効な計測が少なかったこ となどから，潮位計による計測が研究に使われること が多い．潮位計で観測される系列には，天文潮，気象 潮や他の要因などによる変動が記録される．これらは 通常，津波よりも周期が長く，変動が緩やかであるた め，適切な時系列解析手法を適用することにより，こ れらを除くことが可能である．従って，何らかの方法で これらの成分を除いた残差系列を観測系列として用い ることになる．今回は，各観測点毎に１次元ランダム ウォークモデルをあてはめ，その残差系列を使用した． 以上により得られた１次元観測時系列は，津波によ る海水面変動を表現していることになる．また，津波 未到達時の系列のぶれが，観測ノイズであるとして観 測ノイズの分散を推定する．さらに，これらの１次元 観測時系列が K 点得られたときに，これを K 次元観 測 ytにまとめることにする．以上の手続きをまとめる と，観測モデルは yt= Hxt+ wt (4) と定式化されることになる．ただし H は観測点と状態ベ クトルを対応づける 0-1 行列，wtは N (0, diag(σ12, · · · , σK2)) に従うノイズである．以上から，データ同化における シミュレーションモデルならびに観測データについて， 非線形非ガウス状態空間モデルで表現可能であること がわかる．

2.2

逐次データ同化

逐次データ同化とは，データを得るたびに数値シミュ レーション内の変数を修正していく手法である（図１）． これは非線形非ガウス状態空間モデルの状態推定の問 題に他ならない．この推定を与えるアルゴリズムとし て，地球物理学の分野では拡張カルマンフィルタやア ンサンブルカルマンフィルタなどが用いられている．し かしながら，これらのアルゴリズムでは非線形性の扱 いに問題があり [11] ，また Rtが対角の場合でも密な

x

0

x

1

x

2

x

3

x

T

y

1

y

2

y

3

y

T

θ

図 2: もとの物理モデルに対応するグラフィカルモデル

x

0

x

1

x

2

x

3

x

T

y

1

y

2

y

3

y

T

θ

0

θ

1

θ

2

θ

3

θ

T 図 3: 自己組織化モデルに対応するグラフィカルモデル 逆行列の計算が必要であることから，計算負荷も小さ くない．本研究では，これらの問題点が少ない粒子フィ ルタを用いることにした．

2.3

自己組織化状態空間モデルとグラフィカ

ルモデル

2.1節での設定は，以下で説明する通り，図２および 図３のグラフィカルモデルに対応させて議論できる．も ともと，物理モデルにもとづいて素直にモデリングを した場合，対応するグラフィカルモデルは図２となる． このとき，水深 d は状態ベクトル x·ではなく，全体の システムを規定するパラメータ θ に含まれると見なす 方が自然である．しかしながら，これは密なグラフィ カルモデルでの推定問題であり，システムモデルによっ て xt|t−1を得る計算において，大規模なシミュレーショ ンを利用するという制約のために，時間発展方向のみ の計算が可能であること，さらにその計算に時間がか かるといった問題により，難問となる． そこで，図３のようにパラメータ θ にも時間発展さ せることとし，初期分布を導入し，かつ状態ベクトル と組にして新たな状態ベクトルとみなす，すなわち水 深を状態ベクトルに加えるようにする．このことで，自 己組織的にパラメータを構成していくことが可能とな 図 4: 人工地形における双子実験の結果 る [12]．図３において，破線で示した枝は本来誤差項 の入らない枝，すなわち P (θt= θt−1) = 1となる枝で ある．しかし，このように設定して粒子フィルタによ る推定を行った場合には，退化の影響によりサンプル に偏りが生じるおそれがある．そのため，テクニカル な手法ではあるが，枝に擾乱を加えることにより対応 する．このことにより，退化の影響によるサンプルの 偏りを抑え込むことが可能となる． 純粋にモデリングの観点から，この擾乱項を加える 場合もありえる．途中のパラメータが，シミュレーショ ン内部に存在する何らかのモードに対応して変化する と考えるのが適当である場合である．この場合，時刻 全体において状態をコントロールするパラメータを導 入するかわりに，局所的に定常なパラメータを仮定し， パラメータの切り替わりとモードの切り替わりが対応 するように設計することになり，システムノイズの分 散パラメータは大きめに設定することになる．しかし ながら，今回の目的の範囲外であるため，このアプロー チを取っていない．

3

双子実験

実問題に対してデータ同化を適用した場合には，当 然真の推定になっているとは限らない上，そもそも真 の状態はわからない．そのため，手法の妥当性を検証 する手段の一つとして，双子実験とよばれる数値実験 が用いられる．双子実験は一般に，以下の手続きに従っ

て行われる．まず，状態やパラメータ等について記録 しながらシミュレーションを行い，その結果と観測モ デルから観測系列を作成する．次に，この観測系列を， 先ほどのシミュレーションとは別の設定から始めたも のとデータ同化を行う．この結果得られたシミュレー ション内の状態やパラメータが，もとの興味のある状 態やパラメータに近づいているかどうかによって，評 価を行う．非常に単純な比較の方法であるものの，シ ミュレーションモデルの計算が容易ではないことから， この方法を用いる．すなわち，事前分布を固定した状 態での各時点での分布と尤度の評価を前向きにする程 度が限界であり，なおかつその際に動かせるパラメー タは限られるという制約から，このような方法で評価 を行うことになる． 津波シミュレーションモデルの場合の双子実験の流 れは以下の通りである．まずある海底地形を仮定し，そ れにあわせてシミュレーションを行う．このとき，あ らかじめ設定した観測点の近傍における海水面変動を 記録し，潮位計データとする．次に，別の海底地形を もとにしたシミュレーションモデルと潮位計データを もとにデータ同化を行い，最初のシミュレーションに 用いられた地形が復元されるかどうかを調べる．この 手順の意味合いは単純で，事前分布にバイアスを仮定 し，事後分布が真値周辺に絞られるかどうかを確認し ていることになっている． 図４は，人工地形に対して上記の双子実験を適用し た結果である．各絵の左側に示した図が津波の伝播の 様子，右側は左側の白線に沿った断面について描いた 図である．まず，中央部には海山がないという設定のも とでシミュレーションを走らせて観測系列を作成する． このときの観測点は，白点で示した４点である．次に 中央部に海山があるという設定をし，その周りで分布 させたものを初期分布としたもとで，逐次データ同化 を行っている．図のプロットは，この同化の様子を示 したものである．右側の緑線が真の地形であり，青線 が推定地形である．また，橙線が粒子フィルタの粒子 のうちで最も浅いもの，紫が最も深いものである．津 波伝播に従って，分布の範囲が狭まると同時に，推定 が真値に近づいていることが確認できる．

4

北海道南西沖地震津波による解析

これまでの枠組に従い，実データである潮位計デー タを用いて，日本海底地形の修正を行った．今回使用 したデータは，北海道南西沖地震津波の潮位計データ である．観測点は日本海沿岸ならびに韓国沿岸の計 18 点である．またシミュレーションモデルは，格子点が ５分刻みのもので，経度方向に 192 点，緯度方向に 240 点設定している． 10000 8000 6000 4000 2000 0 0 20 40 60 80 Depth(m) SD/Depth(%) 図 5: ４地形データセットの平均水深（左）とデータ セット間の誤差割合（右） 0 1000 2000 3000 4000 5000 6000 40 45 50 55 60 average shallowest deepest estimation 図 6: 推定結果．大和堆の東経 134 度での南北方向断 面を表示したもので，左側が北側．紫線が平均，緑線 が推定，赤線と青線がそれぞれ最浅と最深の粒子が表 す水深である． 日本海海域において，既述の４データを比較し誤差 の大きい所を調べた結果，沿岸部に加え，日本海中央 部の大和堆と呼ばれる海山においてデータセット間の 差が大きいことが確認された（図５）．そのため，この 部分についてデータ同化により修正を行うこととした． 具体的には，当該領域が４地形データセットの線形和， その他の領域が４データセットの平均であると仮定し， そのもとで初期分布を構成し，さらに微小なシステム ノイズを，線形和の係数に加えながら推定を行った． 推定の結果を表したのが図６である．これは，全て の時系列データを用いた後の d についてプロットした ものである．この図より，以下の２点が確認できる．１ 点目は，もとの４地形の平均と比較して，与えられた 推定が浅くなっていることである．２点目は，大和堆 南斜面について，その他の部分と異なり，与えられた 推定が深くなっていることである．より高低差が大き

く，急峻な地形が与えられているということができる． この結果について，他のモデルでの検証や調査を行っ ている段階である．

5

結論

津波シミュレーションモデルと沿岸潮位計のデータ 同化を行う枠組を示した．特に，海底地形補正の枠組 と，そのグラフィカルモデルでの位置づけを確認した． 数値実験により手法の妥当性を確認するとともに，日 本海底地形の補正の結果，従来のデータセットよりも 浅いという結果を得た．以上の結果についての，他の モデルでの検証などにより補強していくとともに，よ り精密なシミュレーションモデルの導入と予測実験等 が，今後の課題として挙げられる．

謝辞

本研究においては，海上保安庁，国土地理院ならび に韓国成均館大学の Choi 教授の潮位計データを使用さ せていただきました．ここに感謝の意を表します．

参考文献

[1] Hirose, N.: Least-squares estimation of bottom to-pography using horizontal velocity measurements in the Tsushima/Korea straits. Journal of

Oceanogra-phy, Vol. 61, pp. 789–794 (2005).

[2] Courtier, P., Thepaut, T., and Hollingsworth, A.: A strategy for operational implementation of 4DVAR, using an incremental approach. Quarterly Journal of

the Royal Meteorological Society, Vol. 120, pp. 1367–

1387 (1994).

[3] Evensen, G.: Sequential data assimilation with a non-linear quasi-geostrophic model using Monte Carlo methods to forecast error statistics. Journal

of Geophysical Research, Vol. 99, No. C5, pp. 10143–

10162 (1994).

[4] Evensen, G.: The ensemble Kalman filter: Theoreti-cal formulation and practiTheoreti-cal implementation. Ocean

Dynamics, Vol. 53, pp. 343–367 (2003).

[5] Kitagawa, G.: Monte Carlo filter and smoother for non-Gaussian nonlinear state space model. Journal

of Computational and Graphical Statistics, Vol. 5,

No. 1, pp. 1–25 (1996).

[6] Gordon, N. J., Salmond, D. J., and Smith, A. F. M.: Novel approach to nonlinear/non-Gaussian Bayesian state estimation. IEE Proceedings-F, Vol. 140, No. 2, pp. 107–113 (1993).

[7] Satake, K.: Inversion of tsunami waveforms for the estimation of heterogeneous fault motion of large earthquakes: The 1968 Tokachi-oki and the 1983 Japan Sea earthquakes. Journal of Geophysical

Re-search, Vol. 94, pp. 5627–5636 (1989).

[8] Tanioka, Y., and Satake, K.: Coseismic slip distri-bution of the 1946 Nankai earthquake and aseismic slips caused by the earthquake. Earth Planets Space, Vol. 53, pp. 235–241 (2001).

[9] Choi, B. H., and Hong, S. J.: Simulation of prog-nostic tunamis on the Korean coast. Geophysical

Re-search Letters, Vol. 28, pp. 2013–2016 (2001).

[10] Manshinha, L., and Smylie, D.: The displacement field of inclined faults. Bulletin of the Seismological

Society of America, Vol. 61, No. 5, pp. 1433–1440

(1971).

[11] Nakamura, K., Ueno, G., and Higuchi, T.: Data assimilation : Concept and algorithm. Proceedings

of the Institute of Statistical Mathematics, Vol. 53,

No. 2, pp. 201–219 (2005).

[12] Kitagawa, G.: Self-organizing state space model.

Journal of the Americal Statistical Association,