WGEN

模擬地震動の作成―WGEN
プログラムWGE N（S i m ul at ed W ave G enerat i on）は，地震のマグニチュードと震央距離が 与えられたとき ，減衰定数 5％の大 崎スペクト ルを目標とし ，これに適 合する模擬 地震動の加 速度時刻歴を作成するサブルーチン副プログラムである．

WGE N（模擬地震動の作成）【目 的】
地震のマグニチュ−ドと震央距離が与えられたとき，減衰定数5％の大崎スペクトルに適合 するよう，解放基盤面上における模擬地震動の加速度時刻歴を作成する．
【使用法】（1）接続方法
CA L L WGE N（E M，R，NN，IR，A CC，ND，DT，A MA X，V MA X， MXCYCL ，
E RR，V W1 ， V W2 ）引 数型

プログラムを呼ぶときの内容

EM

R

地震のマグニチュード

不 変

R

R

震央距離（単位 km）

不 変

NN

I

加速度時刻歴のデータ数

不 変

IR

I

乱数を発生する関数副プログラムに初期値として与える整数

R
ACC

1次元配列

プログラムから戻ったときの内容

変化している

何も入れなくてよい

模擬地震動の加速度時刻歴（単位 gal）

（ND）
ND

I

配列 ACC，VW1，VW2 の整合寸法
ND≧NN

不 変

DT

R

何も入れなくてよい

加速度時刻歴の時間間隔（単位 sec）

AMAX

R

何も入れなくてよい

同 上 最大加速度（単位 gal）

VMAX

R

何も入れなくてよい

同 上 最大速度（単位 kine）

MXCYCLI

逐次計算の最大繰り返し回数

実際の繰り返し回数

ERR

R

許容誤差（無次元小数）

実際の誤差（無次元小数）

-1-

WGEN

R
VW1

1次元配列

何も入れなくてよい（作業領域）

何も入れなくてよい（作業領域）

（ND）
R
VW2

1次元配列
（ND）

（2）必要なサブル−チンおよび関数副プログラム
ENVL OHS P VELK F AS T C R AC IAC C ER ES R AND（3）注意事項
NN は，4096を越えない2の累乗数でなければならない．
【計算法】
プログラムの構成は，およそ次のようになっている．
ａ）引数 引数のうち ， NN は作成 する模擬地 震動時刻歴 におけるデ −タ数であ って，4096 を越えない２の 累乗数とする ．目的にもよ るが，地震の マグニチュ− ドに応じて， およそ次の よ
うな NN が適当であろう．２の累乗数でないデ −タ個数 NN を希望する場合は，後続のゼロを 付け加えた個数を，２の累乗数とすればよい．
マグニチュード M

678

NN51210242048

継続時間 Td ( sec )

12.2

24.9

50.8

時間間隔 ∆ t ( sec )

0.024

0.024

0.025

引数 IR は，使 用者が自分のラ イブラリから選 び，プログラム 中で引用する 0〜1 間の一様 分
布乱数発生 用関数副 プログラ ム（ここで は，仮に その名前 を RAND として いる）を 起動する ため， RAND ( IR ) の形で与え る初期値である ． IR は単なる整数 ではなく奇数， 非負の整数ある いは負 の整数など，プログラムによって，制限の設けられている場合が多い．初期値 IR の値によって，同じ与条件に対して，異なった地震動が作成できる．
引数 ERR は，繰り返し計算を打ち切るための許容誤差であり，実用的な見地からいって，5％より 小さくする必要はない ． ERR = 0.05 の場 合，計算は大体４〜６ 回の繰り返しで 収束するので，通常の場合，繰 り返し計算の 最大回数を与 える引数 MXCYCL は， 10程度にして おけば十分で あろう．
ｂ）位相差分 サブルーチンENVLを呼び，継続時間 Td を求めて，時刻歴の時間間隔 ∆ t を計 算し，包絡曲線 E ( x ) を，その両端つまり時刻０ および Td に対応する点（ x = 0 および x = 1 ）を 含む33点，いい換えれば ０〜１間のx軸を32等 分する点で定める．次 に，包絡曲線の値 E ( x ) を ，０から順次加算 した上，累積 確率密度分布 EE ( x ) を求め，０〜 １間に一様に分 布する任意の 乱数を p としたとき， EE ( x ) = p となる x の値を求める．
このこ とを， ( NN / 2 − 2 ) 個の 乱数につ いて繰 り返せ ば， ０ 〜１間 に分布 する ( NN / 2 − 2 ) 個の x の値の標本が得ら れるので，これらを０ 〜 ( − 2π ) 間の分布に換算す れば，１と組の位相 差分 PDIF すなわち ∆φ k （ k = 1 , 2 , L , NN / 2 − 2 ）が定まる．このような 位相差分の分布形は， 包絡曲線の形と 似たものにな り，したがっ て，作成され る模擬地震動 は，近似的に サブルーチ ン
ENVLによって与えられた包絡形を有するものになると期待してよい．
ｃ） 位 相 角 位 相 差分 が 決 ま れ ば ， 位相 角 PHI は ， PHI ( 1 ) は 無視 し て も よ く ， PHI ( 2 ) を

-2-

WGEN

φ 2 = 0 と仮定すれば，以下順次に
φ k +1 = φ k + ∆ φ kk = 2 , 3 , L , NN / 2 − 2で計算される．
ｄ） 目標スペクトルとフ−リエ第1 次近 似各成分の振動周期
T k = Td / k

NFOLD は，折曲げ点の番号である．ここではまず，

k = 2 , 3 , L , NFOLDを計算する． T1 は実際は∞であるが，プログラム上は一応 T1 = 2 T2 と仮定している．
目標スペク トルである 大崎スペク トルは，外 挿できない こともない が，本来は 周期 0. 02〜 2s ec の間で 定義さ れたも のであ るから ，周期 がこの範 囲に入 る最小 の次数 KMIN と最 大の次 数
KMAX ，範囲内に入る周期の個数 NE を求める．
次に，サブル ーチンOHS P を呼 んで，各周期点 Tk における減 衰定数０％の速 度応答スペク トル ( S v ) kh = 0 を求め（ 上記の範囲外 は０とする）， これとフ−リ エ振幅スペクト ルが，近似的 に等しいこ とを利用 して，振 幅 Fk の第 1次近似 値とする ．ついで ，OHS P を再 び呼んで ，以後 の目標スペクトルとなる減衰定数５％の速度応答スペクトル ( S v ) ktarget を求めておく．
ｅ）繰り返し計算

まず１回目の繰り返し計算（ NCYCL = 1 ）に入り， Fk の第1次近似値と位相角 φ k から，複素フ−リ エ係数 C k を求め，これ をフ−リエ逆変換して ，加速度時刻歴 ACC の 第１次近似とする．
ｆ）基線補 正 ここで ，サブルー チンC R AC を呼 んで，加速 度時刻歴の 基線補正を しておく のが望ましい．補正後の時刻歴をフ−リエ変換し，改めて振幅 Fk の第1次近似値を求める．
ｇ）フ −リ エ振 幅の 修正

基線 補正 した 第１ 次近似 の加 速度 時刻 歴につ いて ，サ ブル ーチ ン

ER ES を 呼 び， 減 衰定 数 ５％ の 速度 応 答ス ペ クト ル ( S v ) k を 求 めて み ると ， 目標 ス ペク ト ル( S v ) ktarget とは一致しない．そこで，両者の比r ( k ) = ( S k ) ktarget / ( S v ) kを計算して ，振幅 Fk を Fk → r ( k ) ⋅ Fk と修正し ，これを Fk の第２次近 似値として ，２回目の 繰り返し計算（ NCYCL = 2 ）に入る．
ｈ）収束誤差 以下，同 様な計算を繰 り返し，繰 り返し回数が 所定の値（ MXCYCL ）に達し た
場合 ，あるいは誤差を表わす比率 r ( k ) の２乗平均が，許容誤差 ERR より小さくなったところで，計算を打ち切る．この 時点で，配列 ACC に格納 されているデータが， 求める模擬地震動の加 速度時刻歴である．
プログラム の算法は，以 上のとおりで あるが，ここ で，このプロ グラムに関連 して，なお 若干の事項を付け加えておく．
i ）引数 IR に 与える数値によ り，地震の同一 のマグニチュ− ド，同一の震央 距離に対して， 異なった地震動 が作成できる ．このプログ ラムによれば ，作成された 地震動の最大 加速度・最 大速度は， IR の値にかか わらず，それぞれ別掲 のプログラム OHAC に よっ て計算される最大加 速度， VELKによって計算される最大速度と，それほどかけ離れたものとはならない筈であるが，できればごく 近い値になる ことが望まし い．また，客 観的な評価は 難しいが，で きれば，作 成された地震動の形が，自然の地震動らしさをもっていることも，模擬地震動に要求される1つ の条件であろう．
実際には，引数 IR に適宜の値を与えながら， 何回かの試行を繰り返し，これらの条件に合 ったものを，目指 していく訳であ る．よほど運が 悪くないかぎり ，おそらく10回 以内の試行で ，一応満足のいく結果に，到達できるだろう．
i i ）前にｄ）項で述べたとおり， このプログラムの算法では，ま ず各周期点 T ｋ における減衰 定-3-WGEN

数0％の速 度応答ス ペクトル S v hk = 0 を求 め，これと フ−リエ 振幅スペク トルとの 近似性に 着目して，振幅 F k の第1次近似値を定める．つい で，以後の繰り返し計算では，減衰定数5％の 速度応答スペクトル S v ktarget を目標スペクトルとして， F の逐次近似値を求めていく．
しかし， F k の第 1次近似値を定め るだけの目的で， まず S v hk = 0 を求め る段階を省略して ，target代わりに，最 初から減衰定 数 5％ の S v Kを， F k の第１次近似 値としても， 結果的には ，ほとんど差異を生じない．
【プログラム】
CCCCCC* * * * * * * * * * * * * * * * * * * * * * * *SUBROUTINE FOR SIMULATED WAVE GENERATION* * * * * * * * * * * * * * * * * * * * * * * *CODED BY Y.OHSAKISUBROUTINE WGEN(EM,R,NN,IR,ACC,ND,DT,AMAX,VMAX,MXCYCL,ERR,VW1,*VW2)CCOMPLEXC(4096)
DIMENSION ACC(ND),VW1(ND),VW2(ND)
DIMENSION E(33),X(33),EE(33)DIMENSION PDIF(2046),PHI(2049),F(2049),T(2049),SV(2049),H(1),*RES(2049,1),RR(2049)
PARAMETER (PI2=6.283185)
DATA
DX/0.03125/,H0/0./,H/0.05/
CCC

PHASE DIFFERENCESCALL ENVL(EM,TB,TC,TD,33,E,33)
DT=TD/REAL(NN)
NN2=NN/2
NFOLD=NN2+1
X(1)=0.
EE(1)=0.
DO 110 M=2,33
X(M)=REAL(M-1)*DX
EE(M)=EE(M-1)+E(M)110 CONTINUE
DO 120 M=2,33
EE(M)=EE(M)/EE(33)120 CONTINUE
DO 150 K=1,NN2-2
P=RAND(IR)DO 130 J=2,33
IF(P.LE.EE(J)) GO TO 140
130 CONTINUE
140 PDIF(K)=-(X(J-1)+(P-EE(J-1))/(EE(J)-EE(J-1))*DX)*PI2
150 CONTINUECCC

PHASE ANGLESPHI(2)=0.
DO 160 K=1,NN2-2
PHI(K+2)=AMOD(PHI(K+1)+PDIF(K),PI2)
160 CONTINUE

CCC

TARGET SPECTRUM AND FIRST FOURIER APPROXIMATION
T(1)=TD*2.
DO 170 K=2,NFOLD
T(K)=TD/REAL(K-1)

-4-

WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN
WGEN

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152WGEN

170 CONTINUE
DO 180 K=2,NFOLD
IF(T(K).LE.2.) GO TO 190
180 CONTINUE
190 KMIN=K
DO 200 K=KMIN,NFOLD
IF(T(K).LT.0.02) GO TO 210
200 CONTINUEKMAX=NFOLDGO TO 220
210 KMAX=K-1
220 NE=KMAX-KMIN+1
DO 230 K=2,NFOLDCALL OHSP(EM,R,H0,T(K),SV0,K-2)
IF(K.LT.KMIN.OR.K.GT.KMAX) SV0=0.
F(K)=SV0/TD
230 CONTINUE
DO 240 K=2,NFOLDCALL OHSP(EM,R,H(1),T(K),SV(K),K-2)
IF(K.LT.KMIN.OR.K.GT.KMAX) SV(K)=0.
240 CONTINUECCCITERATIVE COMPUTATION
ENN=1./REAL(NN)NCYCL=0250 NCYCL=NCYCL+1C(1)=(0.,0.)
DO 260 K=2,NN2C(K)=F(K)*CMPLX(COS(PHI(K)),SIN(PHI(K)))
C(NN+2-K)=CONJG(C(K))
260 CONTINUE
C(NFOLD)=F(NFOLD)*(1.,0.)
CALL FAST(NN,C,4096,+1)
AMAX=0.
DO 270 M=1,NN
ACC(M)=REAL(C(M))
AMAX=AMAX1(AMAX,ABS(ACC(M)))
270 CONTINUE

CCC

BASE LINE CORRECTIONCALL CRAC(DT,NN,AMAX,ACC,ND,VW1,VW2)
DO 280 M=1,NN
C(M)=CMPLX(ACC(M),0.)
280 CONTINUECALL FAST(NN,C,4096,-1)
DO 290 K=2,NFOLD
F(K)=CABS(C(K))*ENN
290 CONTINUECCC

MODIFICATION OF FOURIER AMPLITUDESCALL ERES(1,H,1,NFOLD,T,2049,DT,NN,ACC,ND,2,VMAX,RES)
DO 300 K=2,NFOLDRR(K)=SV(K)/RES(K,1)
F(K)=F(K)*RR(K)300 CONTINUECCC

ERROR FOR CONVERGENCE
EPS=0.
DO 310 K=KMIN,KMAX-5-WGEN 53
WGEN 54WGEN 55
WGEN 56
WGEN 57
WGEN 58
WGEN 59
WGEN 60
WGEN 61
WGEN 62
WGEN 63
WGEN 64WGEN 65
WGEN 66
WGEN 67
WGEN 68
WGEN 69
WGEN 70
WGEN 71
WGEN 72
WGEN 73
WGEN 74
WGEN 75
WGEN 76
WGEN 77
WGEN 78
WGEN 79
WGEN 80
WGEN 81
WGEN 82
WGEN 83
WGEN 84
WGEN 85
WGEN 86
WGEN 87
WGEN 88
WGEN 89
WGEN 90
WGEN 91
WGEN 92
WGEN 93
WGEN 94
WGEN 95
WGEN 96
WGEN 97
WGEN 98
WGEN 99
WGEN 100
WGEN 101
WGEN 102
WGEN 103
WGEN 104
WGEN 105
WGEN 106WGEN 107
WGEN 108
WGEN 109
WGEN 110WGEN 111
WGEN 112
WGEN 113
WGEN 114
WGEN 115WGEN

EPS=EPS+(1.-RR(K))**2310 CONTINUEEPS=SQRT(EPS/REAL(NE))

WGEN 116
WGEN 117WGEN 118WGEN 119
WGEN 120
WGEN 121
WGEN 122
WGEN 123
WGEN 124
WGEN 125
WGEN 126

C
IF(EPS.LE.ERR) GO TO 320
IF(NCYCL.EQ.MXCYCL) GO TO 330
GO TO 250
320 MXCYCL=NCYCL
330 ERR=EPSRETURN
END【使用例】 マグニチ ュ−ド 7. 3，震央距離 25. 0 km の地震を想定し ，解放基盤面上におけ る
模擬地震動（加速度時刻歴）を策定せよ．
［解］プログラム：DIMENSION ACC(1024),VEL(1024),VW1(1024),VW2(1024),T(513),H(1),*RES(513,1),VRES(513)DATA
EM/7.3/,R/25.0/,NN/1024/,H/0.05/
DATA
IR/101/,MCYCL/10/,ERR/0.05/
C
CALL WGEN(EM,R,NN,IR,ACC,1024,DT,AMAX,VMAX,MCYCL,ERR,VW1,VW2)
CALL IACC(DT,NN,ACC,VEL,VW1,1024,VMAX,DMAX)
TD=10.**(0.31*EM-0.774)
NFOLD=NN/2+1
DO 110 K=1,NFOLD
T(K)=TD/REAL(K)110 CONTINUE
DO 120 K=1,NFOLD
IF(T(K).LE.2.) GO TO 130
120 CONTINUE
130 KMIN=K
DO 140 K=KMIN,NFOLD
IF(T(K).LE.0.02) GO TO 150
140 CONTINUEKMAX=NFOLDGO TO 160
150 KMAX=K-1
160 NE=KMAX-KMIN+1
DO 170 K=1,NE
T(K)=T(K+KMIN-1)
170 CONTINUECALL ERES(1,H,1,NE,T,513,DT,NN,ACC,1024,2,VMAX,RES)
DO 180 K=1,NE
VRES(K)=RES(K,1)
180 CONTINUESTOPEND

1234567891011121314151617181920212223242526272829303132

アウトプット：結果 は，模擬地震動の加速 度時刻歴が，配列 ACC に 格納されており，これ を積分した 速度時刻歴が 配列 VEL に，ま た減衰定数 5％の速度応 答スペクト ル
は，配列 VRES に格納されていて，プロットすれば，それぞれ下図（ａ），（ｂ），（ｃ）のようになる．
図（ｃ）に記入してある細線は，目標とした大崎スペクトルである．

-6-WGEN

マグニチュード M = 7.3

震央距離

R = 25.0 km
最大加速度178.5 gal（ａ）

最大速度 14.93 kine（ｂ）（ｃ）速度応答スペクトル200
減衰定数 5%100応答速度（kine）5020
105

2
1.0
0.5

0.2
0.1
0.01

0.02

0.05

0.1

0.2

0.5

周期（sec）-7-

125