月刊 地球/号外 No.37,‐最近の強震動予測研究‐、海洋出版社、pp179-186, 2002

k-2 モデルによる強震動評価
久田嘉章

本論文では k-2（k-square）モデルによる理論的な強震動評価法を説明する。すべり速度関数には動的震源モデルの結果を元にした中村・宮武（2000）を用い、かつ断層面上でのすべりと破壊開始時間の分布に k-2 モデルを用いた場合、震源近傍における強震動の諸特性（長周期速度パルス波や加速度のランダム性など）を満足し、かつω-2 モデルとも整合する広帯域な強震動の評価が可能となることを示す。
１．はじめにまず簡単な震源モデルを用いて理論的強震動による評価手法における現状を説明し、k-2モデルを導入する必要性を示したい。運動力学的震源モデルによる強震動評価を行う場合、(1)式の表示定理が用いられる。
LW00U i (Y ; ω ) = ∫ ∫{µD (e n + e n )Ukjjk*ik , j}e iω ⋅t r dxdy……

(1)

ここで Y は観測点、ωは円振動数、L と W は断層の長さと幅、μは地盤のせん断剛性、D はすべり（食い違い変位）、 N j は断層面方線方向の単位ベクトル成分、 e k は滑り方向の単位ベクトル成分、 U ik* , j はグリーン関数の j 方向微分であり、i, j, k に関しては総和規約を用いている。また t r は破壊開始時間であり、次式を用いる（Hisada, 2000a）。

t r ( x, y ) =r ( x, y )+ ∆t r ( x , y )
Vr……

(2)

ここで、r は破壊開始点から破壊フロント上の点（x,y）までの距離、 V r は平均破壊伝播速度、 ∆tr は破壊開始時間の平均値からの揺らぎである。
(1)式おけるすべり量やすべり関数、破壊開始時間の分布などは当然ながら断層面上で連続関数でなければならない。これは当たり前のようであるが、意外に無視されている場合が多い。例題として図１に示す一辺 10 km の矩形の横ずれ断層の震源近傍の強震動を計算する。
図１(a)に示すように破壊伝播速度は一定とし（(2)式の ∆tr = 0）、すべり分布は図１(b)に示ように３×３小断層による不連続な矩形分布とする。一方、すべり速度関数は後述する中村・宮武(2000)による動力学的断層モデルに基づいた関数を用いる（図３(b)）。観測点は断層面の中央点から 1 km 離れた点をとり、その点から上方に 15 km にある地表面とする。媒質は等方均質な全無限弾性体（ρ=2.8 t/m3、Vp = 6 km/s、Vs = 3.5 km/s、Qp = 500、Qs =250）とし、自由表面の効果は振幅を単純に２倍して評価する。
(1)式の断層面積分を行う際、小断層ごとの積分点数（N）を 1, 4, 16, 64 点と変化させ、10 Hz までの波形及びスペクトルを計算する。図２(a)は、積分点数を変えた断層面直交成分の加速度フーリエ振幅スペクトルを示す。積分点数 N が 1 や４と少ない場合、破壊フロントの連続性が崩れる。このため大きな高振動数成分を励起しており、加速度スペクトルは振動数とともに増加し、変位スペクトルではωにほぼ逆比例して減少するω-1（ω-inverse）
モデルとなる。それに対して積分点を 16、64 と増やし、破壊フロントの連続性を保つと、加速度スペクトルはほぼ一定の振幅となる。すなわち変位スペクトルではωの２乗に逆比例するω-2（ω-square）モデルを得る。但し、図２(b)に N=64 の加速度波形を示すが、すべりの不連続分布に起因する人工的なパルス列が現れており、現実的な波形とは言えない。
以上のことから理論手法を用いて広帯域な強震動評価を行う場合、十分な積分点数（Vrによる最小波長に対し 5～10 点以上）をとらないと破壊フロントが不連続になり、人工的な高振動数波を励起してしまうこと、一方、破壊フロントを連続させた場合でもの小断層によ1月刊 地球/号外 No.37,‐最近の強震動予測研究‐、海洋出版社、pp179-186, 2002
るすべり分布では、すべりの不連続に起因する人工的なパルス列を生じさせること、などが分かった。従って、すべり分布は断層面上で連続関数でとするのが理想的である。それではすべりや破壊開始時間を断層面上で分布させる場合、どのような関数を用いたら現実に近い震源モデルとなるのであろうか？ 経験的に震源スペクトルはω-2 モデルになることが知られているため、ω-2 モデルとも整合する必要がある。すべり分布に関しては、これまでフラクタル理論などを用いて様々な研究が行われているが、破壊開始時間やすべり速度関数まで含め、強震動評価を行っている研究はその諸についたばかりである（例えば、圓・河野,2001 など）。k-2（k-square）モデルは、すべりや破壊開始時間に関し断層面上で自己相似に基づく分布関数を与え、強震動評価法を行うための震源モデルである。このモデルはHerrero and Bernard(1994)により提案され、Hisada(2000a, 2001)により物理的により妥当なモデルに改良された。このモデルに関しての詳細はすでに Hisada(2000a, 2001)などがあるため、ここでは要点のみを解説したい。

(a)(b)

図１：（a）一定の破壊伝播速度（ V r = 2.8 km/s）による滑らかな破壊開始時間の分布、および、（b）３×３小断層による不連続な矩形のすべり分布（平均すべり量は 1 m）

Fourier Amplitude (gal*sec)

1000
N= 1
N= 4
N=16 100
N=6410

1
0.1

1(a)

10(b)

0.1
frequency (Hz)

図２：(a)小断層の積分点数 N と断層直交成分の加速度フーリエ振幅スペクトルの関係、(b)N=64 の断層直交成分の加速度波形

2月刊 地球/号外 No.37,‐最近の強震動予測研究‐、海洋出版社、pp179-186, 2002
２．k-2 モデルとω-2 モデル
k-2 モデルとω-2 モデルの関係を示すため、(1)式から震源スペクトルを導く。すなわち地盤を均質な等方弾性体と仮定し、グリーン関数の遠方項のみを用い、S 波・P 波に分離すると次式を得る。
Ri……… (3)U i (Y ; ω ) =M (ω )4πρrβ 3
ここで、Ri はラディエーションパターン、ρは密度、βはせん断波速度、ｒは震源から観測点までの距離である。Ｍは震源スペクトルであり、次式で表せる。
WM (ω ) = µ ∫0

∫ D( x, y ) F ( x, y;ω ) exp{iω (t β + t )}dxdy
L0r………

(4)

ここで D は最終すべり分布、F はすべり速度関数（最終すべり量を１とした時の速度関数）、tβ は破壊先端部のシグナルが観測点 Y まで到達する時間である。
震源スペクトルはω-2 モデルになることが経験的に知られている。そこで(4)式のすべり、すべり速度関数、及び破壊開始時間の分布のそれぞれが、どのようにω-2 の震源スペクトルに寄与しているのかを示し、k-2 モデルとの関係を説明する。
1) すべり速度関数とフーリエ振幅スペクトル
ここでは様々なすべり速度関数とその振幅スペクトルを調べる。代表的なすべり速度関数
として、表１にあるδ関数、矩形関数、三角形関数、指数関数（２種類）、Hisada(2000a、2001)、及び中村・宮武(2000)による関数を用いる。図３にその時刻歴波形を、図４にすべり加速度のフーリエ振幅スペクトルを、それぞれ示す。
δ関数は Herrero and Bernard(1994)によるオリジナルの k-2 モデルに用いられたすべり速度関数である。物理的に妥当性はなく、速度スペクトルの振幅１となり、加速度スペクトルは図４(b)に示すようにωに比例して増大し、非常に大きな高振動数成分を生じる。なおBernard 他（1996）ではδ関数の代わりに、空間スケール依存の継続時間を持つすべり速度関数を導入しているが、やはり速度の振幅スペクトルは１と仮定している。
一方、矩形関数と三角形関数の加速度スペクトル（ここではτ=1 sec とする）は、図４(a)に示すようにコーナー周波数（矩形では約 0.5/τ、三角形では約 1/τ）以上の振動数で、前者は一定に、後者はωに比例して減少する。矩形関数の大きな高振動数成分は、図３(a)に示すように立ち上がり時と終了時におけるすべり速度の不連続により生じている。
同様に指数関数１と指数関数２のスペクトルは、図４(a)に示すように高振動数にて前者は一定、後者はωに比例して減少している。図３(a)に示すように指数関数１には立ち上がりに不連続であるため、大きな高振動数成分を生じている。
一方、中村・宮武(2000) や Hisada(2000a、2001)の関数は、動力学震源モデルの解析結果をもとに構築したすべり速度関数であり、図３(b)に示すように急激な立ち上がりと滑らかな尾部を持つ特徴がある。中村・宮武(2000)の関数の主なパラメータは、立ち上がり部分の形状とスペクトルの高振動数におけるコーナー振動数をコントロールする fmax、すべりの継続時間の ts、最終変位量に対する最大速度の比の Rvd である。この関数は時刻歴で与えられており、周波数領域の関数形は数値的にフーリエ変換する必要がある。一方、久田（2001）の関数は、様々な継続時間を持つ三角形の重ね合わせですべり速度関数を表現している。主なパラメータは、fmax（最小の継続時間の逆数）、要素三角形の継続時間の比 Trと面積の比 Ar、三角形の重ね合わせ数 n、である。三角形要素の重ね合わせであるため、周波数領域での解析解が使用可能である。
図３(b)と図４(b)には、Hisada(2001)による M7 地震のパラメータ設定例を参考にして、中村・宮武(2000)では fmax= 5 Hz、ts= 3.2 秒、Rvd =1.12 を、Hisada(2001)では fmax=5 Hz、Tr=1.74、Ar=1.4、n=6（継続時間は 3.2 秒）とした例を示している。両図より二つの関数は良く似た特性を示しており、特に加速度スペクトルでは高振動数で fmax までほぼ一定の振幅を示している。fmax より高振動数では Hisada(2001)の方が急激な振幅の低下を示しているが、これは最小の継続時間の要素三角形におけるフーリエ振幅の落ち込みに対応する。但3月刊 地球/号外 No.37,‐最近の強震動予測研究‐、海洋出版社、pp179-186, 2002
し、実際の震源モデルでは、すべり速度関数が断層面上で変化するため、このような急激な振幅の落ち込みは回避可能である。
以下の解析では断層面上で一定のすべり速度関数を用いるため、中村・宮武(2000)の関数を使用する。

表１：代表的なすべり速度関数の時刻歴波形とフーリエスペクトル
時間領域
周波数領域

種類
δ関数

f (t ) = δ (t )
1 τ ,f (t ) =  0,矩形関数

F (ω ) = 1(0 ≤ t ≤ τ )(τ < t )F (ω ) =sin( χ 1 )χ1exp(iχ 1 ) ,三角形関数

 4t τ 2 ,(0 ≤ t ≤ τ 2)
2f (t ) = 4(τ − t ) τ , (τ 2 ≤ t < τ )0,(τ < t )

指数関数1f (t ) = exp(−t / τ ) τF (ω ) = 1 (1 − iωτ )指数関数2f (t ) = t ⋅ exp(−t / τ ) τ 2

F (ω ) = 1 (1 − iωτ )三角形の重ね合わせ（Hisada, 2000a, 2001）
中村・宮武 中村・宮武(2000)

0.5

τ=1 sec

1.5

12
3time (sec))2

三角形の重ね合わせによる解析解（Hisada, 2000a, 2001）
中村・宮武(2000)の数値解

0.002 sin( χ 2 )ωτF (ω ) = exp(iχ 2 ) , ( χ 2 =)4
 χ2

slip velocity

slip velocity

1.0

rectangulartriangle
exponential 1
exponential 2

ωτ2

久田

2.0

(χ1 =

1.2
1.0
0.8
0.6
0.4
0.20.04

HisadaNakamura & Miyatake
fmax=5 Hz
τ=1 sec0

1

2
3
time (sec)4

図３：表１の各種すべり速度関数の時刻歴波形

τ=1 sec

10.01

0.1
Rectangular
Triangle
exponential 1
exponential 2

1
0.1
0.01

frequency (Hz)10

Fourier Amplitude

Fourier Amplitude

10

fmax=5 Hz
τ=1 sec

10
10.01

0.1

110

δ-function 0.1
HisadaNakamura & Miyatake

0.01frequency (Hz)

図４：表１の各種すべり速度関数の加速度フーリエ振幅スペクトル4

月刊 地球/号外 No.37,‐最近の強震動予測研究‐、海洋出版社、pp179-186, 2002
以上のことから、実際の震源に近いと思われるすべり速度関数は図３(b)、図４(b)に示す特性を持ち、その加速度スペクトルの振幅は fmax までの高振動数側でほぼ一定、すなわち速度スペクトルの振幅は振動数ωの逆数にほぼ比例すると考えられる（ω-1 モデル）。
2) すべり分布と k-2 モデルの導入
前述したように図１(b)に示す矩形のすべり分布は実際には存在しないばかりか、高振動数での強震動を評価する際、すべりの不連続による人工的なパルス波を発生する、などの問題も生じる。そこで、すべり分布の連続性を保証し、現実的なすべり分布を発生させる方法
として k-2 モデルが提案されている。k-2 モデルでは空間（L、W 方向）に関するすべり分布のフーリエ振幅スペクトルが空間波数 k の２乗に逆比例すると仮定する。近年の震源逆解析で得られたすべり分布も、その適用波数の範囲で k-2 モデルの妥当性が確認されている。位相スペクトルに関しては、Herrero and Bernard(1994)によるオリジナルの k-2 モデルではランダム位相を仮定している。この場合、断層境界で大きなすべりが発生したり、変位の不連続が生じたり、またアスペリティーの位置や大きさ、個数などを調整することができないなどの問題が生じる。そこで Hisada(2001)では小断層サイズで決まるナイキスト波数（＝1/2⊿L、1/2⊿W）以下の波数では、震源逆解析の結果などで決まるすべり分布を決定論的に用い、一方、ナイキスト波数以上では k-2 モデルを導入するハイブリッドは手法を提案した。
k-2 モデルとして、Hisada(2001)では２次元 Butterworth 関数の次式を用いている。
NMD( x, y ) = ∑∑n =1 m =1yx
 cos 2π ⋅ m ⋅ + θ mn  cos 2π ⋅ n ⋅ + θ n LW 
1 + (m + n )
D2……

(5)

2 2

ここで、 D は平均すべり量、ｍとｎはそれぞれ L、W 方向の無次元化波数、θmn、θn は 0 から２πまでのランダム数である。
図１のすべり分布をもとに Hisada(2001)によるハイブリッド手法を用いて、k-2 モデルの構築法を紹介する。まず図１のすべり分布で断層周辺のすべりを 0 と拘束し、２次元 Cubic
Spline 補間を用い、滑らかで連続したすべり分布を生成する（図５(a)）
。上述したように、
この分布は小断層サイズの解像度であるナイキスト波数（k=0.15 1/km）まで有効である。
それ以上の高波数には(5)式の k-2 モデルによる自己相似則を仮定してすべり分布を発生させ、最終的に先の滑らかなすべり分布と重ね合わせる（図５(b)）。ここで(5)式の重ね合わせ数 M、N には、それぞれ Lf / V r 、 Wf / V r を用いる（ f は対象とする最大周波数で、ここでは 10 Hz のため、M=N=36 とした）。

(a)

(b)

図５：(a)図１のすべり分布から Cubic Spline 補間による滑らかなすべり分布の生成、及び、(b)(5)式を重ねた k-2 モデル

5

月刊 地球/号外 No.37,‐最近の強震動予測研究‐、海洋出版社、pp179-186, 2002
図５(b)の k-2 すべり分布を用い、図２で計算した同じ観測点（断層面中央 1 km 離れた点から上方に 15 km の地表面）で(1)式より強震動を計算する。破壊伝播速度は一定とし、破壊開始時間は図１(a)の滑らかな分布とする。すべり速度関数はオリジナルの k-2 モデルであるδ関数と、より現実的な図３(b)に示した中村・宮武（2000）の関数を用いる。
図６に２つのすべり速度を用いた断層直交成分の加速度フーリエ振幅スペクトルを示す。
δ関数を用いたオリジナルの k-2 モデルは加速度スペクトルも高振動数で一定となるω-2モデルとなっている。それに対して中村・宮武（2000）の関数を用いた場合、加速度スペクトルが高振動数でωと共に減少するω-3（ω-cube）モデルとなってしまう。従って、現実的なすべり速度関数を用いた場合、すべり分布に k-2 モデルを導入しただけでは高振動数を十分に励起できないことが分かる。

Fourier Amplitude (gal*s)1001010.01

0.1

110

δ-function
Nakamura & Miyatake

0.1
frequency (Hz)

図６：破壊伝播速度を一定、かつ図５(a)の k-2 すべりモデルを用いた震源近傍の加速度フーリエ振幅スペクトル（すべり速度関数にはδ関数と中村・宮武の関数を使用）

3) 破壊開始時間の分布と k-2 モデルの導入
図６の例より、すべり分布に k-2 モデルを用いても、すべり速度関数に現実的なHisada(2000a, 2001)や中村・宮武（2000）の関数を用いた場合、震源スペクトルはω-3 モ
デルになることが分かった。この理由は破壊伝播速度を一定とし、図１(a)に示した滑らかな破壊開始時間を仮定したためである。現実の震源でも破壊フロントは複雑に変化し、高振動数が励起されていると考えられている。このため Hisada(2000a, 2001)は、(2)式の破壊開始時間の揺らぎである ∆tr に、k-2 モデルを導入することを提案した。
∆tr に関する k-2 分布の与え方は現状では詳しくは分からないが、すべり分布と破壊伝播速度には相関関係があると報告されている（例えば Day, 1982）。そこで、すべり分布の大きい場所では破壊開始時間も小さいとし、図５のすべりの k-2 分布に逆比例する ∆tr の k-2分布与える。(5)式の D の代わりに 0.8 秒を与えた破壊開始時間の分布を図７に示す。
図５のすべり分布、及び図７の破壊開始時間の k-2 分布、かつ図３(b)の中村・宮武のすべり速度関数を用いた震源近傍の加速度フーリエ振幅スペクトルと、速度・加速度波形と、を図８、図９に示す。加速度スペクトルは 0.4 Hz 程度の第一コーナー振動数から fmax である 5 Hz までがほぼフラットな振幅を持つω-2 モデルとなっている。但し、断層直交成分は断層平行成分に比べ、低振動数だけでなく高振動数でも数倍大きな振幅を示しており、約 1 Hz 以上では両成分に差が無くなるという経験的モデルと矛盾している。これを改善するために、周波数依存のラディエーションパターン導入などの検討が必要である。
図９の波形では、速度におけるパルス波や、加速度におけるランダム波など震源近傍の強震動波形の諸特性が現れている。一方、S 波の立ち上がり非常に大きな振幅が見られるが、中村・宮武（2000）によれば、破壊開始時のすべり速度関数は比較的滑らかな形状である
6

月刊 地球/号外 No.37,‐最近の強震動予測研究‐、海洋出版社、pp179-186, 2002

Fourier Amplitude (gal*s)

知られている。従って、今後は破壊開始点からの距離に依存したすべり速度関数を用いるべきだと思われる。

100
Fault Normal
Fault Parallel
Uo Down
1010.010.1

1100.1
frequency (Hz)

図７：破壊開始時間の k-2 分布モデル(a)

図８：すべりと破壊開始時間に k-2 分布モデルを用いた震源近傍の加速度フーリエ振幅スペクトル（すべり速度関数には中村・宮武の関数を使用）(b)

図９：すべりと破壊開始時間に k-2 分布モデルを用いた震源近傍の速度及び加速度波形（断層直交成分、すべり速度関数には中村・宮武の関数を使用）

３．おわりに（強震動評価から強震動予測へ）以上のことから、すべり速度関数に動的震源モデルを元にした中村・宮武（2000）やHisada(2000a, 2001)の関数を用い、断層面上でのすべりと破壊開始時間の分布に k-2 モデ
ルを用いれば、震源近傍における強震動の諸特性（長周期速度パルス波や加速度のランダム性など）を満足し、かつ経験的震源モデルであるω-2 モデルとも整合する広帯域な強震動評価モデルが構築できることが分かった。ちなみに Hisada(2000a, 2001)では、すべり速度関数のω-1 モデルと、すべりと破壊開始時間の k-2 分布によるω-1 モデルとの合積からω

7

月刊 地球/号外 No.37,‐最近の強震動予測研究‐、海洋出版社、pp179-186, 2002
-2 モデルが構成されているため、このモデルをω-12（ω-inverse-squared）モデルと呼んでいる。
本モデルを強震動予測に応用する場合の利点は、理論的手法による強震動評価を任意の短周期領域までの拡張を可能にすることである。すなわち、現在使用可能な震源逆解析から得られている震源パラメータや、アスペリティーサイズなどに関する各種震源スケーリング則は長周期（長波長）に限定されているが、k-2 モデルによる自己相似則を用いることで、より短周期（短波長）へのパラメータの拡張が可能になる（例として図１(b)から図５(b) のすべり分布の構築）。特に M8 クラスの巨大地震を対象とした場合、逆解析で用いられる小断層サイズは 20～30 km と大きく、使用可能な震源パラメータの有効周期範囲は 3～4 秒程度以上の長周期である。従って規模の大きな地震を対象とするほど本手法はより有効になる。
他方、本手法の欠点は、破壊フロントやすべり分布の連続性を保証するため、高振動数では膨大な数の断層積分点が必要となることである（用いる積分則によるが、対象波長に対して５～１０点は必要）。特に現実的な多層地盤のグリーン関数を用いた場合、計算時間も飛躍的に増大する。従って理論的にはどんな高振動数でも計算は可能であるが、実用上は 1～2 Hz 程度が適用限界になると思われる。それ以上の高振動数では経験的手法を用いるなど、ハイブリッドは手法を用いた方が賢明である（例えば、Hisada(2000b)による 1985 年メキシコ地震の震源近傍の強震動評価）。一方、このモデルで用いられているパラメータの構築方法（すべり速度関数や破壊開始時間の断層面上での分布のさせ方、など）も、まだ確立された訳ではなく、試行錯誤の段階である。今後は様々な地震の観測記録に適用することで、強震動予測に向けて本モデルの信頼性を高めていきたいと思っている。
謝辞
本論文では、すべり速度関数の計算に中村・宮武(2000)によるプログラムを使用させて頂きました。
参考文献[1] Bernard, P., A. Herrero, and C, Berge (1996): Modeling directivity of heterogeneous
earthquake ruptures, Bull. Seism. Soc. Am., 86, 1149-1160.
[2] Herrero, A. and P. Bernard (1994): A kinematic self-similar rupture process forearthquake, Bull. Seism. Soc. Am., 84, 1216-1228.
[3] Hisada, Y. (2000a): A theotetical ω -squared model considering the spatial
variation in slip and rupture velocity, Bull. Seism. Soc. Am., 90, 387-400.
[4] Hisada, Y. (2001): A theotetical ω -squared model considering the spatialvariation in slip and rupture velocity. Part 2: Case for a Two-Dimensional Source
Model, Bull. Seism. Soc. Am., 91, 651-666.
[5] 中村洋光、宮武 隆(2000): 断層近傍強震動シミュレーションのための滑り速度関数の近似式、地震、第５３巻、第１号、1-11[6] 圓 幸史朗・河野允宏(2001): 断層面の不均質フラクタル性を考慮した震源のモデル化と地震動、日本建築学会構造系論文集、第 543 号、51-60[7] Day, S. M. (1982): Three-dimensional simulation of spontaneous rupture: The effectof nonuniform prestress, Bull. Seism. Soc. Am., 72, 1991-1902.
[8] Hisada, Y. (2000b): A hybrid method for predinctiong strong ground motion at broadfrequencies near M8 earthquakes in subduction zone, Proc 12 World Conf. on Earthq.
Engng., Auckland, New Zealand, CD-ROM.

8