特集「環境リスクと統計解析―データ基盤構築と解析―」［原著論文］

統計数理（2013）
第 61 巻 第 2 号 257–270
c 2013 統計数理研究所


地震音波データ同化システムの開発
双子実験による検証

1

長尾 大道 ・樋口 知之2（受付 2013 年 2 月 20 日；改訂 9 月 9 日；採択 10 月 15 日）

要旨2011 年 3 月 11 日に発生した東北地方太平洋沖地震（東日本大震災）は，東北地方太平洋沿岸部を中心に，特に巨大津波による未曾有の大被害をもたらした．近い将来起こるとされる東海・
東南海・南海地震の際にも巨大津波が発生することが確実視されており，今回の経験を教訓にして被害を最小限に食い止めるために，あらゆる科学分野において様々な視点からの努力がなされている．大地震発生直後に，津波の規模と到来時刻をいち早く予測するための一つの手段として，地震によって励起される地震音波を利用する方法が考えられている．東日本大震災の際には，津波起源と考えられる音波によって電離層を含む地球の大気が大規模に揺り動かされている様が，陸域の高感度気圧計や GNSS 受信機によって捉えられた．音波が津波よりも速く伝搬する特性を利用し，津波が到達する前にその規模および到達時刻を予測する津波早期警戒システムの実現に期待が寄せられている．本稿では，地震音波伝搬の数値シミュレーションと微気圧観測データの融合によって, 震源に関する物理パラメータならびに各時刻における音波
伝搬の状態を推定する地震音波データ同化システムについて紹介し，東南海地震を想定した仮想地震について，震源パラメータを再現する双子実験を実施することにより，本システムを検証する．

キーワード： データ同化，マルコフ連鎖モンテカルロ法，東北地方太平洋沖地震，微気圧，地震音波，津波．
1. はじめに2011 年 3 月 11 日午後 2 時 46 分に発生した，観測史上 4 番目の大きさとなるマグニチュー
ド 9.0 の東北地方太平洋沖地震（東日本大震災）は，東北地方太平洋沿岸部を中心に，特に巨大津波による未曾有の大被害をもたらした．ここ 20 年間の地球観測ネットワークの急速な整備や，数値シミュレーション技術の大幅な向上により，特に東日本大震災のような海洋プレート
境界型地震やそれに伴う津波については，その物理メカニズムの理解が進んだにも関わらず，発生の予測はおろか，津波からの被害を防ぐことができなかった．近い将来，確実に起こるとされる東海・東南海・南海地震に向けて，自然科学をはじめ，工学や社会科学といった様々な科学分野の研究者が，それぞれの視点から知恵を絞り努力を積み重ねている．地震や津波を解
1

統計数理研究所：〒190–8562 東京都立川市緑町 10–3（現 東京大学地震研究所：〒113–0032 東京都文京区弥生 1–1–1）
2

統計数理研究所：〒190–8562 東京都立川市緑町 10–3

258

統計数理

第 61 巻

第 2 号 2013

析する上で，従来の地震波解析とは異なる切り口の一つとして，それらによって大気中に励起される「地震音波」がある．東日本大震災の際には，津波起源の音波によって励起されたと考えられる地震音波によって，地球の大気ならびに高度数十ないし数百 km に存在する電離層が大規模に揺り動かされている様が，日本およびその周辺国の微気圧計および GNSS 受信機によって観測された（Arai et al., 2011; Tsugawa et al., 2011; Nagao and Higuchi, 2012）．大気中における音波の伝搬速度は秒速約 335 m であるのに対し，津波の伝搬速度は海洋プレート境界型地震
発生域の典型的な深度 4,000 m において秒速約 200 m 程度であり，深度 1,000 m 程度の大陸棚に達すると秒速約 100 m に減速し, さらに深度が浅くなる沿岸部に近づくにつれて急速に遅くなる．音波が津波よりずっと速く伝搬するというこの性質を利用し，地震音波を詳しく解析することにより，津波が到来する前にその規模と到達時刻をいち早く知らせることのできる津波
早期警戒システムの実現に期待が寄せられている．もっとも，海洋で発生する全ての地震が津波を引き起こすわけではないが，数値計算との比較から，この時に観測された微気圧変動や電離層変動には，震源域近傍の海面上昇によって励起されたとしか考えられない Lamb 波と思われる変動が含まれており（Matsumura et al., 2011），微気圧計および GNSS による観測でこのLamb 波に着目することにより，津波が発生したかどうかを判断できるのではないかと注目を浴びている．また，国家予算級の十分な資金を確保できるのであれば，海洋研究開発機構の地
震・津波観測監視システム（DONET）のように，沖合に海底地震計や水圧計等の観測センサー
を稠密に展開し，観測データをリアルタイムで回収して解析した方が，津波の早期検知という目的を達成するには優れていることは議論の余地がないものの，微気圧観測は地球物理観測の中でもかなり安価に実施することが可能であり，空間的に稠密な観測を展開することが比較的容易と考えられることも，地震音波研究を実施する利点の一つである．
津波早期警戒システムの実現のためには，地震音波伝搬に関するシミュレーションと観測
データの解析を同時に実施することが必要であることは自明であり，そのための基盤技術であるデータ同化の役割が非常に大きくなるであろうことは想像に難くない．本稿では，第 2 節で地震音波の特徴を実際の観測事例に基づいて概説した後，地震発生後に地震音波伝搬に関する
数値シミュレーションと微気圧観測データ解析を実施し，震源に関するモデルパラメータおよび各時刻における地震音波の状態を推定する「地震音波データ同化システム」について，その原
理および概要を第 3 節において述べる．第 4 節では，東南海地震を想定した仮想的な地震による地震音波について，実際の微気圧観測点において観測されるべき気圧変動を数値シミュレー
ションによって算出し，逆にそれを擬似観測データとして本システムに入力して震源パラメー
タを再現する双子実験により，本システムの性能を検証する．
2. 地震によって励起される地震音波

東日本大震災は，本震後 1 年以上経った現在でも余震活動が非常に活発であり，それに伴う気圧変動が高感度の微気圧計によって観測されている．著者らは，固体地球と大気の物理的相互作用，特に地震音波の伝搬特性を明らかにすることを目的に，2011 年 12 月には電気通信大
学 菅平宇宙電波観測所に，2013 年 3 月には京都大学 潮岬風力実験所に，それぞれ微気圧観測点を設置した．図 1（a）は，本震の約 1 年後となる 2012 年 3 月 14 日に千葉県沖で発生したマ
グニチュード 6 クラスの余震に伴う地震音波を，震央から見て西北西方向約 250 km に位置す
る菅平微気圧観測点の微気圧計によって捉えた観測データである．地球大気の背景気圧および気象変動に伴う短周期の気圧変動を取り除くため，1 ∼ 20 秒周期の帯域通過フィルタを適用した．地震発生から約 1 分後に現れる振幅約 1 パスカル（Pa）の波群は Rayleigh 波と呼ばれる表面波であり，秒速約 1 km で伝搬する地震波によって，直上の大気が振動していることを示す．

地震音波データ同化システムの開発259

図 1. （a）菅平微気圧観測点において観測された 2012 年 3 月 14 日に千葉県沖で発生した東日
本大震災の余震に伴う気圧変動．（b）
ノーマルモード重合による数値シミュレーションで算出した理論波形．

約 15 分後に現れる振幅約 0.5 Pa の波群は, 震源域から大気中を伝搬した地震音波である（長尾他, 2012a）．この地震音波を詳しく解析すれば，この余震のような震源の浅い地震については，
震源の深さを精度良く決定できることが指摘されている（長尾 他, 2012b）．
与えられた固体地球と大気から成る地球構造モデルおよび断層モデルの下で, 観測点における気圧変動を再現するための数値シミュレーション法として, ノーマルモード重合と呼ばれる手法がある（Kobayashi, 2007）．計算手順を簡単に述べると, まずは固体地球の標準モデルである
PREM（Dziewonski and Anderson, 1981）および大気の標準モデルである NRLMSISE-00（Picone
et al., 2002）を組み合わせることにより，地球の構造モデルを構築する．本稿では海洋で発生した地震によって励起された音波が，海洋上を伝搬する距離が長いことを考慮し，PREM に海洋層を含めたままにしておく．地球の構造の空間変化は，水平方向よりも半径方向の方が圧倒的に大きいため, ここでは地球中心からの距離のみに依存する 1 次元モデルを採用し，実際の観測データには含まれていると考えられる水平方向の不均質性の効果をなるべく取り除くた
めに，30 秒以上の長周期変動のみに着目する．この構造モデルについての固有モード（ノーマ
ルモード）およびそれに対応する固有関数を求めることにより，震源-観測点間の Green 関数を構成することができる．次に，地震の震源となる断層モデルを構築する．地震は，地球表層である地殻内において岩盤がある地点（震源）から割れ進むことに端を発し, それによる弾性波が周囲に伝搬していくという現象であり，岩盤の割れ目が断層として残される．したがって，断層の破壊伝搬現象は，断層に沿って複数個の震源を並べ, それらを次々に発震させるというモデルで近似することができる．個別の地震の震源における力源モデルとしては，気象庁や各大学・研究機関が公表している地震のモーメントテンソル（centroid moment tensor; CMT）解⎞⎛Mrr Mrθ Mrφ
⎟
⎜
(2.1)
M = ⎝ Mθr Mθθ Mθφ ⎠Mφr Mφθ Mφφが利用可能である．ただし，(r, θ, φ) は震源を原点とする球座標である．Mαβ（α，β は，r ，θ，
φ のいずれか）は, α 方向に働く力が β 方向に与える影響を表し, モーメントの釣り合い条件から Mαβ = Mβα が成り立つ．すなわち，CMT 解は 2 階 3 次の対称テンソルであり，自由度は6 である．前述した Green 関数と，この CMT 解の畳み込み積分を実行することにより, その

260

統計数理

第 61 巻

第 2 号 2013

応答となる任意の観測点における気圧変動の理論波形が得られる（Lognonné et al., 1998）．
千葉県沖の余震について，気象庁が発表した CMT 解を単一震源とする地震と仮定し，ノー
マルモード重合によって算出された理論波形を図 1（b）に示す．観測波形（図 1（a））と比較すると，表面波および地震音波の走時（地震発生から出現するまでの時間）はほぼ一致しており，観測された波群が本余震に伴う現象であることが確認できる．ただし，両者の波形そのものはそれほど合致しているとは言えないことから，単一震源による観測波形の説明は不十分であり，複数個の震源を並べた断層モデルの導入ならびに各震源に関するパラメータのより高精度の推定が必要であることを示唆している．第 3 節で概説する地震音波データ同化システムにおいては，複雑な断層モデルの取り扱いや，マルコフ連鎖モンテカルロ法（MCMC）による高精度なモ
デルパラメータの推定が可能となっている．
3. 地震音波データ同化システム
3.1 状態空間モデル
数値シミュレーションと観測・計測データを融合するための基盤技術であるデータ同化（正確には，逐次データ同化）では，ある時刻までの「状態」を記憶している状態変数から，シミュレー
ションモデルに基づいて次の時刻における状態を予測し，それと実際に入力された観測データを突き合わせることにより，その時刻における状態を推定するという作業を繰り返す．データ
同化は，モデルパラメータを推定するという従来のオフライン解析方法を包含している上に，シミュレーションモデルに含まれるモデルパラメータが時間的に変化しても構わないように自然に拡張されている．ここで最も重要なことは，状態やモデルパラメータは単一の値ではなく，全て分布関数として表現されていることにある．本節では，所与の観測データに対するシミュ
レーションモデルパラメータの推定というオフライン解析の問題を取り扱うが，実際の地震や微気圧の観測データはオンライン的にどんどんと入ってくるため，近い将来の解析システム構
築に向けた準備として，現段階で逐次データ同化による解析の土台を作っておく．データ同化に関する包括的な解説については，上野 他（2011）を参照されたい．
データ同化では，システムモデルと観測モデルから成る一般状態空間モデルを出発点とする．
(3.1)

xt = ft (xt−1 , vt ),v t ∼ q(·)(3.2)

y t = ht (xt , wt ),w t ∼ r(·)（3.1）式で与えられるシステムモデルは，時刻 t − 1 における状態 xt−1 に対してシミュレーショ
ン ft を実施し，一期先時刻 における状態 xt を予測することを意味している．ここで v t は，ある分布関数 q(·) に従うシステムノイズであり，シミュレーションモデルが実際の現象を表現しきれない「モデリングエラー」を意味する．
（3.2）式で与えられる観測モデルは，状態変数 xt
から観測演算子 ht によって観測可能な量を抽出および算出し，ある分布関数 r(·) に従う観測
ノイズ wt を考慮した上で観測データと比較することを意味している．さて，地震音波データ
同化システムで採用しているノーマルモード重合による数値シミュレーションは線形計算であるため，（3.1）式および（3.2）式は次式で与えられる線形の状態空間モデルに書き直すことがで
きる．
(3.3)

xt = Ft xt−1 + Gv t

(3.4)

y t = Ht xt + wt

ここでは，観測ノイズ w t は平均 0, 共分散行列 diag(σk2 ) の正規分布に従うものとし，σk2 は k
番目の微気圧観測点における観測ノイズの分散を表す．状態ベクトル xt および観測ベクトル

地震音波データ同化システムの開発

図 2.

261

地震音波データ同化システムにおけるデータ同化プロセスのフローチャート.

y t は，具体的に次式で与える．
(3.5)

sim
sim
sim

xt = (δpsim
1,1,t , . . . , δpN,1,t , δp1,2,t , . . . , δpN,K,t )

(3.6)

obs 
y t = (δpobs
1,t , . . . , δpK,t )ただし， は転置を意味し， δpsim
n,k,t (n = 1, . . . , N ; k = 1, . . . , K) は，断層上に配置した n 番目の震源（サブイベント）が励起した地震音波によって，k 番目の微気圧観測点において観測される
べき気圧変動の理論波形を，また δpobs
k,t は，k 番目の微気圧観測点において得られた実際の観測波形を表す．また，シミュレーションモデルとなる行列 Ft には，ノーマルモードや固有関数あるいは CMT 解に関する情報が含まれている．なおここでは，モデルパラメータの推定に最大限に注力するため，システムノイズ v t は考慮しないこととする．サブイベントは数値シ
ミュレーションの都合上，便宜的に設定したものであり，各サブイベントに対する気圧応答である状態ベクトル xt の各成分は，個別に観測することが不可能な潜在変数となっている．実際に観測されるのは，これらを各サブイベントのモーメントマグニチュードによって重み付け
した線形和であるため，観測行列 Ht は次式のようにデザインする．

Ht = m1 B1 (tw ) · · · mN BN (tw )(3.7)
ここに，mn (n = 1, . . . , N ) は各サブイベントのモーメントマグニチュードである．また Bn (tw )は，長さ L の断層上に配置したサブイベントが，断層の破壊伝搬速度 V に従って次々と発震することを表現するために導入した一種の時間遅延演算子であり，任意の時系列 z(t) に対して
以下のように作用するものと定義する．
(3.8)

Bn (tw )z(t) = z t − tw −nL
(N − 1)V

地震音波が伝搬する際には，大気中を吹く背景風の影響を受けるため，それに伴う時間遅延を表現するために tw というパラメータも導入している．
3.2 パラメータ推定
図 2 に，地震音波データ同化システム内で実施されているプロシージャの概要を示す．本シ
ステムによって推定するモデルパラメータは，次のパラメータベクトルにまとめられる．
(3.9)



θ = (C  , m , V, tw , σ 2 )

統計数理

262

第 61 巻

第 2 号 2013

ここに，C = (Cθ , Cφ , Cz ) は震源の位置であり，震源の緯度 Cθ ，経度 Cφ および深さ Cz を成分に2 持つ．また m = (m1 , . . . , mN ) はサブイベントのモーメントマグニチュード，σ 2 = (σ12 , . . . , σK)は微気圧観測点における観測ノイズの分散である．したがって，モデルパラメータの個数はdim θ = N + K + 5 となる．
本システムでは，まず予め用意された過去の地震リストから選択された地震イベントについて，気象庁が発表した震源に関するパラメータ，即ち，震源の位置 C model ，CMT 解 M model ，
断層の長さ Lmodel および破壊伝搬速度 Vmodel を得る．添字 model は，地震波解析によって公表された値であることを示す． MCMC 法によるパラメータ推定の際には，初期値として次のものを設定する．
θ model = (C model , mmodel , Vmodel , 0, 12 )

(3.10)

ただし，12 はすべての成分が 12 であるベクトルを表す．また，断層に沿って同一の CMT 解で表されるサブイベントを等間隔に N 個並べ, それぞれのモーメントマグニチュードは M modelの大きさ Mmodel を等分し, mn,model = Mmodel /N (n = 1, . . . , N ) としておく．
（3.3）∼（3.8）式をまとめ, かつパラメータの事前分布 p(θ) を適当に定めることにより, 理論波形と観測波形の適合度を表す事後分布を，次式のように定義することができる．
(3.11)p(θ|δpobs
1:K,1:T ) ∝ p(θ)
TK×
t=1 k=1

 
2 

N

1
obs
sim
2
exp − δpk,t −mn δpn,k,t−tw − nL
2σk
(N −1)V
2πσk2
n=1

obs
右辺に現れる δpsim
n,k,t は，ここでは時刻 1 から t − 1 までの観測データ δp1:K,1:t−1 が所与の下での一期先予測分布の平均を意味し, 状態ベクトルの定義（（3.5）式）に現れた確率変数とは，本来は表記を変えるべきものである．しかしながら，本稿ではシステムノイズを考慮していないた
め，各時刻において観測データ δpobs
1:K,t に基づくフィルタリング操作を省略することができ, シ
ミュレーションで得られた値と観測データとを直接比較すれば良い．そこで，記号が煩雑になることを防ぐために，ここでは表記を同じものにした．
（3.11）式は，各サブイベントに対応す
る理論波形 δpsim
をサブイベントのモーメントマグニチュード
mn で重みづけし，断層の破壊
n,k,t
obs
伝搬速度 V を考慮した上で重合した波形と，観測波形 δpk,t とを比較することを示している．
4. 双子実験によるシステムの検証

本システムのアルゴリズムを検証するために，仮想的な震源モデルから数値シミュレーショ
ンによって生成した気圧変動を擬似観測データとし，本システムによって震源パラメータを再現する「双子実験」を試みる．図 3 に示す通り，この双子実験においては近い将来起こるとされる東南海地震を想定し，遠州灘沖においてマグニチュード 8 クラスの大地震が起こったと仮定する．ここでは真の震源の深さを Cz,0 = 10 km，断層の長さを L0 = 40 km，断層の破壊伝搬速度を V0 = 2 km/sec，そして地震全体のモーメントマグニチュード M0 = 1021 ニュートン・メー
トル（Nm）とし, 断層に沿った 5 個のサブイベントのモーメントマグニチュードを，それぞれ
m1,0 = 0.15M0 ，m2,0 = 0.25M0 ，m3,0 = 0.35M0 ，m4,0 = 0.20M0 , m5,0 = 0.05M0 であるとする．
添字に 0 を付加した場合は，モデルパラメータの真の値を意味するものとする．実際の東海・東南海・南海地震においては，もっと長い断層が形成されることが想定されているが，あまりに長大な断層を仮定した場合には，サブイベントの個数も大幅に増やす必要が生じ，双子実験
においてそのような大規模計算は検証実験の目的から不要であるため, ここでは比較的小規模
な断層を仮定した．このように断層の長さを短くとると，観測点における気圧変動の振幅が増地震音波データ同化システムの開発

図 3.

263

地震音波データ同化システムを検証するための双子実験において仮想する地震. 長さ40 km の断層上に 5 個のサブイベントを配置し，北東端のサブイベントから順に，断層の破壊伝搬速度に従って発震するものと仮定する. 菅平観測点および潮岬観測点における理論波形を擬似観測データとし，それを再現するようなモデルパラメータを探索する．

大するため，モーメントマグニチュードを小さくした方がより現実的な振幅に近づくが，これは単純にスケーリングだけの問題であり，データ同化による推定には本質的に何ら影響を及ぼさないので、そのままにしておく．震源域から励起された地震音波による大気変動は，断層の中心から 400 km 離れた菅平観測点，ならびに 100 km 離れた潮岬観測点において気圧変動として観測されるものとする．なお，この仮想地震のようないわゆる海洋プレート境界型地震の場
合，震源の水平位置に関しては，地震波解析の段階で大きな推定誤差が生じることは考えにく
いため，本実験においては推定パラメータから除外する．また，大気中の背景風による地震音
波伝搬への影響をノーマルモード重合計算によって再現することは難しいため，遅延パラメー
タ tw についても推定パラメータから除外する．したがって，本実験において推定すべきパラメータは，震源の深さ Cz ，破壊伝搬速度 V ，サブイベントのモーメントマグニチュード m, お
よび観測ノイズの分散 σ 2 となる．
図 4 に，数値シミュレーションによって算出した各サブイベントに対する気圧応答波形 δpsimn,k,tを示す．実際の微気圧観測（図 1（a））において捉えられた特徴，すなわち地震発生後まもなく現れる表面波，およびしばらく経過した後に現れる地震音波が再現されていることが分かる．
これらを，上述の各サブイベントのモーメントマグニチュード mn,0 によって重み付けした上2で，破壊伝搬速度 V0 を考慮して重ね合わせ，平均 0，分散 σk,0の正規分布に従う観測ノイズ
wk,t を加えた

統計数理

264

図 4.

第 61 巻

第 2 号 2013

数値シミュレーションによって理論的に予測された（a）菅平観測点および（b）潮岬観測点における仮想地震によって励起される気圧変動. 1 段目から 5 段目 δpsim
k,n,t (k =1, 2; n = 1, . . ., 5) は各サブイベントに対する理論波形を示し，最下段 δpobs
k,t はそれらを各サブイベントのモーメントマグニチュード mn,0 で重み付けし，断層の破壊伝搬速度 V0 を考慮して重ね合わせた上で，正規分布 N (0, σk2 ) に従う観測ノイズを加えた擬
似観測データ δpobs
k,t である．

(4.1)

δpobs
k,t =

N


mn,0 δpsim
n,k,t + wk,tn=1を，本実験で用いる擬似観測データであるとする（図 4 最下段）．ここでは，海岸付近にある潮岬観測点における観測ノイズは，山間部にある菅平観測点のそれよりも大きいことを考慮し，
2
2
観測ノイズが従う正規分布の分散を σ1,0= 202 Pa2 ，σ2,0= 502 Pa2 であるとする．
図 5 に，MCMC による各モデルパラメータを推定する過程, ならびに得られたサンプリングから推定した各パラメータ ξ ∈ θ についての周辺化事後分布 p̂(ξ|δpobs) MCMC では, 目
1:K,1:T （的分布の規格化定数を推定することが困難であるため，p̂ と記述した）を示し, 表 1 に各パラ
メータの真の値 ξ0 ，初期値 ξmodel , MCMC の際に与えた提案分布（対角共分散行列を持つ正規分布）の分散値, および周辺化事後分布の平均と標準偏差ならびにモードをまとめた．時系列の長さとしては，地震音波に伴う気圧変動が両観測点で顕著となる A: 2,400 秒間の他，津波が海岸に襲来するまでに緊急予測することを想定した B: 300 秒間の 2 つの場合についてデー
タ同化を実施した．なお実際には，気象庁およびそれ以外の研究機関が発表した震源パラメー
タの統計を取り，その結果を基に事前分布 p(θ) を構築することが必要であるが, 本実験においてはもちろんそのようなことは不可能であるため，ここでは事前分布は一様であることを仮定している．MCMC のアルゴリズムとして Metropolis 法（例えば，伊庭 他, 2005）を採用した
ところ, 約 25,000 ステップ目あたりから定常状態に達したため（図 5（a））, 30,005 ステップ目から 80,000 ステップ目まで 5 ステップごとにサンプリングを収集することによって得られた

地震音波データ同化システムの開発

図 5.

地震発生後，A: 2,400 秒間および B: 300 秒間の擬似観測データを用いてデータ同化を実施した場合の（a）事後分布からの MCMC サンプリングの過程，および（b）得られた
サンプリングから推定した各パラメータの周辺化事後分布.

265

統計数理

266
表 1.

第 61 巻

第 2 号 2013

地震発生後，A: 2,400 秒間および B: 300 秒間の擬似観測データを用いてデータ同化を実施した場合の，双子実験における各モデルパラメータの真の値，MCMC による事後分布推定の際の初期値および提案分布の標準偏差，推定された周辺化事後分布の平均と標準偏差ならびにモード．

10,000 個のサンプリングを用いて事後分布を推定した（図 5（b））．本実験においては，MCMCのためのパラメータの初期値として, あえて最悪の値を与えた上で実施したものの, A，B のどちらの場合についても，各パラメータの周辺化事後分布は，真の値の周辺に分布していることが分かる．なお，ノーマルモードの計算アルゴリズムの制約上, 震源の深さ Cz については
1.25 km おきの離散値に対してのみ計算可能であり, また他のパラメータと比較して非常に早く収束したため，得られた周辺化事後分布は Dirac のデルタ関数の形状となっている．気圧変動は，震源の深さに対して非常に敏感であるという特徴を持つため（長尾 他, 2012b），このよ
うな早い収束をもたらしたと考えられる．最後に，MCMC によって得られた事後分布のモードを与えるパラメータ（maximum-a-posteriori 解; MAP 解）から再現された理論波形 δpMAPと,k,t
擬似観測波形 δpobs
k,t を比較したものを，図 6 に示す．理論波形は観測波形の細部に至るまで説明できていることは明らかであり, このことからも地震音波データ同化システムのアルゴリズ
ムが適切に動作していることが分かる．特に注目すべきは，B の場合は地震音波に伴う気圧変動が起こっている途中であるにも関わらず，波形がきちんと再現されている点であり，本シス
テムが早期警戒システムへと発展する可能性を十分に秘めていることを示唆している．
5. 終わりに

本稿では，地震の際に励起される地震音波に伴う気圧変動について，数値シミュレーションと微気圧観測データを融合することにより，震源パラメータを推定することを目的とした地震音波データ同化システムについて紹介した上で，近い将来起こるとされる東南海地震を想定した仮想地震を用いた双子実験を実施することにより，本システムのアルゴリズムを検証した．
双子実験においては，仮定した真の震源パラメータを MCMC によって正しく推定することに成功し，本システムのアルゴリズムの妥当性を示すことができた．提案手法は CMT 解の存在を前提としているため，地震波がほとんど励起しないような津波地震の場合に直接適用することは，現時点では困難が生じる．ただし，第 1 節で述べたように，東北地方太平洋沖地震の際には，津波波源付近の海面上昇によって励起されたと考えられる Lamb 波に伴う微気圧変動が地震音波データ同化システムの開発

図 6.

267

地震発生後，A: 2,400 秒間および B: 300 秒間の擬似観測データを用いてデータ同化を実施した場合の（a）菅平観測点 (k = 1) および（b）潮岬観測点 (k = 2) における擬似観測
MAP
波形 δpobs
k,t （黒線）と MAP 解から再現された理論波形 δpk,t （白線）の比較．

観測されたことが報告されており，また Kobayashi（2007）のコードは Lamb 波についての計算が可能であることからも，津波地震の場合であっても微気圧観測や提案手法の有効性を損なうものではない．
将来的な早期津波警戒システムの構築に向けては，微気圧観測網の整備，大地震による地動や停電等に対して頑健なシステムの構築，さらには数値シミュレーションモデルの精緻化やデータ同化計算の高速化等，数多くのハードルを乗り越える必要がある．本稿の双子実験においては，PC クラスタを用いた 120 並列の計算によって約 1 時間の計算時間を要したが，観測データの回収やその初期処理に要する時間を考慮すれば, 計算時間は少なくとも 10 分の 1以下に抑える必要がある．そのためには，ノーマルモード重合計算の並列化や，レプリカ交換
MCMC 法等の採用による分布関数の推定の効率化を図る必要があるものの，計算資源の単純

統計数理

268

第 61 巻

第 2 号 2013

な大規模化によって十分に達成可能な目標である．地震音波データ同化システムの開発によって，その礎を築くことはできたものの，これより先はあらゆる分野の研究者が知恵を出し合ってプロジェクトを大規模化し，共通の目標に向かって動いていく必要がある．東海・東南海・南海地震においては，地震発生から大津波が襲来するまで，早いところでわずか数分という予測もなされていることからも，津波防災システムには現代における最新鋭の技術を結集させることが求められる．
謝辞菅平および潮岬における微気圧観測点の設置にあたり, 電気通信大学 宇宙・電磁環境研究センターの冨澤一郎准教授, および京都大学防災研究所 流域災害研究センターの林泰一准教授から, 多大なるご協力を賜りました．また，ノーマルモード計算のコードは，宇宙航空研究開発機構の小林直樹助教からご提供頂きました．これらの皆様に，この場を借りて厚く御礼申し上げます．本研究は，情報・システム研究機構 新領域融合研究センターの育成融合プロジェクト
より，ご支援を賜っております．

参 考 文 献Arai, N., Iwakuni, M., Watada, S., Imanishi, Y., Murayama, T. and Nogami, M.（2011）. Atmosphericboundary waves excited by the tsunami generation related to the 2011 great Tohoku-Oki
earthquake, Geophysical Research Letters, 38, doi:10.1029/2011GL049146.
Dziewonski, A. M. and Anderson, D. L.（1981）. Preliminary reference Earth model, Physics of the
Earth and Planetary Interiors, 25, 297–356.
樋口知之 編著，上野玄太，中野慎也，中村和幸，吉田 亮（2011）. 『データ同化入門 — 次世代のシミュ
レーション技術 —』，朝倉書店，東京．
伊庭幸人，種村正美，大森裕浩，和合 肇，佐藤整尚，高橋明彦（2005）. 『計算統計 II マルコフ連鎖モンテカルロ法とその周辺』，岩波書店，東京．
Kobayashi, N.（2007）. A new method to calculate normal modes, Geophysical Journal International,
168, 315–331.
Lognonné, P., Clévédé, E. and Kanamori, H.（1998）. Computation of seismograms and atmosphericoscillations by normal-mode summation for a spherical earth model with realistic atmosphere,
Geophysical Journal International, 135, 388–406.
Matsumura, M., Saito, A., Iyemori, T., Shinagawa, H., Tsugawa, T., Otsuka, Y., Nishioka, M. andChen, C. H.（2011）. Numerical simulations of atmospheric waves excited by the 2011 oﬀ thePaciﬁc coast of Tohoku Earthquake, Earth, Planets and Space, 63（7）, 885–889.
Nagao, H. and Higuchi, T.（2012）. Data assimilation of the Earth’s atmospheric and ionosphericoscillations excited by large earthquakes, The Proceedings of 15th International Conferenceon Information Fusion, 2472–2479.
長尾大道，冨澤一郎，家森俊彦，金尾政紀，徳永旭将，樋口知之（2012a）. 菅平微気圧観測とその初期
解析，月刊地球，34（10），580–584．
長尾大道，小林直樹，深尾良夫，樋口知之，石原 靖，坪井誠司（2012b）. 微気圧変動データ同化による 2008 年岩手・宮城内陸地震の断層パラメータの推定，月刊地球，34（10），566–573．
Picone, J. M., Hedin, A. E., Drob, D. P. and Aikin, A. C.（2002）. NRLMSISE-00 empirical model ofthe atmosphere: Statistical comparisons and scientiﬁc issues, Journal of Geophysical Research,
107, doi:10.1029/2002JA009430.
Tsugawa, T., Saito, A., Otsuka, Y., Nishioka, M., Maruyama, T., Kato, H., Nagatsuma, T. and Mu-

地震音波データ同化システムの開発

269

rata, T.（2011）. Ionospheric disturbances detected by GPS total electron content observation
after the 2011 oﬀ the Paciﬁc coast of Tohoku Earthquake, Earth, Planets and Space, 63,
875–879.

270

Proceedings of the Institute of Statistical Mathematics Vol. 61, No. 2, 257–270(2013)

Data Assimilation System for Seismoacoustic Waves— Verification by a Twin Experiment —Hiromichi Nagao1 and Tomoyuki Higuchi2
1 The Institute of Statistical Mathematics;
Now at Earthquake Research Institute, The University of Tokyo
2 The Institute of Statistical Mathematics

The oﬀ-Tohoku Paciﬁc Ocean Earthquake and especially the accompanying giant
tsunamis, which occurred on March 11, 2011, seriously damaged the Paciﬁc coast in the
Tohoku region. Since a Tokai-Tonankai-Nankai Earthquake and giant tsunamis are expected to occur in the future, researchers in various ﬁelds are making eﬀorts to minimize
damage by putting our experience to use. Large-scale oscillations in the Earth’s atmosphere including the ionosphere probably excited by the tsunamis were observed by highly
sensitive barometers and GNSS receivers. This is expected to realize a tsunami early warning system, which will predict the magnitude of an approaching tsunami and its arrival
time by utilizing the feature that acoustic waves propagate faster than tsunamis. The
present paper introduces a data assimilation system that infers model parameters related
to the hypocenter and states of propagating seismoacoustic waves at every time point byintegrating numerical simulation and observed microbarometer data, and veriﬁes the system by a twin experiment, in which the hypocentral parameters of a synthetic earthquakethat models the forthcoming Tonankai Earthquake are properly estimated.

Key words: Data assimilation, MCMC, the oﬀ-Tohoku Paciﬁc Ocean Earthquake, infrasound, seismoacoustic wave, tsunami.