Journal of Natural Disaster Science, Volume 25, Number 2, 2003, pp93-101

Simulation of Earthquake Motion from Phase Information
Tadanobu SATO
Disaster Prevention Research Institute, Kyoto University, Uji Kyoto 611-0011, Japan
Yoshitaka MURONO
Rail Way Research Institute, Hikari-Cho, Kokubunji 185-8540, Tokyo Japan(Received for 10. May, 2003 and in revised form 1 Dec., 2003)

ABSTRACTWe developed a method based on the concept of wavelet transformation to simulate earthquake motion thatuses phase spectra of earthquake motions. By using Mayer’s analyzing wavelet a simultaneous linear equation
was derived from which the relative values of wavelet coefficients can be determined from the phase spectrum oneach compact support of a scale function. The relationship between the power of earthquake motions and thewavelet coefficients on each compact support was used to determine the absolute values of wavelet coefficients.
The efficiency of the proposed method was investigated by comparing resimulated with observed earthquake
motion. A stochastic model with which to simulate the phase spectrum is proposed that is based on the concept ofgroup delay time. A sample phase spectrum simulated by this stochastic model was used to simulate an artificial
earthquake motion based on the proposed method.

1. INTRODUCTION
Estimation of the responses of structural systems based on
dynamic analyses is essential for the seismic design of civil structures. How design earthquake motions should be defined has beena main concern of civil engineers since the 1995 Hyogoken Nambu
earthquake in Japan (JSCE 1996).
Methods used to simulate design earthquake motion fall into
three categories: (i) Theoretical methods based on the source rapture mechanism and elastic wave propagation theories, (ii)
Methods using the empirical Green’s function, (iii) Stochastic simulation methods. Methods in (i) and (ii) can be used to simulate
precise earthquake motion at a design point if detailed information
can be provided to define the fault mechanism, propagating pathand local site effects. Earthquake motions that are used for establishing measures against future earthquake damage sometimes are
simulated theoretically using scenario earthquakes and taking into
account the seismic environment in the area of concern (Irikura
2000). But to design conventional civil structures, such as small
bridges, this method is not practical because the position of theactive fault and its source parameter, which control the rupture
mechanism, can not be easily determined.
Methods in class (iii) often are used to define design earthquake motions because they have simple forms and only small
number of parameters are needed to define the models. A primitive
method developed during the early stage studies of this category
was to generate a stationary time history by summing up the harmonic waves with random amplitudes and phases (Yang 1972).
The nonstationary power spectrum was modeled by a regression
equation as a function of earthquake magnitude and epicentral distance (Kameda 1977), and based on the results, a method to simulate nonstationary strong motion was developed (Goto et al. 1979).

Most often used method simulates earthquake motion from amplitude spectra expressed by the attenuation relationship and a random phase. To guarantee nonstationarity, an envelope function ismultiplied to the stationary time history simulated by means of arandom phase (Jennings et al. 1968).
The modeling of amplitude spectra has claimed the attentionof many researchers. The attenuation relationship of the Fourier
amplitude was derived in the U.S.A. in the late 1970s (Trifunac1976) and that of response spectra was proposed in Japan in late1980s (Kamiyama and Ynagisawa 1986, Kawashima et al. 1984).
Moreover to account for the fault extent the concept of equivalentsource distance was introduced (Ohno et al. 1993). Based on the
two-stage least square method (Joyner and Boore 1981) that eliminates the dependency between magnitude and distance, Fukushima(1994) derived an attenuation relationship by using shortest distance to earthquake faults. Rapture directivity also was introduced
into the attenuation of response spectrum in the case of using theequivalent source distance (Ohno et al. 1998).
Several pioneer works clarified the nonstationary characteristics of earthquake motions through the analyses of their phasecharacteristics. Osaki et al. (1984) showed a similarity between the
distribution width of the phase difference and the duration of earthquake motion. Izumi and Katsukura (1980) who used the conceptof group delay time, showed that the average arrival time of earthquake energy and the duration of earthquake motion can be evaluated by the mean and standard deviation of the group delay time.
Sawada et al. (1986), Ishi et al. (1987), and Satoh et al. (1996) analyzed the phase difference or the group delay time of earthquake
motions and derived attenuation relationships related to the duration, mean and standard deviation of the group delay time. We alsostudied the phase characteristic of earthquake motions and developed two methods (Sato et al. 1999a, 1999b) to simulate phase

KEY WORDS: Phase spectra, simulation of earthquake motion, wavelet coefficient, wavelet transform, group delay time

T. SATO AND Y. MURONO

94

spectra. One is the modeling of phase spectra near the source
region, taking into account the source rupture mechanism, propagating path and local soil condition (Sato et al. 1999a). The other isthe modeling of phase spectra using the concept of group delay
time and wavelet analyses, in which we proposed regression relationships that express the mean group and the standard deviation of
the group delay time as functions of epicentral distance and earthquake magnitude (Sato et al. 1999b). Using the models developed
to simulate phase spectra, we simulated earthquake motions compatible with design acceleration response spectra.
The simulated motions based on these methods, however, do
not satisfy causality because the amplitude and phase spectra mustbe given independently. The purpose of the study reported here
was to overcome this deficiency in simulating design earthquake
motions and to develop a simulation method compatible with given
phase spectra.

(1) Algorithm for Simulating Earthquake Motions
Compatible with Phase Spectra
A method to simulate a time history on the scale factor ofMayer’s analyzing wavelet is developed when the group delay
time or phase spectrum is given on this scale factor. If the wavelet
coefficients a jk of the time history f (t) are given, the discrete
inverse wavelet transformation is defined by (Sasaki et al. 1992)
N -1

MÂ Â a ◊ j (t )
jk

j =1

(1)

jk

k =0

in which N is given by 2j ( j=1,2,…,M ), and the total numberof digitized points in time is 2M, and ψjk(t) is the analyzing wavelet
defined by(

j jk (t ) = j 2 j T -1t - k)(2)

The decomposed time history fj(t) of an earthquake motion on
the jth scale factor is expressed by
f j (t ) =t jk = Tkj
2

(7)

(

)
T)

Yj (w ) = T 1 2 2 - j 2 Y 2 - j Tw
=T

12

2

-j 2

◊

(

◊ exp - iw 2

( F( 2- j -1

Tw- j -1(8))) + (F( 2 Tw))
2

-j

2

where Φ(ω) is a scaling function. If the amplitude and phase spectra of the summation term in equation (6) are expressed by
Ajp (w ) and f pj (w ) thenAjp (w )e

2. THEORYf (t ) =in which Aj (ω) and φj (ω) respectively are the amplitude andphase spectra of Fj (ω) . The values of tjk and Ψj (ω) have the following form when Meyer’s formulation of the analyzing wavelet(Mayer 1989) is used;

- if

p
w
j( )

=

N -1

Âa e

- iwt jk

(9)

jk

k =0

in which Ajp (w ) and f pj (w ) are defined by2

ÏÔ N -1¸Ô ÏÔ N -1¸ÔA (w ) = Ì
a jk sin wt jk ˝ + Ìa jk cos wt jk ˝ÔÓ k =0Ô˛ ÔÓ k =0
Ô˛
p
j

Â

( )

Â

( )

¸
Ï N -1
a jk sin wt jk ÔÔÔÔ
f p (w ) = tan -1 Ì- Nk =-01˝
Ôa jk cos wt jk ÔÔÔ˛
Ó k =0

Â

( )

Â

( )

2

(10)

(11)

The relationship between φj (ω) and φp (ω) is obtained fromequations (6), (8) and (11)f j (w ) = f pj (w ) - w Ê Tj +1 ˆ
Ë2 ¯

(12)

N -1

Â a ◊ j (t )
jk

(3)

jk

If the value ofφj (ω) is given, the following equation is
obtained by substituting equation (11) into equation (12)k =0

If we compose an analyzing wavelet by the method of Mayer(1989), this wavelet has a compact support in the frequency
domain, and the Fourier transform of equation (3) is

N -1

Â {cos(wt ) ◊ b (w) + sin(wt )} ◊ a = 0
jk

j

jk

jk

(13)

k =0

Fj (w ) =N -1

Âa jk ◊ T 1 2 ◊ 2 - j 2 ◊ e- iw Tk
2j(◊ Y 2 - j Tw

k =0

)

(4)

in which T is the duration of earthquake motion and Ψ(ω)
the Fourier transform of the analyzing wavelet ψ(t) expressed by
Y (w ) =•Ú j(t ) ◊ e dt
- iwt

(5)

-•

- f j (w )

N -1

Âa ◊e

= Yj (w )

jk

k =0

- iwt jk(b j (w ) = tan f j (w ) + w T 2 j +1

)

(14)

To solve equation (13) with respect to ajk, we consider m discrete circular frequency points (i=1,2,…,m) on the compact support of the jth scale factor and define a new variable

(

)

(

a jk (w i ) = b j (w i ) ◊ cos w it jk + sin w it jk

Furthermore, we can rewrite equation (4):Fj (w ) = Aj (w )ein which(6))

k=0,1,…, 2 j-1, i=1,2,…,mEquation (14) is then rewritten(15)

SIMULATION OF EARTHQUAKE MOTION FROM PHASE INFORMATION

È a j 0 (w 1 ) a j 1 (w 1 ) L a jN -1 (w 1 ) ˘ Ï a j 0 ¸ Ï0 ¸
˙Ô
Í
Ô Ô Ô
a jN -1 (w 2 ) ˙ Ô a j 1 Ô Ô0 Ô
Í a j 0 (w 2 ) a j 1 (w 2 )=Ì˙ M ˝ ÌM ˝
Í
M
O
M
Ô Ô Ô
˙Ô
Í
ÍÎa j 0 (w m ) a j 1 (w m ) L a jN -1 (w m )˙˚ ÔÓa jN -1 Ô˛ ÔÓ0 Ô˛(16)

È a j 1 (w 1 ) a j 2 (w 1 ) L a jN -1 (w 1 ) ˘ Ï a j 1 ¸˙ÔÍÔ
a jN -1 (w 2 ) ˙ Ô a j 2 ÔÍ a j 1 (w 2 ) a j 2 (w 2 )
Ì
˙ M ˝
Í M
O
M
Ô
˙Ô
Í
ÍÎa j 1 (w m ) a j 2 (w m ) L a jN -1 (w m )˙˚ ÔÓa jN -1 Ô˛
Ï a j 0 (w 1 ) ¸ÔÔÔ a (w ) Ô= -Ì j 0 2 ˝ÔÔ MÔa j 0 (w m )Ô
˛
Ó

Ú

2
f j (t ) dt = 1
2p
-•

Ú

•

-•

Fj (w ) dw
2

-•

Ú

(20)

Taking into account the fact that Mayer’s analyzing wavelet comprises a completely orthonormal system;N -1Êf j (t ) dt =
a Á 1
Ë 2p
-•
k =0

Ú

•

Â

2

2
jk

ˆ
Yj (w ) dw˜
¯
-•

Ú

•

2

N -1
Ê
ˆ •
2
2
Fj (w ) dw = a j20 ◊ Á 1.0 +a jk2 ˜ ◊Yj (w ) dw
-•
-•
Ë
¯
k =1

Ú

Â

Ú

(21)

(22)

If Mayer’s analyzing wavelet is used, the integral appearing
on the right side of equation (22) isÚ Y (w) dw = 2 p
-•

2

j

Âa

2
jk

(2) Algorithm for Simulating Earthquake Motion Compatible
with the Amplitude Spectrum
The wavelet coefficients can be determined by using theamplitude spectrum of earthquake motion. If the amplitude spectrum of the time history on the jth scale factor Aj (ω) is given thefollowing equation is obtained from equation (6)Aj (w ) = Ajp (w ) ◊ Yj (w )(26)

in which Ajp (w ) is defined in equation (10). Substituting equation(10) into equation (26) gives
Aj2 (w ) =

N -1

Â
k =0

2

a 2jk Yj (w ) + 2{(

N -1 N -1

Â Â a Y (w ) ◊jkjk =0 l = k +1

a jl Yj (w ) ◊ cos w t jl - t jk

)}

(27)

Because this is a nonlinear equation with respect to ajk, it is
not so easy to construct an algorithm to obtain ajk. In the following
analyses a method to determine an ajk compatible with the phase
spectrum therefore was used.

This was made using the observe ground motion at Morioka
during the 1994 Off Sanriku earthquake. Fig.1 shows the wavelet
coefficients, ajk, on the scale factor of j=5〜12, and Fig.2 thedecomposed time history fj (t) on each scale factor. The Fourier
amplitude spectrum is shown in Fig.3. The frequency range ofMeyer’s analyzing wavelet on the compact support of the jth scale
factor is{ f | 2 3T £ f £ 2
j

•

2p

(25)

n-1

3. VERIFICATION OF THE METHOD DEVELOPED

Combining equations (19) and (21) and taking into accounta jk = a j 0 ◊ a jk , the equation (19) becomes•ÂOnce aj0 is given, the value of ajk compatible with the phase
spectrum φj(ω) is determined, after which the decomposed time
history of earthquake motion on the jth scale factor fj (t) can be resimulated.

2

Ú Â

n-1
Ê
ˆ
2 pÁ 1.0 +a jk2 ˜
Ë
¯
k =1

(17)

(19)

Ê N -1ˆf j (t ) dt =a jk j jk (t )˜ dt
Á
-•
-• Ë
¯
k =0
•

2

(24)

Substituting equations (23) and (24) into equation (22) gives theexpression for aj0;

Using equation (3), the left hand side integral is expressed by
•

2

j

k =0

in which the dimension of matrix [A] is m×(N-1). The number of
unknowns is N-1 but the order of equation (18) is m. The leastsquares method is used to obtain a jk by means of the QR deconvolution technique of matrix [A].
To determine the absolute value of ajk, a relation is required todefine the value of aj0. We assume that the power of earthquake
motion is exerted on the compact support of jth scale factor.
Perseval’s equality gives•Ú F (w) dw

= lj

(18)

jk

•

aj 0 = l j

For simplicity equation (17) is rewritten[ A] ◊ {a } = -{B}The left side of equation (22) is equal to the power of the earthquake motion on the compact support of jth scale factor. We
assume that this is expressed by l2j as
l2j =Because equation (16) is a homogeneous simultaneous equation, we assume the condition aj0=1 to obtain the ratio of ajk to aj0,
as expressed by a jk (k=1,2,…,N-1)

95

(23)

j +2

3T

}

(28)

in which T is the duration of earthquake motion. Each support isoverlapped by neighboring supports (Fig.3). The central frequancy
range of each support (Sato et al. 1999b) is

T. SATO AND Y. MURONO

96

Fig. 1{f | 2

Wavelet coefficients for each scale factor j (j=5〜12)

j -1

T £ f £ 2j T}

Fig. 2

Wavelet decomposed wave for each scale factor j (j =5〜12)

Fig. 3

Fourier amplitude of the wavelet decomposed wave for eachcompact support(29)

shown by broken lines in Fig.3.
To determine the ajk wavelet coefficients as being able toexpress precisely the phase spectrum for the frequency range
defined by equation (29) we selected m discrete circular frequencies, ωi (i=1, 2,…, m), in this frequency range and calculatedφj (ωi) after which we detemine the ajk wavelet coefficients using
equations (17) and (25).
For the case of j=5 Figs.4 and 5 show the time history ofwavelet coefficients a5k and the recomposed time history fj=5(t).
Three m values (32, 64, 128) were considered to determine theeffect of number of discrete points, m, on the identified wavelet
coefficients. The case of m=32 is equivalent to 2 j ( j=5), and is one
value larger than the number of unkown values a jk (k=1,…,31).
As the number of discrete points increases, the duration of thecomposed time signal increases. We therefore used m=2j for thefollowing caluculations.
For the case of j=5〜12 the time history of the ajk wavelet
coefficients are shown in Figs.6 and the recomposed time history
fj (t) in Fig.7. A comparison of Figs. 2 and 3 shows that the recomposed time histories agree well with the decomposed time histories
from the observed earthquake motion for all the j values.
Resimulated earthquake motion summing up the recomposed time
history for j=5〜12 is shown in Fig.8. Nonlinear response spectra

SIMULATION OF EARTHQUAKE MOTION FROM PHASE INFORMATION

Fig. 4

Fig. 5

Simulated wavelet coefficients for j=5

Fig. 7

Simulated decomposed waves for the scale j (j=5〜11)

Fig. 8Comparison of observed and simulated seismograms

Simulated decomposed wavelet for j=5

Fig. 9Fig. 6

Simulated wavelet coefficients for each scale factor j (j=5〜11)

Comparison of non-linear response spectra calculated fromobserved and simulated seismograms97

T. SATO AND Y. MURONO

98(ductility demand spectra) calculated with this resimulated earthquake motion are shown in Fig.9. For comparison those calculatedusing the band-passed (correspond to the frequency range j=5〜12) observed earthquake motion also are shown. Not only wave
forms, but the nonlinear response characteristics of the structural
system obtained using resimulated earthquake motion agree wellwith those of the observed earthquake motion.
4.

MODELING PHASE SPECRA USING GROUPDELAY TIME

We developed a simple model to simulate phase spectra using
the concept of group delay time and wavelet analysis. The groupdelay time of each earthquake motion was calculated on each compact support of Mayer’s analyzing wavelet using existing data setsof observed earthquake motions. If the phase spectrum of a time
history of the jth scaleφ( j ) (ω) is given, the group delay time on thejth compact support at the circular frequency of ωi is defined byt (grj ) =f
df( j )=dw

( j)(w ) - f ( ) (w )
j

i

Dw

i +1

(30)

Because of the fluctuating nature of the group delay time oneach compact support the mean value of the group delay time,js tgr(w ) , and standard deviation, stgrj (w ) , were calculated;( j)m tgr =jN( )t (grj ) (w i )Â N( )
j

i =1

( j)

s tgr =

Fig. 10

j
N( )

Â(

1
( j)
t (grj ) (w i ) - m tgr
N ( j ) i =1

)

2

(31)

Regression equations as functions of epicentral distance, Δ,
and the earthquake magnitudes, M, for these values were derived(Sato et al. 1999b).
( j )M( j)m tgr= a 1( ) ¥ 10 1bj( j )M( j)s tgr= a 2 ( ) ¥ 10 2
b

j

g

( j)

¥D1
g

(32)

( j)

¥D2

(33)

The power on the j th support, l2j , also is expressed by the
regression equation( j )Ml2j = a 3 ( ) ¥ 10 3
j

b

g

( j)

¥D3

(34)

The observed earthquake motion data sets used are from the1993 Off South West Hokkaido earthquake (M7.8), the 1994 Off
East Hokkaido earthquake (M8.1), the 1995 Far-off Sanriku earthquake (M7.5), the 1995 Hyogoken Nambu earthquake (M7.2) andthe 1997 East Kagoshimaken earthquake (M6.3).
Results of the regression analyses are shown in Table 1 and
2
Fig.10. Except for j=7, the value of l j increases as M increases
and attenuates with the epicentral distance. The effect of M and Δ
2on the value of l j becomes less remarkable for small j values. The
2value of l j is not solely affected by earthquake magnitude andepicentral distance, but by soil conditions of the site as well.
Classifying the value of a (3j ) from the soil or site classification(Kamiyama and Yanagisawa 1986, Kawashima et al. 1984)
improves the accuracy of the regression equation, but for simplicity those effects are neglected. The saturation effect of attenuation
relationship at the near source region is also not considered.

Attenuation relationships for the power of earthquake motion on the j th compact support

SIMULATION OF EARTHQUAKE MOTION FROM PHASE INFORMATION

Table 1

5.

Regression analyses results

SIMULATION OF EARTHQUAKE MOTIONSFROM SIMULATED PHASE SPECTRA

Phase spectra of earthquake motions can be simulated by the
regression equations given in section 4. Appling the method developed in section 3, earthquake motions compatible with the phase
spectra can be simulated as follows:(1) Procedure
The proposed simulation method is composed of next nine
steps.
(i) Decide an earthquake magnitude M and an epicentral distance Δ.
j(ii) Calculate m tgr(w ) and stgrj (w ) using the regression
equations.
(iii) Generate random numbers based on normal distributionjj
N m tgr, s tgrthen calculate t grj (w ) . Obtain a phase spectrum φj (ω) by integrating this value with respect to ω.
(iv) Simulate a time history f j¢(t ) on the jth scale factor
using a simulated phase spectrum φj (ω) and amplitude
spectrum Fj¢(w ) with the value of 1.0 within the frequency range {f | 2 j−1/T≦f ≦2 j/T} and 0.0 outside of that
range.
(v) Take Fourier transform on f j¢(t ) and obtain a new phase
spectrum f¢j (w )(vi) Calculate relative values of wavelet coefficient a jk . using
this phase spectrum f¢j (w ) .
(vii) Calculate aj0 from the regression equation of power l2j
and obtain wavelet coefficient ajk by multiplying aj0 by
a jk .
(viii) Obtain a time history fj (t) on the jth scale factor using theinverse wavelet transformation.
(ix) Summing up fj (t) for a certain range of j values give theartificial earthquake motion.
The wavelet coefficient ajk can be obtained directly form thephase spectrum on the jth scale factor by the method developed in
section 2. Steps (iv) and (v) seem redundant, but without them atime history, fj (t) can not be properly simulated. The phase spectrum and Fourier amplitude of fj (t) overlap the neighboring scale
factors in the frequency ranges {2 j /3T<f <2 j -1 /T} and{2j/T<f <2j +1/3T}. This effect is avoided by assigning 0.0 to the

(

)

99

Fourier amplitude within these frequency ranges in steps (iv) and(v).
(2) Examples of simulated earthquake motions
Earthquake motions of magnitude M=8 are simulated for three
epicentral distances: 50(km), 100(km) and 200(km). Because a setof random values is used to simulate group delay time on the scale
factors, concerned the simulated earthquake motion is a sample
from the mother set of time histories. Examples of simulated
results are shown in Figs.11 and 12. Time histories for several
scale factors that cover important frequency ranges for the seismic
design of civil structures are shown in Fig.11. Artificial earthquake
motions obtained by summing up all the time histories related tothe scale factors concerned are shown in Fig.12.
Fig.11 shows that the main arrival time of the time history
delays as a scale factor decreases as the epicentral distance increases. Fig.12 shows a decrease in maximum acceleration and elongation of duration as the epicentral distance increases and that longer
period motion predominates in the later part of simulated earthquake motions, especially for a long epicentral distance.
(3) Verification of simulated earthquake motions
Values of the maximum accelerations and velocities obtained
from simulated earthquake motions were compared with an existing attenuation relationship of maximum acceleration and velocity
to show the efficiency of the proposed simulation method. Several
attenuation relationships are proposed based on the concept ofequivalent source distance (Ohno et al. 1993) and shortest distanceto the earthquake fault (Fukushima 1994). But the attenuation relationship used here is the simple one proposed by Kawashima et al.
(1984) because the earthquake magnitude M and epicentral distance Δ were selected as parameters of the regression relationshipof group delay time. Results are shown in Fig.13 in which dataobtained for M=8.0 with is plotted for epicentral distances of 50,
100 and 200 km. Several data points are at the same epicentral distance because the sample earthquake motions could be simulated
from different sets of random values. Simulated maximum accelerations have small values for long epicentral distances when compared with Kawashima’s attenuation relationship but the maximum
velocities fit the attenuation relationship well. In our regression
relationship there is no term to account for the saturation effect in
the near source region; therefore, the estimation of large maximum
accelerations and velocities in the near source region is possible.

100

T. SATO AND Y. MURONO

Fig. 12

Simulated seismograms for M=8, Δ=50, 100 and 200(km)

Fig. 13

Maximum acceleration and velocity values of simulated
waves as comparedwith M-Δ relation (Kawashima et al )6. CONCLUSION

Fig. 11

Simulated wave for scale factors j =7−11 based on attenuation relationships (M=8)

Using the wavelet transformation technique we developed amethod with which to simulate earthquake motions compatible
with given phase spectra. The following results were obtained.
(1) Because, using Mayer’s analyzing wavelet, a decomposed
time history is defined on a compact support in the frequency domain (the amplitude spectrum of the decom-

SIMULATION OF EARTHQUAKE MOTION FROM PHASE INFORMATION

(2)

(3)

(4)

(5)

(6)

posed time history has values in a certain frequency
range), the relative values of the wavelet coefficients
could be obtained from the phase spectrum of this compact support.
If the power of the decomposed time history on a compact
support is given, the absolute values of the wavelet coefficients on this support can be calculated. Therefore theoriginal time history was shown to be recomposed basedon the inverse wavelet transform.
The efficiency of simulating earthquake motion by the
proposed method was evaluated by resimulating theobserved earthquake motion recorded during the 1995
Far-off Sanriku earthquake.
A method to simulate phase spectra that uses the conceptof group delay time was introduced, as were the regression relationships of the mean value and standard deviation of the group delay time on the compact support asfunctions of earthquake magnitude and epicentral distance. The regression relationship of the power of earthquake motion on the compact support was derived.
A phase spectrum was obtained by integrating the simulated group delay time using proposed regression relationship. Earthquake motion compatible with the phase spectrum then was simulated.
The appropriateness of the earthquake motions simulated
by our novel method to for seismic design purpose wasverified by a comparison with existing attenuation characteristics of maximum acceleration and velocity
REFERENCES

Arakawa, T., Kawashima, K. and Aizawa, K. (1984). “Input groundmotions for step-by-step integration dynamic response analysis bymodifying response spectral characteristics”, Civil Engineering
Journal, Vol.26, No.7 (in Japanese).
Fukushima,Y. (1994), “Empirical prediction for strong ground motion
reflected on theoretical backgrounds of source and propagation ofseismic wave”, ORI Report 93-07, Ohsaki Research Institute, (inJapanese).
Goto, H., Kameda, H. and Sugito, M. (1979), “Prediction of strong earthquake motions by the evolutionary precess model”, J. Struct. Mech.
Earthquake Eng., JSCE, No.286, 37-51 (in Japanese).
Irikura, K. (2000), “Prediction of strong motions from future earthquakes
caused by active faults: case of the Osaka basin”, Proc. the 12th
World Conference on Earthquake Engineering, Paper No.2687.
Ishii, T. and Watanabe, T. (1987). “Relation between phase characteristics
of ground motions and magnitudes, focal distances and focal depths
of earthquakes”, Summaries of Technical Papers of Annual Meetingsof AIJ, Structures I, 385-386 (in Japanese).
Izumi, M. and Katukura, H. (1980). "Representation of nonstationary
nature of seismic waves based on phase (Part I)”, Proc. of theInternational Research Conference on Earthquake Engineering,199212.
Japan Society Civil Enginners (1996), “Proposal on earthquake resistancefor civil engineering structures”, Special report of the committee of101

earthquake resistance on civil engineering structures, Japan.
Jennings, P. C., Housner, G. W. and Tstai, N. C. (1968). “Simulated earthquake motions for design purposes”, Proc. of the 4-th World
Conference on Earthquake Engineering, A-1, 145-160.
Joyner, W. B. and Boore, D. M. (1981): “Peak horizontal acceleration andvelocity from strong-motion records, including records from 1979
Imperial Valley, California, earthquake”, Bulletin of SeismologicalSociety of America, 71, 2011-2038.
Kameda, H. (1977), “On a method of computing evolutionary power spectra of strong motion seismograms”, J. Struct. Mech. Earthquake Eng.,
JSCE, No.235, 55-62 (in Japanese).
Kamiyama, M. and Yanagisawa, E. (1986), “A statistical model for estimating response spectra of strong earthquake motions with emphasis
on local soil conditions”, Soil and Foundations, 26, 16-32.
Kawashima, K., Aizawa, K. and Takahashi, K. (1984), “Attenuation ofpeak ground motions and absolute acceleration response spectra”,
Proc. the 8th World Conference on Earthquake Engineering, 257264.
Meyer, Y. (1989), “Orthonormal wavelets in wavelets”, Springer, 21-27.
Ohno, S., Ohta, T. Ikeura, T. and Takemura M. (1993), “Revision of attenuation formula considering the effects of fault size to evaluate strong
motion spectra in the near field”, Tectonophysics, 218, 69-81.
Ohno, S., Takemura, M. and Kobayashi, Y. (1998), “Near-fault rupture
directivity effects on strong-motion records”, Proc. of the 10 th Japan
Earthquake Engineering Symposium, 133-138 (in Japanese).
Ohsaki, Y., Kanda, J., Iwasaki, R., Masao, T., Kitasa, T. and Sakata, K.
(1984). “Improved methods for generation of simulated earthquake
ground motions”, Proc. the 8th World Conference on Earthquake
Engineering, 573-580.
Sasaki, F., Maeda, T. and Yamada, M. (1992). “Study of time history datausing the wavelet transform”, J. Struc. Eng., AIJ, Vol.38BÅC9-20 (inJapanese).
Sato, T., Murono, Y. and Nishimura A. (1999a). “Modeling of phase characteristics of strong earthquake motion”, J. Struct. Mech. Earthquake
Eng., JSCE, No.612/I-46, 201-213 (in Japanese).
Sato, T., Murono, Y. and Nishimura A. (1999b). “Modeling of phase spectrum to simulate design earthquake motion”, Optimizing post-earthquake lifeline system reliability, Technical Council on Lifeline
Earthquake Engineering, Monograph No.16, 804-813.
Satoh, T., Uetake, T. and Sugawara, Y. (1996). “A study on envelope characteristics of strong motions in a period range of 1 to 15 seconds byusing group delay time”, Proc. of 11th World Conference on
Earthquake Engineering, Paper No.149.
Sawada, T., Nagae, M. and Hirao, K. (1986), “A definition of the durationof earthquake ground motion based on phase differences and its statistical analysis”, J. Struct. Mech. Earthquake Eng., JSCE, No.368/I-5,
373-382 (in Japanese).
Trifunac, M. D. (1976). “Preliminary empirical model for scaling Fourier
amplitude spectra of strong ground accleration in terms of earthquake
magnitude, source-to-station distance, and recording site conditions”,
Bulletin of Sesimological Soc. of Am., Vol.66, No.4, 1343-1373.
Yang, J. N. (1972), “Simulation of random envelope process”, Journal ofSound and Vibration, Vol.21, No.73.