学際大規模情報基盤共同利用・共同研究拠点 2021 年度共同研究 最終報告書

jh210054-NAH

大規模地震波シミュレーションによる沈み込み帯の波形トモグラフィー:
2011 年東北地震震源域と南⻄諸島域
岡元 太郎（東京工業大学）
概要 本課題研究の最終目標は波形トモグラフィー手法を用いて 2011 年東北地方太平洋沖地震震源域および南西諸島域の３次元速度構造モデルの推定と改良を行うことである。この目標に向けて今年度は以下の研究を行なった。東北沖パートでは東京工業大学の TSUBAME-3.0 の GPU を利用した大規模計算結果をもとにして、①波形トモグラフィーで用いる感度カーネルの時間領域における特徴を検討した。②波形トモグラフィーの逐次近似１回目に相当する逆問題計算を試行して、地震波速度としては現れにくいパラメータ個別の修正量を得られる可能性などを示すことができた（§5.1）。南西諸島域パートでは、名古屋大学の「不老」スーパーコンピュータ Type-I サブシステムによる大規模計算結果をもとにした計算波形と観測グリーン関数波形との比較により、設定した構造モデルの特徴を検討することができた。特に、① 南西諸島域のほぼ全域に相当する約1000 km にわたる伝播距離に関する理論波形の計算、② その理論波形と観測グリーン関数との比較検討により観測グリーン関数の特徴をある程度まで再現できる場合があることを見出したこと、③ 比較的に近距離（≲ 200 km）かつ長周期（周期 10 〜20 秒）では波形を良好に再現できることの確認、などの成果を得た。
（§5.2）。

1. 共同研究に関する情報

象領域は東北地方太平洋沖地震の震源域と南西諸(1) 共同研究を実施した拠点名島域の２地域とする。改良した３次元構造モデル

東京工業大学 名古屋大学

を用いて短周期の波形再現性を向上させて震源パラメータや破壊過程の解析を高精度化し、短波長(2) 共同研究分野域において地震発生場と地球内部構造との相関関超大規模数値計算系応用分野係を探ることが本課題研究の最終的な目標である。
地震の震源パラメータや断層破壊過程などの(3) 参加研究者の役割分担

研究では、計算地震波波形を用いた解析が欠かせ

東京工業大学・理学院・研究代表者・TSUBAME-ない。このとき計算波形の精度が不十分であると、

岡元太郎3.0 での地震波計算と波形解析結果にバイアスを生じたり分解能に制約を受インバージョンけたりする可能性がある。我々は、プレート沈み

岡山大学・学術研究院副代表者・不老（TypeI）で自然科学学域・竹中博の地震波計算と構造モデ込み型の巨大地震である 2011 年東北地方太平洋

士

ルの検討

沖地震の震源域において、たとえ３次元構造モデ

東京工業大学・学術国際

GPU 計算最適化・大規模情報センター・青木尊之並列化・可視化への助言

岡山大学・大学院自然科

不老（TypeI）での地震波計学研究科・小松正直

算と構造モデルの検討

ルを用いて理論波形を計算しても、周期 10 秒程度以下の短周期地震波波形の再現性が良くないことを見出した（Okamoto et al. 2018）。このことは、波形計算に用いた既存の３次元構造モデルが数2.

10 km 程度の波長域で十分ではないために、短周研究の目的と意義

本研究では大規模地震波シミュレーションの結期（短波長）の計算波形の精度が低下しているこ果を用いた全波形トモグラフィーによって、沈みとを示唆する。この計算波形の精度限界のために、込み帯の 3 次元不均質構造モデルを改良する。対短波長現象である破壊現象の考察や、地震発生場
1

学際大規模情報基盤共同利用・共同研究拠点 2021 年度共同研究 最終報告書の短波長構造の考察においても現状では精度に限

4. 前年度までに得られた研究成果の概要界があるのではないか、と考えられる。
そこで本研究では全波形トモグラフィーによる本課題研究の目標は波形トモグラフィー手法を沈み込み帯の３次元速度構造モデルの改良をテー用いて 2011 年東北地方太平洋沖地震震源域の３次マとする。本研究の意義は、短波長側が改良され元速度構造モデルの改良を行うことである。前年度た３次元構造モデルを提案できること、さらにそには、東京工業大学の TSUBAME3.0 の GPU を利用したの改良したモデルにより震源パラメータや破壊過大規模計算結果をもとにして、物性パラメータ摂動程の解析において短周期（短波長）側の精度を向に対する波形の変化量を定量的に表した「感度カー上させて、地震発生場と地球内部構造との相関関ネル」の試行計算と検討を行った。この感度カーネ係を探る研究の短波長域への拡張を目指すことがルは波形トモグラフィーで用いるものであり、それできる、というものである。

自体に波動伝播の情報が含まれたものになっている。

本年度は前年度に引き続いて、東北地方太平洋計算の結果、浅い地震に対しては深さ 0-12km 程度の沖地震震源域についての陸上観測点データを用い範囲の浅い領域の感度カーネルの振幅が強くなるこた解析を進め、トモグラフィー解析の試行検討をと（つまり構造モデルの浅い部分が観測波形に大き行った。また、南西諸島域を新たな対象として、な影響を及ぼすこと）のような深さ方向の変化や、この地域における地球内部構造モデルの検討を実剛性率（S 波弾性率）の感度カーネルが最も強く P 波施した。

弾性率や密度の効果はやや小さくなることなどのパラメータごとの振る舞いの違い、感度カーネルのパ3.

当拠点公募型研究として実施した意義

ターンが屈曲するなどの空間的な振る舞いなど、基沈み込み帯は、海水層・不規則海底地形・厚い本的な性質を明らかにして波形トモグラフィー実施

堆積物層・沈み込む海洋プレートなどから構成さに向けた検討材料を得ることができた。

れており、不均質性が非常に強い。そのため地震また我々は、GPU とは異なるアーキテクチャのも波動場には 3 次元不均質構造による伝播経路の屈とでの計算手法開発も進めている。前年度には新規曲などの効果が強く表れる。しかも周期数秒程度導入された名古屋大学の「不老」スーパーコンピュの地震波計算では格子間隔を 100－200 m 程度にータ Type-I サブシステム（富岳と同じアーキテクチする必要があることから、数百 km スケールの対ャの富士通 FX1000）において大規模地震波シミュレ

象領域での地震波シミュレーションは必然的に大ーションを実施し、計算結果の確認や性能の考察な規模計算となる。そのため、東京工業大学のど、当システムを利用していくうえでの基本的な検TSUBAME-3.0 および名古屋大学の FX1000 という討を行うことができた。

スーパーコンピュータの利用が必要となる。本研
究では具体的には、東京工業大学の TSUBAME-3.0

5. 今年度の研究成果の詳細

で格子間隔 0.15 km、格子サイズ4896 × 3360 ×

今年度は次の２つのパートで研究を進めた。

800（734km × 504km × 120km）という大規模計［1］2011 年東北地方太平洋沖地震の震源域につ算を実施してトモグラフィーの試行解析を行うこいて、陸上観測点データを用いた解析を進め、構とができた。また名古屋大学の FX1000 では格子造モデル改良のためのトモグラフィー解析の試行

間隔が約 0.2 km、格子サイズ 2201 × 5501 × 1001的な検討を行った。
（§5.1）（ 約 483km × 約 1206km × 200km） と い う 大 規［2］南西諸島域を新たな対象として、大規模地震模計算により 1000 km を超える南西諸島全域に波伝播シミュレーション結果と、地震波干渉法に渡る地震波伝播をシミュレートできて、この地域よって得られる観測グリーン関数との比較によっの地球内部構造モデルの検討を行うことができた。 て、この地域における地球内部構造モデルを検討2

学際大規模情報基盤共同利用・共同研究拠点 2021 年度共同研究 最終報告書

した。
（§5.2）のパラメータを表 1 に示す。

5.1 東北沖での陸上観測点データによる解析
時間領域カーネルの検討

我々は、構造パラメータ摂動量に対する地震波波形の変化量（感度）を示す「感度カーネル」をもとにして構造パラメータ推定の研究を進めている。例えば剛性率𝜇! (𝒙)に関する周波数領域の感度カーネル𝑲" (𝒙、𝜔) は、シンボリックには下式のように物性パラメータ変化
分布𝛿𝜇! (𝒙)に関する体積積分として波形の変化量𝛿𝑢(𝜔)を表す際の積分核として使われる。
𝛿𝑢(𝜔) = 8 𝑲" 9𝒙、𝜔: ;𝛿𝜇! (𝒙)< 𝑑𝒙𝜇! (𝒙)①

もしも波動場の経路上に不均質𝛿𝜇! (𝒙)があれば、感度カーネルと不均質の積に応じた変化を波形に

図 2 解析対象地域。本研究でこれまでに震源パラメータを解析した地震（赤い記号でプロットしたもの）、およびデータを用いた観測点（赤い三角）を示す。

もたらす。つまり感度カーネルは構造モデル中の各部分が波形に対して及ぼす影響を定量的に表す

表 1 TSUBAME-3.0 での計算パラメータ。

ものとなっている（図1）。

𝑢 + 𝛿𝑢感度カーネル 𝑲" (𝒙、𝜔)
震源

波形の変化
不均質 𝛿𝜇! (𝒙)

図 1 感度カーネルと波形の変化量の関係を表す概念図。

全格子サイズ4896×3360×800

副領域サイズ544×480×400

格子間隔150 m

時間ステップ数

35000

時間間隔

0.007143 s

ノード数

32GPU 数126

シミュレーション回数

57 回（=17×3+6）

1 回あたり計算時間

3886 s (平均値)

これまでに、東北沖の3次元構造モデルのもとで周期10秒以下の波動場について周波数領域の感度カ

図 3 には短周期の時間領域感度カーネルの一例

ーネルを計算し、直線経路に関して非対称で屈曲を示す。周波数領域感度カーネルでは、構造モデしたパターンになることなどを見出してきた（4.

ルの浅い部分（深さ 0 km〜数 km）に感度が非常に研究成果の概要）。

強い部分が表れていた。これは時間領域感度カー

今年度はまず、これまでに周波数領域で得られネルでは波形の後半に対応する 140 秒付近以降のた特徴がどのような波動に由来しているのかを把時間帯に現れていることが見て取れる。つまり、握するために、時間領域の感度カーネルを計算し強い感度を持つ浅い領域は、主に波形の主要部分て考察した。感度カーネルは各観測点から広がるが通り過ぎた後に到来する群速度が遅い後続波やアジョイント波動場と震源から広がるフォワード散乱波に対応していると見なされる。このことか波動場の積として求められる。そこで、東工大ら、東北沖の浅部構造をこれらの短周期データでTSUBAME-3.0 を用いた大規模地震波シミュレーシ探るには後続波のモデリングが重要になることがョンによってこれらの波動場を計算した。本研究示される。一方で、後続波は複雑な経路を通っての対象地域や対象地震を図 2 に、また大規模計算きた散乱波も含まれていると考えられ、直逹波よ3

学際大規模情報基盤共同利用・共同研究拠点 2021 年度共同研究 最終報告書

りもモデル化することは難しくなることが考えら

ここでは逆問題解析の試行として、図 1 に示しれる。また、浅い部分の感度が極端に強いと、深た 17 観測点のうち、広帯域地震計が使われているい部分のパラメータ摂動量を求めにくくなる場合

F-net 観測網（7 地点）の観測データを使った。まもあると考えられる。

た、地震は 15 個を用いた。逆問題解析は対象地域

そこで今回はまず、波形の前半部分（振幅の大をブロックに分割し、各ブロック内の平均摂動量きい部分）を解析対象データとして、トモグラフを求めるように設定する（図 4）。

ィー解析のための逆問題計算を行うこととした。

逆問題計算では周波数領域で①式を逆に解いて摂動量を推定する。ここでは 15 地震・7 観測点（各3 成分）の波形データから周期 9 秒−41 秒のフーリエ係数を利用したので、周波数ポイント数は20,868 となる。逆問題計算には平滑化の拘束条件を付けた最小二乗法を用いた。この解析のための

100 s

120 s

感度カーネル生成も TSUBAME-3.0 で実施した。
パラメータがどのような方向に修正されていくかを見るため、図 5 に逐次近似１回目の剛性率（S波弾性率）に関する修正量を示す。まず、浅い層（第 1 層）で修正量が最も大きく、層が深くなる

140 sと修正量が小さくなることがわかる。これは波形

160 sにおいて振幅が大きい部分が主に表面波に対応し
図 3 東北沖の浅い地震（2003/11/01 Mw5.8）と観測点
TYS（上下動成分）に関する、周期 6-8 秒の波動に対する時間領域の感度カーネル（深さ 0-6 km のブロック）。
感度が強い部分が色の濃い領域で示される。ていることが理由になっていると考えられる。また浅い層（第 1 層）に顕著に表れているように、剛性率に関しては全体的に負の修正量（初期値よりも S 波速度を低下させる方向）が卓越している。
図 6 には P 波弾性率、剛性率（S 波弾性率）、密度逆問題解析の試行

上記のように波形の主要部分について、それぞれの第２層（深さ 8.4-16.8km）を対象とするために、フィルターをかけた観測波での修正量を比較した。この図から剛性率（S 波弾形データについて、振幅が最大値の 30 %に低下す性率）で負の修正量が大きい部分では、P 波弾性率るまでの時間帯を解析対象とした。と密度は正の修正量となっており、大まかにはこれらのパラメータが相補的になっている傾向が示唆される。地震波速度は（弾性率／密度）の平方根なので、剛性率と密度の傾向は S 波速度を初期値よりも遅くする方向になっている。P 波については弾性率も密度も同程度の正の修正量なので、P 速度変化は小さいことになる。
このように本研究では、波形を使った解析では密度をパラメータとすることができて、地震波速

図 4 ブロック分割の概念図。本研究での実際のブロック度の変化が小さくてもパラメータには変化があっ数は水平方向 39×28、鉛直方向が 5 個である。ブロックサて、それを検出できる可能性が示された。しかし

イズは 12.6km×12.6km×8.4km とした。

各パラメータの寄与度の精度についてはさらに検討が必要となる。また、現時点では逐次近似の最4学際大規模情報基盤共同利用・共同研究拠点 2021 年度共同研究 最終報告書

初の段階であるため、今後、逐次近似を繰り返し計算波形とを比較して、南西諸島の初期速度構造て最終的な修正量を推定する必要がある。

モデル（小松、2018）の検討を実施した。
「地震波干渉法」とは、ある観測点で観測され第１層第2層た地動波形と、別の観測点で観測された地動波形との相互相関関数を計算することによって、これらの２地点の一方を震源、他方を観測点とするグリーン関数を抽出するものである（Campillo &Paul, 2003 など）。この観測グリーン関数には現実の不均質構造の効果が含まれているので、速度構造モデルの妥当性を検証する上で非常に有用な情報となる。本研究では、観測グリーン関数は観測微動波形をもとにして抽出したもの（邱・ほか、第３層第4層2017）を用いた。

図 5 逐次近似１回目における各層の剛性率（S 波弾性率）の修正量（初期値に対する百分率）。青色は正の修正量（≒密度一定なら S 波速度が初期値よりも速い）、橙色は負の修正量（≒S 波速度が初期値よりも遅い）ことを示す。
図 7 緑色破線ならびに青色破線が計算領域。黒三角形が剛性率
P 波弾性率（S 波弾性率）

F-net 観測点を示す。赤色線とピンク色線は各領域におけ密度る観測点間を結ぶ線。

表 2 沖縄本島周辺・八重山諸島の計算パラメータ。沖縄
本島周辺の計算時間は３回の平均値である。
格子サイズ

沖縄本島周辺

八重山諸島

3601×3401×801

2501×3501×601

格子間隔

0.1 km ×0.1 km ×0.1 km

図 6 逐次近似１回目における第２層（深さ 8.4-16.8 km）

時間ステップ

50001の P 波弾性率、剛性率（S 波弾性率）、密度の各修正量。

時間間隔

5.2

南西諸島の波動場の検討

0.005秒

ノード数

432288

プロセス数

17281152
12スレッド数計算時間

本年度は、地震波干渉法によって抽出した地震95分

63分

観測点の間を伝播する観測グリーン関数と、大規
本年度はまず、2 つの計算領域（図 7 の A と B）

模地震波伝播シミュレーションによって生成した5

学際大規模情報基盤共同利用・共同研究拠点 2021 年度共同研究 最終報告書を設定して、比較的短い距離を伝播する理論波形以上の距離でも、位相に若干のずれが見られる程を計算し、観測グリーン関数と比較した。計算に度であり（図 8 下段の AMM-ZMM）、この例ではは我々が開発した差分法（Nakamura et al. 2012）構造モデルの上部数 10 km のパラメータを改良すと名古屋大学情報基盤センターの FX1000 を使用ることによって波形の一致が改善されることが示した（表 2）。この地震波干渉法で抽出されるのは、唆される。

主に鉛直加振で生成される表面波成分に関するグ

次に、邱・ほか（2017）で抽出したすべての観リーン関数であると考えられるので、シミュレー測グリーン関数との比較を行うために、与那国島

ションの震源も鉛直加振とした。

から鹿児島県を含む全長 1000 km の領域（図 7 のC）を伝播する理論波形を大規模シミュレーションによって生成した。伝播経路がおよそ 1000 kmにわたることから、この計算では地球の曲率の効果を考慮するために我々が開発した準球座標による差分法の計算スキーム（Takenaka et al. 2017）を利用した。計算は前のケースと同様に、名古屋大学情報基盤センターの FX1000 を使用して実行した（表 3）。

表 3 南西諸島全域の計算パラメータ。計算時間は７回の平均値である。

格子サイズ

2201×5501×1001

格子間隔0.002度×0.002度×0.2 km

時間ステップ

60001

時間間隔

0.01秒

ノード数

432

プロセス数

1728

スレッド数

12

計算時間130分

すべての観測点ペアにおける観測グリーン関
図 8 沖縄本島周辺・八重山諸島（表 2）を対象とした比数と理論波形を比較すると（図 9）、領域 A と B の較的に近距離（〜300km まで）の計算結果。周期 5～10 秒

計算結果と同様に長周期かつ観測点間距離（上段）と 10～20 秒（下段）における観測グリーン関数

200 km 以下では、理論波形が観測グリーン関数を波形（黒線）と理論波形（赤線）との比較。横軸は時間、良好に再現することを確認した（図 10 上段）。さ

縦軸は観測点間距離。波形は上下動成分の速度波形である。 らに、それ以上の長距離を伝播する理論波形についても一部のフェーズの再現性が良い観測点ペアが存在した（図 10 下段）。今後、波形再現性に地比較の結果、長周期（10〜20 秒）かつ観測点間距域性があるか否かなどについて、図８を例とする

離 200 km 以下では、理論波形が観測グリーン関ような局所的な計算結果と、南西諸島全域に渡る数を良好に再現することが示された（図 8）。それ大規模な計算結果（図 10）とを相互に比較しなが6

学際大規模情報基盤共同利用・共同研究拠点 2021 年度共同研究 最終報告書

ら引き続き分析する必要がある。の検討も行い、波動場の特徴を把握することができた。② さらにその結果を反映させたデータ処理を施したのちにトモグラフィーにおける逆問題計算の最初の部分を進めることができた。構造モデ
ルパラメータ相互の関係から、地震波速度についての初期モデルからの修正量が小さくても個別のパラメータには一定の修正量がある場合が考えられる。今年度の結果から、本研究の手法ではそのような個別のパラメータの修正量を検出できる可能性が示された。

図 10 伝播距離の異なる 2 つの観測点ペアにおける周期
図 9 南西諸島全域（表 3）を対象とした遠距離（〜1000km） 10～20 秒の観測グリーン関数（黒線）と理論波形（赤線）を含む計算結果。周期 5～10 秒（上段）と周期 10～20 秒の比較例。図の上部に観測点間距離を明記している。波形（下段）の観測グリーン関数波形（黒線）と理論波形（赤は上下動成分の速度波形である。なお、下段で 0〜150 秒線）との比較。横軸は時間、縦軸は観測点間距離。波形は付近の観測波形（黒線）に現れている震動はノイズである

上下動成分の速度波形である。と考えられるので、比較対象にはならない。

6. 今年度の進捗状況と今後の展望

ただし、各パラメータの寄与度の精度（パラメ

東北地方太平洋沖地震震源域の構造推定のパータ間の分解能）についてはトレードオフがありートでは、① 今年度は新たに時間領域のカーネル得るのでさらに検討が必要となる。また、現時点
7

学際大規模情報基盤共同利用・共同研究拠点 2021 年度共同研究 最終報告書

では逐次近似の最初の段階であるため、今後、逐(5) 公開したライブラリなど次近似を繰り返して最終的な修正量を推定する必該当なし。

要がある。これらの点や、海底地震計データの追(6) その他（特許、プレスリリース、著書等）

加を検討することについては次年度（2022 年度）該当なし。の課題とする。
南西諸島域のパートでは観測グリーン関数波

参考文献

1.

形と計算波形の比較により、設定した構造モデル

Campillo, M. and A. Paul, Long-Rangeの特徴を検討することができた。特に、名古屋大

correlations in the diffuse seismic coda,学の「不老」Type-I サブシステム（FX1000）の性

Science, 299, 547-549 (2003).
2.

能を活用した大規模計算によって、① 南西諸島域

小松正直、南西諸島における地震波減衰構造のほぼ全域に相当する約 1000 km にわたる伝播距の推定とそれに基づく地震動シミュレーショ

離に関する理論波形計算ができたこと、② その理ンのための広域地下構造モデルの構築、岡山

論波形と観測グリーン関数との比較検討により観大学博士学位論文（2018）．
3.

測グリーン関数の特徴をある程度まで再現性でき邱 冬晨・小松正直・豊国源知・竹中博士、南る場合があることを観察できたことは大きな成果

西諸島における連続地動記録を用いた雑微動であると考えられる。また、③ 検討した地域では、の相互相関解析、日本地震学会 2017 年秋季

比較的に近距離（≲ 200 km）かつ長周期（周期 10大会、S06-P16（2017）．
4.

〜20 秒）では波形の再現性が良好であることも確

Nakamura T, Takenaka H, Okamoto T,
Kaneda Y, FDM simulation of seismic-wave

認できた。
これらの成果を受けて、次年度（2022 年度）にpropagation for an aftershock of the 2009

は、波形再現性の地域性や、どの程度の改良が必Suruga bay earthquake: effects of ocean-

要であるかなどについて地域的な計算結果と、南bottom topography and seawater layer, Bull.

西諸島全域に渡る大規模な計算結果とを相互に比Seismol. Soc. Am., 102, 2420–2435 (2012)．
5.

較しながら引き続き分析する予定である。

HiroshiTakenaka,

Masanao

Komatsu,

Genti Toyokuni, Takeshi Nakamura, TaroOkamoto,

7. 研究業績一覧

finite-difference computation of seismic wave

（発表予定も含む。投稿中・投稿予定は含まない）(1)

Quasi‑Cartesian

学術論文 （査読あり）

propagation for a three‑dimensional sub-該当なし。

global model, Earth, Planets and Space, 69,67 (2017).

(2) 国際会議プロシーディングス （査読あり）
該当なし。
(3) 国際会議発表 （査読なし）
該当なし。
(4) 国内会議発表 （査読なし）
岡元 太郎・竹中 博士・中村 武史、東北地方太平洋沖地震震源域の波形トモグ
ラフィーに向けた波動場の検討、日本地震学会 2021 年度秋季大会、
S01P-05、2021
年 10 月 15 日。
8