学際大規模情報基盤共同利用・共同研究拠点 平成 30 年度共同研究 最終報告書 2019 年 5 月

jh180052-NAJ

地盤の非線形応答を考慮した
地震波シミュレーション：強震動予測と震源過程解析
竹中博士（岡山大学）
概要 本研究では、浅部地盤の非線形性を数値的に計算する順解析及び逆解析の手法の開発と、非線形計算手法の強震動予測や破壊過程解析への応用を進める。今年度は順解析について、南西諸島域を対象として研究を進めた。非線形波形生成のプロセスは大きく分けると、(i) 非線形物性を有する薄い浅部地盤に入射する地震波（基盤波）を大規模線形差分法によって計算する部分と、(ii) その基盤波を入射させたときの浅部地盤の非線形挙動を計算する部分の 2 つに分けられる。今年度の前半では (i) に相当する部分の研究を実施した。対象とした南西諸島では大規模線形差分法により基盤波の計算を行い、モデルの有効性を検討した。その結果、同じ地域でも観測点により観測波形の再現性が異なることが分かり、さらなる構造モデルの改良の必要性と方向性に関する知見を得た（§5.1）。今年度の後半には、南西諸島域の八重山諸島域を対象として(ii)に相当する部分の研究を実施した。この部分では、短周期波動まで含む基盤波を生成して、大きさの異なる地震を想定した複数の非線形シミュレーションを行い、八重山地域の観測点の土質によっては中規模の地震でも非線形効果が大きく表れる可能性があるなどの知見を得た（§5.2）。また、北海道では今年 9 月 6 日に 2018 年北海道胆振東部地震（M6.7）が発生し、強い地震動による液状化が起こるなどの広範な地震災害がもたらされた。そのため急遽、この地域も対象とすることとして、震源パラメータ推定に用いるグリーンテンソル波形生成のための地震波シミュレーションを実施した（§5.3）。

防災科学技術研究所・地震津波防災研
究部門・中村武史
東北大学・地震・噴火予知研究観測セン
ター・豊国源知
岡山大学・大学院自然科学研究科・小松正直

1. 共同研究に関する情報(1) 共同研究を実施した拠点名
東京工業大学 学術国際情報センター
名古屋大学 情報基盤センター(2) 共同研究分野■

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

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

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

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

共同研究者・FX100 でのCPU 版地震波プログラムの開発
共同研究者・CPU 版地震波プログラムを用いた
FX100 での計算
CPU 版地震波プログラムを用いた FX100 での計算

2. 研究の目的と意義
日本列島では、1995 年兵庫県南部地震（マグニ(3) 参加研究者の役割分担

チュード（M）7.3；阪神淡路大震災）、2011 年東北地方太平洋沖地震（M9.0；東日本大震災）、2016 年

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

研究代表者・全テーマの統括

熊本地震（M7.3）、2016 年鳥取県中部地震（M6.6）、副代表者・GPU 計算の最適化、大規模並列化・計算結果の可視化への助言
副代表者・TSUBAME でのGPU 版プログラムによる計算・連絡責任者

地震とそれに伴う地震災害が発生してきた。なか2018 年北海道胆振東部地震（M6.7）などの大きなでも 1995 年兵庫県南部地震では住家全壊が 10 万棟を越え、2016 年熊本地震でも震度７に達する強い揺れが 2 回発生して最新の耐震基準による住家にも倒壊が発生するなど、強震動による大きな地1

学際大規模情報基盤共同利用・共同研究拠点 平成 30 年度共同研究 最終報告書 2019 年 5 月

震災害がもたらされてきた。このような強震動のもとでは地表付近の浅部地盤での歪が大きくなり、 5. 今年度の研究成果の詳細線形限界を超えて非線形挙動を示すと考えられる。 5.1 沖縄・奄美諸島の大規模シミュレーション
したがってこれらの強震動を理解し予測するため（中間報告書の概要）には、地盤の非線形挙動をモデル化し定量的に計将来発生が予想される地震の強震動予測につい算することが必要となる。また地震断層の浅部でて、地震動を再現するためには、差分法などを用の破壊過程を探るうえでも、地盤の非線形挙動をいた大規模計算が必要となる。その際、高精度な把握しておくことは重要である。そこで本研究で3 次元地下構造モデルの導入が必要である。これは、浅部地盤の非線形性を数値的に計算する順解は、本課題における非線形応答計算のための基盤

析及び逆解析の手法の開発と、非線形計算手法の波を得るために重要である。

強震動予測や破壊過程解析への応用を進める。計

今年度前半は鹿児島県南方沖から沖縄県にかけ算手法としては、これまでに我々が研究してきたての南西諸島を対象とし、3 次元地震波速度構造線形の時間領域差分法（Finite-Difference Timeおよび非弾性減衰（Q）モデルの改良を行った。地Domain、 FDTD）による大規模地震波シミュレーシ

殻、マントル、プレートについては、最新の反射ョン手法との接続を考慮して、時間領域の差分法

法地震探査［1］やフィリピン海プレートの最新深をもとにして研究を行う。

度モデル［2］などを参考にモデルの修正を施した。

今年度は、既存の線形３次元差分法との接続をさらに、南西諸島のうち、沖縄・奄美諸島を対考慮することから、非線形地盤と線形の基盤とを象とし、地震波伝播の大規模シミュレーションを分離し、基盤からの地震波（S 波、横波）が地盤に行い、モデルの有効性を検討した。その結果、同鉛直に入射するものとして地盤を１次元の非線形じ地域でも観測点により観測波形の再現性が異な媒質として取り扱って計算をする。

ることが分かり、さらなる構造モデルの改良が必要であるなどの知見を得た。以上の結果について
は本課題の中間報告にて報告した。

3. 当拠点公募型共同研究として実施した意義
5.2 八重山諸島の非線形シミュレーション

本課題では基盤に入射する地震波を、3 次元構

5.2.1 対象領域と計算パラメータ造モデルのもとで大規模な線形地震波伝播シミュレーションを実施して計算する。最高周波数が数今年度後半は、短周期波動まで含む地震波伝播Hz 程度の短周期地震波を計算するためには数シミュレーションによる非線形応答の計算を試み

10 m ～ 100 m 程度の格子間隔を用いる必要があた。対象領域については、本年度前半に奄美・沖る。計算対象領域は数百 km×数百 km のサイズで縄本島域の計算と合わせて実施した予備的な計算あるため、大規模な計算格子サイズが必要となる。

では八重山諸島域での短周期波動での再現性が比

そのため、名古屋大学情報基盤センターおよび東較的に良好であったことと、八重山諸島域では浅京工業大学学術国際情報センターのスーパーコンい地震が実際に発生していることなどから、八重ピューターが不可欠であることから、本課題では山諸島域を対象として非線形計算を試みることとこれらの拠点を利用させていただくこととした。

した。計算領域を図 1 および図 2 に示す。本計算
では非線形応答を高周波数まで評価するため、格子間隔および時間ステップを細かく取る必要があ

4. 前年度までに得られた研究成果の概要

る。今回は表 1 に示すように、格子間隔を 20 m と本研究課題は今年度新規課題である。

した。これは今年度前半に行った沖縄・奄美諸島の広域地震波伝播シミュレーションの際に設定し2

学際大規模情報基盤共同利用・共同研究拠点 平成 30 年度共同研究 最終報告書 2019 年 5 月

図 1 南西諸島のシミュレーション範囲の広域図。
赤い枠が計算領域、黄色い星印は設定した地震の震央、赤い三角は観測点である。シミュレーションに用いた地震のメカニズムも示した。

格子サイズ

4001×5001×1051

格子間隔20 m

時間ステップ数

30001

時間間隔0.001 s

最高周波数10 Hz

非弾性要素数3

ノード数432

プロセス数864

スレッド数

16

計算時間

2 時間 12 分

表 1
八重山諸島のシミュレーションに用いた差分計算のパラメータと計算時間。計算には名古屋大学情報基盤センターの FX-100 を使用した。

大地震の波形は、この計算波形に震源の効果を畳み込むことによって合成できる。
5.2.2 弱震（Mw3.2）のシミュレーション
まず、上記のシミュレーションによって得られた計算波形をそのまま利用して、Mw 3.2 相当の弱い地震の場合を考察する。図 3 に計算波形を示す。
図 2 南西諸島のシミュレーション範囲の拡大図。
赤い枠が計算領域、黄色い星印はシミュレーションで設定した地震の震央、赤い三角は観測点である。

これらの 2 つの波形を非線形応答の計算に用いるた格子間隔の 5 分の 1 という細かなものである。

点と OKN011 観測点（図 2）の基盤でそれぞれ観この格子間隔のもとで、最高周波数として 10 Hz測される波形をシミュレーションしたものに相当までの波形モデリングが可能となる。

する。波形は水平成分のうちのトランスバース成基盤波として用いる。これらの 2 つの波形は、防災科学技術研究所の強震観測網の OKN009 観測分（震央と観測点とを結ぶ大円経路に垂直な方向シミュレーションには我々が開発したスキームの震動成分）である。

による差分法（HOT-FDM［3］）を利用した。計算には名古屋大学情報基盤センターの FX100 を使

本研究では鉛直入射を仮定した SH 波の非線形用し、MPI と OpenMP によるハイブリッド並列応答を計算するために、弾塑性レオロジーモデル計算を行った。震源パラメータには実際に発生し［5, 6］に基づくスキーム［4］を採用する。このた地震（2018 年 3 月 1 日 22 時 46 分、気象庁マ

モデルはバネとスライダーからなる複数の要素をグニチュード 5.6）のものを用いた。この地震によ直列に連結することで構成されている。これによる地震動を計算するための震源パラメータを表 2

り砂地盤や粘土地盤における地震動の非線形挙動に示す。なお、計算では実際よりも小さなモーメを計算することができる。今回はこれらのモデル

ント・マグニチュード（Mw 3.2）を使って近似的を 1 次元地盤構造モデルにおける線形 SH 平面波に「インパルス応答」に相当する波形を生成した。の時間領域差分法計算プログラム［7］に組み込み、3

学際大規模情報基盤共同利用・共同研究拠点 平成 30 年度共同研究 最終報告書 2019 年 5 月

緯度24.2400°

OKN009

S 波速度

層厚経度

123.8267°

第1層

190 m/s

3.6 m

深さ15 km

第2層

320 m/s

15 m

地震モーメント

7.94×1013 Nm

第3層

600 m/s

18.8 m―モーメント・マグニチュード

3.2 相当第4層

1100 m/s――

震源パルス幅0.1 s

震源メカニズム

気象庁 CMT 解

OKN011

S 波速度

層厚

土質

第1層150 m/s

3m

第2層

310 m/s

2.8 m

第3層

680 m/s

4.2 m

―

第4層

1100 m/s――

表2

シミュレーションに用いた震源パラメータ。

土質粘土

砂

表 3 非線形応答の計算に用いた構造。両観測点ともに第 1 層と第 2 層の非線形応答を計算した。

図 3 観測点 OKN009 と OKN011（図 2）とを対象として Mw3.2 の地震を想定して合成した基盤波波形。

鉛直入射を仮定して計算を行った。
今回計算を行った 2 観測点の表層地盤構造は、

図 4 OKN009 と OKN011 における Mw3.2 の地震を想定した非線形応答の計算結果（波形の全体図）。黒線が線形応答、赤線が非線形応答の計算結果。

観測点直下のボーリング調査による検層結果を用いて設定した。非線形の計算に用いた表層地盤の土質は、OKN009 は粘土層、OKN011 は砂層としーク振幅の減少が見られる。

た。設定した構造を表 3 に示す。

次に S 波部分の拡大図を図 5 に示す。上述のよ設定した構造について、線形と非線形のそれぞれの応答を計算した。まず波形の全体図を図 4 にうに、砂地盤を設定した OKN011 では S 波到達示す。これらの時系列波形では、大きな S 波パル付近では非線形計算波形（赤）の方が線形計算波スに引き続く高周波の後続波形について、減衰は形（黒）よりも高周波数の震動が減衰している。

顕著ではない。これは震動の振幅があまり大きくこのことから、このような弱い震動でも非線形効なく、非線形効果が小さいことによるものと考え果がある程度現れることが示唆される。

られる。一方、砂地盤を設定した OKN011 ではピ4学際大規模情報基盤共同利用・共同研究拠点 平成 30 年度共同研究 最終報告書 2019 年 5 月

図 5 OKN009 と OKN011 における Mw3.2 の地震を想定した非線形応答の計算結果（S 波到達付近の波形の拡大図）。それぞれ黒線が線形応答、赤線が非線形応答の計算結果。

図 7 観測点 OKN009 と OKN011（図 2）とを対象として Mw6.2 の地震を想定して合成した基盤波波形。

次にこれらの波形の S 波部分のフーリエ振幅スペクトルを図 6 に示す。特に砂層で計算した化のために震源距離（断層面上の各点と観測点OKN011 については非線形挙動として、スペクトまでの距離）の違いによる効果は省略した。合成

ル振幅の減少がある程度見られる。一方、卓越周波形（トランスバース成分）を図 7 に示す。

波数がシフトする効果は小さい。

この中規模地震の合成波形では、有限大きさの断層での破壊伝播モデルを考慮したことにより、図 3 よりも長い周期の震動が含まれており、S 波主要部分の幅も長くなっていることがわかる。またこの例では図 3 の場合よりも振幅も大きくなっている（ただし合成波形の振幅や卓越周期は設定した断層モデルや破壊伝播モデルにも依存する）。
まず非線形計算波形の全体図を図 8 に示す。これらの時系列波形では、高周波の後続波形について、非線形計算の場合に振幅が減衰する効果が顕

図 6 OKN009 と OKN011 における Mw3.2 の地震を想定した非線形応答の計算結果のスペクトル。
それぞれ黒線が線形応答、赤線が非線形応答。

著に表れている。さらに S 波主要部分を拡大してみると（図 9）、主要部分では高周波数成分の減少が顕著である。さらにもう一つの非線形効果として、長周期成分の位相遅れも大きいことがわかる。

5.2.3 強震（Mw6.2）のシミュレーション
次に強い震動について考察するために Mw6.2

次にこれらの波形の S 波部分のフーリエ振幅スの地震の場合を考える。ここでは有限な大きさをペクトルを図 10 に示す。図 6 の場合と同様に、砂持つ断層を前提として、文献［8］の方法に従って層で計算した OKN011 について高周波数でのス差分法波形から合成波形を生成した。ただし簡略

ペクトル振幅の減少が大きく、非線形挙動が顕著5

学際大規模情報基盤共同利用・共同研究拠点 平成 30 年度共同研究 最終報告書 2019 年 5 月

図 8 OKN009 と OKN011 における Mw6.2 の地震を想定した非線形応答の計算結果（波形の全体図）。黒線が線形応答、赤線が非線形応答の計算結果。

図 9 OKN009 と OKN011 における Mw6.2 の地震を想定した非線形応答の計算結果（S 波到達付近の波形の拡大図）。それぞれ黒線が線形応答、赤線が非線形応答の計算結果。に現れていることがわかる。
このように中規模地震においても、3 次元地球内部構造モデルを考慮して計算した基盤波波形をもとにして、南西諸島域の地震波観測地点においても土質によっては地盤の非線形効果が大きく現れる可能性があることが推定された。高周波数の減衰や位相の遅れなどの効果は、震源スペクトルや震源位置などの物理量を推定するうえでバイア
図 10 OKN009 と OKN011 における Mw6.2 の地震を想定した非線形応答の計算結果のスペクトル。
それぞれ黒線が線形応答、赤線が非線形応答。

スともなりえるものであるので、今後も検討を継続する予定である。
5.3 北海道胆振東部地震のシミュレーション
2018 年 9 月 6 日の未明（3 時 7 分頃）に北海道厚ら、急遽この地震も本課題の対象とした。

真町付近で発生した 2018 年北海道胆振東部地震

この地震では、破壊開始時の断層メカニズムと（M6.7）では、最大震度 7 に達する強い地震動に破壊過程全体の平均的な断層メカニズムが異なっ

よって 200 件に上る土砂災害（がけ崩れ・山崩れ）ていることが報告されている（防災科学技術研究が引き起こされ、数多くの人命が失われた（内閣所 WEB: http://www.hinet.bosai.go.jp/topics/

府調べ）。また北海道全域で長時間の停電が発生すishikari180906/?LANG=ja）。これは断層破壊過程るなど、ライフラインにも大きな被害があった。

が複雑なものであったことを示唆するため、まずこのように、北海道全域に大きな影響をもたらし、本震の破壊過程や余震の震源パラメータを推定す

強い震動とそれに伴う液状化等も発生したことかることに取り組む必要があると考えられる。そこ
6学際大規模情報基盤共同利用・共同研究拠点 平成 30 年度共同研究 最終報告書 2019 年 5 月

格子サイズ

1280×1280×1100

格子間隔50 m

計算領域サイズ

64 km×64 km×55 km

時間ステップ数

16000

時間間隔0.0025 s

最高周波数1.333 Hz

非弾性要素数6

ノード数8

副領域サイズ（1GPU あたり） 320×320×550
GPU（NVIDIA P100）の個数

32

1 シミュレーションあたり計算時間

1558.3 s（平均値）

シミュレーション回数33（11 観測点）

図 11 北海道胆振東部地震の計算領域（実線で囲まれた四角）。星印は気象庁の震央。黄色い丸は 1 カ月間の余震（M ≥ 3）。三角は計算に用いた観測点位置。波線で囲まれた範囲が、グリーンテンソル波形を出力した領域を示す。出力領域は深さ方向に10 km から 45 km までの範囲となる。

表 4 北海道胆振東部地震のシミュレーションに用いた差分法計算のパラメータと計算時間。計算には東工大のTSUBAME を使用した。

で本研究ではまず、これらを推定するためのグリーンテンソル波形計算を実施することとした。
石狩平野・勇払平野の長周期地震動の計算を対象として既存のボーリングや反射法・屈折法地震た。今後、この計算結果を用いて震源解析を試み探査・微動探査などの調査結果をまとめて作成さる予定である。

れた 3 次元構造モデル[9]を利用して、差分法計算用の構造モデルを作成した。シミュレーションの

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

計算パラメータは表 5 に示す。計算領域は図 11 に実線で示した四角領域である。このシミュレーシ［§5.1 および中間報告書］ 南西諸島を対象とし、ョンでは我々が開発した差分法スキーム（HOTFDM [3]）に基づいて、複数 GPU 演算を導入した

3 次元地震波速度構造および非弾性減衰（Q）モデプログラム[10]を利用した。計算には東京工業大

ルの改良を行い、さらに地震波伝播の大規模シミ学学術国際情報センターの TSUBAME を使用し、ュレーションによってモデルの有効性を検討した。

MPI による並列計算を行った。そして、図 11 で波その結果、同じ地域でも観測点により観測波形の線で囲んだ領域内で格子状に候補点震源を設定し、 再現性が異なることが分かり、さらなる構造モデ
ルの改良の方向性に関する知見を得た。
それらの格子点に点震源を置いたときの観測点でのインパルス応答に相当するグリーンテンソル波［§5.2］ 3 次元地球内部構造モデルを考慮して形を出力した。この計算では水平方向 1 km 間隔、計算した基盤波波形をもとにして、南西諸島域の深さ方向 0.5 km 間隔に、全 90,241 個の点震源に地震波観測地点における非線形応答の順解析を行対応する波形を出力した。

い、弱い震動と強い震動のそれぞれに対する非線形効果に関する知見を得た。

本課題ではこのような計算を図 11 に示した 11観測点の 3 成分に対して計算し、合計で約 300 万［§5.3］ 北海道胆振東部地震（2018 年 9 月 6 日）個のグリーンテンソル波形を生成することができの発生に伴い緊急に北海道地域の 3 次元構造モデ
7

学際大規模情報基盤共同利用・共同研究拠点 平成 30 年度共同研究 最終報告書 2019 年 5 月

ルを作成し、震源パラメータ推定に用いるグリー

Earthquake: Effects of Ocean-Bottom Topography andンテンソル波形生成のための地震波シミュレーシ

Seawater Layer, Bulletin of the Seismological Society of

ョンを実施した。

America,

102,

No.

6,2420-―2435,

doi:

10.1785/0120110356, 2012.
[4] Bardet JP, Tobita T, NERA: A computer program for7.

nonlinear earthquake site response analyses of layered

研究成果リスト(1) 学術論文

soil deposits, Department of Civil Engineering,なし。

University of Southern California, Los Angeles, 2001.
https://sites.google.com/site/tt60898/home/software.

(2) 国際会議プロシーディングスなし。

[5] Iwan, W. D., On a class of models for the yielding(3) 国際会議発表

behavior of continuous and composite systems, Journalof Applied Mechanics, ASME, Vol. 34, pp.612-617.

なし。
(4) 国内会議発表[6] Mróz, Z., On the description of anisotropic

なし。

workhardening, J. Mech. Phys. Solids, Vol.15, pp.163175, 1967.

(5) その他（特許、プレス発表、著書等）[7] 田中宏樹・竹中博士, 鉛直方向任意不均質弾性なし。

媒質における平面波入射問題の時間領域差分解法,謝辞

地震 第 2 輯, 57, 343-354, 2005.

本研究では（独）産業技術総合研究所活断層研究センターによる「石狩平野・勇払平野 3 次[8] Irikura K., Semi-empirical estimation of strong

元地盤構造モデル」を利用させていただきました。

ground motions during large earthquakes, Bull. Disas.

また、このモデルについて堀川春央氏から多くのPrev. Res. Inst., 33, 63-104, 1983.

助言をいただきました。気象庁と防災科学技術研[9] 吉田邦一・吉見雅行・鈴木晴彦・森野道夫・究所からは震源パラメータや地震波形データを提滝沢文教・関口春子・堀川晴央,長周期地震動計供していただきました。記して感謝いたします。

算のための石狩平野および勇払平野の 3 次元堆積
盆地構造モデル, 活断層・古地震研究報告, 第 7号, 1-29, 2007.

参考文献[10] Okamoto, T., H. Takenaka, T. Nakamura, and T.
[1] Nishizawa A., K Kaneda, M. Oikawa, D. Horiuchi,Aoki, Accelerating Large-Scale Simulation of Seismic

Y. Fujioka, C. Okada, Variations in seismic velocity

Wave Propagation by Multi-GPUs and Three-distribution along the Ryukyu (Nansei-Shoto) Trench

Dimensional

subduction zone at the northwestern end of theSolutions to Multi-scale Problems in Science andPhilippine Sea plate, Earth Planets Space, 69, 86, 2017.

Engineering, D. A. Yuen et al. (eds.), Chapter 24,[2] Iwasaki, T., Sato, H., Shinohara, M., Ishiyama, T.,

DOI:10.1007/978-3-642-16405-7_23, Springer-Verlag

Hashima, A., Fundamental structure model of island

Berlin Heidelberg, 2013.

arcs and subducted plates in and around Japan, 2015
Fall Meeting, American Geophysical Union, San
Francisco, Dec. 14-18, T31B-2878, 2015.
[3] Nakamura, T., H. Takenaka, T. Okamoto, and Y.
Kaneda,

FDM

SimulationofSeismic-Wave

Propagation for an Aftershock of the 2009 Suruga Bay8Domain

Decomposition,inGPU