活断層・古地震研究報告，No. 2, p. 341-357, 2002

大阪平野における地震動予測の試み−動的地震破壊シナリオと 3 次元地下構造モデルによる地震動計算−
Ground motion prediction in Osaka Plain based on dynamic rupture
scenario and three-dimensional subsurface structure model
関口春子１・加瀬祐子２・堀川晴央 3・佐竹健治 4・杉山雄一 5・Arben Pitarka6
Haruko Sekiguchi1, Yuko Kase2, Haruo Horikawa3,
Kenji Satake4, Yuichi Sugiyama5 and Arben Pitarka6
1, 2, 3, 4, 5

活断層研究センター（Active Fault Research Center, GSJ/AIST, haruko.Sekiguchi@aist.go.jp,
kase.yuko@aist.go.jp, h.horikawa@aist.go.jp, kenji.satake@aist.go.jp, sugiyama-y@aist.go.jp）
6
URS Corporation（URS Corporation, Arben_Pitarka@URSCorp.com）

Abstract: We present ground motion simulations in the Osaka sedimentary basin from a rupture on the
Uemachi fault system. This is the first result of new project of the Active Fault Research Center, AIST, topredict ground motion for expected future large earthquakes based on realistic earthquake rupture
modeling and three-dimensional velocity structure constrained by geological and geophysical data. The
Uemachi Fault system is a 45 km-long reverse fault underlying the center of Osaka City with an uplift
rate of 0.4 m/ky (Osaka Prefecture, 1999). Rupture of the whole fault plane would cause an M 7earthquake and heavily shake the whole Osaka plain. We simulated ground motions, using
three-dimensional finite difference method (Pitarka, 1999), from the dynamic rupture scenarios (Kase etal., 2002, in this issue) in the three-dimensional structure model (Horikawa et al., 2002, in this issue). The
simulation assumes linear site response and the frequency range up to 1 Hz. Modeling space includes the
Uemachi fault at its center, and expands 38 km in EW direction and 58 km in NS direction. Different
rupture propagation directions in each scenario resulted in different ground motion distributions. The
distribution of strong motion seems to be mainly controlled by the rupture process, and the strength ofshaking is affected by the subsurface structure. Rupture propagations from different hypocenter locations
and segment-gap distances cause variable forward directivity effect. Seismic waves are diffracted at thebasin edge or basement step buried beneath the Neogene-Quaternary deposits and are amplified by theyounger deposits with low velocity. Ground motions become largest in a 7 to 8 km-wide zone along the
updip projection of the fault plane. The computed maximum peak ground velocities, about 200~300 cm/s
or more, would be reduced by an appropriate heterogeneity in the dynamic rupture, which would weaken
the forward directivity effect.
キーワード：上町断層系，大阪平野，3 次元差分法，動的地震破壊シナリオ，セグメント
Keywords: Uemachi fault (system), Osaka basin, three-dimensional finite difference method, dynamic
rupture scenario, segment

１．はじめに活断層研究センターでは，平成 13 年度のセンター発足と同時に，活断層情報から，アスペリティなどの不均質震源特性を予測する可能性を検証するため，
3 次元地下構造モデルと起震断層モデルによる，将来の地震と地震動の予測研究を開始した．地質調査所時代からの活断層および地下構造調査の蓄積の多い大阪平野をモデル地域とし，平成 16 年度には強震
動予測地図を作成する予定である．本報告では，上町断層系を想定震源断層とし，大阪平野における 1
Hz までの強震動のシミュレーション結果を報告す

る．なお，一連の研究のうち，地下構造のモデリン
グについては，堀川ほか（2002, 本報告書），破壊シ
ナリオの想定については加瀬ほか（2002, 本報告書）
を参照していただきたい．
大阪平野は人口が集中した堆積平野であり，大きな地震被害を受ける恐れのある環境にある．大阪平野に強震動を生じさせ得る起震断層のうち，今回は，上町断層系を取り上げた．上町断層系は，千里丘陵の佛念寺山断層，大阪市街地直下に存在する上町断層，泉北丘陵の坂本断層などからなる全長約 45 km，
平均上下変位速度 0.4 m/千年に達する逆断層系である（大阪府，1999）．全体が破壊すれば M7 前後の規

関口春子・加瀬祐子・堀川晴央・佐竹健治・杉山雄一・Arben Pitarka

模になると考えられる．そして，平野直下に伏在する断層と柔らかく厚い堆積層により，大阪平野全体に大きな地震動が生じるものと予想される．三浦ほか（2002）の反射法地震探査とボーリングによる複合調査によれば，新淀川付近の上町断層の最新活動
時期は 9,000 年以上前であり，このときの断層の上下変位量は 1.8 m 以上と推定されている．仮に現在まで，0.4 m/千年程度の平均上下変位速度を保っているとすると，最新の活動で解放されただけのひずみが，現在までに再び蓄積されている可能性がある．
地震破壊シナリオは，動的破壊過程の数値シミュ
レーションに基づいて作成した（加瀬ほか，2002）．
反射法地震探査やボーリング，地形地質調査結果など，これまでに得られている情報をもとに，2 つのセグメントからなる上町断層系モデルを設定した．
破壊開始点位置，セグメント間のジョグの構造，断層破壊時の応力降下量などの未知要素は，地下構造，
地形などから考え得る範囲内で様々な可能性を検討した．
大阪盆地の 3 次元地下構造モデルは，重力異常の分布，反射法地震探査やボーリングの結果，地層の堆積過程を総合してモデル化し，3 次元地下構造モ
デルのメッシュ数値データを作成した（堀川ほか，
2002）．
地震動のシミュレーションには，3 次元差分法（Pitarka, 1999）を用い，周波数 1 Hz 程度まで，線形応答領域のみの計算を行った．計算空間は上町断層系を中心とする東西 38 km，南北 58 km，深さ 25 kmの領域である．本稿では，3 次元差分法のプログラ
ムの概略，動的震源モデルおよび 3 次元地下構造モデルの導入方法，地震動シミュレーション結果について述べる．

テンソル，λ , µ はラメ定数，∂ t と ∂ m は時間および空間 m 方向の微分， δ lm はクロネッカーのデルタ（ δ lm = 1

for l = m ， δ lm = 0

for l ≠ m ）．

(1)式の両辺を ρ で割り， ∂ t u l = v l とする．
∂ t vl = b(∂ mτ lm + f l )(3)

ここで， b = 1 ρ と置き直している．(2)の両辺を時間で微分すると，∂ tτ lm = λδ lm∂ l vm + µ (∂ l vm + ∂ m vl ) .

( 4)

次に，(3)，(4)式中の時間微分を 2 次の前進差分，
空間微分を 4 次の中央差分で近似する． x に関する任意関数 g (x ) を考える．差分を取る値の存在する間隔が一定値 h であれば， g ( x ) の 2 次の前進差分は，
1(g (x + h ) − g (x ) )h
4 次の中央差分は，
1 9∂ x g ( x) ≅ (g ( x + h ) − g ( x − h ) )2h  8∂ x g ( x) ≅

−

(5)

1
(g ( x + 3h) − g ( x − 3h) )
24


( 6)

となる．
4 次の中央差分について，x の離散間隔が一定でない場合を考える．このとき，(6)式右辺の各項の係数は格子間隔に依存した値を取る．Fig. 2 に示されるよ
うなグリッド間隔の分布を考え，∂ x g ( x ) ≅ Dx g (x ) = c1g ( x + ∆1 ) + c2 g ( x − ∆ 2 )( 7)と置く． g ( x ) = g' exp(ikx ) という単一波長の平面波を

地震動のシミュレーションは，Pitarka (1999) による不等間隔格子 3 次元差分法を用いた．Pitarka (1999)の差分法は，不等間隔のスタッガードグリッドを使い，速度‐応力によって弾性波の場が定式化されている（Fig. 1）．空間に関する差分は 4 次精度である．
グリッド間隔は不等間隔に指定できるため，等間隔グリッドの計算に比べ，計算の精度を保ちつつ数値
計算の規模を小さくすることができる．以下にプロ
グラムについて概略を説明する．
2.1 波動方程式の差分表現
連続体中の要素の運動方程式は，等方弾性体中の応力‐ひずみ関係式は，
τ lm = c lmpq e pq = λδ lm ∂ l u m + µ (∂ l u m + ∂ m u l )分の q 方向への移動に関する変化），c lmpq は弾性定数+ c3 g ( x + ∆ 3 ) + c4 g ( x − ∆ 4 )

２．数値計算法

ρ∂ t 2ul = ∂ mτ lm + f l ,は m 方向ベクトルに垂直な面に働く応力の l 成分，
f l は要素に働く体積力， e pq はひずみ（変位の p 成

(1)

( 2)

である．ここで， ρ は密度， ul は変位の l 成分，τ lm

考え，これを(7)式の両辺に代入すると，（左辺）＝ ik ⋅ g ' exp(ikx )（右辺）＝ c1g ' exp[ik ( x + ∆1 )] + c2 g ' exp[ik ( x − ∆ 2 )]+ c3 g ' exp[ik ( x + ∆ 3 )] + c4 g ' exp[ik ( x − ∆ 4 )]よって，
ik ≅ c1 exp(ik∆ 1 ) + c 2 exp( −ik∆ 2 ) + c3 exp(ik∆ 3 )+ c 4 exp( −ik∆ 4 )(8)

これを，4 次精度で等しくする．指数関数の Taylor
展開の 4 次精度
1
1
exp(ik∆ ) = 1 + ik∆ − k 2 ∆2 − ik 3∆3 + O ( ∆4 )2
6を(8)式の右辺に当てはめると，(8)式右辺＝

大阪平野における地震動予測の試み−動的地震破壊シナリオと 3 次元地下構造モデルによる地震動計算−( c1 + c2 + c3 + c4 ) + ik ( c1∆1 − c2 ∆ 2 + c3 ∆ 3 − c4 ∆ 4 )1 2
k ( c1∆ 12 + c2 ∆ 2 2 + c3 ∆ 3 2 + c4 ∆ 4 2 )2
1− ik 3 ( c1∆ 13 − c2 ∆ 2 3 + c3 ∆ 33 − c4 ∆ 4 3 )
6−+ O ( ∆ 14 , ∆ 2 4 , ∆ 3 4 , ∆ 4 4 )( 9)いかなる波数の波に対しても，(8)式が 4 次精度で満たされる必要から，次のような連立方程式が立てられる．
1
11   c1   0 
 1
   

∆
−
∆
∆
−
∆ 4   c2   1 
 1
2
3
(10)
 ∆ 2 ∆ 2 ∆ 2 ∆ 2   c  =  0
3
1
2
3
4
   
 3
3
3
3 
  
∆
 1 − ∆ 2 ∆ 3 − ∆ 4   c4   0 
これを解くことにより，差分オペレータ D x が得られる．スタッガードグリッドで半グリッドずれると，グリッド間隔の関数である ∆ i の表現が変わるので，1 組の応力と速度のセットについて，(7)式の差分オ
ペレータは 2 種類必要になる（Fig. 2）．
(3)，(4)式の時間微分を(5)式の 2 次差分，空間微分を(6)式の 4 次差分で置き換えると，
v x n +1i , j ,k = v x n −1 2i , j ,k(+ ∆t bx Dxτ xx + D yτ xy + Dzτ xz + f x) n i , j ,k(+ ∆t b y D xτ yx + D yτ yy + D zτ yz + f yv x n +1i +1 2, j ,k +1 2 = v x n −1 2 i +1 2, j ,k +1 2(

+ ∆t b x D xτ zx + D yτ zy + D zτ zz + f z

τ xx

n

i +1 2, j ,k = τ xx i +1 2, j ,k + ∆t)]n +1 2+ λ (D x v x + D z v z )]n +1 2

(

τ yy
τ zz

n +1

n

i +1 2, j ,k = τ yy i +1 2, j ,k + ∆tni +1 2, j ,k = τ zz i +1 2, j ,k + ∆t(

) n i +1 2, j +1 2,k) n i +1 2, j,k +1 2[(λ + 2µ )D x v x+ λ D y v y + Dz v z
n +1

i +1 2, j ,k[(λ + 2µ )D y v y(11)

される 1 タイムステップ分の減衰係数，
 − π f 0 ∆t A( x, y, z ) = exp (12) Qs( x, y, z ) の 2 乗したものを毎タイムステップ，速度が更新される度に掛けている．なお，Graves (1996)では(11)式の減衰係数を，毎タイムステップ，速度および応力値に掛ける．ここでは，速度と応力それぞれに掛ける代わりに，速度のみに(12)式の 2 乗を掛けている．
これは，初期応力場がゼロである場合に同じ効果になると考えられ，数値実験でもほぼ同様の効果が確認されている．ある周波数 f における非弾性減衰の)] n+1 2 i +1 2, j,kτ xy n +1i , j +1 2,k = τ xy n i , j +1 2,k + ∆t [µ xy ( D y v x+ Dx v y )] n +1 2i, j +1 2,kτ xz n +1i , j ,k +1 2 = τ xz n i , j ,k +1 2 + ∆t [µ xz ( Dz v x+ Dx v z ) ]n +1 2i , j , k +1 2[

τ yz n +1i +1 2, j +1 2,k +1 2 = τ yz n i +1 2, j +1 2,k +1 2 + ∆t µ yz ( D y v z+ Dz v y ) − πft 
A = exp  Q であるから，(12)式による減衰は，f
Q = Q0f0(13)というように Q 値が周波数依存性を持つことになる．
また，P 波と S 波に異なる Q 値を与えることはできないため，シミュレーションの主たるターゲットである S 波の Q 値を与える．

i +1 2, j ,k[(λ + 2µ )D z v z+ λ D y v y + Dx v x

2.2 非弾性減衰
非弾性減衰は，媒質の各点にはたらくフィルターのようなものであり，時間領域の数値計算には導入しにくい．ここでは，Graves (1996) による近似的な方法を用いる．計算空間内の各点，各タイムステッ
プにおいて，その場所で定義されたターゲット周波数 f 0 に対するＳ波の Q 値， Qs( x, y, z ) を用いて計算

定義は，

v y n +1i +1 2, j +1 2,k = v y n −1 2 i +1 2, j +1 2,k

n +1に応じて周囲のグリッドの値の調和平均を取り，実
効物質定数（effective material parameter）としている．
実効物質定数を導入することにより，物質境界面における応力の連続条件がより正確に近似される（Randall et al., 1991 など）．

] n +1 2i +1 2, j +1 2,k +1 2

ただし，上付きの添え字は時間ステップのインデックスを示し，下付きの i, j, k などはグリッドのインデックスである．なお，bx , by , bz , µ xy , µ xz , µ yz は差分式

2.3 震源の表現
震源は，面での食い違いと等価なトラクションの変化を，対応するグリッドの応力値に付加する方法で与える（Coutant et al., 1995）．任意の向きの面での食い違いをモーメントテンソルで表示し，差分グリッドの座標系におけるモーメントテンソル要素の重みを計算する．モーメントテンソルの要素は，次のように応力値に変換できる．
モーメントテンソル各成分の偶力の腕の長さを微小長さ a と置き，これを 1 対の力（力の向き l ，大きさ ∆f l ）に変換する．腕の方向を m とする．偶力の腕のところに小さな直方体を考える．この直方体は，偶力の腕と平行で同じ長さの 4 つの辺と，これらの辺に垂直で偶力それぞれの作用点を含む面積 A の 2
つの面を持つ．この小直方体内部の応力は，偶力がある場合とない場合で，∆f l / A だけ応力が変化する．

関口春子・加瀬祐子・堀川晴央・佐竹健治・杉山雄一・Arben Pitarka

逆にいえば，微小直方体内に ∆f l / A の大きさで応力の増分を与えることにより，対応するメカニズムのモーメントテンソル震源を導入することになる．
∆τ lm = ∆f l / A = ∆M lm Aa = ∆M lm V
Coutant et al. (1995) は，差分法の計算においては，a , A の値として，それぞれ，グリッド間隔とその 2乗を用いるのが適当であることを数値実験により示した．
差分スキームの中では，通常の差分の計算により応力値が 1 タイムステップ進んで更新された後，震源にあたるグリッドに，震源に起因するトラクショ
ンの各タイムステップでの増加分を与える．震源の中心は τ ii を計算するグリッドに置く． M xy , M yz ,M zx はそれぞれ， τ xy , τ yz , τ zx を用いるので，単独のグリッドを使うと重心がずれてしまう．そこで，τ ii
を囲む 4 つのグリッドに重心が揃うように重みを付けて応力増分を配分する（Fig. 3）．また，小直方体の体積 V には，不規則格子間隔を考慮する．
τ n lm = τ n lm − ∆tM&lm (t ) V(l , m = x, y , z )(14)
2.4 領域端部の境界条件
直方体の計算領域に対し，地表面にあたる境界面には自由表面の境界条件，残りの 5 つの面には，吸収境界条件とダンピングを課している．
吸収境界は，Clayton and Enquist (1977) による．計
算領域の境界において，１）境界面に垂直に P 波速度で伝播し，境界面に垂直な方向に偏向する波の成分と，２）境界面に垂直に S 波速度で伝播し，境界面に平行に偏向する波の成分について，領域の外向きに伝播することのみを許す，つまり，領域の内側に伝播することを禁止する，という条件を課す方法である．
この方法では，境界面に斜めに入射する波の一部
は内側へ反射されてしまう．そのため，合わせて，ダンピングゾーンを境界に沿って設け，境界へ入射する波，および，境界で反射されて内側へ向かう波の振幅を抑えている．自由表面を除く境界から 10 グ
リッド程度の領域を，ダンピングゾーンとし，この中を伝播する波の振幅に１以下の定数を掛ける方法
である．ダンピング定数は境界面からの距離にしたがって変化し，最も内側で１，境界に近づくにしたがって指数関数的に小さくなるよう設定している．
地表面には，Graves (1996)の zero-stress 自由表面境
界条件を用いる． k = k 0 を自由表面とした場合を考える．自由表面に働くトラクションをゼロという条件，
τ xz k = k 0 = τ yz
= τ zz k = k 0 = 0(15)
k =k 0および，鏡像を考えた場合に自由表面を境に応力が点対称となる関係，

τ xz k = k 0 − 1 = −τ xz k = k 0 + 1
2

τ yz

k =k 0−2

1 = −τ yz
2

k =k 0+1
2

τ zz k = k 0 −1 = −τ zz k = k 0 +1

(16)

τ xz k = k 0 − 3 = −τ xz k = k 0 + 3
2τ yz

k =k 0−2

31 = −τ yz
2

k =k 0+32を(4)式に導入することにより，自由表面での差分式
はDz v z =[]−λD v +D v(λ + 2µ ) x z y y k = k 0

Dz v x = − [Dz v x + Dx v z ] k = k 0 + 1 − Dx v z k = k 0 − 1[

Dz v y = − D z v y + D y v z

2

2

2

2

(17)

] k = k 0+ 1 − D y vz k = k 0− 1となる．また，自由表面近傍の差分も 4 次精度で計算できる．
2.5 安定条件，有効周波数範囲
差分計算の安定条件は，Pitarka (1999) の中で数値実験により求められている．
∆t < 0.495h Vmax(18)
これは，個々の差分グリッドにおいて満たされるべき条件である． h および Vmax は各差分グリッドにおけるグリッド間隔と最大の波動伝播速度である．また，計算波形の有効周波数範囲は，
f < Vmin / 5h(19)で与えられる． Vmin は各差分グリッドにおける最小
波動伝播速度である．
３．計算領域と地下構造モデル
地震動シミュレーションの領域は，Fig. 4 に示した
東西方向 38 km，南北方向 58 km の領域である．この領域は，平面的には堀川ほか（2002）で作成した大阪平野の地下構造モデル領域と同じであり，深さ25 km まで延長している．
この領域の深さ 3 km までの部分には，堀川ほか（2002）モデルを当てはめた．堀川ほか（2002）の基盤深度分布を Fig. 5 に示す．それより深い部分は，
京都大学防災研究所地震予知研究センター阿武山観測所で，微小地震の震源決定に用いている水平成層
構造モデル（前田・渡辺，1984）を当てはめた．ただし，以下のような調整をしている．１）深度 3 kmの基盤岩上部における堀川ほか （2002）モデルと阿武山観測所モデルの物性値が若干異なるため，堀川ほか（2002）の物性値を採用し，阿武山観測所モデルの対応する層の物性値をこれで置き換えてつないでいる．２）堀川ほか（2002）では地形もモデル化されているが，地震動シミュレーションにおいては，

大阪平野における地震動予測の試み−動的地震破壊シナリオと 3 次元地下構造モデルによる地震動計算−

海抜 0 m 以上の部分は削除し，逆に，大阪湾内の地表面が海抜以下になっている部分に関しては，海底面部分の媒質を海抜 0 m まで延長している．３）堀川ほか（2002）のモデルで S 波速度 550 m/s 未満の部分に対応する領域では，一律に S 波速度を 550 m/sとし，P 波速度，密度，Qs 値も S 波速度 550 m/s に対応する値にしている．
計算領域内の浅い部分では，ほぼ全域にわたって，
水平方向に S 波速度の小さい堆積層が広がっている
ため，差分のグリッドは，深さ方向（ z 軸方向）にのみ不連続格子とし，東西方向（ x 軸方向）および南北方向（ y 軸方向）には一様とした．不連続格子のメリットを最大限に活かし，できるだけ効率のよい計算をするため，以下の方法を採った．まず，領域中の最も小さい S 波速度および計算したい周波数帯域から，(19)式より最小のグリッド間隔 h を決める．
h ≈ VS min / 5 f max(20)
この h を差分グリッドの z 軸方向の初めのグリッ
ド間隔 dz (1) とする．次に，堀川ほか（2002）モデルのメッシュデータの中で深さ dz (1) のところの S 波速度の最小値（ VS min ）を調べる．VS min と f max から(20)式を使って dz (2) を決める．次に深さ dz (1) + dz ( 2) のS 波速度の最小値を調べる．以上を繰り返すことによって，無駄のない不連続格子取りができる．本研
究では，1 Hz までの周波数帯をターゲットとし，最小の S 波速度が 550 m/s であることから，dz (1) = 100m とした．不連続格子を用いることにより， z 軸方向のグリッド数は不連続格子を用いない場合の 25％以下となった．
上記のように z 軸方向のグリッド間隔を決めると，最も S 波速度の大きい層での z 軸方向のグリッド間
隔は約 700 m となる．水平 2 方向へは 100 m のグリ
ッド間隔に固定しているため，非常に扁平なグリッドになる．この影響を調べるため，大阪平野 3 次元構造モデルの平野中心部付近の堆積層構造に近い，
水平成層構造モデル（Table 1）中に，点震源を置き，不連続格子を取った場合と取らない場合の地震動を計算し，結果を比較した（Fig. 6）．領域は東西・南
北・鉛直方向とも 25 km で，東から 12.5 km，南から12.5 km，深さ 15 km に，走向 0°, 東落ち 60°, す
べり角 90°の逆断層メカニズムの点震源を置いた．
点震源が全く同じ位置に来るようにするため，不連続格子の間隔は上記の方法による値からずれている．
波形は，計算領域の南端から 8 km の東西のラインに沿って計算したものである．両者は非常によく一致しており，この程度の扁平不連続格子の影響は，ごく小さいものと言える．
４．破壊シナリオの導入
加瀬ほか（2002）で計算された破壊シナリオは，断層面上に細分されたグリッド点におけるすべり破

壊のメカニズムとすべり時間関数で表現される．各グリッド点でのすべりの時刻歴を(14)式を用いて導入する．そのためには，１）動的破壊シミュレーシ
ョンの震源のグリッド点を地震動シミュレーションの差分グリッドにおいて，最も近いグリッド点に配置し，２）動的破壊シミュレーションで得られたすべり時間関数にフィルターをかけ，地震動シミュレーションの時間ステップの増分にあわせてサンプリ
ングしなおし，モーメントレートの時刻歴に置き換えて与える必要がある．Fig. 7 に断層面上のある点における破壊シナリオのすべり時間関数，これに0.1-1.0 Hz のバンドパスフィルターをかけたもの，さらに微分してすべり速度関数にしたもの，それをリサンプルしたもの，を示した．地震動シミュレーシ
ョンのグリッドの 3 軸が，水平 2 方向と鉛直方向であるのに対し，動的破壊シミュレーションのグリッ
ドは，傾斜角 60°の断層面内に 2 軸が入るように設定されているため，たとえ，2 つのシミュレーショ
ンで同じグリッド間隔を用いたとしても，震源のグリッド点をすべて同地点に置くことはできない．さらに，地震動シミュレーションでは，不連続格子を用い，岩盤内部では，動的破壊シミュレーションのグリッド間隔 200 m より大きなグリッド間隔を与えるため，位置の誤差は大きくなっている．しかし，不連続格子のグリッド間隔はターゲットとする周波数帯域の計算に十分なように与えるため，震源の位置の誤差も，この周波数帯では無視できる．
動的破壊計算も，地震動計算もどちらも差分法を用いているが，それぞれ別に計算を行っている．その理由は，動的破壊計算においては，破壊条件を判定しながら計算を行うため，波動計算のみの差分法に比べ，変数が多く計算容量がかさむ．このため，
地下構造モデルを単純にして容量を押さえている．
これに対して地震動計算では，計算結果の信頼性を高めるため，詳細な地下構造モデルを導入する必要があるからである．
５．地震動シミュレーション結果

Table 2 に，地震動シミュレーションを行った破壊シナリオの一覧を示す．すべてのシナリオを，堀川ほか（2002）の地下構造モデルを用いて計算した．3
次元地下構造モデルの影響を調べると共に，破壊シ
ナリオの違いによる地震動分布の違いを見やすくす
るため，すべてのシナリオに対して水平成層構造でも計算を行った．ここで用いた水平成層構造は，計算領域の南西端から東へ約 13 km，北へ約 42 km，基盤岩上面深度約 1.2 km 地点の堀川ほか（2002）モデ
ルの鉛直方向の地下構造を計算領域全体に広げたものである．
Fig. 8 に，すべての地震破壊シナリオについて，堀川ほか（2002）モデルおよび水平成層構造に基づく
地震動シミュレーションの，地表面最大速度値の分

関口春子・加瀬祐子・堀川晴央・佐竹健治・杉山雄一・Arben Pitarka

布を列挙した．これらを見ると，強震動域の分布は，第一義的には震源過程で決まり，第二義的にその強さを地下構造モデルが左右しているように見える．
地震破壊シナリオの違いによる強震動域の分布の影響のみを比較するためには，Fig. 8 の右側に並べている水平成層構造モデルによる結果を比較するのが良い．どの場合においても，地震動が最も強くなるのは，傾斜した断層面の延長が地表面と交わる線に沿ったゾーンである．このゾーンは，上盤側に広くなっている．このゾーンの内部での地震動の強さ分
布は，地震シナリオによって大きく変化している．
すべりの方向は純粋に傾斜方向成分のみであり，断層面はほぼ南北の走向を持っている．すべりの方向に沿って破壊が伝播する先では，断層面の延長と地表との交線を中心に，断層面に直交する方向の地震
動が破壊の前方指向性によって強められる．これが東西成分と上下成分へ分配されていると考えられる（地震シナリオ s1d801，水平成層構造の計算の東西成分と上下成分，北セグメントの断層上端付近など）．
すべりの方向と直交する方向へ破壊が伝播する先では，断層面と 45 度をなす方向を中心に前方指向性による効果が強く現われている．断層の幅方向に破壊
フロントが広がるため，この方向では破壊フロントを一斉に発した波は東西成分が打ち消し合いやすく，結局，南北成分と上下成分に現われていると考えられる（地震破壊シナリオ s1d801，水平成層構造の計
算の南北成分，南セグメントの断層の上端東側に沿った部分など）．Fig. 9 にそれぞれの場所のパーティ
クルモーションを示す．
また，破壊開始点が同じであっても，セグメント
間のステップ幅の違いにより，破壊伝播の様相が異なると，地震動分布にも大きな違いが現われた．地震破壊シナリオ s1d801 と s1d802 では，南セグメン
ト上の破壊伝播が異なり，前者ではほぼ水平に，後者では浅いほうから深いほうへ広がるように伝播している．地震動分布を比べると，前者では南セグメント沿いでの前方指向性が強く出ているが，後者ではそれが弱くなっている．地震破壊シナリオ s2d801と s2d802 のように，2 つ目のセグメントに乗り移る
場合と乗り移らない場合では，当然，地震動分布は大きく変わる．
地下構造モデルの違いによる強震動域の分布への影響は，概観では，地震破壊シナリオの違いほど大きくはないように見える．これは，小さな Q 値を持
つ堆積層中を伝播する地震波が急速に減衰し，断層面上端の近くに強震動域が集中するためと考えられる（Fig. 10）．しかし，堀川ほか（2002）モデルで計算された地動最大速度値と水平成層構造での地動最大速度値の比（3 次元／水平成層比）の分布を見ると，堀川ほか（2002）モデルの基盤岩上面形状が浮かび上がってくる（Fig. 11）．水平動の 3 次元／水平
成層比は，空間的な変動が大きく，大阪平野内部は周囲の岩盤上よりも 2〜3 倍の揺れになっている．上下動の 3 次元／水平成層比は，変動が小さいが，上町断層の西側近傍では，上下動成分の比も大きくなっている．これは，上町断層による地下の基盤岩上面の段差（高度差）により，S 波が回折を受け，上下動成分を大きくしているためと考えられる．
Fig. 12 に堀川ほか（2002）の地下構造モデルで計算した 9 つの地震破壊シナリオによる最大地動速度分布の平均と，最大値の分布を示す．地震破壊シナ
リオの違いにより，強震動のピークの位置が大きく変わるため，平均をとることにより，全体に小さな値になってしまう．地震動予測地図の作成には，複数の地震シナリオの予測値を融合させる必要がある
が，平均を取る方法は，地震動を過小評価する危険があるといえよう．各地震シナリオの予測値の最大値を取った場合には，断層面の延長と地表面との交線の東側（断層上盤側）の，幅 7〜8 km のゾーンで200〜300 cm/s，ところによっては 300 cm/s を超える
地動速度となった．この値は過去の地震で記録された地震動に比べても大きな値である．ただし，断層面上の破壊過程を推測する際に，今回の計算には取り込んでいない動的パラメータの不均質分布を考慮することにより，破壊の前方指向性が弱められ，地震動の予測値は小さくなると考えられる．
６．まとめ
活断層研究センターが取り組み始めた「3 次元地
下構造モデルと活断層情報に基づく，起震断層モデ
ルによる将来の地震と地震動の予測研究」の初のモ
デルケースとして，大阪平野における上町断層系の地震と地震動の予測を行った．
全長 45 km の上町断層系を 2 つのセグメントからなる断層としてモデル化し，すべり弱化の摩擦構成則にしたがって動的破壊過程を計算した．その結果，
2 つのセグメントが連動する地震（Mw〜7.2）と連動しない地震（Mw〜7.0）が得られた．
強震動域の分布は，第一義的には震源過程で決まり，第二義的にその強さを地下構造モデルが左右しているように見える．これは，破壊開始点の位置やセグメント間のステップ幅の違いにより，破壊伝播の様相が異なると，破壊の前方指向性によって地震動の強くなる場所が変わるためである．
断層面より発せられた地震波は，堆積層で増幅され，盆地端部や堆積層下の基盤岩上面の形状によって干渉し合う．断層面の延長と地表との交線の東側（断層の上盤側）に沿って，大きいところで 200〜
300 cm/s，ところによっては 300 cm/s を超える地動速度となった．ただし，断層面上の破壊過程を推測する際に，今回の計算には取り込んでいない動的パラメータの不均質分布を考慮することにより，破壊の前方指向性が弱められ，地震動の予測値は小さくなると考えられる
複数のシナリオによる地震動予測の分布を，平均

大阪平野における地震動予測の試み−動的地震破壊シナリオと 3 次元地下構造モデルによる地震動計算−を取る方法と最大値を取る方法で求め，両者を比較した．その結果，平均を取る方法では，せっかくのシナリオ計算の予測値が生かされないことがわかった．
今後は，今回の計算よりも高周波数部分の予測，
地震破壊シナリオのさらなる拘束，他の活断層の破壊シナリオの検討などを進めたいと考えている．
謝辞 図の作成には，Generic Mapping Tools (Wessel
and Smith, 1995) を使用した．以上，記して謝意を表
します．
文

献

Clayton, R. W. and B. Engquist (1977) Absorbing
boundary conditions for acoustic and elastic wae
equations. Bull. Seism. Soc. Am., 67, 1529-1540.
Coutant, O., J. Virieux, and A. Zollo (1995) Numericalsource implementation in a 2D finite difference
scheme for wave propagation. Bull. Seism. Soc. Am.,85, 1507-1512.
Graves, R. W., (1996) Simulating seismic wave
propagation in 3D elastic media using staggered-gridfinite differences. Bull. Seism. Soc. Am., 86,
1091-1106.
堀川晴央・水野清秀・佐竹健治・関口春子・加瀬祐子・杉山雄一・横田 裕・末廣匡基・Arben Pitarka（2002）大阪平野の 3 次元地盤構造モデルの作成．活断層・古地震研究報告，no. 2（本報告書），

291-324．
加瀬祐子・堀川晴央・関口春子・佐竹健治・杉山雄
一（2002）上町断層系の動的破壊過程の推定．
活断層・古地震研究報告，no. 2（本報告書），
325-340．
前田直樹・渡辺 晃（1984）微小地震の破壊様式−
近畿地方中北部の微小地震活動について−．地震 2, 37, 579-598.
三浦健一郎・七山 太・内海 実・杉山雄一・安原
盛明・横山芳春・北田奈緒子・竹村恵二（2002）
新淀川群列ボーリングコアの高精度解析に基づ
く上町断層の完新世活動性評価．活断層・古地
震研究報告，no. 2（本報告書），109-123．
大阪府（1999）平成 10 年度地震関係基礎調査交付
金・上町断層帯に関する調査成果報告書．
Pitarka, A. (1999) 3D elastic finite-difference modelingof seismic motion using staggered grids withnonuniform spacing. Bull. Seism. Soc. Am., 89,
54-68.
Randall, C. J., D. J. Scheibner, and P. T. Wu (1991)
Multiple borehole acoustic waveforms: synthetic
logs with beds and borehole washouts. Geophysics,
56, 1757-1769.
Wessel, P., and W. H. F. Smith (1995) New version of theGeneric Mapping Tools released. EOS Trans. AGU,
76, 329.
（受付：2002 年 8 月 23 日，受理：2002 年 9 月 17 日）

関口春子・加瀬祐子・堀川晴央・佐竹健治・杉山雄一・Arben Pitarka

第１表．不連続格子の影響を調べる為に用いた水平成層構造モデルと，各層でのグリッド間隔．
Table 1. Velocity structure and grid spacings for numerical tests to check the influence of
inhomogeneous grid spacing.

layer top depth
(km）
0.0
0.6
1.4

Vp
(km/s)
1.8
2.5
5.3

Vs
(km/s)
0.55
1.0
3.0

5.0
18.0

6.0
6.7

3.46
3.87

variable grid spacing
（km）
0.1
0.1
0.1 (shallower than 3km)
0.5 (deeper than 3km)
0.6
0.7

homogeneous grid spacing
(km)
0.1
0.1
0.1
0.1
0.1

第 2 表．地震動シミュレーションに用いた破壊シナリオ．
Table 2. Earthquake rupture scenarios on the Uemachi fault system used for ground motion prediction.

rupturescenario
s1d801s1d802s2d801s2d802s3d801s3d802s4d801s4d802s1d8f1

hypocenter
north segment / north end
north segment / north end
north segment / south end
north segment / south end
south segment / north end
south segment / north end
south segment / south end
south segment / south end
north end bottom

fault-parallel
distance (km)
0.
0.
0.
0.
0.
0.
0.
0.
-

fault-normal
distance (km)
0.4
0.8
0.4
0.8
0.4
0.8
0.4
0.8
-

Mw
7.2
7.2
7.2
6.9
7.2
7.0
7.2
7.2
7.2

大阪平野における地震動予測の試み−動的地震破壊シナリオと 3 次元地下構造モデルによる地震動計算−

第 1 図．Pitarka (1999) の 3 次元差分法コードにおける不等間隔格子のスタッガードグリッド.
Fig. 1. Schematic representation of unit cells with variable grid spacing adopted in Pitarka (1999).

第 2 図．不等間隔格子における 4 次差分オペレーターの空間増分．
Fig. 2. Spatial increments of the fourth-order finite-difference operator in staggered grids with nonuniform grid spacing.

関口春子・加瀬祐子・堀川晴央・佐竹健治・杉山雄一・Arben Pitarka

第 3 図．モーメントテンソルの各要素を表現するために与える力の分配．
Fig. 3. Forces distributed to represent each moment tensor element.

大阪平野における地震動予測の試み−動的地震破壊シナリオと 3 次元地下構造モデルによる地震動計算−

第 4 図．地震動シミュレーションの計算空間．
Fig. 4. Space for ground motion simulation.

第 5 図．堀川ほか（2002）モデルの基盤岩上面深度．
Fig. 5. Distribution of depth to the top of basement in the Horikawa et al.’s (2002)
three-dimensional velocity structure model.

関口春子・加瀬祐子・堀川晴央・佐竹健治・杉山雄一・Arben Pitarka

第 6 図．不連続格子の影響を見るための比較計算結果．東西・南北・上下それぞれ 25km の計算空間中，東端から 12.5km，南端から 12.5km，深さ 15km に，走向 0°傾斜角 60°すべり角 90°の逆断層メカニズムの点震源を置いた．波
形は，計算空間南端から 8km の地表面上の東西測線に沿ったもの．波形の上には測線沿いの最大振幅分布を示した．赤は，不連続格子を用いた場合，青は領域全体で均質な格子間隔を用いた場合の結果．
Fig. 6. Comparison of waveforms synthesized using FDM mesh with variable grid spacing (red lines) and that with homogeneous
grid spacing (blue lines). A point source with (strike, dip, rake)=(0., 60., 90) at 12.5km east and 12.5km north to thesouthwestern corner of the simulation space and at 15.0km depth. Waveforms are taken from an east-west line at 8km northto the southern edge. Above the waveforms, the peak amplitude distribution along the line is shown.

第 7 図．断層破壊のシミュレーションより得られたすべり関数(黒線)，すべり時間関数に 0.1〜1.0 Hz の範囲のバンド
パスフィルターを掛けたもの（緑線），フィルターをかけたすべり時間関数を微分してすべり速度時間関数にしたもの（青線），さらにデータ末尾をテーパー処理したもの（赤線）．(左)すべり量の大きい，南セグメント中央付近の点のすべり時間関数．
(右)最もノイズの大きい点の 1 つである，南セグメント南端部のすべり時間関数．
フィルターにより高周波ノイズは充分落とされている．
Fig. 7. Examples of slip time functions (black line) obtained by the dynamic rupture process simultation, band-pass filtered slip
time functions (green line) in 0.1-1.0 Hz, slip velocity functions (blue line) differentitated from the band-pass filtered slip
time functions, and tapered slip velocity functions (red line) which are introduced into the ground motion simulation. (Left)
A slip time function at a point near the center of the southern segment. (Right) A slip time function at a point near thesouthern edge of the southern segment where numerical noise is biggest over the fault plane. The high frequency noise iswell eliminated by the bandpass filter.

大阪平野における地震動予測の試み−動的地震破壊シナリオと 3 次元地下構造モデルによる地震動計算−

第 8 図．各破壊シナリオに対する地震動シミュレーションの最大地動速度の分布．左の列は 3 次元地下構造モデル，右の列は水平成層構造モデルの結果．上から，破壊シナリオ s1d801，s1d802，s2d801，s2d802（このページ），s3d801，
s3d802，s4d801，s4d802，s1d8f1（次ページ）（第 2 表参照）による．各破壊シナリオ，各地下構造モデルごとに左から，地動速度の東西成分，南北成分，上下成分，断層面上の最大すべり量分布，破壊時刻の分布（細線コンタ
ーは 0.5 秒間隔，太線コンターは 2 秒間隔）が示されている．地図中の緑色の点線は，断層面の地表面上への投影
位置を示す．
Fig. 8. Peak ground velocity distribution of simulated ground motions. Frames in the left-hand side are of Horikawa et al.’s (2002)
velocity structure model and those in the right-hand side are of the laterally homogeneous model. From the top to the bottom,
rupture scenario of each result is in the same order with Table 2. The velocity structure model and the rupture scenario used
are indicated above the illustrations of each result. Maximum slip distribution on the fault plane and the rupture time for the
scenario (thin contours at 0.5s interval, thick contours at 2.0 s interval) are also illustrated. Green dashed lines in the maps
indicate the fault plane projected onto the Earth surface.

関口春子・加瀬祐子・堀川晴央・佐竹健治・杉山雄一・Arben Pitarka

第 8 図．続き．
Fig. 8 (continued).

大阪平野における地震動予測の試み−動的地震破壊シナリオと 3 次元地下構造モデルによる地震動計算−

第 9 図．地震破壊シナリオ s1d801，水平成層構造モデルでの地震動シミュレーション（第 8 図左の列最上部）の，（上）領域南西端から東へ 18.6 km，北へ 46 km の地点のパーティクルモーションと，（下）領域南西端から東へ 16 km，
北へ 8 km の地点のパーティクルモーション．両地点の位置は，第 8 図左の列最上部の地図中に示されている．
Fig. 9. Particle motions of simultated ground motion assuming rupture scenario s1d801 and the laterally homogeneous velocity
structure, (top) at 18.6 km east and 46 km north to the southwestern corner of the simulation area and (bottom) at 16 km
east and 8 km north to the southwestern corner. These points are indicated in Fig. 8.

関口春子・加瀬祐子・堀川晴央・佐竹健治・杉山雄一・Arben Pitarka

第 10 図．（上）破壊シナリオ s1d801，堀川ほか（2002）の地下構造モデルによる最大地動速度の東西成分，南北成分，
上下成分，断層面上の最大すべり量分布，破壊時刻の分布（細いコンターは 0.5 秒間隔，太いコンターは 2.0 秒間隔）が示されている．地図中の緑色の点線は，断層面の地表面上への投影位置を示す．（下）上記地図中の黄色線の測線上の点で計算された最大地動速度振幅分布，地震動速度波形と，この測線直下の地下構造断面．左から東西
成分，南北成分，上下成分．
Fig. 10. (Top) Peak ground velocity distribution of simulated ground motion from rupture scenario s1d801 and the three-dimensional
velocity structure of Horikawa et al. (2002). Maximum slip distribution on the fault planes and the rupture time (thin contours
at 0.5s interval, thick contours at 2.0 s interval) for the scenario are also illustrated. Green dashed lines in the map indicate thefault planes projected onto the Earth surface. (Bottom) Simulated ground motion waveforms of EW, NS and UD components
along the lines A and B indicated in the map above. The peak ground velocity distribution along the lines and undergroundstructure beneath the lines are also illustrated.

大阪平野における地震動予測の試み−動的地震破壊シナリオと 3 次元地下構造モデルによる地震動計算−

第 11 図．3 次元地下構造モデルと水平成層構造モデルの最大地動速度の比．震源モデルは，左上が s1d801，右上が s2d801，
左下が s3d801，右下が s4d801．
Fig. 11. The distributions of ratio of peak ground velocities from Horikawa et al.’s (2002) model to those from the laterally
homogeneous velocity structure model for rupture scenarios, s1d801, s2d801, s3d803 and s4d801.

第 12 図．（左）全シナリオの最大地動速度値の平均．（右）全シナリオの最大地動速度の各点での最大値．
Fig. 12. (Left) Distribution of the average peak ground velocities of all the rupture scenarios listed in Table 2 with Horikawa et al.’s(2002) velocity model. (Right) Distribution of the maximum of all the peak ground velocities.