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

jh130008-NA04

GPGPU による地震ハザード評価
青井真（独立行政法人防災科学技術研究所）

概要 平成 24 度までに、TSUBAME2.0 上に実装した GPU 版の地震動シミュレータ
GMS により大規模かつ現実的なモデルによる計算を行う環境が構築された。平成 25 年
度は、南海トラフ・相模トラフで発生が懸念される巨大海溝型地震を対象として数 100ケースの地震波伝播シミュレーションを実施し、長周期地震動の予測とそのばらつき評価を行った。また、TSUBAME2.5 での GMS の性能評価を行い、220 億格子のモデルの 1024GPU による性能が 107 TFLOPS になることを確認した。

1. 研究の目的と意義ュータシステムに積極的に取り入れている東京工兵庫県南部地震を契機に、地震に関する調査研

業大学の TSUBAME 上で大規模な地震ハザード評究の成果が国民や防災を担当する機関に十分に伝価を行う。これにより、これまでは主にリソース達され活用される体制になっていなかったというの制約から「決めうち」的に限られた想定震源モ課題意識の下に、行政施策に直結すべき地震に関デル（シナリオ）に対してのみ行われてきた従来する調査研究の責任体制を明らかにし、これを政のハザード評価の枠組みを超え、想定外を可能な府として一元的に推進するため、地震防災対策特限り生じさせないハザード評価のあり方を検討す

別措置法に基づき総理府に（現・文部科学省に）る。さらに、南海トラフなどで想定される巨大海地震調査研究推進本部が設置されている。地震調

溝型地震による長周期地震動のばらつきを含めた査研究推進本部の地震調査委員会において、地震定量的な評価を行うことを目的としている。

に関する観測、測量、調査または研究を行う関係行政機関、大学等の調査結果を収集、整理、分析

2. 当拠点公募型共同研究として実施した意義し、並びにこれに基づき総合的な評価を行ってお地震の被害を軽減するためには、個々人の地震り、この中で地震動予測地図の高度化を進めていへの防災意識を高め、地震に対する備えを促すこる。日本周辺で発生する全ての大地震に関して、とが不可欠である。このため、防災科学技術研究

地震のリスク評価の基礎となり得る精度で地震ハ所では、日本全国で発生する地震を対象として、ザードを予測できるよう、手法・モデルの高度化地震調査研究の成果の集大成である地震動予測地を目指している。そのために、必要な精度、分解図を高度化し、地震ハザード・リスク評価に関す能を持つ地盤構造の開発を行うとともに、高精度る研究を行うとともに、WebGIS 等の技術を用いかつ汎用性のある強震動シミュレーション手法のて、地震ハザード・リスク情報、地下構造データ開発に関する研究を行っている。

等の関連情報を網羅的に提供可能な地震ハザー

詳細な 3 次元地下構造を十分な精度で離散化しド・リスク情報ステーションを構築している。

高精度な計算を行い、短周期まで計算するために本研究では東京工業大学との超大規模数値計算は細かな格子が必要であるため、実用的な計算に応用分野での共同研究を行なっており、これによおいては格子数が数億から数百億に及ぶ極めて規り地震ハザード・リスク評価に用いられる地震波模の大きなモデルを扱う必要がある。我々はこれ

伝播シミュレーションに対する GPGPU の適用までに GPU で地震ハザード評価を行うシステムについて、実務等に使用される汎用計算機環境かを構築し、従来は困難であった規模のシミュレーら TSUBAME のような大規模なスーパーコンピションを大量に行うことが可能となった。本研究ュータシステムまでの様々な規模の計算機環境にでは、GPU コンピューティングをスーパーコンピおける評価が可能となっている。また、防災科学
1

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

技術研究所と京都大学防災研究所では、我々がこれまでに TSUBAME 上で構築してきた GPU による地震ハザード評価を行うシステムを最大限活用し、従来は極めて困難だった大規模モデルによるばらつきを含めた定量評価を実施し、日本の防災力向上に大きく貢献することが期待される。

3. 研究成果の詳細と当初計画の達成状況(1) 研究成果の詳細について
3.1.1 GMS による波動場の計算
本研究では、不連続な食い違い格子（Aoi and

図 1 （中央）計算に用いる不連続格子。（右）Fujiwara, 1999）を用いた、空間四次・時間二次不連続格子の垂直断面。領域 I と領域 II の接続

精度の差分演算子による実用コードである GMS部分で、内挿のために格子が重なっている。
（左）（Ground Motion Simulator、青井・他、2004）
をベースに GPGPU の適用性の評価を行なった。
GMS は、防災科学技術研究所によってパッケージ

領域 II の最上面（A 面）と領域 I の最下面（D面）における不連続格子の水平断面。

化された、3 次元有限差分法(FDM)により地震波3.1.2 GPU での差分計算処理

伝播シミュレーションを行うためのツール群でおGMS の差分計算処理については、ほぼ全てり、主に Fortran90 で書かれた差分計算ソルバはGPU 上で処理するように実装されている。開発環

ソースコードも公開されている。
差分法による地震波伝播シミュレーションを行

境には、NVIDIA 社から提供されている CUDA

う場合、格子サイズは計算すべき最短波長により（Compute Unified Device Architecture）を使用決定されるため、モデルのごく一部のみが低速度している。の媒質である場合でも計算領域全体を小さな格子

平成 24 年度までのプログラム開発により、に分割せざるを得ず、大規模なモデル計算の大きGMS 差分計算ソルバの単体 GPU 版プログラムは、な 障 害 と な っ て い た 。 こ れ に 対 し 、 Aoi andCPU 版プログラムに比べて 20.4 倍の性能となっFujiwara（1999）や青井・他（2004）は、大きた。また、複数 GPU 版プログラムについては弱さの異なる格子を組み合わせることにより効率的スケーリングに関しては、ほぼ線形の性能が得らかつ高精度に計算を行うことの出来る不連続格子

れ、約 220 億格子の 1024GPU による性能はによる差分法の定式化を提案した。GMS では、軟87TFLOPS になることが確認されている。

弱（＝地震波速度が低い）な表層が存在する浅い部分（領域 I）の格子点間隔は細かく、深い部分

3.1.3 TSUBAME2.5 におけるベンチマークテスト
2013（領域 II）の格子点間隔は領域 I の 3 倍の粗い格

年9月 にTSUBAME2.0

か ら子点間隔を有する格子モデルを用いている（図 1）。 TSUBAME2.5 にバージョンアップしたことに伴2 つの領域はオーバーラップしており、波動場のい、GMS 差分計算ソルバについて複数 GPU を用連続性が保たれるよう内挿される。典型的な盆地

いた並列性能の評価を実施した。計算モデルには、構造モデルの計算において、均質な格子による場21,560,000 格子のモデル(Unit420 と呼ぶ)を単位合と比較し数倍から十数倍程度効率が良いことがモデルとし、これを水平方向に 2×2、3×3、4×4、分かっている。

8×8、16×16、32×32 個並べて、全体の格子数がそれぞれ 4 倍、9 倍、16 倍、64 倍、256 倍、10242

学際大規模情報基盤共同利用・共同研究拠点 平成 25 年度共同研究 最終報告書 2014 年 5 月倍となるようなモデルを使用した（それぞれ、1.4

Unit840 、 Unit1260 、 Unit1680 、 Unit3360 、1.35

Unit6720、Unit13440 と呼ぶ）。浮動小数点演算

1.3は単精度で行っており、コンパイラには CUDA4.1性 1.25
能向 1.2上比 1.15および Intel Compiler 11.1 を用いた。なお、中間報告では Unit3360 までの評価を行っていたが、今回はさらに大規模なモデルに対して 1024 並列1.1までの評価を行っている。

1.05

まず、TSUBAME2.0 において同条件で取得し11.E+00ていた性能評価結果と比較したところ、性能向上
比 は 最 大 で 1.38 倍 で あ っ た （ 図 2 ）。

Unit420
Unit3360

TSUBAME2.0 から 2.5 へのアップグレードによ
図2

り単精度演算性能は約 3.5 倍、メモリバンド幅は約1.71.E+011.E+02

並列数
Unit840
Unit1260
Unit6720
Unit13440

1.E+03

Unit1680

TSUBAME2.0 に対する TSUBAME2.5 の性能向上比と並列数との関係。

倍 に 向 上 し た と さ れ て い る（http://tsubame.gsic.titech.ac.jp/node/966）。これまでの検討により、GMS 差分計算ソルバの律速

1.E+03

972.44

要因は演算処理ではなくメモリ転送速度であることが分かっており、性能向上比の最大値はメモリバンド幅の向上に規定されていると考えられる。
1.E+02

弱スケーリングに関しては、TSUBAME2.0 と高速化
率

同 様 に ほ ぼ 線 形 の 性 能 が 得 ら れ た （ 図 3 ）。
Unit13440（約 220 億格子）の 1024GPU による性能は 107TFLOPS になることが確認された。

1.E+01

なお、GMS 差分計算ソルバを CUDA5.0 にてコンパイルし上記と同様の性能評価を実施したが、性能は CUDA4.1 の方が高いという結果であった。
今回は CUDA5.0 対応のために最小限のコードチ

1.E+00
1.E+00

1.E+01

1.E+02

1.E+03

並列数ューニングを行ったが、今後は CUDA5.0 や 6.0の新機能を活かしたコードチューニングを行うことにより、更なる性能向上も期待される。

図 3

Unit420

Unit840

Unit1260

Unit3360

Unit6720

Unit13440

Unit1680

TSUBAME2.5(CUDA4.1)における並列数と高速化率の関係。

3.1.4 南海トラフの地震の長周期地震動評価
日本周辺では千島海溝、日本海溝、相模トラフ、く存在している大都市の大多数は平野上に発達し南海トラフなどのプレート境界でマグニチュードており、海溝型地震に対しては強震動と津波の評

8 を超える巨大海溝型地震が発生し、強震動や津価とともに長周期地震動の評価も重要である。平波によって甚大な被害が引き起こされてきた。こ成 25 年度は、南海トラフ・相模トラフで発生するれまでに発生した巨大海溝型地震では震源から遠巨大海溝型地震による長周期地震動についての検く離れた平野や盆地においても高層ビルや石油タ討を行ったが、ここでは代表として南海トラフにンクなどの長大構造物が長周期地震動による被害ついての結果を示す。
2013 年 5 月に公表された「南海トラフの地震活を受けたケースが知られている。長大構造物が多3

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

動の長期評価（第二版）について」［地震調査委員会（2013）；以下、長期評価と略す］では、これまで考えられてきた「固有地震モデル」に基づいた評価ではなく、2011 年東北地方太平洋沖地震の経験を踏まえ、南海トラフにおいて過去に発生したことが確認されていないものの、現在の科学的知見に基づいて推定されうる最大クラスの地震も含めた地震の多様性が考慮されている。次の地震として多様な想定震源モデル（シナリオ）が想定される中では、少数の限られたシナリオのみを考慮した地震動の評価では不十分であり、多数のシナリオを設定して、各シナリオに対する地震動シミュレーションを行うことに加え、それら予測結果を用いてばらつきを評価することが重要となる。
しかしながら、これまでは主にハードウェア－リソースが不十分であったため、大規模なシミュレーションを伴う長周期地震動の評価は極めて限られた数の想定震源に絞らざるを得なかった。そこで、本研究では大規模な GPU 環境が使用可能な図 4 南海トラフの地震として想定した 15 通り

TSUBAME をプラットフォームとし、数百ケースの震源断層面。灰線で示した 18 の小領域の組みの長周期地震動シミュレーションに基づいた評価合わせにより、赤線で示した震源域を設定していを試みる。

る。

地震波伝播シミュレーションに用いる地下構造モデルと震源モデルは、次のように設定した。まず、地下構造モデルは、全国１次地下構造モデル
から、東西 1150 ㎞、南北 950 ㎞、深さ 100 ㎞までの領域を切り出して、水平方向 200m、深さ方向 100m（深さ 8km 以深ではその 3 倍の大きさ）の格子サイズでモデル化した。格子数は約 32 億である。
震源モデルについては、長期評価を参考として15 通りの震源域（図 4）を設定し、さらに、震源図 5 アスペリティの配置（i～iv）と破壊開始点域内で最初に破壊が始まる場所（破壊開始点）やの配置（v）。

地震波を強く励起する領域（アスペリティ）についても数通りを設定した（図 5）。これらの中から関数を変えた 90 個の震源モデルを作成し、計 369

計算結果のばらつきが大きくなると考えられる組個の震源モデルに対する計算を行った。

合せを選んで 279 個の震源モデルを作成した。さ地震波伝播シミュレーションでは、地下構造モらに特定の震源域を対象として、他のばらつき要デルを水平方向に 9×9 の部分領域に分割して、そ因となる断層面上での破壊伝播速度やすべり時間れぞれ１つの GPU を割り付けて計算処理を行っている。計算結果は、GMS の 1D Dump 機能を4

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

図 6 計算時間の頻度分布。
利用し、地動速度の 1 成分を 1 ファイルとして東西、南北、上下の 3 成分について出力している。1ファイルには、陸域に約 2km 間隔で設定した約 8万観測点の毎タイムステップ（全部で 72000 ステップ=600 秒）の地動速度値が記録されており、1ファイルで約 90GB である。以上の計算では、TSUBAME2.5 の S キュー（X キュー）を利用し、81GPU（27 ノード）を使用した実行時間は概ね2.5～2.7 時間程度であったが、最大で 11 時間を要したケースもあった（図 6）。防災科学技術研究所の計算機（Altix4700）を用いた計算では異なるデータ出力方式を採用しているため厳密な比較はできないが、同じ地下構造モデルと震源モデルに対して 252core（Intel Itanium 1.66GHz）を用いた場 合 の 実 行 時 間 は 約 43 時 間 で あ り 、TSUBAME2.5 を使用することで数百ケースの計

図 7

算を現実的な時間で実行することが可能となって南海トラフの地震を想定した地震波伝播シミュレーションによる、最大地動速度値（PGV）いる。の空間分布。黒線で囲んだ領域が震源断層面（図出力された 1D Dump ファイルに対するポスト

4 の No.4 に対応）
、震源断層面内の小領域がアス処理として、フィルタ処理、ダウンサンプリングペリティ（図 5 の i と iv の n1 との組合せ）、赤星処理を行い、各観測点の速度時刻歴を抽出すると印は破壊開始点を示している。

ともに、最大地動速度（PGV）や速度応答スペクトル（Sv）の空間分布を求める面的処理を行った。

解できる。

計算結果の一例を図 7 に示す。これら 10 ケース

多様な震源モデルを考慮して計算された長周期では、震源域（No.4）とアスペリティ配置は同じ

地震動がどの程度のばらつきを有しているかを見（地震規模が同じ）であるが破壊開始点の違いにるために、大阪府庁を代表地点として最大地動速より長周期地震動の分布が大きく異なることが理5

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

ケースに対するヒストグラム（最上段）は数 10～100 倍程度の大きなばらつきを有しており、その分布形状はほぼ対数正規分布で表されている。これらを 15 個の震源域毎に分けて作成したヒストグラムからは、震源域が広いほど（地震規模が大きいほど）振幅レベルが高い傾向が認められ、さらに、すべりの大きな領域が西側（△）で、破壊開始点が西側（赤いバー）の場合に振幅が大きくなる傾向もみられる。
図 4 に示した震源域には、最大クラスの超巨大地震から、南海地震や東南海地震といった南海トラフの一部が破壊するものまで様々な規模の地震が含まれており、それらによる長周期地震動の大きさは図 8 に示すように数 10 倍から 100 倍近いばらつきを有している。一方、地震には規模の大きなものほど発生頻度が低いという一般的な関係があり、図 8 に含まれるひとつひとつの計算結果が同じ頻度で発生するわけではない。
発生頻度を考慮した評価として、確率論的地震ハザード評価がある。そこでは、ある地点における地震動強さとその超過確率の関係がハザードカーブとして表現される。地震調査研究推進本部から公表されている確率論的地震動予測地図（確率地図）は震度に関するハザード評価結果である。
確率地図での震度の評価は経験式（距離減衰式）に基づいており、それは地震規模と震源からの距離で決まる平均値とそのばらつきで表現される。
震度は本研究で対象としているものよりも短周期帯域の地震動を反映した指標であり、非常に多くの観測記録を基にその経験式が構築されている。

図 8 大阪府庁における、最大速度値と速度応答

一方で、本研究で対象としている長周期地震動に値（固有周期 5 秒、10 秒、減衰 5％）の分布。

ついては特に巨大地震による記録が少ないという最上段は 279 ケースについてのヒストグラム。2

課題がある。そこで、ここでは経験式ではなく、段目以下は各震源域についてのヒストグラムと本研究で得られた多数のシミュレーションによる累積頻度分布。No. 1～15 は、図 4 に示した震源地震動評価に基づいた長周期地震動の確率論的地域番号に対応している。ヒストグラムは破壊開震ハザード評価を試行する。

始点と最大値の関係、累積頻度分布はアスペリ

本研究で行った 369 ケースの計算結果に、本研ティ配置と最大値の関係を示している。

究以前に行っていた 100 ケースの計算結果を加え度値と相対速度応答値（周期 5 秒、10 秒；減衰 5%）た 469 ケースの結果を用いる。また、確率地図のの頻度分布（ヒストグラム）を示す（図 8）。279

中で設定されている南海トラフの地震の発生パタ
6学際大規模情報基盤共同利用・共同研究拠点 平成 25 年度共同研究 最終報告書 2014 年 5 月(2) 当初計画の達成状況について大規模計算を要する長周期地震動のシミュレー
ション研究では、主として過去に発生した地震の再現を目的として数～10 ケース程度の計算に基づいた研究が行われているが、本研究では巨大海溝型地震による長周期地震動の予測を目指し、数百ケースのシナリオを設定し、各シナリオに対する予測に加え、シナリオによる結果のばらつきの評価を行うことを計画していた。我々はこれまで防災科学技術研究所のスーパーコンピュータを使用して、2 年間で 100 ケース程度の長周期地震動シミュレーションを実施してきたが、平成 25 年度
は最終的に南海トラフでは 369 ケース、相模トラ
フでは 246 ケースのシナリオに対する長周期地震動シミュレーションと、その結果をもとにしたばらつき評価を行った。このような数百ケースの計算結果に基づいたばらつきの評価は、TSUBAMEを用いたことではじめて実現できたものといえる。
図 9

南海トラフの巨大地震に対する長周期地震

さらに、当初計画には含まれていなかったTSUBAME2.5 におけるベンチマークテストを実動の条件付きハザードマップ。

施した。平成 25 年度は、当初配分された計算資源ーンと各パターンの重み（発生頻度に対応する数量が課題申請時の想定を下回っていたが、研究期値）を利用する。これにより、南海トラフの巨大間中に計算資源を追加で配分していただいたこと地震が発生したという条件付きのハザードカーブで、数百ケースにおよぶ計算を実施し、さらに高を任意の地点に対して求めることができる。この並列環境における性能評価を行うことが可能とな作業を面的に行い、ある超過確率での長周期地震った。

動強さ（最大速度、速度応答）を地図として示し4.

たのが図 9 である。条件付き超過確率 50%（2 回今後の展望に 1 度は地図に示した値を超える地震動強さとな

GMS ソルバの高度化に向けて、CUDA5.0 や 6.0る）でも平野部では長周期地震動が大きくなる傾の新機能を用いたコードチューニングを行うこと向が見られるのに対し、条件付き超過確率 10%により、TSUBAME2.5 の性能を活かした計算を（10 回に 1 度に相当）であっても平野部以外では行うことを目指す。
地震動シミュレーションに関して、本研究では長周期地震動はそれほど大きくなっておらず、長周期地震動の特徴が反映された地図となっている。 震源モデルの不確実性を考慮して多数の長周期地
図 9 は確率論的ハザード評価のある一面を切り出震動シミュレーションに基づいた確率論的ハザーして示したものであり、個別地点のハザードカード評価を試行した。しかしながら、今回は現状のブや、様々な超過確率や地震動強さにおけるハザ

強震動予測で用いられている標準的なパラメーードマップをもとにハザード評価の結果をより詳

タ設定の範囲内での不確実性の考慮となってい細に検討することが重要である。

る。今後は、断層面上のすべりや破壊伝播速度のランダム不均質性やすべり時間関数などにつ7

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

いての更なる不確実性を考慮したシミュレーショ

SE04-A020.

ンにより、ハザード評価に資するデータセットをT. Maeda, N. Morikawa, S. Aoi, and H.

拡充していくことが重要である。また、これまで

Fujiwara (2013) Hazard assessment of

周期 3～20 秒程度の長周期地震動を対象とした検

long-period ground motions for the Nankai

討を行ってきたが、これは 30 階建て程度以上の高Trough

層建物や石油タンクなどの長大構造物に影響を及Geophysical Union 2013 Fall Meeting, in

ぼす周期帯であり、より数多く存在しているよりSan Francisco, S43A-2482.

earthquakes,

American(4) 国内会議発表

小規模な構造物に対しては短周期地震動の評価が必要である。周期 1 秒よりも短周期帯域に対して

青井真・前田宜浩・青木尊之（2013）強震動シ

は地震動の統計的性質に基づいた簡易な地震波合ミュレータ GMS による超大規模計算，物理成法が確立されていることから、今後は周期 1 秒探査学会 第 129 回学術講演会（高知市）．

以上の長周期帯域を対象としたシミュレーション

前田宜浩・森川信之・青井真・藤原広行（2013）

が重要な課題となる。そのためには、より精緻に南海トラフの海溝型地震による長周期地震動

モデル化された地下構造モデルや震源モデルを用ハザード評価，日本地震学会秋季大会（横浜いる必要がある。現在、関東など一部の地域では市），P3-61．

地下浅部の不均質さが反映された地下構造モデル

前田宜浩・森川信之・青井真・藤原広行（2014）が構築されており、その成果を取り入れることに南海トラフの巨大海溝型地震を対象とした長より周期 1 秒程度までを対象とした地震動シミュ

周期地震動ハザード評価，日本地球惑星科学レーションを目指す。

連合 2014 年大会（横浜市），SSS23-13．
森川信之・前田宜浩・岩城麻子・今井隆太・青5.

研究成果リスト

井真・藤原広行（2014）相模トラフの巨大海(1) 学術論文（投稿中のものは「投稿中」と明記）

溝型地震を対象とした長周期地震動ハザード

T. Maeda, N. Morikawa, A. Iwaki, S. Aoi, and評価へ向けた検討，日本地球惑星科学連合

H.

Finite-Difference

2014 年大会（横浜市），SSS23-14．

Simulation of Long-Period Ground Motion(5) その他（特許，プレス発表，著書等）

for

Fujiwara
the(2013)

NankaiTroughMegathrust

本研究は文部科学省の「長周期地震動ハザード

Earthquakes, JDR, Vol.8, No.5, p912-925.

マップ作成等支援事業」と連携して進められておA. Iwaki, N. Morikawa, T. Maeda, S. Aoi, andり、本研究で得られたシミュレーション結果の一

H.

Finite-Difference

部は地震調査研究推進本部（文部科学省）のもとSimulation of Long-Period Ground Motionに設置された強震動評価部会等の委員会に資料提for供され、審議されており、長周期地震動予測など

Fujiwara
the(2013)

SagamiTroughMegathrust

Earthquakes, JDR, Vol.8, No.5, p926-940.のハザード評価に活用される予定である。

(2) 国際会議プロシーディングス(3) 国際会議発表
T. Maeda, W. Suzuki, S. Aoi, and H. Fujiwara(2013)

Finite-Difference

SimulationofLong-Period Ground Motions of the 2011
Tohoku-Oki,

Japan

Earthquake,

Asia

Oceania Geosciences Society 10th Annual
Meeting,inBrisbane,

Australia,8