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

jh200059-NAH

大規模地震波シミュレーションを用いた 2011 年東北地方太平洋沖地震の震源域構造モデルの波形トモグラフィー
岡元 太郎（東京工業大学）
概要 本課題研究の目標は波形トモグラフィー手法を用いて 2011 年東北地方太平洋沖地震震源域の３次元速度構造モデルの推定と改良を行うことである。この目標に向けて今年度は、東京工業大学の TSUBAME3.0 の GPU を利用した大規模計算結果をもとにして、波形トモグラフィーで用いる「感度カーネル」の試行計算と検討を行った（§5.1）。感度カーネルは物性パラメータ摂動に対する波形の変化量を定量的に表したものである。
計算の結果、感度カーネルの深さ方向の変化やパラメータごとの振る舞いの違い、周波数依存性などの基本的な性質を明らかにして、波形トモグラフィー実施に向けた検討材料を得ることができた。また我々は、GPU とは異なるアーキテクチャのもとでの計算手法開発も進めている。本課題研究では、今年度新規に導入された名古屋大学の「不老」
スーパーコンピュータ Type-I サブシステムにおいて大規模地震波シミュレーションを実施し、計算結果の確認や性能の考察など、今後本システムを利用していくうえでの基本的な検討を行うことができた（§5.2）。

1. 共同研究に関する情報
2. 研究の目的と意義(1) 共同研究を実施している拠点名
 東京工業大学 学術国際情報センター

地震の震源物理や強震動生成メカニズムを探 名古屋大学情報 基盤センター

るためには、計算地震波波形を用いた解析が必須となっている。このとき計算地震波波形の精度が不十分であると、解析結果にバイアスを生じたり、(2) 共同研究分野■

超大規模数値計算系応用分野

分解能に制約を受けたりする可能性がある。我々

超大規模データ処理系応用分野

は、プレート沈み込み型の巨大地震である 2011



超大容量ネットワーク技術分野

年東北地方太平洋沖地震の震源域において、たと

超大規模情報システム関連研究分野

え３次元構造モデルを用いて理論波形を計算しても、周期 10 秒程度以下の観測地震波波形の再(3) 参加研究者の役割分担現性が良くないことを見出した。これは波形計算

東京工業大学・理学院・岡元太郎に用いた３次元構造モデルが十分ではないこと

岡山大学・大学院自然科学研究科・竹中博士
東京工業大学・学術国際情報センター・青木尊之
岡山大学・大学院自然科学研究科・小松正直

研究代表者・TSUBAME-3.0 での地震波シミュレーションと波形インバージョン
副代表者・名大新システ
ムでの地震波プログラムの開発・チューニング・検証計算
GPU 計算の最適化、大規模並列化・計算結果の可視化への助言
名大新システムでの地震波プログラムの開発・チューニング・検証計算を示唆する。そこで本研究では全波形トモグラフ
ィーによる東北地方太平洋沖地震震源域の３次元速度構造モデルの推定・改良を行う。本研究の特徴は、①初期震源パラメータの不確定さを最小化 す る た め に 独 自 の FAMT （ First-motionAugmented Moment Tensor）解析を用いること、②構造パラメータ摂動に対する波形の感度カーネ
ルを高周波近似の波線理論等によらず波動論に基づいて計算すること、である。そして、最終的にはこの改良した３次元構造モデルを用いて短1

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

周期の波形再現性を向上させて震源パラメータ度を向上させて、地震発生場と地球内部構造とや破壊過程の解析を高精度化し、地震発生場と地の相関関係を探る研究の短波長域への拡張を目

球内部構造との相関関係を探る研究を短波長域指すことができる、というものである。に拡張することを目指す。

4. 前年度までに得られた研究成果の概要
3. 当拠点公募型研究として実施した意義

新規課題のため該当せず。

沈み込み帯に位置する日本列島では、2011 年東5. 今年度の研究成果の詳細

北地方太平洋沖地震のようなプレート境界型の巨大地震がしばしば発生する。このような巨大地本課題研究の目標は波形トモグラフィー手法

震の断層破壊過程や地震の震源分布・発震機構解を用いて 2011 年東北地方太平洋沖地震震源域のは、地震発生場の応力状態・物性や大規模構造と３次元速度構造モデルの推定と改良を行うこと関係すると考えられることから、地震の各種の解である。この目標に向けて今年度は、東京工業大

析結果と不均質構造との相関関係を考察するこ学の TSUBAME3.0 スーパーコンピュータの GPUとが地震発生場を理解する上で重要な手がかり（NVIDIA P100）を活用した大規模計算により、となる。

波形トモグラフィーで用いる「感度カーネル」の試行計算と検討を行った（§5.1）。

これらの研究では地震波波形が重要なデータとなる。しかしプレート境界型地震が発生する

また我々は、GPU とは異なるアーキテクチャの沈み込み帯は前述のように強い３次元不均質性もとでの計算手法開発も進めている。本課題研究を持っており、地震波は伝播経路が歪められるでは、今年度に新規導入された名古屋大学のスーなどの大きな影響を受ける。この影響を考慮しパーコンピュータ「不老」を利用することとした。

ない場合には地震の解析結果にバイアスが発生「不老」は新たなプロセッサが用いられたシステする可能性がある。我々は、この計算波形によるムである。このシステムを利用するのは初めてで観測波形の再現性の問題に取り組んできた。そあることから、今年度は地震波シミュレーションの結果、既存の 3 次元構造モデルと最適震源モプログラムの実行結果の確認や性能の考察などデルを用いて計算波形を生成しても波形全体にの基本的な検討を実施する（§5.2）。これは今後ついては周期 10 秒程度以下の短周期帯での再現

このシステムで大規模計算を用いた研究を進め性は低いことを見出した。これは現状の 3 次元ていくための準備段階と位置付けられる。

構造モデルの精度が数 10km 程度の波長域で充5.1 波形トモグラフィーのための感度カーネルの分ではないことを示唆する。この計算波形の精度限界のために、短波長現象である破壊現象の考察考察や、地震発生場の短波長構造の考察におい感度カーネルとはモデルパラメータの摂動にても現状では精度に限界があるのではないか、関する波形の変化量を示すものであり、下記のと考えられる。

ような体積積分の積分核𝑲𝑲𝑀𝑀 (𝑿𝑿, 𝒙𝒙, 𝜔𝜔)として表現される。

そこで、本研究では波形トモグラフィーによって 3 次元構造モデルを改良し、周期数秒程度𝛿𝛿𝛿𝛿 (𝒙𝒙)� 𝑑𝑑𝒙𝒙𝛿𝛿𝛿𝛿(𝑿𝑿, 𝜔𝜔) = � 𝑲𝑲𝑀𝑀 (𝑿𝑿, 𝒙𝒙, 𝜔𝜔) �𝑀𝑀 (𝒙𝒙)の短周期帯での波形再現性を向上させることをここで𝛿𝛿𝛿𝛿(𝑿𝑿, 𝜔𝜔)は位置𝑿𝑿の観測点における波目的とする。本研究の意義は、短波長側が改良さ形の周波数領域での変化量、𝑀𝑀 (𝒙𝒙)は物性パラメれた 3 次元構造モデルを提案できること、さらにその改良したモデルにより震源パラメータやータ、𝛿𝛿𝛿𝛿 (𝒙𝒙)はその摂動量をそれぞれ表す。この

感度カーネルが波形トモグラフィーにおける分破壊過程の解析において短周期（短波長）側の精2

学際大規模情報基盤共同利用・共同研究拠点 2020 年度共同研究 最終報告書解能（解像度）に関係する。そこで、感度カーネ
ルの試行計算によって、その性質を詳細に検討した。
5.1.1 対象地震と観測点
本課題研究では 2011 年東北地方太平洋沖地震の震源域を対象とする。この対象地域に発生した 19 個のマグニチュード 6 前後の中規模地震を解析対象とした（表 1 および図 1）。また陸上の観測点網である KiK-net、F-net から合計 17 個
表 1：解析対象地震のリスト（図１の赤いマーク）。の観測点を選択してその観測データを利用した（図 1）。これらの地震に対して FAMT 解析（「2.
研究の目的と意義」を参照；Okamoto et al., 2017,5.1.2 地震波シミュレーション2018）を適用して震源位置・震源メカニズム・震本研究では波動場のボルン近似に基づいて、源時間関数の各パラメータを推定し、感度カー

震源位置から広がるフォワード波動場と観測点ネル生成の際にこれらのパラメータを利用したから広がる相反波動場（アジョイント波動場）と（§5.1.3）。

を掛け合わせる方法により感度カーネルを生成する。アジョイント波動場については 51 回、フォワード波動場については 19 地震について
TSUBAME-3.0 を用いた大規模計算を実施した（表 2）。
全格子サイズ4896×3456×900

副領域サイズ544×576×454

格子間隔

150 m

時間ステップ数

35000

時間間隔

0.007143 s

ノード数27

GPU 数108

シミュレーション回数

70 回（=17×3+19）

1 回あたり計算時間

5795 s (平均値)

表 2

アジョイント波動場およびフォワード波動場計算のための TSUBAME-3.0 での計算パラメータ。

図 1 対象地域で 1996 年 1 月～2019 年 10 月に
発 生 し た 浅 い 地 震 の Global-CMT 解 （ www.

計算にあたっては TSUAME3.0 の NVIDIA

globalcmt.org）。赤でプロットしたものが解析対

P100 GPU に演算ブロックサイズを最適化したプ

象地震（19 個）。赤い三角は観測点の位置を示すログラムを利用した。また本課題の計算では（17 点）。
3

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

GPU 1 基に一つの副領域を割り当てる。その際、ル、下段は周期 7.1 秒の波に対する感度カーネル

副領域サイズとして GPU のグローバルメモリである。カラースケールは図ごとに異なる因子（16 GiB）のほぼ全てである 15.6 GiB を利用すで規格化したものであり、また、スケールは図ごるように設定した。とに異なっている。

5.1.3 感度カーネルの考察① ３次元不均質構造の効果

感度カーネルは地震の震源と観測点のペアにまず３次元不均質構造が波動場に及ぼす効果対して定義される（5.1.1 の積分式を参照）。ここ

を見るために、半無限構造（平坦な自由表面を持では 2017 年のアウターライズの地震 EV2017 とつ一様均質な媒質）の場合と比較する。図２に示観測点 TYS を例として、代表的ないくつかの場すように半無限構造モデルにもとづいて計算し合における感度カーネルの特性を観察する（震た場合は感度カーネルは震源と観測点を結ぶ直源位置と観測点位置は図 1 を参照）。

線経路に関してほぼ対称な空間パターンを示す

今回の検討では、対象領域を水平方向 10.8 kmが、３次元不均質構造モデルに基づいて計算し×10.8 km、深さ方向 6 km のブロックで分割し、た場合には直線経路に関して非対称で、やや屈各ブロック内でのカーネルの平均値をそのブロ曲したパターンを示す。これらは３次元付近異ックでのカーネルの値として図示する。波形ト質構造によって波動場や波の伝播経路が歪めら

モグラフィーの逆問題計算においても、このブれていることに対応する。このような３次元不ロック分割をそのまま利用することができる。

均質性の効果は短周期（短波長）の波ほど顕著になってくることもわかる（図 2 の下段）。
」② 深さ方向の変化
次に感度カーネルの深さ方向の変化を見る。
図 3 には①と同じ観測点・震源ペアに関する感度カーネルを深さごとに描いた。感度が強い（色が濃い）部分は 12 km よりも浅い部分に集中しており、とくに最も浅い 0-6 km のブロックの感度が強くなっている。これは、この例で震源として選択した EV2017 が浅い地震（深さ 13.1 km）であることを反映していると考えられる。浅い地震の場合には表面波が卓越する。本解析では感度カーネルの計算において全波形を用いているために、表面波が卓越する場合には結果的に

図 2 ブロック分割の第 2 層目（深さ方向 6-12 km）

表面波の影響が大きく現れる。そのため浅いブの観測点 TYS と震源 EV2017 のペアに関する P

ロックにおける感度が強く表れていると考えら

波弾性率の感度カーネルの実部。水色は感度がれる。このことから、深い層の構造モデルの改良ない（ゼロ）ことを示し、青または赤の色が濃いを効率的に進めるには、深い地震を用いたり、深ほど感度が強いことを示す。左列は３次元構造

い場所を主に通過する波（実体波）だけを抽出しモデルに基づいて計算した感度カーネル実部、て感度カーネルの計算に利用したりするなどの右列は半無限構造に基づいて計算したものを示対応が必要になると考えられる。

す。上段は周期 10.2 秒の波に対する感度カーネ4

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

次元構造モデルに基づいて計算した感度カーネル実部。右列は S 波の弾性率（剛性率）、左列は密度の感度カーネルである。上段は周期 10.2 秒の波に対する感度カーネル、下段は周期 7.1 秒の波に対する感度カーネルである。カラースケー
ルは図ごとに異なる因子で規格化したものであり、また、スケールは図ごとに異なっている。
③ パラメータによる違い
前と同じ観測点・震源ペアを用いた場合の、P 波弾性率以外のパラメータに関する感度カーネルを図 4 に示す。剛性率は３次元不均質性の
図 3 ブロック分割の第 1~4 層目（深さ方向 0-影響を受けやすく、直線経路に関する非対称性

6 km、6-12 km、12-18 km、18-24 km）の各層ごとが強く現れる。一方、密度は３次元不均質性の影の周期 7.1 秒の波に対する P 波弾性率の感度カ響がやや弱くなるようであるが、それでもややーネルの実部（規格化したもの）。3 次元構想モ短周期（7.1 秒）ではパターンの屈曲や非対称性デルに基づいて計算した。図 2 と同様に水色はが顕著になってくることがわかる。

感度がないことを示し、青または赤の色が濃い④ RMS 振幅の変化ほど感度が強いことを示す。0-6 km の図で白い部分は、該当ブロック全体が海水で構成されて波形に対する各パラメータの影響度が深さにいることを示す。この検討では、海水は物性の不よってどう変化するかを総合的に見るために、均質性が小さいとみなして感度カーネル計算か深さ方向のブロック層ごとに感度カーネルのら除外している。カラースケールは４つの図にRMS 振幅値を求めた。図 4 と 5 にブロックの深おいて共通の因子で規格化したものである。

さに対する各パラメータの RMS 振幅値をプロットした。図 4 の RMS 値はブロック層が深くなると急激に振幅が小さくなり、12 km よりも深いブロックでは RMS 値が浅いブロックの約 1/10 程度以下になる。この現象は短周期ほど顕著である（図 4 の下図）。これらは②で述べた内容に対応している。前述のようにこのケースでは表面波が卓越し、今回の計算ではその表面波が主な対象となることから浅いブロックの感度が卓越するようになると考えられる。
一方、陸寄りの深い地震である EV2007 では
深いブロック層でも波形への影響度がある程度保たれる（図 5）。特に震源（深さ 32.3 km）が位置する層（上から 7 番目、30-36 km）は RMS 値が大きくなっている。これは震源周辺の媒質に図 4 ブロック分割の第 2 層目（深さ方向 6-12 km）よる波形への影響が大きいことを示している。の観測点 TYS と震源 EV2017 のペアに関する３

このことから、深い場所に関する分解能を向上
5

学際大規模情報基盤共同利用・共同研究拠点 2020 年度共同研究 最終報告書させるためには、このような深い震源を解析にとでの計算手法開発も進めている。本課題研究

積極的に導入する必要がある。

では新規に導入された名古屋大学のスーパーコ

また、図 4、5 のいずれの場合でも S 波弾性率

ンピュータ「不老」の Type-I サブシステム（Fujitsu（つまり剛性率）の波形への効果が、P 波弾性率

PRIMEHPC FX1000）において地震波シミュレーや密度の効果よりも大きいこともわかる。これ

ションプログラムの大規模計算を実施した。まも、この計算では表面波が主な対象になっていたその結果を、過去に Fujitsu PRIMEHPC FX100ることと関係していると考えられる。また P 波で計算した結果（2019 年度 JHPCN 課題）と比較弾性率と密度の効果は同程度である。これらの検討した。

結果は、逆問題解析によってパラメータ摂動を

本報告書では大規模計算の例として、鹿児島推定するうえでは S 波弾性率（剛性率）が主要な県から沖縄県にかけての南西諸島を中心とする対象になることを意味する。さらにパラメータ

全長 1000 km 以上におよぶ範囲を設定したものを加える場合には P 波弾性率と密度を同時に追を説明する（図 6）。

加しなければならない。
計算マシン

図 4

深さ方向の各層ごとに感度カーネルのRMS 振幅値を深さに対してプロットしたもの。
観測点：TYS、震源：EV2017。MP： P 波弾性率

FX100

FX1000

格子サイズ

2201×5501×1001

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

時間ステップ

30001

時間間隔0.01 秒

ノード数432432

プロセス数8641728

スレッド数1612

計算時間7766 秒4309 秒表 3 FX100 と FX1000 の計算パラメータ（紫色）。MS：S 波弾性率（オレンジ色）。RHO：密度（グレー）。上図：周期 10.2 秒の波に関するグラフ。下図：周期 7.1 秒の波に関するグラフ。

計算方法は地球の曲率の影響を考慮した
quasi-Cartesian 座 標 に お け る 計 算 ス キ ー ム（Takenaka et al., 2017）を採用し，FX100 で行っ。
た過去の計算と同一のノード数とした（表 3）
計算対象とした地震は過去の計算と同じもので、2008 年 7 月 8 日 16 時 42 分に沖縄本島近海で発生した深さ 53 km，MJMA 6.1 の地震である。計算には名古屋大学情報基盤センターの不老 Type-Iシステム（FX1000）を使用し，MPI と OpenMPによるハイブリッド並列計算を行った。

図 5 震源を変更して図 4 と同様の RMS 振幅値をプロットしたもの。観測点：TYS、震源：EV2007。

計算の際、プロセス数は 1 ノードあたり 4 個（総計 1728 個）とした。これは FX1000 のプロセッサ（A64FX）が内部的には 12 個の演算コア

5.2 FX-1000 におけるプログラムの検討（＋1 個のアシスタントコア）及び HBM2 メモ

我々は GPU とは異なるアーキテクチャのも6

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

リからなる Core Memory Group（CMG）で構成さ1 ノード（＝1 プロセッサ）あたりのプロセス数れていることによる。1 つのプロセッサには 4 個も 4 個とした。また 1 プロセスあたりのスレッの CMG が内蔵されており、CMG 間はリングバ

ド数は 1 つの CMG 内の演算コア数と同じ 12 個スで結合されている。このようにプロセッサ内とした。

部が 4 つのグループに分けられていることから、

表 3 に示したように、計算に要した時間を比べると、FX1000 では FX100 の約 55%の計算時間で計算を行うことができた。この計算で用いているのは時間領域有限差分法（Finite DifferenceTime Domain: FDTD）であり、メモリ律速の手法となっている。そのためこの性能差はメモリバンド幅の差（FX100：ノードあたり 480 GB/s、FX1000：ノードあたり 1024 GB/s）が主要な原因ではないかと考えている。また、計算結果（出力波形）も FX100 の場合と FX1000 の場合とで一致していることを確認できた（図 7）。

6. 今年度の進捗状況と今後の展望
東北地方太平洋沖地震震源域の構造推定のパートでは、ほぼ予定通りの数の地震と観測点（それぞれ 19 個と 17 点）を選択して計算を実行す
図 6 計算領域。灰色のハッチがかかった領域がることができた。しかし今年度は年度当初から計算領域。逆三角形は観測点。震央とメカニズム

新型コロナ感染症が世界的に拡大し、この感染を震源球で示す。破線は琉球海溝。

症に対応するための業務等が大幅に増加した。
その結果、計算から得られる大規模な出力の整理や、その出力にもとづく解析などの、時間を要

図 7 N.YNGF における FX100（青線）と FX1000（赤破線）による計算波形の比較。フィルタはかけていない。左が波形全体図、(a)と(b)は左図の四角枠内をそれぞれ拡大表示したもの。
7

学際大規模情報基盤共同利用・共同研究拠点 2020 年度共同研究 最終報告書する作業に遅延を生じた。そのためこのパート

initial model, 日本地震学会 2020 年秋では、解析に用いる地震の震源パラメータを推季大会, 2020 年 10 月 29 日．

定する FAMT 解析の実施と、その結果にもとづ(5) 公開したライブラリなど該当なし。

いて波形トモグラフィーのための感度カーネル(6) その他（特許，プレスリリース，著書等）を生成し詳細に検討するステージまでを実施す該当なし。

ることとした。波形トモグラフィーによって目標とした構造モデル改良を試みることは継続課
参考文献

題として次年度（2021 年度）に実施する。その際には今年度得られた知見に基づいて、解析対象地震として深い地震を今よりも増やして深い

Hiroshi Takenaka, Masanao Komatsu, Genti

部分の物性パラメータをカバーできるようにす

Toyokuni, Takeshi Nakamura and Taroるなどの対応を追加する。

Okamoto, Quasi‑Cartesian finite‑difference

名古屋大学の「不老」Type-I サブシステム

computation of seismic wave propagation for（FX1000）ではほぼ予定通りに大規模地震波シa three‑dimensional sub‑global model, Earth,

ミュレーションを実施して、その結果の確認やPlanets and Space (2017) 69:67. DOI性能の考察を行うなど、基本的な検討を進める

10.1186/s40623-017-0651-1

ことができた。
Taro Okamoto, Hiroshi Takenaka, Takeshi

7. 研究業績一覧（発表予定も含む．投稿中・投稿予定は含まない）(1)

Nakamura

andTatsuhikoHara,simulationofearthquakes

offFDMwestern

学術論文 （査読あり）

Kyushu, Japan, using a land–ocean unified

該当なし。

3D structure model, Earth, Planets andSpace (2017) 69:88. DOI 10.1186/s40623-017-(2) 国際会議プロシーディングス （査読あ0672-9り）
該当なし。

Taro Okamoto, Hiroshi Takenaka and Takeshi(3) 国際会議発表 （査読なし）
T.

Okamoto,H.

Takenaka,Nakamura,T.

EvaluationofaccuracyofNakamura, Toward short-period (< 0 s)

synthetic waveforms for subduction‑zone

full waveform tomography in andearthquakes by using a land–ocean unified

around the 2011 Tohoku-Oki source

3D structure model, Earth, Planets andarea using land-ocean unified 3D

Space (2018) 70:98. DOI 10.1186/s40623-018-

initial model, American Geophysical0871-zUinon 2020 Fall Meeting, 2020 年 12 月
17 日．
(4) 国内会議発表 （査読なし）
岡元太郎・竹中博士・中村武史, Toward
short-period(<0s)

fullwaveform

tomography in and around the 2011 TohokuOki source area using land-ocean unified 3D8