EQGEN86 AN ARTIFICIAL EARTHQUAKE
ACCELEROGRAM SIMULATION PROGRAM
by
BonHsiang Lien
B.A., National Taiwan University, Taiwan, R.O.C.
A thesis submitted to the
Faculty of the Graduate School of the
Universtiy of Colorado at Denver in partial fulfillment
of the requirements for the degree of
Master of Science
Department of Civil Engineering
1986
This thesis for the Master of Science degree by
BonHsiang Lien
has been approved for the
Department of Civil Engineering
by
Lien, BonHsiang (M.S., Civil Engineering)
EQGEN86 An Artificial Earthquake Accelerogram Simulation Program
Thesis directed by Professor NienYin Chang
For the purpose of simulating earthquake acceleration time
histories for use in a dynamic analysis or dynamic testing, a user
oriented artificial earthquake simulation program, EQGEN85, was
coded by Huang and Chang (1985). At the present time, the EQGEN85
is developed on PRIME minicomputer system and written in FORTRAN77.
According to userspecified earthquake magnitude, sourcesite
distance and local site condition, and a design spectrum, many
spectrumcompatible accelerograms can be generated in just a few
iterations each. Three earthquake strong motion data by Chang,
Seed, and Trifunac, respectively, were implemented in EQGEN85. A
peak ground acceleration database, a duration database, and an
attenuation database were implemented in the revised computer
program, EQGEN86, in this study.
In order to simulate an actual ground motion record, an
envelope function, which is magnitudedependent, is used in EQGEN85.
In EQGEN86, the envelope function was modified as a function of
magnitude and sourcesite distance. The EQGEN86 can be used to
obtain a simulated accelerogram for preliminary designs.
ACKNOWLEDGEMENT
I wish to express my sincere gratitude to my advisor
Professor NienYin Chang for his support and guidance during my
graduate study at the University of Colorado at Denver. Gratitude
is also extended to Professor Tzong H. Wu and Professor John R. Mays
for serving in the examination committee.
This project was partially sponsored by the Waterways
Experiment Station (WES), the U.S. Army Corps of Engineers,
Vicksburg, Mississippi under a Contract No. DAC3984M4780 with the
University of Colorado at Denver. The project manager from WES was
Dr. Ellis L. Krinitzsky of the Geotechnical Laboratory. The above
sponsorship was gratefully appreciated.
Appreciation also goes to my family in Taiwan, Republic of
China for their support and understanding.
I would like to express special thanks to Pauline LeBlanc
and HsienHsiang Chiang for their help in the preparation of this
thesis.
Finally and most importantly, I would like to express my
thanks and love to my wife, ChuanChuan Chen, for typing my thesis
and for her encouragement during the past one and half years.
CONTENTS
ACKNOWLEDGEMENT................................................. iv
TABLES........................................................1 ix
FIGURES......!................................................... x
CHAPTER
I. INTRODUCTION............................................... 1
1.1 Earthquake Ground Motions........................... 1
1.2 Factors Affecting Ground Motion
Characteristics.................................. 3
1.3 Selection of Design Earthquakes.............. ... 7
1.4 Purpose of the Development of Artificial
Accelerograms...................................... 9
1.4.1 Dynamic Testing and Dynamic Analysis.... 9
1.4.2 Application of Artificial Earthquake
Records.................................. 9
1.5 Objectives and Scope of Study..................... 11
II. DEVELOPMENT OF EQGEN..................................... 13
2.1 Flow Chart of the EQGEN86 Computer Program...... 14
2.2 Generation of the Simulated Accelerogram........ 19
2.2.1 Fourier Transform........................... 20
2.2.2 Convolution............................. 23
2.2.3 Discrete Fourier Transform.............. 25
2.2.4 Fast Fourier Transform..................... 31
vi
CHAPTER
2.3 Baseline Correction in EQGEN.................... 39
2.4 Calculations of the Velocity and Displacement
According to the Generated Accelerogram....... 43
2.5 Generation of Response Spectra.................. 45
2.6 Improvements to the execution of EQGEN.......... 56
III. REQUIRED INPUT PARAMETERS OF EQGEN..................... 57
3.1 Design Earthquake Magnitude......................... 58
3.1.1 Definition of Earthquake
Local Magnitude, Ml.................... 58
3.1.2 Revised Curve of the Amplitude
Attenuation Function Proposed
by Richter, 1958. .................... 63
3.1.3 Relation Between Magnitude
and MM Intensity........................... 64
3.1.4 Relation Between Magnitude
and Slipped Fault Size................. 71
3.2 SourceSite Distance................................ 71
3.3 Local Site Condition............................... 73
3.4 Duration of Strong Motion....................... 76
3.4.1 Definition of Strong Motion Duration... 77
3.4.2 Estimation of Strong Motion Duration... 77
3.4.3 Strong Motion Duration Database in EQGEN.. 83
3.5 Peak Ground Acceleration
(For Chang's model only).................;. ... 86
3.5.1 Estimation of Peak Ground Acceleration... 87
3.5.2 Epicentral Peak Ground Acceleration
Database in EQGEN........................ 103
3.6 Choice of the Type of Design Spectrum Input..... 108
3.7 Choice of the Type of the Accelerograms............ 109
3.8 Damping Values of Designed Structures............. 109
Vll
CHAPTER
3.9 Choice of Different Graphic Output................. 109
3.10 Further Comment of the Required
Input Parameters.............................. 110
IV. EMPIRICAL MODELS IN THE EQGEN STRONG MOTION DATABASE... Ill
4.1 Chang's Model (1981)............................. 112
4.1.1 Improvement of the Chang's Model (1981)
Implemented in EQGEN86................... 119
4.2 Seed's Model (1982)................................ 119
4.3 Trifunac's Model (1976).......................... 121
4.4 Limitations of Each Empirical Model................ 125
4.4.1 Chang's Model.............................. 126
4.4.2 Seed's Model............................. 126
4.4.3 Trifunac's Model........................... 128
V. ENVELOPE FUNCTION IN EQGEN............................. 131
5.1 The Envelope Function Proposed
by Huang and Chang (1985)........................ 131
5.2 Development of a Revised Envelope Function...... 135
5.2.1 The Revised Envelope Function for
Nearfield (R < 10 km),
Hard Site Condition...................... 138
5.2.2 The Revised Envelope Function
for Farfield (R > 10 km,
Hard Site Condition...................... 141
5.2.3 The Revised Envelope Function
for Nearfield (R < 10 km),
Soft Site Condition...................... 141
5.2.4 The Revised Envelope Function
for Farfield (R > 10 km),
Soft Size Condition.....:............. 142
viii
CHAPTER
VI. FIDELITY OF EQGEN COMPUTER PROGRAM..................... 144
6.1 Discussion of EQGEN Simulation Results............ 144
6.1.1 Simulation Results by Using Chang (1981)
Strong Motion Database................... 144
6.1.2 Simulation Results by Using Seed (1982)
Strong Motion Database...................... 148
6.1.3 Simulation Results by Using Trifunac (1976)
Strong Motion Database...................... 148
6.1.4 Result Comparisons of the Implemented
Earthquake Strong Motion Models............. 148
6.2 Consistency of Earthquake Energy................... 157
6.3 Effect of the Selection of Earthquake
Ground Motion Duration............................ 157
6.4 Compatibility of Peak Ground Acceleration
(for Chang's Model Only)........................ 168
6.5 Judgments of Spectrum Compatible
Simulation Accelerograms.......................... 169
6.5.1 Uniformity Between the Calculated Local
Magnitude and the Design Magnitude.......... 170
6.5.2 Comparison of Predominant Frequency........... 172
6.5.3 Compatibility of the Generated Spectrum
and the Design Spectrum..................... 172
VII. SUMMARY, CONCLUSION AND RECOMMENDATIONS
FOR FURTHER RESEARCH.................................... 174
7.1 Summary and Conclusions......................... 174
7.2 Recommendations for Future Research................ 176
BIBLIOGRAPHY................................................... 179
TABLES
Table
1 Properties of the Fourier Transform (from Brigham,
1974, p.46)....................................... 22
2 .............................................................. 52
3(a) Log10Ao(R) Versus Epicentral Distance R (from
Richter, 1958)............................................ 59
3(b) Revised value of Log10Ao(R) Versus Epicentral
Distance R (From Trifunac, 1976)...................... 60
4 Abridge ModifiedMercalli Intensity Scale (from
Wiegel, 1983, p. 90).................................... 67
5 Related Coefficients of Equations (3.3) and (3.4)......... 83
6 Coefficient r of Equation (3.5).......................... 86
7 Comparison of numerical values of some suggested
relationships between peak ground acceleration
(in cm/sec2) and MM intensity, or equivalent
(from Trifunac and Brady, 1975).......................... 89
8 Values of Mmax and Mm^n used in Equation (4.17)
as a function of selected period. (from
Trifunac, 1976).......................................... 125
9 Data Distribution of Chang's Model (1981)............... 127
10 Data Distribution of Trifunac's Model (1976)............ 129
11 Parameters of the Envelope Function (defined
by Huang and Chang, 1985)................................ 133
12 Parameters of the Revised Envelope Function
(for nearfield, hard site condition)................... 140
13 Parameters of the Revised Envelope Function
(for nearfield, soft site condition).................... 143
6
15
24
26
27
28
30
37
38
47
55
61
FIGURES
Comparison of the accelerograms xinder different
local site conditions (September 19, 1985,
Mexico (Michoacan) Earthquake) INSTITUTO DE
INGENIERIA UNAM, IPS10A, 200985, and IPS10B,
210985).........................................
Flow chart of EQGEN86 computer program..........
Graphical Example of Convolution (from Brigham,
1974, p. 55)....................................
Graphical example of the convolution theorem
(from Brigham, 1974, p. 60).....................
Graphical example of the frequency convolution
theorem (from Brigham, 1974, p. 63).............
Graphical development of the discrete Fourier
transform (from Brigham, 1974, p. 92)...........
Properties of discrete Fourier transform (from
Brigham, 1974, p. 129)..........................
Comparison of multiplications required by direct
calculation and FFT algorithm (from Brigham,
1974, p. 152)...................................
FFT signal flow graph, N=4 (refer to Brigham,
1974, p. 153)...................................
A onedegreeoffreedom system under an
earthquake excitation..........................
Misrepresentation of a highfrequency wave form
by a lowfrequency sampling (from Hovanessian,
1976, p. 173)...................................
Nomogram for determining earthquake magnitudes
from trace amplitudes in millimeters of a
standard torsion seismogram (from Gutenburg and
Richter, 1942)..................................
xi
Figure
3.2 Modified curve of log10Ao versus sourcesite
distance for computing local magnitude from
a generated accelerogram (from Jennings and
Kanamori, 1983)................................. 65
3.3 Correlation of three intensity scale (from
Trifunac and Brady, 1975)........................... 66
3.4 Relationship between magnitude scale and MM
intensity scale (from Krinitzsky and Chang, 1977)... 69
3.5 Relation of magnitude, MM intensity and epicentral
distance (from Krinitzsky and Chang, 1977).......... 70
3.6 Relationships between earthquake magnitude and
the length of main fault surface rupture (from
Bureau, 1976)....................................... 72
3.7 Comparison of significant durations of a rock
site and a soil site (1971, San Francisco
earthquake) (from Dobry), Idriss, and NG, 1978)..... 75
3.8 Typical available relations between magnitude
and significant duration (from Dobry, idriss,
and Ng, 1978)....................................... 79
3.9 Comparison of nearfield, earthquake ground
motion duration (from Chang and Krinitzsky, 1977)... 80
3.10 Duration versus epicentral distance and
magnitude for rock sites (from Chang and
Krinitzsky, 1977)....................................... 81
3.11 Duration versus epicentral distance and
magnitude for soil sites (from Change and
Krinitzsky, 1977)....................................... 82
3.12 Curve for estimating the earthquake ground
motion duration under nearfield (R<10 km),
hard site condition.................................... 85
3.13 Relationships between peak ground acceleration
an MM intensity scale, or equivalent (from
Trifunac and Brady, 1975)............................. 88
3.14 Correlations of peak ground acceleration
(horizontal component) and epicentral distance .
under MM intensity V (from Gupta, 1980)................. 91
Figure
xii
3.15 Correlation of peak ground acceleration
(horizontal component) and epicentral distance
under MM intensity VI (from Gupta, 1980)............. 92
3.16 Correlation of peak ground acceleration
(horizontal component) arid epicentral distance
under MM intensity scales IV and VII (from
Gupta, 1980)............................................ 93
3.17 Determination of peak ground acceleration
as a function of magnitude and epicentral
distance (from Gutenburg and Richter, 1956).......... 97
3.18 Seed's empirical earthquake strong motion
model (from Seed and Idriss, 1982)................... 101
3.19 Available correlations of epicentral peak
ground acceleration and magnitude...................... 102
3.20 Available attenuation curves for peak
ground acceleration. .................................. 104
3.21 Curves for determining epicentral peak ground
acceleration; implemented in the peak ground
acceleration database in EQGEN86.............. . .. 106
3.22 A representative attenuation curve for peak
ground acceleration................................... 107
4.1 Chang's empirical earthquake strong motion
model (from Chang, 1981)............................. 113
4.2 Relationships between peak ground acceleration
and "base" average acceleration under different
local site conditions (from Chang, 1981)............. 114
4.3 Seed's empirical earthquake strong motion model;
average acceleration spectra for different local
site conditions (from Seed and Idriss, 1982)......... 120
4.4 Effect of shaking intensity on normalized
acceleration spectra for deep cohesionless soil
sites (from Seed, Ugas, and Lysmer, 1976)............ 122
4.5 Trifunac's empirical earthquake strong motion model;
coefficients for scaling Fourier amplitude spectra
in terms of earthquake magnitude, sourcesite dis
tance, and localsite condition (from Trifunac, 1976). 124
5.1 The envelope function proposed by Huang and
Chang (1985)......................................... 132
Figure
xiii
5.2 Available relationships of earthquake strong
motion duration under nearfield, rock site
condition, and rupture duration as functions
of magnitude......................................... 137
5.3 Curves for defining an envelope functions
(CURVE A) and determining earthquake ground
motion duration (CURVE B) under nearfield
(R<10km) hard site condition....................... 139
6.1(a) Simulation results by using Chang's model
(1981) (Ml = 6.5, R = 10 km, rock site
condition, duration = 22.46 sec., and input
Amax 0.39 g); zero iteration's results..................... 145
6.1(b) Simulation results by using Chang's model
(1981) (Ml =6.5, R = 10 km, rock site
condition, duration 22.46 sec., and input
Amax = 0.39 g) ; three iteration's results................... 146
6.1(c) Simulation results by using Chang's model (1981)
(Ml =6.5, R 10 km, rock site condition,
duration = 22.46 sec., and input Amax = 0.39 g);
generated accelerograms after three iterations....... 147
6.2(a) Simulation results by using Seed's model (1982)
(Ml = 6.5, R 10 km, rock site condition, and
duration 22.46 sec.); zero iteration's results.... 149
6.2(b) Simulation results by using Seed's model (1982)
(Ml = 6.5, R 10 km, rock site condition, and
duration =22.46 sec.); three iteration's results... 150
6.2(c) Simulation results by using Seed's model (1982)
(Ml = 6.5, R 10 km, rock site condition, and
duration = 22.46 sec.);generated accelerograms
after three iterations................................. 151
6.3(a) Simulation results by using Seed's model (1982)
(Ml = 6.5, R 10 km, deep cohesionless soil
site condition; and duration = 22.46 sec.);
one iteration's results................................ 152
6.3(b) Simulation results by using Seed's model (1982)
(Ml = 6.5, R 10 km, deep cohesionless soil
site condition; and duration = 22.46 sec.);
one iteration's results................................ 153
XIV
Figure
6.4(a) Simulation results by using Trifunac's model (1976)
(Ml 6.5, R = 10 km, rock site condition, and
duration = 22.46 sec.); zero iteration's results.... 154
6.4(b) Simulation results by using Trifunac's model(1976)
(Ml =6.5, R = 10 km, rock site condition, and
duration 22.46 sec.); threeiteration's results.... 155
6.4(c) Simulation results by using Trifunac's model
(1976) (Ml 6.5, R 10 km, rock site
condition, and duration 22.46 sec.);
generated accelerograms after three iterations....... 156
6.5(a) Comparison of the simulation results by using
different empirical models; generated
accelerograms after two times iterations
(Ml 6.5, R 10 km, rock site conditions;
and duration 22.46 sec.)................................... 158
6.5(b) Generated ground acceleration timehistory,
velocity timehistory, and displacement
timehistory........................................ 159
6.5(c) Generated pseudovelocity response spectra............. 160
6.5(d) Generated absolute acceleration response spectra.... 161
6.5(e) Generated relative velocity response spectra........... 162
6.5(f) Generated relative displacement response spectra.... 163
6.5(g) Generated Fourier amplitude spectra.................... 164
6.6(a) Comparison of simulation results by selecting
different earthquake ground motion durations;
applying Seed's model. (M 6.5, R = 10 km,
and rock site condition)...................................... 165
6.6(b) Generated absolute acceleration response spectra.... 166
6.6(c) Generated Fourier amplitude spectra.................... 167
CHAPTER I
INTRODUCTION
Due to the paucity of actual ground motion records from
earthquakes with magnitude 6 or larger, frequently, an artificial
accelerogram is needed in seismic response analyses. For the pur
pose of generating artificial accelerograms, a computer program,
EQGEN85, was coded by Huang and Chang (1985). The EQGEN85, as a
useroriented, manualdriven computer program is developed on PRIME
minicomputer and written in FORTRAN 77. Based on a design spectrum,
using a stationary stochastic process with a deterministic envelope
function, and the fast Fourier transform algorithm, an artificial
spectrumcompatible accelerogram can be generated in just a few
iterations. In this study, a revised EQGEN computer program,
EQGEN86 is accomplished to enhance the capability of EQGEN85. In
this study, the term EQGEN indicates both EQGEN85 and EQGEN86 in a
general way.
1.1 Earthquake Ground Motions
Generally speaking, earthquake ground motions are induced by
stress waves imposed on earth crust due to fault rupture. The
stress waves are generated by the energy released by the rupturing
process along a fault and propagated as a complicated and random
process. The randomness as pointed out by Housner (1947 and 1955)
may be originated by acceleration pulses random in direction,
in magnitude and in time, and the complicated reflection and
refraction behavior when the stress waves propagate through various
soil and rock strata. Nowadays, more understanding of the
mechanism of earthquake ground motions has been achieved. Most
earthquakes are resulted from the release of the accumulated strain
energy in an active fault or faulting system.
An earthquake ground motion can be recorded in three dimens
ions; two horizontal components and one vertical component. Usually
an accelerometer is used to record the acceleration timehistory,
i.e. the accelerogram. In most cases, the two horizontal components
of an accelerogram have about the same shaking intensity. The
severity of an earthquake ground motion is basically characterized
by the maximum acceleration, the strong ground motion duration and
the frequency content of the shaking. According to the past
earthquake records, the vertical component frequently has a peak
ground acceleration which is equal to about one third to two thirds
of the peak ground acceleration of the horizontal component and has
a higher frequency content. This also implies that the vertical
ground motion attenuates more rapidly than the horizontal components
do (Wiegel, 1983, p.79).
Generally speaking, by inspecting accelerograms alone,
earthquake ground motions may be characterized by the following five
features (Hudson, 1979, EERI); i.e. (1) peak ground acceleration,
(2) duration of earthquake ground motion, (3) rough idea of the pre
dominant frequency and the frequency content, (4) comparison of the
characteristics of vertical and horizontal ground motions, e.g.the
shaking intensity and the frequency range, and (5) approximate
3
hypocentral distance. In engineering practice, it is, in general,
not sufficient to discern the characteristics of an earthquake
ground motion by inspecting a recorded accelerogram alone. To
better understand the characteristics of an earthquake ground
motion, a spectral analysis becomes a necessary tool. The spectrum
of an earthquake, as defined by Housner, Martel, and Alford (1953),
is a plot of the maximum response of a simple onedegreeoffreedom
oscillator under the earthquake excitation versus the period (or
frequency) of the oscillator. In other words, the spectral analysis
shows the effect of an earthquake ground motion to simplified,
viscoelastic, onedegreeoffreedom structures with different
damping values and nature frequencies. It is also noted that
different accelerograms can be considered similar if they have
similar spectral shapes. This provides a reasonable comparison
between different earthquake ground motions, which usually can not
be achieved by just inspecting an accelerogram. Detailed discussion
about spectral analysis is presented in Section 2.5.
1.2 Factors Affecting Ground Motion Characteristics
It is very important to select appropriate earthquake ground
motions for seismicresistant designs. Before selecting an appro
priate design earthquake ground motion, it is, however, necessary to
understand the factors affecting the characteristics of earthquake
ground motions.
The factors affecting earthquake ground motion character
istics are presented as follows:
4
(1) Earthquake magnitude. When describing the size of an
earthquake, the Richter's magnitude, which is used in
this study, is a measurement related to the released
energy of an earthquake (Gutenburg and Richter, 1942,
and 1956). This means that the earthquake magnitude
would reflect some characteristics of ground motions,
e.g. the ground shaking intensity, the peak ground
acceleration, and the strong motion duration (see
Sections 3.4.2, and 3.5.1).
(2) Sourcesite distance. Due to material damping and
geometric damping, the released energy of an earthquake
would attenuate with increasing sourcesite distance.
The attenuation characteristic is usually expressed by
the relationship between the peak ground acceleration
and sourcesite distance (see Section 3.5.1).
(3) Local site condition. Effects of local site conditions
have been discussed by many researchers (Trifunac and
Brady,.1975; Trifunac, 1976; Seed, et al, 1976; Boore,
Joyner, Oliver, and Page, 1980; Joyner and Boore, 1981;
Campbell, 1981; and Dobry et al. 1978, Chang, 1981).
The affected earthquake ground motion characteristics
include the peak ground acceleration, velocity,
displacement and frequency content etc. Section 3.3,
shows the situ condition effect on the shape of
accelerograms from soil sites and rock sites. The
accelerogram from soil site has a moderate intensity
with low frequency tail part. The geometry and the
5
depth of soil deposits may also affect the character
istics of surface ground shaking. The Michoacan
earthquake of September 19, 1985 induced ground motion
recorded in Mexico City (see Figure 1.1), where a very
deep clay deposit is contained in a bowlshaped
bedrock. The bedrock motion propagated upwards,
amplified and resulted in a significant surface ground
shaking. Besides, the multiple reflection might have
resulted in an extended duration of strong motion.
More detail discussions on the local site condition
effect are presented in Section 3.3.
(4) Site topography. Because seismic waves propagate in a
three dimensional space, site topography may affect the
characteristics, of earthquake ground motions, e.g. a
steep topography may result in a higher peak ground
acceleration (Campbell, 1981).
(5) Errors involved in data processing. In standard earth
quake ground motion data processing procedures, a digi
tized earthquake accelerogram is applied in further
ground motion characteristics analysis, e.g. response
spectral analyses, and calculations of ground velocity
and displacement. On many occasions, when digitizing a
photographtype accelerogram, some manmade errors
are inevitable; e.g. high frequency noise may be
induced by an operator (Hudson, 1979, EERI). The above
errors may be reduced if digital accelerographs are
used. The determination of the baseline of a recorded
ACCELERATION (gals) ACCELERATION (gals)
6
 NS component 
50 gals
 EW component 
 Vertical component 
______i_________ _______ TIME, s...
__i
11
21 31 41
(a) Stiff soil site
51
120 
 NS component 
120 
(b) Soft soil site
Figure 1.1 Comparison of the accelerorams under different local
site conditions (September 19, 1985, Mexico (Michoacan)
Earthquake) (INSTITUTO DE INGENIERIA UNAM. IPS10A,
200985, and IPS10B, 210985)
7
accelerogram may also produce some errors (Schiff and
Bogdanoff, 1967). The baseline correction is discussed
in Section 2.2.
1.3 Selection of Design Earthquakes
In a seismicresistant design problem, the selection of
design earthquakes is very important, especially for the critical
structures built in seismically active areas. The following
discussions about the selection of design earthquakes is based on
the paper published by Seed (1982a).
Due to the controversy oh the procedure of selecting design
earthquakes, Seed (1982a) pointed out that three different aspects
should be taken into consideration; i.e. the Maximum Credible
Earthquake, the Seismic Safety Evaluation Earthquake, and the
Seismic Engineering Design Earthquake.
(1) The Maximum Credible Earthquake. The Maximum Credible
Earthquake is defined as "the largest rationally
conceivable event that could occur in the tectonic
environment in which the project is located." Its
determination is a function of the recurrent periods of
major earthquake events. The selection of a proper
recurrent period is related to the importance of a
designed structure, and should be selected only by
expert engineering seismologists.
(2) The Seismic Safety Evaluation Earthquake. The Seismic
Safety Evaluation Earthquake is "the largest event and
the strongest motions which a structure should be
designed to withstand." It may be equal to or less
8
than the magnitude of the Maximum Credible Earthquake.
The magnitude of the Seismic Safety Evaluation Earth
quake depends upon the factors: the design life of a
structure; the seismic risk in engineering and societal
points of view; and the consequences of underestimating
and overestimating the seismic risk. Also due to the
variety of the considerations, the determination of the
Seismic Safety Evaluation Earthquake should be made by
a group of experts in the fields of seismology, seismic
geology, engineering, seismic risk, social science,
economics, and public policy, and government.
(3) The Seismic Engineering Design Earthquake or the
Engineering Analysis Earthquake. This is the earth
quake motion to be used in dynamic analysis or testing
to evaluate the performance and safety of a structure.
The following factors may affect the determination of
Seismic Engineering Design Earthquake and the associ
ated ground motion:
(a) The analysis method used and its conservatism.
(b) The material and system properties of structure.
(c) The soilstructure interaction effects.
The consideration in selecting a Seismic Engineering
Design Earthquake for an earth structure and for a
concrete structure may be different. For an example,
the effect of high frequency ground shaking is more
pronounce in concrete structures than in earth struc
tures. Thus the Seismic Engineering Design Earthquake
9
and the associated ground motion should be assigned by
experienced engineers.
1:4 Purpose of the Development of Artificial Accelerograms
' 1.4.1 Dynamic Testing and Dynamic Analysis
Because a strong earthquake can impose severe damages to
structures, it is important to understand the performance of
structures under earthquakeinduced ground shaking, so as to achieve
an appropriate seismicresistant design. Performances of
structures subjected to earthquake ground shaking can be determined
through dynamic testing or dynamic analysis. Using dynamic testing,
the earthquake excitation to be imposed onto a structure model is
either generated by a real time signal generator or signals stored
in a magnetic tape. Different structure responses are monitored
throughout a test to evaluate the structure performance. For
instance, shaking table tests and centrifuge tests are in this
category.
When using dynamic analysis, an accelerogram, and the
dynamic material properties of a structure and its foundation soils
are required as input data. An accelerogram is usually expressed in
terms of mathematical functions or in digitized time history.
Theoretical analyses are performed to determine the response of the
structure under a prescribed earthquake excitation to examine the
effectiveness of a seismic resistant design.
1.4.2 Application of Artificial Earthquake Records
As discussed in Section .1.2, when performing a seismic
resistant design of a structure under earthquake excitations, it
10
requires appropriate earthquake accelerograms as inputs. Either
recorded accelerograms or artificially simulated accelerograms can
be used.
(1) Recorded Earthquake Accelerograms. Due to the appli
cation of digital recording accelerographs, it becomes
more and more convenient to use recorded earthquake
accelerograms in seismic analysis. In a design pro
ject, according to the local site condition, source
site distance and the design earthquake magnitude,
existing recorded earthquake accelerograms obtained
from sites with similar conditions can be used in
seismic analysis and/or dynamic testing. During the
process of selecting a recorded accelerogram for design
purpose, it is very important to understand the factors
that affect earthquake ground motion characteristics,
as discussed in Section 1.3, so that an appropriate
accelerogram can be. chosen. Because earthquake ground
motions are complicated random events, it is
recommended to use a number of accelerograms in the
analysis to find a range of performance.
(2) Artificially Simulated Accelerogram. Although the
number of recorded accelerograms has been rapidly
increasing during the last two decades, the near field
records of events with magnitudes larger than 7.5 are
still scarce. Thus, there exists a strong need to
generate artificial earthquake accelerogram. It is
hoped that enough ground motion records will become
11
available in the future. Until then, the ground motion
simulation technique remains to be a viable tool for
obtaining design accelerograms. Using a desired
response spectrum, an accelerogram can be generated on
a digital computer to represent the ground motion
expected at a selected site. In general, the selection
of a design spectrum is based on the following design
information: earthquake magnitude, sourcesite dis
tance, and local site condition.
1.5 Objectives and Scone of Study
A useroriented computer program, EQGEN85, for generating
spectrumcompatible, artificial earthquake acceleration time
histories was coded by Huang and Chang (1985). Following the
existing main structure of the EQGEN85 computer program, the scope
of this study includes:
(1) Give a complete description of the methodology applied
in EQGEN85 computer program.
(2) Discuss the input parameters required in executing
EQGEN85 computer program.
(3) In addition to the empirical earthquake strong motion
models implemented in EQGEN85 database, in the revised
computer program EQGEN86, an earthquake strong motion
duration database and a peak ground acceleration
database are implemented to assist a user in executing
EQGEN. In the peak ground acceleration database, the
determination of epicentral peak ground acceleration
and an attenuation curve are included.
12
(4) Modify the execution procedure of EQGEN85 to enhance
the procedure of developing a useroriented computer
program.
(5) Discuss the empirical earthquake strong models imple
mented in the database of EQGEN85 and their limit
ations .
(6) Modify the envelope function to more closely simulate
an actual earthquake ground motion record.
(7) Discuss the fidelity of the simulated accelerograms.
CHAPTER II
DEVELOPMENT OF EQGEN
The program EQGEN85 was developed on the PRIME minicomputer
and written in FORTRAN 77 (Huang and Chang, 1985) by using the char
acter manipulation capability of FORTRAN 77 and interactive graph
ics. EQGEN85 is a useroriented and manualdriven computer pro
gram. The manual provides adequate guidance for a user to input
data and to obtain outputs, so as to allow a user to interactively
generate artificial accelerograms in front of a graphic terminal.
It is suitable for installation in a minicomputer or microcomputer
environment.
According to a prescribed earthquake design spectrum, a
stationary stochastic process is used to simulate earthquake ground
motion accelerograms in EQGEN85. The prescribed design spectrum can
be a Fourier amplitude spectrum or a relative velocity response
spectrum. The design spectrum can be read from a file, typed inter
actively by a user, or obtained from a strong motion database. If
the strong motion database is chosen, a user has to specify design
earthquake magnitude, sourcesite distance and local site condition.
An earthquake strong motion duration database and a peak ground
acceleration database are also implemented in EQGEN86. An
interactive graphic system allows a user to view simulation results
14
and to compare the design spectrum with a generated one. Whenever
their difference is intolerable, an iteration procedure can be
applied repeatedly until a spectrumcompatible accelerogram is
obtained.
2.1 Flow Chart of the EOGEN86 Computer Program
Basically, the ground motion simulation procedure used in
the original EQGEN includes six phases: INPUT, GENERATION, PROCESS,
CHECK, PLOT and OUTPUT phases. After Huang and Chang (1985)
developed the artificial earthquake simulation computer program,
EQGEN85, in this research, some modifications have been implemented
in the development of the revised computer program, EQGEN86.
In order to provide an overall understanding of the execu
tion of EQGEN86, a flow chart of the computer program is presented
in Figure 2.1. The six phases involved in the program execution are
described sequentially in the following sections.
(1) The INPUT phase asks a user to specify the design
earthquake magnitude, sourcesite distance and local
site condition. A user also has to provide a design
spectrum to the program. Alternatively, the design
response spectrum may also be obtained from one of the
three empirical models implemented in the strong motion
database, i.e. a design Fourier amplitude spectrum from
Chang's model (1981), or Trifunac's model (1976), or a
5% damping, design relative velocity response spectrum
from Seed's model (1982). It is noted that if the
Chang's model (1981) is applied, an input of peak
ground acceleration is required. The peak ground
Figure 2.1 Flow chart of EQGEN86 computer program
16
acceleration may also be obtained from a peak ground
acceleration database, which is discussed in Section
3.6.2. Then, a user has to assign the number of
acceleration data points and the time increment between
two consecutive data points to specify the duration of
the generated accelerogram. An earthquake strong
motion duration database in EQGEN86, discussed in
Section 3.4.3, can help a user in choosing an
appropriate duration.
(2) The GENERATION phase applies the inverse fast Fourier
transform algorithm to generate a time series X(t)
according to a specified design Fourier amplitude
spectrum and a randomly generated phase spectrum. The
Fourier transform of a time function, X(t) can be
expressed as
F(f) R(f) + j.I(f) F(f).eJ^(f) (2.1)
where R(f) and 1(f) are the real part and the imaginary
part of F(f), respectively; j= J1; F(f) the Fourier
amplitude spectrum of X(t); ^(f) the phase spectrum of
X(t); f frequency and t time (Hsu, 1967, p.85 and
p.104). Both F(f) and ^(f) are functions of fre
quency f. When applying the inverse fast Fourier
transform algorithm, the input design spectrum is
interpolated to have sufficient data points for
computing the time series X(t). Either linear inter
r.
polation or logarithmic interpolation can be applied.
17
When using the strong motion database, the linear
interpolation is. applied if the Chang's model (1981) or
the Seed's model (1982) is used; and the logarithmic
one is applied if the Trifunac's model (1976) is used.
The phase spectrum with phase angle values ranging from
0 to 27T is generated by a pseudo random number genera
tion subroutine with a uniform phase angle distribu
tion. It is noted that if a design spectrum is
retrieved from the Seed's model (1982) in the strong
motion database, the retrieved values of the 5% damp
ing, relative velocity response spectrum are set equal
to the values of the initial design Fourier amplitude
spectrum. This assignment is reasonably made according
to the close relationship between the Four.ier amplitude
spectrum and the undamped, zero damping relative
velocity spectrum (see the discussion in Section 2.5).
The generated time series X(t) is then multiplied
by an envelope function, ENV(t) to obtain a new time
series Z(t), i.e.
Z(t) X(t) ENV(t) (2.2)
The purpose of applying the envelope function is to
produce a realistic accelerogram. Detail discussions
on the envelope function are presented in Chapter 5.
The multiplication of the envelope function would
distort the generated Fourier amplitude spectrum and
change the energy content specified by the design
18
spectrum. According to Parseval's theorem (Hsu, 1967,
p.129), the energy content that specified by the design
spectrum is defined as
^design
r
X(t)2dt
if
F(o>) \zdo> J"F(w) 2df
d cd
(2.3)
where F(o>)  is the Fourier amplitude of X(t); F(w)2
is the energy spectrum or energy spectral density
function of X(t). In order to assure that the gene
rated accelerogram has the same energy content as that
specified by the design spectrum, an adjustment is
performed by the following equation to generate a new
time history:
a(t) SZ(t) (2.4)
The constant S is equal to the square root of the ratio
between the area under the design energy spectrum and
the area tinder the energy spectrum of Z(t).
(3) The PROCESS phase applies a parabolic baseline correct
ion, as discussed in Section 2.3, to the generated
accelerogram, a(t), so as to obtain a corrected
accelerogram. The corrected accelerogram is then used
to obtain other related information, such as the local
magnitude, the ground velocity and displacement
timehistories, and the different response spectra. In
this phase, the input of structure damping ratio is
19
needed if the response spectra are required. Also the
program iteration times can be prescribed at this
stage.
(4) The CHECK phase gives a graphic comparison between the
generated spectrum and the design spectrum. Their
incompatibility is basically resulted from the modific
ation procedures used in the GENERATION phase. If they
are not sufficiently compatible, a user can ask for
further iteration, which will be discussed in Section
2.6, to generate a spectrumcompatible accelerogram.
(5) The PLOT phase allows a user to obtain the graphical
representation of different simulation results, e.g.
the various response spectra, the generated acceler
ogram, ground velocity and displacement time histories
on a graphic terminal and hardcopies by using a graphic
printer.
(6) The OUTPUT phase asks a user to assign a file name to
store the generated spectrumcompatible accelerograms.
The stored data format can be assigned by a user, so
that the data can be conveniently used in seismic
analysis or testing. The user can also ask for another
accelerogram simulation by using the same design
spectrum used in the foregoing simulation.
2.2 Generation of the Simulated Accelerogram
In the EQGEN85, the inverse fast Fourier transform algorithm
is applied to convert a function in frequency domain, which is
defined by a design Fourier amplitude spectrum and a randomly
20
generated phase spectrum (refer to Equation 2.1), to a function in
time domain; i.e. a simulated, acceleration time history. The fast
Fourier transform algorithm is an efficient computational method
based on the theory of Fourier transform. The following sections
give an introduction of the development of fast Fourier transform
algorithm.
2.2.1 Fourier Transform
The Fourier transform of a time function, a(t) is defined as
the Fourier transform is a complex value that can be expressed as
where R(f) and 1(f) are the real part and imaginary part of F(f),
respectively; F(f) is called the Fourier amplitude spectrum of
x(t) and is equal to Jr2(Â£) + I2(f); and ^(f) tan1[I(f)/R(f)] is
the phase spectrum of the Fourier transform.
The inverse Fourier transform is defined as
OD
(2.5)
where F(f) is a function of frequency, f; and j = J~T. In general,
F(f) R(f) + jI(f) F(f).eJ*(f)
(2.6)
a(t) F(f).eJ27rft df
(2.7)
Usually the functions a(t) and F(f) with the above relation
ships are called a Fourier transform pair and expressed as
21
a(t) <=> F(f) (2.8)
On the other hand, the Fourier transform pair can be
expressed as
F()
a(t)e"Jwt dt
(2.9)
a(t) c2J F(w)eJwt dw (2.10)
.00
where w2?rf is the radian frequency; and it is required that
Cj *c2=l/(27r). A summary of the properties of the Fourier transform
is shown in Table 1 (Brigham, 1974, p.46).
Herein, it is important to mention that there is an alter
nate form of the inverse Fourier transform as
0
a(t) ~ [
F*(f)*ej25rft df]*
00
(2.11)
where F*('f) represents the conjugate of F(f) ; i.e. the superscript *
means a "conjugate" term. It is realized that, when comparing
Equations (2.5) and (2.11), there exists a common term,
This is important when discussing the fast Fourier transform in
Section 2.4.
22
Table 1 Properties of the Fourier Transform (from Brigham,
1974, p.46)
Time domain Frequency domain
Linear addition Linear addition
x(t) + y(t) X(f) + Y(f)
Symmetry H(t) symmetry h(f)
Time Scaling h(kt) Inverse Scale Change (lA)'H(fA)
Inverse Scale Change (1A) h(tA) Frequency Scaling H(kf)
Time Shifting h(ttQ) Phase Shift H(f) e"j27rfto
Modulation h(t)eJ2fftfo Frequency Shift H(ff0)
Even Function Real Function
he(t) He(f) Re(f)
Odd Function hQ(t) Imaginary H0(f) 
Real Function h(t) hr(t) Real part even Imaginary part odd H(f) Re(f) + jI0(f)
Imaginary Function h(t) = jhi(t) Real part odd Imaginary part even H(f) = Ro(f) + jIe(f)
Note: x(t) <=> X(f), y(t) <=> Y(f), and h(t) <=> H(f)
23
2.2.2 Convolution
Before developing discrete Fourier transform to enable
computer programming to perform the Fourier transform of the sampled
data from a continuous waveform signal, it is necessary to introduce
the concept of convolution. The convolution integral is defined as
y(t) 
x(r)h(tr)dr x(t)*h(t)
(2.12)
or an another equivalent form as
y(t)  h(r)x(tr)dr h(t)*x(t) (2.13)
' 00
where y(t) is called the convolution of the functions x(t) and
h(t). Brigham (1974) applied a graphical method, as shown in Figure
2.2, to explain the physical meaning of the convolution.
The convolution theorem is expressed by the Fourier trans
form pair as
h(t)*x(t) <=> H(f) X(f) (2.14)
On the other aspect, the frequency convolution theorem is expressed
as
h(t)x(t) <> H(f)*X(f)
(2.15)
24
FOLDING
x(T)
[Take the mirror image of
h(r) about ordinate axis]
(b)
DISPLACEMENT
[Shift h(r) by the amount of t]
(c)
r
r
A x(Tlh(tr)
r
MULTIPLICATION h(T)*(tT)
[Multiply the shifted function
h(tr) by x(r)]
(d)
t.h r
y(t)
INTEGRATION
[Area under the product of h(tr)
and x(r) is the value of the
convolution at time t]
(e)
SHADED
'AREA
t
Note: The explanations are for the lefthand side figures.
Figure 2.2 Graphical Example of Convolution
(from Brigham, 1974, p.55)
25
where H(f) and X(f) are the Fourier transforms of h(t) and x(t),
respectively. Figures 2.3 and 2.4 show the graphical examples of
the above theorems.
2.2.3 Discrete Fourier Transform
Discrete Fourier transform is an approximation of the
continuous Fourier transform discussed in Section 2.2.1, which makes
the Fourier transform performed by using digital computer feasible.
Figure 2.5 shows the development of the discrete Fourier transform
(Brigham, 1974, p.92). It Is noted that both the continuous
functions in time domain and frequency domain, i.e. h(t) and H(f) in
Figure 2.5, should be modified into discrete forms, fi(t) and H(f) by
sampling and extending to be periodic functions.
The discrete Fourier transform of a continuous function h(t)
can be expressed as
n
H() S h(kt)e'j2,rnk/N (2.16)
NT k=0
where n=0, 1, ..., Nl; N is the number of sampled data; T sampling
time interval; h(kt) sampled value; and the notation means
"approximation". Though, from the discrete point of view, the
discrete Fourier transform of a set of discretely sampled periodic
function g(kt) is
n Nl j2jmk/N
G() = S g(kt)e
NT k=0
(2.17)
26
Figure 2.3 Graphical example of the convolution
theorem (from Bringham, 1974, p. 60)
M>
27
T
H(f)
(c)
1 I
T
Id)
Figure 2.4 Graphical example of the frequency
convolution theorem (from Brigham,
1974, p. 63)
Figure 2.5 Grahical development of the discrete
Fourier transform (from Brigham, 1974,
p. 92)
29
n=0, 1, Nl the inverse Fourier transform is
1 Nl n i27mk/N
g(kt) =  2 G() e (2.18)
N k=0 NT
k=0, 1, ...,N1, and an alternative form of the inverse Fourier
transform is
g(kt)
1
N
. n j2imk/N
2 G*() e
k=0 NT
*
(2.19)
In summary, i.e.
1 Nl n i2jmk/N n Nl j2jrnk/N
g(kt) = 2 G()e <> G() = 2 g(kt) e (2.20)
N k=0 NT NT k=0
or
g(kt)
N1 ^ n
S G*()e
k=0 NT
'j2ffnk/N
n
<=> G()
NT
1 i2 7mk/N
S g(kt)e
k=0
(2.21)
and it is required that both the time and frequency domain functions
as shown in Equations (2.20) an (2.21) are periodic functions; i.e.
n (rN + n)
G() G[] rO, 1, 2,..: (2.22)
NT NT
g(kt) = g[(rN + k)T] r=0, 1, 2,... (2.23)
The properties of the discrete Fourier transform are shown
in Figure 2.6.
30
Fourier transform
x(t)+y(t) <=> X(f)+Y(f)
H(t) <=> h(f)
h(ttQ) <=> H(f)eJ2lTfto
h(t)eJ27rtfo <=* H(ff0)
[0*(f)e3^ftdf]*
he(t) <=> Re(f)
h0(t) <=> jl0(f)
h(t)=he(t)+hQ (t)
y (t)=/_^x (x )h (tT ) dx=x (t) *h (t)
y(t)*h(t) <=> Y(f)H(f)
y(t)=/_"x(T)h(t+r)dT
y(t)h(t) <=> Y(f)*H(f)
i2 h2(t)dt=f_2 H(f)i2df
Discrete Fourier transform
(A) x(k)+y (k) <=> X(n)+Y(n)
(B) iH(k) <=> h(n)
(C) h(ki) <=> H(n)ej2rrni/N
(D) h(k)ej2irk:L/N <=Â£> H(ni)
(E) Si* (n) e^ ^Irlcn^^ *
LJNn=o
(F) he(k) <=> Re(n)
(G) h0(k) <=> jl0(n)
h(k)=he(k)+hQ(k)
(H) = h(k) H (Nk) j + jh(k) h(Nk) j
(I) y (k)=NZ^x(i)h(ki)=x(k)*h(k)
i=o
(J) y(k)*h(k) <=> Y (n)H (n)
(K) y (k)= ^xCDMk+i)
i=o
a) y(k)h(k) <=> il(n)*H(n)
(M) Vh2(k)= i^lHCn)!2
k=o N n=o
PROPERTIES:
(A) Linearity
(B) Symmetry
(C) Time shifting
(D) Frequency shifting
(E) Alternative inversion
formula
(F) Even functions
(G) Odd functions
Figure 2.6
(H) Decomposition
(I) Convolution
(J) Time convolution theorem
(K) Correlation
(L) Frequency convolution
theorem
(M) Parsevals theorem
Properties of discrete Fourier transform
(from Bringham, 1974, p. 129)
31
2.2.4 Fast Fourier Transform
Due to the application of digital computer, the fast Fourier
transform (FFT), an algorithm introduced by Cooley and Tukey in
1965, becomes a powerful computational tool to handle the calcu
lation of the discrete Fourier transform (DFT). Due to the high
efficiency of applying the fast Fourier transform, it is adopted in
EQGEN computer program to generate simulated accelerograms from the
prescribed design Fourier amplitude spectrum by using inverse FFT
method, and perform the response spectral analysis of the generated
accelerogram. The following discussion is basically based upon the
introduction by Brigham (1974), which provides a clear understanding
of the fast Fourier transform algorithm.
For convenience of later discussions, we express the
discrete Fourier transform, Equation (2.17) in the following form.
j2jmk/N
X(n) Z xQ(k)* e (2.24)
k=0
where n0, 1, ..., Nl, and simplified notations are applied;
i.e. kt and n/NT are replaced by k and n, if compared with Equation
(2.17). Then let
 J27T/N
W = e
we can rewrite Equation (2.25) as
(2.25)
^1 nk
X(n) = S xQ(k).W
k=0
(2.26)
32
where n=0, 1, ..., Nl. Using N=4 as an example, i.e.
X(0) .1111 x0()
X(1) 1 w1 w2 w3 xQ(l)
X(2) 1 w2 w w2 x0(2)
. X(3) . 1 w3 w2 w1 . Xq(3) .
If there are N=2r data points involved in the calculation, where r
is an integer, we express the integers n and k in a binary form with
r bits. Then with N4, r2, i.e.
k 0, 1, 2, 3 as k ( kx kQ) =* 00, 01, 10, 11
n 0, 1, 2, 3 as n ( nx, no) 00, 01, 10, 11
i.e.
k 2kx + kQ, and n 2nx + iiq (2.28)
where kQ, kx, nQ, and nx are either 0 or 1. Then the Equation
(2.26) can be rewritten as
1 1 (2n1+n0)(2k1+k0)
X(nltn0) > S S x0(k1,k0)W
ko0 kx=0
(2.29)
It is noted that the number of summations is equal to r. Consider
ing the W term in Equation (2.29), we have
(2nx +nQ) (2kx +kQ) ^(2n^+n0)2k1 ^(2n1+n0)k0
= (^1) \j^noki ( ^ni +nG ^
= (1) .w(2ni+no)ko
w2n0ki .w(2ni+n0)k0
(2.30)
33
because for N=4,
w4niki = e(jair/N)(4n1k1) = e(j*27r) (nxkx ) = 1
Then Equation (2.29) can be written as
1
X(nx,nQ) ko=0
Equation (2.31) represents the basic equation of the fast Fourier
transform. We set
k SQ x0(k1(k0)W
2n0kx
W
(2^+noJko
(2.31)
1 2n0k1
Xi(no,k0) S x0(k1(k0)W
kx=0
i.e. in a matrix form as
(2.32)
xx (0,0) ' ' 1 0 W 0 xo(0,0)
xx (0,1) 0 1*0 w xo(0,l)
xx(1,0) o is o r4 xo(l,0)
Xid.l) . 0 10 2 . . Xo(l,l) .
(2.33)
Then we set
1 ^nj+n^ko
2 Xi^.koJW
k0=0
*2(nQini)
(2.34)
34
i
e.
x2(0,0) 1 i* *3 o o 1 xl(0,0)
x2(0,1) 1 W2 0 0 Xx(0,1)
x2(1,0) 0 0 1 W1 Xid.0)
. x2(l,l) . o o t* . Xid.l) .
(2.35)
Compare Equations (2.31) and (2.34), we have
X(nx ,nQ) = XaCno,^) (2.36)
Combine Equations (2.32), (2.34), and (2.36), we have
1 2n0kr
Xi.(n0,k0) = S x0(k1,k0)W
kjo
1 (2^+^)^
x2(n0,ni)  2 x^no.k^W
k0=0
X(nlfn0) = x2(n<5,n1) (2.37)
Equation (2.37) represents the formulations of the FFT algorithm
with N=4 condition. Also combine Equations (2.33), (2.35), and
(2.36), we have
' X(0) ' X(0,0) ' ' x2(0,0) ' x2 (0)
X(2) X(1,0) x2(0,1) x2(l)
X(l) X(0,1) x2(1,0) x2 (2)
. X(3) . . X(l,l) . . x2(l,l) . . M3) .
35
1 w 0 0
1 w2 0 0
0 0 1 w1
0 0 1 w3
1 o o ' xo(0,0)
0 1 0 w xo(0,1)
1 0 W2 0 x0(l,0)
_ 0 1 0 w2 . x0(l,l) .
(2.38)
It is noted that, in Equation (2.38), there exists a bitreverse
relation between X and x2 terms.
The efficiency of the FFT algorithm can be illustrated as
following. Also the calculation with N4, r2 conditions is used as
an example. From Equation (2.33), we have
Xi(0)  xo(0) + W0xo(2) (2.39)
and
(2) = xo(0) + W2*x0(2)
xo(0)  Wx0(2) (2.40)
In Equations (2.39) and (2.40), only one complex multiplication and
two complex additions are needed. The calculations of xx(l) and
Xi(3) are in thej same situation. From Equation (2.35), we have
x2 (0)  Xl (0) + W x, (1) (2.41)
and x2(l)  xx (0) + t^x^l)
 Xj(0) Wx1(l) (2.42)
which show the same situation as Equations (2.39) and (2.40). So,
in total, it needs Nr/2=4 complex multiplications and Nr=8 complex
36
additions when using the FFT algorithm. If the direct discrete
method in Equation (2.26) is applied, there need N2 complex multi
plications and N(Nl) complex additions. If assumed the computing
time is'proportional to the number of multiplications, the ratio of
the computing time of the direct discrete method to the FFT algor
ithm is
N2 2N
Nr/2 r
(2.43)
This comparison is also shown in Figure 2.7, where we can realize
the efficiency of using the fast Fourier transform algorithm,
especially when the number of sampled data is large. Another
advantage of applying the FFT algorithm is the economy of assigned
memory spaces in computer programming. This can be explained by
using a signal flow graph of the FFT algorithm shown in Figure 2.8.
In Figure 2.8, the calculation of the node value in each computa
tional array is the summation of the input array paths from the
foregoing node values; the factor, Wterm shown along the path is a
multiplication factor applied to the foregoing node value. It is
noted that xo(0) and x0(2) only affect xx(0) and xt(2) values; so as
the x1(0) and xx(l) to x2(0) and x2(l). This implies an inplace
computation can be performed so as to save computer memory spaces.
In Figure 2.8, although the final results of Xterm are in a
scrambling order, a bitreverse procedure can be applied to make the
X terms in an unscrambling order. The characteristic of the
bitreverse procedure is shown in Equation (2.38).
NUMBER OF MULTIPLICATIONS (xlOOO)
1024
512
256
128
64
Direct i Calculation i
FFT Algorithm
i
64 128 256 512 1024
N (number of sample points)
Figure 2.7 Comparison of multiplictirons required
by direct calculation and FFT algorithm
(from Brigham, 1974, p. 152)
38
COMPUTATIONAL ARRAYS
Data Array
x0(k)
Array 1
xi (k)
 Bitreverse
Array 2 Procedure
x2 (k)
X(0)
\
X(l)
X(2)
X(3)
Figure 2.8 FFT signal flow graph, N=4 (Refer to Brigham,
1974, p. 153)
39
Examining the Equations (2.17), (2.19), and (2.21), though
the inverse fast Fourier transform algorithm was not pointed out in
the above discussions, it is essentially the same as the fast
Fourier transform algorithm mentioned above, because there exists a
common term, ej27mk/Ni between the Fourier transform and the
inverse Fourier transform.
The above discussions of the fast Fourier transform algo
rithm is in the base of 2, i.e. the sampled data number is the power
of 2, which is adopted in the EQGEN85. More general cases have been
discussed detailly by Brigham (1974).
2.3 Baseline Correction in EOGEN
In the standard data processing of a recorded accelerogram,
due to the characteristic of an accelerograph, it is difficult to
recognize the true axis of an accelerogram with zero acceleration.
The difficulty comes from: (1) the initial part of the earthquake
shaking that happened before an accelerometer is triggered is not
recorded; (2) the final acceleration or velocity is not zero because
of the existence of ground noise; (3) the net ground displacement
after shaking is not known (Schiff and Bogdanoff, 1967). Before any
further processing of an accelerogram, such as spectral analysis, is
conducted to characterize an earthquake ground motion, it is neces
sary to apply baseline correction to locate an appropriate zero
acceleration axis of a recorded accelerogram. The baseline correc
tion method proposed by Berg and Housner (1961) adopts the baseline
with a parabolic form as Cx + C2t + C3t2, where t is time in second,
and the constants Cj^, C2 and C3 are determined with the requirement
of a minimum mean square computed ground velocity. It is also
40
assumed that the initial ground velocity and ground displacement are
zeros. More detailed discussion about the parabolic baseline
correction are discussed below (Schiff and Bogdanoff, 1967).
In a recorded accelerogram, the baseline, i.e. zero accel
eration axis is assumed to have a parabolic form as
where the constant term, Cx implies that the acceleration may not be
zero when t=0. The firstorder term, C2t implies that the assigned
baseline may not be parallel to, the true baseline. As to the
secondorder term, C3t2, Schiff and Bogdanoff (1967) pointed out
that there was no physical justification for it, but suggested it as
a correction for the distortion of recording paper.
ground velocity, the coefficients Cj/s are determined as discussed
below:
Cx + C2t + Cgt2
.(2.44)
Using the principle of minimizing the mean square calculated
3
x(t) xQ(t) 2 (^tCJ1)
j1
0 < t < T
 x0(t) Cx
t = 0
x(t) x0(t) 2 7 Cjt
31 J
x(t) xQ(t)
(j+D
VQt do
(2.45)
j=l 3(3+1)
41
where x0 is the unadjusted acceleration and x is the adjusted
acceleration. With the following assumptions:
V0 ~ io<0) *<) = o
dQ xo(0) x(0) 0 (2.46)
The mean square ground velocity is
T
 x2(t)dt (2.47)
0
where T is the duration of a recorded accelerogram,
value of Equation (2.47), we have
T
a r
x2(t) dt = 0
3C J
0
where i 1,2,3. From Equations (2.45), (2.46), and (2.48), we have
3 T+1> Vo
S CiBi
j=lJj(i+j+D (i+l)T
To minimize the
(2.48)
 Bi
(2.49)
where i = 1,2,3; and
T
tWirft) dt (2.50)
1
T(i+2)
0
42
From Equation (2.49), the coefficients Cj/s are determined as
Cx 300 Bx 900 B2 + 630 B3
C2 (1800 Bx + 5760 B2 4200 B3)/T
C3 ( 1890 Bx 6300 B2 + 4725 B3)/T2 (2.51)
sa that the corrected accelerogram can be obtained by using
Equation (2.45).
The effects of applying the corrected accelerogram in
further analyses (Schiff and Bogdanoff, 1967) are briefly discussed
below:
(1) The attempt to determine permanent earthquakeinduced
ground displacement by integrating the corrected
accelerogram is impractical becauise of the imposed
baseline correction procedure and the lack of the
initial part of the accelerogram preceding the
triggering of an accelerograph. Some uncertainties may
also be induced when digitizing a recorded
accelerogram.
(2) The adjustment of the baseline of an accelerogram would
not affect the calculation of response spectra over the
frequency range. Though the errors involved in data
reading procedure do have significant effect on
response spectrum calculations.
(3) The assumptions made in the parabolic baseline
correction mentioned above, i.e. zero initial ground
velocity and ground displacement, are not strictly
43
correct. Thus, a certain degree of uncertainty is
induced when using a corrected accelerogram.
Although digital accelerographs, which give the acceler
ograms not requiring the baseline correction are currently
available, most recorded accelerograms used in characterizing
earthquake ground motions, e.g. those used in developing the
empirical strong motion models, such as Chang's model (1981), Seed's
Model (1982), and Trifunac's model (1976), are believed to be
phototype accelerograms. So, EQGEN still adopts the parabolic
baseline correction method proposed by Berg and Housner (1961) in
the PROCESS phase to adjust a generated accelerogram, a(t) from the
GENERATION phase. The corrected accelerogram is then used in
further analyses.
2.4 Calculations of the Velocity and Displacement According
to the Generated Accelerogram
Generally speaking, most of the available earthquake ground
motion records are accelerograms. Given an accelerogram, the
associated ground velocity and the ground displacement time
histories can be calculated by integrating and double integrating
the accelerogram, respectively. Theoretically, using the above
procedure is better than using the differentiation of displacement
time history to derive the ground velocity and acceleration time
histories.
Though the calculations of ground velocity and ground
displacement are straight forward, a main difficulty is that the
baseline of an accelerogram is not exactly known. There are many
studies discussing the best way to carry out the integrations
44
(Hudson, 1962, and Schiff and Bogdanoff, 1967). The integration
methods adopted in EQGEN85 is referred to those proposed by Berg and
Housner (1961). After the zero baseline correction, as discussed in
Section 2.3, is performed to generate a corrected accelerogram, the
ground velocity and ground displacement time histories are
calculated as follows:
vi vo + (ao +ax)h/2
xi x0 + h*v0 + (2a0 + ax)h/6 (2.52)
where aQ, vQ and xQ are the ground acceleration, velocity and
displacement at time t, respectively; ax., vx and xx are the values
at time t+h; and h is the time interval between two consecutive data
points. It is also assumed that the initial velocity and
displacement (i.e. at t=0) are zeros. It has been pointed out that,
when processing recorded accelerograms, it is not reasonable to
assume zero initial ground velocity and ground displacement, because
usually the ground shaking happens about 0.1 to 0.2 seconds before
an accelerometer is triggered (Boyce, 1970, and Schiff and
Bogdanoff, 1967). At the present time, however, initial ground
velocity and ground displacement are assumed zeros in EQGEN85.
It is necessary to discuss the accuracy, of the above
calculations imposed on actual recorded accelerograms, especially
the accuracy of the ground displacement time history calculation.
As pointed out by Wiegel (1983, p.110), "the largest error in most
double integration calculations appears to be in the digitization of
the original accelerogram and in the location of the true baseline
of zero acceleration." Besides, it is understood that in the
45
earthquake accelerogram data processing procedure (Hudson, 1979,
EERI, P.34), a digital filtering method is involved. It is,
therefore, impractical to calculate permanent ground deformation by
double integrating a digitized accelerogram.
2.5 Generation of Response Spectra
For structure engineers, it is important to understand the
structural behavior under earthquake excitations. To assist in
assessing the effects of an earthquake on structures, spectral
analyses are usually conducted to determine the response spectra of
an accelerogram. From the spectral analysis, one can discern the
characteristics of an earthquake, such as the predominant frequency,
so that an aseismic structure design can be properly conducted. The
concept of applying spectral analysis in earthquake engineering
problems was first introduced by Benioff (1934). The response
spectrum of a transient loading, e.g. an earthquake ground motion,
is the curve of "the maximum response of a single degree of freedom,
linear, second order system versus the undamped nature frequency of
the system for various fraction of critical damping (Schiff and
Bogdanoff, 1967, p.872)." The system response can be velocity,
displacement or absolute acceleration.
In EQGEN85, a corrected accelerogram, a(t) generated in the
PROCESS phase is' used as an input motion to singledegreeoffreedom
structures of different damping factors and natural frequencies to
obtain different response spectra. At the present time, if the
strong motion database in EQGEN85 is used, only one horizontal
component or vertical component is treated. According to the study
of Mostaghel and Ahmadi (1982), the resultant horizontal
46
acceleration has a peak value about 1.2 times the larger peak
acceleration of the two horizontal components. This implies that
the response spectra generated by EQGEN should be used with care in
aseismic structural design.
In response spectrum calculations, a linear, singledegree
offreedom structure is subjected to the earthquake ground motion.
The associated rheological model is shown in Figure 2.9. In Figure
2.9, m is the mass of a structure, k the spring constant of the
system, C the coefficient of viscous damping, Z(t) the earthquake
induced ground acceleration time history, and x the relative
displacement of the structure subjected to the ground motion. The
equation of motion of the above system is:
mx(t) + cx(t) + kx(t) mZ(t) (2.53)
Equation (2.53) can be rewritten as
x(t) + 2junx(t) + rt2x(t) Z(t) (2.54)
where wn ./k/m undamped circular natural frequency,
(2.55)
f = C/Ccr =,viscous damping factor or damping ratio,
(2.56)
Ccr = 2 ,/km = 2mwn 2k/wn = critical damping coefficient.
(2.57)
47
x
Figure 2.9 A onedegreeoffreedom system under an
earthquake excitation
48
Considering an underdamped system, i.e. 0
conditions x(0)=0 and x(0)=0, the solution of Equation (2.54) is
t
1 f fn(tr)
x(t,wn,D Z(r) e sin dr (2.58)
wd J
0
where d = ^n A_f2 (2.59)
 the damped circular natural frequency of the structure.
Then the relative displacement response spectrum, is defined as:
I x(t,Wn,f) I max I x(tm,Wn,f) I (2.60.)
where tm is the time when the structure has a maximum relative
displacement response. Usually 0 < tm < T, where T is the duration
of the earthquake ground motion.
From Equation (2.58), we may have
x(t,wn,r)
Z(r)e r^cosW(j(tr)
0
dr
+ i
f
Jo
Z(r)*e r^sinti)(j(tr) dr
(2.61)
and then the relative velocity response spectrum, Sv is determined
as
Sv =  x(t,w^f) max = I xC^m^n'S) I
(2.62)
49
and the absolute acceleration, g (t, wn, (*) is determined as
g(t,n,f)
n(l2r2)
fwn(tr)
Z(7)e sinw^Ctr) dr
t
f rn'(tr)
+ 2Ca>n I Z(r)*e cosWjjCtr) dr (2.63)
0
If assume the f value is very small, so that l2f2 1, Vlf2 1,
f/yif2 0, and u,j ti^yif2 c,^, the absolute acceleration,
gCt,^,^) can be expressed as
gCt.Wn.O ~ x(t) + Z(t)
n2x(t,n,f) 2wnrx(t,wn,n (2.64)
Error induced by using the above assumptions will increase with
increasing damping factor value, though the induced error is of
little significance because of the existed random characteristic of
earthquake ground motions (Hudson, 1962). Then the absolute
acceleration response spectrum, Sa is defined as
Sa  x(t) + Z(t) max
I g(t,<*>nf)  max I g(^mÂ£<>hf) I (2.65)
Also., according to Equation (2.61), using the same assumptions, a so
called pseudovelocity response spectrum, SpV can be defined as
50
S
pv
t
fWnCtr)
Z(r)e cosW(j(tr) dr
0
max
t
I
rn(tr)
Z(r)e sinu^Ctr) dr
max '
0
>irsd
(2.66)
and from Equation (2.63), a pseudoacceleration spectrum, Spa can be
defined as
S
pa
t
(fwn(tr)
Z(r)e sinw
0
max
n2'Sd nspv (2.67)
In EQGEN computer program, a numerical method is applied to
calculate the response spectra discussed above (refer to Craig,
1981, Chapter 7). The solution of Equation (2.54) can be expressed
in the following form:
Xi+l All A12 Xi Bn X2 Zi
, . + *
. xi+1 . . A21 ^2 2 . . Xi . . 21 2 2 . . zi+l .
 [ a ] ' Xi + [ b ] ' Zi
. % ; . zi+1 .
where i=0, 1,..., Nl, and N is the number of data points. Known
Zj/s values and assuming Xo=0 and Xo=0, the subsequent relative
displacement and relative velocity can be calculated. The absolute
51
acceleration is then calculated by using Equation (2.64). For an
underdamped single degree of freedom system, the elements of matrix
[A] and [B] are functions of and f, and are shown in Table 2.
Herein it is noted that the damping factor is assumed to be a very
small value.
Using the above numerical method, x(t), x(t) and x(t)+Z(t)
values are evaluated step by step. The maximum value of each of
them is the S Sv, and Sa, respectively. Also the pseudovelocity
spectrum, SpV is calculated by using Equation (2.66).
In a special case, assuming a single degree of freedom, un
damped, linear system, we can relate the Fourier amplitude spectrum
and the undamped relative velocity spectrum as discussed below.
Known f 0, from Equation (2.54), we have the relative
velocity of an undamped system as
t
^(t.Wn)  Z(r)cosa)n(tr) dr
0
t
 cos Wnt J* Z(r)coswndT + sin unt
0
t
Z(r)sinwndr
0
(2.69)
The undamped relative velocity response spectrum, (SV)Q is
defined as the maximum value of the relative velocity, so that
Table 2 Elements of [A] and [B] matrix of Eq. (2.68)
(Refer to Craig, 1981, p. 142)
, Jf
11 4dh 1
coso>dh
,^h
, 2/Joy
2
n
2 a2
<*dP
2
n
 0h
)
s inudh 
+
2
W,
n
12
^h
2 fl2
2
n
cosw^h + Wjjh 
20Wc
2
n
)
sinwdh +
2u>g0
2
L n
An e
./8h
cosw^h +
 sinwdh
l >d
A!2 [ ^ sinw^h
B
21 2 ,
nudh
^h
(0 + unh) sin wdh + udcos wdh
 wd
B
22 2 ,
wnwdb
2
Un
l wd
_ AL
 e p (jj sin w^h + to
e"^1 sinudh
*
a22 e_/3h coswdh  Â§_ d . sinwdh
where 0 a fud and h a Atj
53
(Sv)o
t t
[ ZCrJcosc^dr]2 + [ ^ Z ( r ) s ino>ndr ]2
0 0
(2.70)
Further, it is assigned that the maximum response occurs at
t=T, where T is the duration of the earthquake ground motion. This
assumption is not quite correct, because frequently the maximum
response happens before the end of earthquake shaking. Based upon
this assumption, we may have:
(Sv)0 / [
Z(r)coswndr']2 + [
Z(r)sinwndr]:
(2.71)
Next the Fourier transformation of Z has the following form:
T
I* it^r
F(w) Z(r).e dr (2.72)
0
According to Equation (2.72), we have the Fourier amplitude
spectrum as
I F(W)  = / [
Z(T)cosu>ndr ]2 + [
Z(r)sinwndr ]
(2.73)
54
From Equations (2.71) and (2.73), one can readily observe
that the undamped relative velocity spectrum and the Fourier
amplitude spectrum are equal. When applying the Seed's model (1982)
in the strong motion database in EQGEN, the above expressions
provide a base of equating the retrieved 5% damping, relative
velocity response spectrum values to the initial design Fourier
amplitude spectrum values in the GENERATION phase of EQGEN as
mentioned in the Section 2.1.
When processing photographtype accelerograms, the accuracy
of spectral analysis depends upon the sampling frequency of the
digitized acceleration values. Higher sampling frequency would
produce more accurate result. Usually the sampling frequency should
be at least two times the highest frequency of an accelerogram, so
as to prevent aliasing interpretation: Figure 2.10 shows that a
high frequency wave is taken as a lower frequency wave due to a low
sampling frequency (Hovanessian, 1976). If a digital recording
seismograph is applied, the response spectrum can be easily
calculated according to the digitized accelerogram.
As indicated in Section 2.1, the modification procedures
used to generate the simulated accelerogram a(t) in the GENERATION
phase may result in incompatibility between the generated spectrum
and the design spectrum, and further iteration becomes necessary.
In EQGEN, new Fourier amplitude spectrum values are used in
the iteration procedure to generate a more spectrumcompatible
accelerogram (Huang and Chang, 1985). If the design spectrum is a
Fourier amplitude spectrum, the new Fourier amplitude spectrum,
FAS(f)new is determined as
FAS(f)new FAS(f)old + 0.2[FAS(f)design FAS(f)old] (2.74)
55
Time
Figure 2.10 Misrepresentation of a highfrequency wave
form by a lowfrequency sampling (from
Hovanessian, 1976, p. 173)
56
and the phase spectrum, ^(f) remains the same as the old one. The
constant 0.2 used in Equation (2.64) is an arbitrarily assigned
value.
If the design spectrum is a relative velocity response
spectrum, then the new Fourier amplitude spectrum, FAS(f)new is
computed as
FAS(f)new FAS(f)old[ SV(f, f)design/SV(f, f)Qld 1 (275)
where SV(f, f) is the relative velocity response spectrum. Also the
old phase spectrum is used in the iteration procedure. The above
iteration procedure can be executed repeatedly until a spectrum
compatible accelerogram is generated.
2.6 Improvements to the execution of EOGEN
The EQGEN computer program coded by Huang and Chang (1985)
is a useroriented and manualdriven program. The following
modifications are further performed in a revised computer program
called EQGEN86 to facilitate the execution of EQGEN:
(1) During the execution of EQGEN, a beep from the terminal
will remind a user to input data.
(2) For the convenience of a user, in the PLOT phase of
EQGEN, the user can ask for a complete set of plots and
hard copies with just one command.
(3) Generally speaking, in order to generate a spectrum
compatible accelerogram, at least, two to three inter
ations are needed. When executing EQGEN, it takes a
while for each iteration procedure. So, in the PROCESS
phase of the revised EQGEN, the iteration times can be
predefined for user's convenience.
CHAPTER III
REQUIRED INPUT PARAMETERS OF EQGEN
For the purpose of properly applying EQGEN to simulate
earthquake ground motions, the user has to understand the required
input data for executing the EQGEN. A list of the required input
data is shown below:
(1) Design earthquake magnitude;
(2) Sourcesite distance;
(3) Local site condition;
(4) Duration of strong motion;
(5) Peak ground acceleration;
(6) Choice of the type of accelerograms;
(7) Damping values of design structures; and
(8) Choice of different graphic output.
Details of each required input are discussed in the follow
ing sections. Before further discussions, one should realize that
due to the randomness and complexity of earthquake ground motion,
with the above parameters, it is still not quite enough to perfectly
describe the characteristics of an earthquake ground motion. As in
dicated by Trifunac (1976), under the situation with limited strong
motion records and investigated information, using generally avail
able parameters to determine earthquake ground motions, a certain
degree of uncertainty of the determined ground motion is inevitable.
58
3.1 Design Earthquake Magnitude
The importance of the selection of design earthquakes has
been discussed in Section 1.2. In the process of the selection, the
design earthquake magnitude is the most fundamental parameter that
has to be assigned. The following sections discussed the determi
nation of the earthquake magnitude.
3.1.1 Definition of Earthquake Local Magnitude. Mj_
At the present time, the magnitude scale proposed by Richter
(1935) is still the most popular and useful indicator of the
strength of an earthquake shaking. The Richter's magnitude scale,
Ml, also called local magnitude, is defined as
MLlog10(A) logio(Aq)
or
Ml log10A(R) log10Ao(R)
In Equation (3.1), A represents the maximum amplitude recorded by a
WoodAnderson torsion seismograph at an epicentral distance of 100
km, and Aq is equal to one thousandth of a millimeter. A(R), in
Equation (3.2), is the traced amplitude of the WoodAnderson
seismograph in millimeters, and log10Ao(R), the distance correction
coefficient, is shown in Table 3(a). Both A(R) and log10Ao(R) are
function of the epicentral distance, R, in km, of the recording
site. To determine the magnitude of an earthquake from the known
traced amplitude in millimeters of a standard torsion seismogram, a
nomogram was proposed by Gutenberg and Ritcher (1942) as shown in
Figure 3.1. As a standard torsion seismometer, the WoodAnderson
(3.1)
(3.2)
59
Table 3(a) Log10Ao(R) Versus Epicentral Distance R
(from Richter, 1958)
R(km) R(km) R(km) l6io Aq(R)
0 1.4 140 3.2 370 4.3
5 1.5 150 3.3 380 4.4
10 1.5 160 3.3 390 4.4
15 1.6 170 3.4 400 4.5
20 1.7 180 3.4 410 4.5
25 1.9 190 35 420 4.5
30 2.1 200 3.5 430 4.6
35 2.3 210 3.6 440 4.6
40 2.4 220 3.65 450 4.6
45 2.5 230 3.7 460 4.6
50 2.6 240 3.7 470 4.7
55 2.7 250 3.8 480 4.7
60 2.8 260 3.8 490 4.7
65 2.8 270 3.9 500 4.7
70 2.8 280 3.9 510 4.8
80 2.9 290 4.0 520 4.8
85 2.9 300 4.0 530 4.8
90 3.0 310 4.1 540 4.8
95 3.0 320 4.1 550 4.8
100 3.0 330 4.2 560 4.9
110 3.1 . 340 4.2 570 4.9
120 3.1 350 4.3 580 4.9
130 3.2 360 4.3 590 4.9
60
Table 3(b) Revised value of log^A^R) Versus Epicentral Distance
R (from Trifunac, 1976)
R(km) logioA0(R) R(km) *l6i o R(km)
0 1.400 140 3.230 370 4.336
5 1.500 150 3.279 380 4.376
10 1.605 160 3.328 390 4.414
15 1.716 170 3.378 400 4,451
20 1.833 180 3.429 410 4.485
25 1.955 190 3.480 420 4.518
30 2.078 200 3.530 430 4.549
35 2.199 210 3.581 440 4.579
40 2.314 220 3.631 450 4.607
45 2.421 230 3.680 460 4.634
50 2.517 240 3.729 470 4.660
55 2.603 250 3.779 480 4.685
60 2.679 260 3.827 490 4.709
65 2.746 270 3.877 500 4.732
70 2.805 280 3,926 510 4.755
80 2.920 290 3.975 520 4.776
85 2.958 300 4.024 530 4.797
90 2.989 310 4.072 540 4.817
95 3.020 320 4.119 . 550 4.835
100 3.044 330 4.164 560 4:853
110 3.089 340 4.209 570 4.869
120 3.135 350 4.253 580 4.885
130 3.182 360 4.295 590 4.900
61
B = Maximum trace amplitude, measured from
Figure 3.1 Nomogram for determining earthquake magnitudes
from trace amplitudes in millimeters of a
standard torsion seismogram (from Gutenburg
and Richter, 1942)
62
instrument has a nature period of 0.8 second, a damping ratio of
0. 8, and nominal static magnification of 2800. The standard
seismograph can be treated as .a onedegreeoffreedom oscillator
with an equation of motion as shown in Equation (2.53), and its
response under an earthquake excitation is the same as the response
spectral analysis discussed in Section 2.5. As a result, the
maximum response of the WoodAnderson seismograph does not necessary
happen simultaneously with the maximum ground velocity or the
maximum ground acceleration (Kanamori and Jennings, 1978).
In engineering applications, it is very important to realize
that the proposed Richter's earthquake magnitude scale is basically
applied for describing the earthquakes with focal depth about 20 km,
1. e. shallow focus earthquakes, happened in the Western United
States, especially in California. For the earthquakes with focal
depth deeper than about 40 km, the Richter's magnitude scale is not
applicable (Gutenburg and Richter, 1942). Also cautions should be
exercised when applying the Richter's magnitude scale in the areas
other than the Western United States.
Although there are several other earthquake magnitude
scales, e.g. bodywave magnitude scale, and surfacewave magnitude
scale, the local magnitude, Ml is most popular in engineering
applications. The WoodAnderson instrument can sensitively record
the ground motion with shorter period (typically 0.2 to 3 sec.),
which is of the most interest to engineers. One of the main
deficiencies of using the local magnitude scale is raised due to
frequent offrange records during strong earthquakes of the
WoodAnderson torsion seismograph. As a result, records are,
sometimes, incomplete and the standard measurements of Ml are not
available. Kanamori and Jennings (1978) proposed an alternate
method to determine the local magnitude. That is to use the
accelerogram recorded by an accelerograph as a strong motion input
to the equation of motion of the WoodAnderson torsion seismograph,
and a synthetic seismogram is generated. This synthetic seismogram
is then used in determining the local magnitude of an earthquake.
63
The methodology of generating the synthetic seismogram Is the same
as the spectral analysis of an accelerogram and is implemented in
the EQGEN computer program to calculate the local magnitude of the
generated accelerogram, and the calculated local magnitude is then
compared with the design earthquake magnitude.
Due to different site conditions, orientations to fault, and
epicentral distances, it may be possible to have different magnitude
measurements of an earthquake at different seismological stations.
It is suggested that the magnitude of an earthquake is determined by
i
averaging the magnitude values measured by several seismological
stations. The station measurement correction was discussed by
Gutenburg and Richter (1942), and Richter (1958).
3.1.2 Revised Curve of the Amplitude Attenuation Function
Proposed by Richter. 1958
Ideally, the measurement of local magnitude Ml, describing
the size of an earthquake, should be independent of the epicentral
distance of a recording station. As indicated by Kanamori and
Jennings (1978), basically the original amplitude attenuation
function proposed by Richter (1958) (see Table 3(a)), does not have
to be modified for calculating the local magnitude of a shallow
focus, nearfield earthquake with magnitude up to about 6.5
except in some particular cases. For instance, in Borrego Mountain
Earthquake, 1968, a large amount of energy is released in a single
large impulse, and the response WoodAnderson seismograph would be
larger in the nearfield.
After the method proposed by Kanamori and Jennings (1978)
that calculates local magnitude Ml by using synthetic WoodAnderson
64
seismograms as discussed in the Section 3.1.1, more magnitude
measurements become available. A detail study about the uniformity
of the. local magnitude measurements was conducted by Jennings and
Kanamori (1983). Their conclusions showed that if the standard
attenuation curve, log10Ao(R) (Richter, 1958) was applied, the
values were consistent between farfield and zero epicentral
distance measurements. At about .40. km and 50 to 60 km away,
however, there were about 1/4 magnitude unit decrease and increase,
respectively. Based upon the differences, a revised curve of the
distance correction coefficient, log10Ao(R) was proposed as shown in
Figure 3.2. In EQGEN, the revised form of the distance correction
coefficient is adopted in the calculation of the local magnitude of
the generated accelerogram.
3.1.3 Relation Between Magnitude and MM Intensity
When performing the seismicity study of a design site, it is
necessary to review historic earthquake records. In the early time,
earthquakes are described by an intensity scale.. There are three
major intensity scales in the world, i.e. the Modified Mercalli
(MM) intensity scale used in the United States, the JMA (Japan
Meteorological Agency) intensity scale used in Japan, and the
GEOFIAN (Geophysics Institute of the Academy of Sciences) scale used
in Russia. The correlation of the above intensity scales are shown
in Figure 3.3. Herein the Modified Mercalli intensity scale, which
is described in Table 4, is used in this thesis. The disadvantage
of the MM intensity scale, as pointed out by Housner (1947), is that
the determination is strongly related to building damages, which may
differ in different sites subjected to the same shaking intensity.
R, km
Figure 3.2 Modified curve of log^QAo versus sourcesite
distance for computing local magnitude from
a generated accelerogram (from Jennings and
Kanamori, 1983)
66
(A)
(S)
(C)
I II in IV V VI VII
I II III IV . V VI VII VIII IX X XI XII
I II III IV V ' vr, VII VIII IX X XI XII
NOTE: (A) JAPAN (Kawasumi, 1951)
(B) RUSSIA GEOFIAN (Medvedev, 1953)
(C) UNITED STATES MODIFIED MERCALLI
(Wood and Newman, 1931)
Figure 3.3 Correlation of three intensity scale
(from Trifunac and Brady, 1975)
67
Table 4 Abridged ModifiedMercalli Intensity Scale (from
Wiegel, 1983, p.90)
I. Detected only by sensitive instruments.
II. Felt by a few persons at rest, especially on upper
floors; delicate, suspended objects may swing.
III. Felt noticeably indoors, but not always recognized as a
quake; standing autos rock slightly, vibration like
passing truck.
IV. Felt indoors by many, outdoors by a few; at night some
awaken; dishes, windows, doors disturbed; motor cars
rock noticeably.
V. Felt by most people; some breakage of dishes, windows
and plaster; disturbance of tall objects.
VI. Felt by all; many are frightened and run outdoors;
falling plaster and chimneys; damage small.
VII. Everybody runs outdoors; damage to building varies,
depending on quality of construction; noticed by drivers
of autos.
VIII. Panel walls thrown out of frames; fall of walls, monu
ments, chimneys, sand and mud ejected; drivers of autos
disturbed.
IX. Buildings shifted off foundations, cracked, thrown out
of plumb; ground cracked; underground pipes broken.
X; Host masonry and frame structures destroyed; ground
cracked; rail bent; landslides.
XI. New structures remain standing; bridges destroyed;
fissures in ground; pipes broken; landslides; rail bent.
XII. Damage total; waves seen on ground surface; lines of
sight and level destroyed; object thrown up into air.
68
This means the MM intensity can not accurately portray the size and
severity of an earthquake. Nevertheless, it is still necessary to
have the correlation between the magnitude and the MM intensity.
When the Richter's earthquake magnitude was proposed in
1958, a correlation between the Modified Mercalli Intensity scale
and the magnitude scale was also proposed as shown below:
Magnitude 2 3 4 5 6 7 8
Maximum
Intensity III III V VIVII VIIVIII IXX XI
Cautions should be made when using the above correlation,
because it is "for ordinary ground conditions in metropolitan cen
ters in California" (Richter, 1958, P.353.). Chang and Krinitzsky
(1977) proposed the relationship between magnitude and epicentral MM
intensity as shown in Figure 3.4. Also the relation among magni
tude, epicentral distance and MM intensity was proposed as shown in
Figure 3.5.
When interpreting the Modified Mercalli intensity scale from
the magnitude scale, some limitations should be taken into
consideration (Wiegel, 1983, pp. 169170), i.e.
(1) Effect of Different Site Conditions. Under a specific
bedrock excitation, different site conditions may result
in different surface ground motions with different pre
dominant frequency or frequency content (Dobry and
Idriss and Ng, 1978). In other words, site conditions
may affect the severity of structure damages on which
the intensity scale is based.
MM INTENSITY, IQ
Figure 3.4 Relationships between magnitude scale and MM
intensity scale (from Krinitzsky and Chang,
1977)
EPICENTRAL DISTANCE,km EPICENTRAL DISTANCE
70
I
1000 
SYMBOL EARTHQUAKE M MM
o NEW MADRID, MO., 1811 7.5 XXI
_ . CHARLESTON, S.C.,1886 7.0 IXX
X NEW MADRID, MO., 1843 6.5 VII
 NEW MADRID, MO., 1895 VII
VIRGINIA, 1887 VIII
' A SOUTHERN ILL., 1968 5.6 VII
 V SOUTHERN ILL., 1857 VI
* k NEW MADRID, MO.,1963 5.1 VI
XVv
100 r
\ \ \X
Ox \\ X
A
X
X
II IV VI VIII
MM INTENSITY
500
100
10
II IV VI VIII X
MM INTENSITY
Figure 3.5' Relation of magnitude, MM intensity and epicentral
distance (from Krinitzsky and Chang, 1977)
71
(2) Focal Depth of an Earthquake. The distribution of the
released energy of an earthquake is a function of focal
depth. The deeper the focal depth, the less energy is
available to each unit ground surface area, i.e. less
. structure damages will result.
(3) Duration of shaking. At a specific released energy, the
duration of earthquake shaking may affect the damage
degree of structures. For instance, the liquefaction
potential of a soil deposit may be higher when subjected
to a moderate, longer shaking rather than subjected to a
stronger, but shorter shaking.
3.1.4 Relation Between Magnitude and Slipped Fault Size
Figure 3.6 (USBR, 1976) shows the relationships between
earthquake magnitude and the length of surface rupture of fault.
Although there exists good agreement among the proposed correlations
shown in Figure 3.6, one needs to notice that the type of an
active fault has significant effect on the magnitude of induced
earthquakes. For instance, with the same slipped fault size, a
reverse fault would cause a larger earthquake than a normal fault.
3.2 SourceSite Distance
Many terms are currently available to specify sourcesite
distance, e.g. hypocentral distance, epicentral distance and the
closest distance to a slipped fault. In EQGEN, the recommendation
of Jennings and Kanamori (1983) was adopted. The sourcedistance is
measured from the site to the closest point of the surface trace
of a slipped fault if the site is located within the circle with
72
LENGTH OF MAIN FAULT SURFACE RUPTURE, km
NOTE: (1) Bonilla (1967)
(2) Tocher (1958)
(3) Curve of. the result of leastsquare fit
(4) Greensfelder (1974)
(5) Housner (1969)
Figure 3.6 Relationships between earthquake magnitude and
the length of main fault surface rupture (from
Bureau, 1976)
73
diameter equal to the longitudinal dimension of the fault and with a
center at the center of faulting. Outside this circle, the distance
is measured to the center of the circle. As they pointed out, using
this definition would result in more consistent local magnitude
calculations.
3.3 Local Site Condition
When characterizing earthquake ground motions, it is very
important to consider the effect of local site conditions. It
should be realized, however, that there are still many uncertainties
resulted from the complicated effects of different geological
conditions and soil dynamic properties. Based on the past research,
the effects of local site conditions are summarized in the following
qualitative descriptions.
(1) Earthquake shakings with the same MM intensity produce a
higher peak ground velocity and peak ground displace
ment, and lower peak ground acceleration at an alluvium
site than those at rock sites (Trifunac and Brady,
1975). Some studies indicate that local site conditions
do not have significant effects on peak ground accelera
tions (Trifunac, 1976; Boore, Joyner, Oliver and Page,
1980; Joyner and Boore, 1981; and Campbell, 1981), and
have more important effects on peak ground displacements
and velocities (Trifunac, 1976). Though in some partic
ular cases, for instance, the earthquake ground motion
happened in Mexico City, Mexico (see Figure 1.1), the
shaking intensity seems to amplify at soft soil site.
This behavior, although further research is needed, may
74
be resulted from the existence of bowlshape boundary
between the soft soil deposits and the base rock at
Mexico City.
(2) As pointed out by Dobry, Idriss and Ng (1978), under a
given earthquake shaking, the frequency contents of the
acceleration timehistories, measured at rock and soil
sites, are quite different (see Figure 3.7). At a soil
site, the existence of tail shakings with a lower fre
quency content and a moderate intensity is explained as
a consequence of the amplification effect due to multi
ple reflection and refraction of the bodywave propaga
tion and the surface wave effects (Dobry et al., 1978).
The above effects also result in longer strong motion
duration observed in soil sites. Knowing the above
characteristics, it should be pointed out that, at the
present time, the EQGEN computer program does not have
the capability of simulating the time variation of the
frequency content of an earthquake ground motion. This
insufficiency should be taken into consideration, espe
cially under soil site conditions, where the frequency
content of an earthquake shaking may shift to a lower
frequency range as time goes.
In EQGEN, one of the following iocal site condition can be
assigned:
(1) Rock sites.
(2) Stiff soil sites.
(3) Deep cohesionless soil sites.
bo
Â§
3
w
ij
w
o
0.20
0.10
0.00
0.10
0.20.........
0 8 16 24 32 40 48 56
TIME, seconds
 Rock site 
Duration = 11.3 seconds
Probable arrival of direct
S wave from the source
(90 0.20 r 0.10
A a o H H 0.00
3
W 
w CJ 0.10
0.20
Duration = 26.5 seconds
Probable arrival of direct
S wave from the source
i 11 i i i i i i
0 8 16 24 32 40 48 56
TIME, seconds
 Soil site
Figure 3.7 Comparison of significant durations of a rock
site and a soil site (1971, San Fernando
earthquake) (from Dobry, Idriss, and Ng, 1978)
76
(4) Soft to medium clay and sand sites.
3.4 Duration of Strong Motion
In addition to the shaking intensity and frequency content,
the duration of an earthquake ground motion is another important
factor that greatly affects the severity of earthquakeinduced
structure damages, especially for a lightly damped linear system or
a nonlinear system (Dobry, Idriss and Ng, 1978; Vanmarcke and Lai,
1980). For an example, it is well known that the cumulative
characteristic of the pore water pressure generation in a soil
deposit under cyclic loading may cause liquefaction failure,
especially in prolonged shakings. In the past research, when
characterizing earthquake ground motions, simple parameters, such as
earthquake magnitude, sourcesite distance, local site condition and
peak ground acceleration etc. are usually used, and the duration of
an earthquake ground motion is usually considered as an implicit
parameter. Conversely, when specifying design earthquake ground
motion or performing artificial earthquake simulation, e.g. using
the GQGEN computer program, the duration of an earthquake ground
motion becomes one of the critical input parameters (Krinitzsky and
.Marcuson, 1982).
In EQGEN, the duration of a simulated accelerogram is
expressed by two parameters input by a user; i.e. (i) the number of
total acceleration data points; and (ii) the time internal. At
between two consecutive data points. It is required that At equal
to or less than 0.04 second; i.e. the sampling frequency should be
at least 25 Hz. This requirement is to prevent aliasing from
happening when analyzing a waveform function by using discretely
77
sampled data. The sampling frequency should be at least twice the
highest frequency involved in the waveform (Hovanessian, 1976).
Generally speaking, in engineering practice, the dynamic loading
with frequency higher than 10 Hz has little effect on structures
(Chang and Krinitzsky, 1977). Thus the time interval from 0.01 to
0.04 seconds can be used.
3.4.1 Definition of Strong Motion Duration
Before discussing the characteristic of the duration of an
earthquake ground motion, it is necessary to understand the defini
tion of earthquake strong motion duration, which is sometimes called,
the significant duration. According to the past research, there are
many different definitions of the strong motion duration. Some
studies define the strong motion duration by inspecting the degree
of contribution to the total seismic energy of ah earthquake
excitation; herein Parseval's theorem is applied (Trifunac and
Brady, 1975; Dobry, Idriss and Ng, 1978; and Vanmarcke and Lai,
1980). Some other studies apply the characteristic of rootmean
square acceleration to. define the strong motion duration (McCann.and
Shan, 1979). The definition proposed by Bolt (1973) is, however,
widely applied. It defines a so called "bracketed" duration as the
time interval of an accelerogram within the first and the last
occurrences of a prescribed shaking level, usually 0.05 g.
,3.4.2 Estimation of Strong Motion Duration
When executing EQGEN computer program, for the purpose of
selecting an appropriate strong motion duration, it is suggested
that the proposed correlation based on "bracketed" duration (with
78
acceleration greater than 0.05 g) be used, if Chang's model (1981)
or Seed's model (1982) in the strong motion database is used. When
using Trifunac's model (1976), it is recommended that the corre
lation proposed by Trifunac and Brady (1975) be referred to.
Concerning the correlation for estimating the "bracketed"
duration of an earthquake, a summary of typically, available
relations between magnitude and strong motion duration is given by
Dobry, Idriss and Ng (1978) as shown in Figure 3.8. Nevertheless,
the correlation shown in Figure 3.8 does not reflect the effects of
source.site distance and local site condition. As indicated by some
studies, under a given earthquake, the duration at a soil site is
about two times of that at a rock site, and the duration decreases
with increasing epicentral distance. The correlations including the
site condition effect were proposed by Chang and Krinitzsky (1977)
as shown in Figures 3.9, 3.10 and 3.11.
When applying Trifunac (1976) strong motion database to
generate artificial accelerograms, it is recommended to use the
e
correlation discussed below to estimate the strong motion duration.
Defining the strong motion duration based on the degree of contri
bution to the total seismic energy, Trifunac and Brady (1975)
proposed a correlation as shown in Equation (3.3) and Table 5.
SIGNIFICANT DURATION, seconds
79
NOTATION
# Housner (1965)
H Donovan (1972)
Bolt (1973); 0.05 g
A Bolt (1973); 0.10 g
Figure 3.8 Typical available relations between magnitude
and significant duration (from Dobry, Idriss,
and Ng, 1978)
DURATION, seconds (acceleration Â£ 0.05 g)
80
________________NOTATION__________________
Soil site; Chang and Krinitzsky (1977)
A Rock site; Chang and Krinitzsky (1977)
Soil site; Page et al. (1972)
O Rock site; Bolt (1973)
Figure 3.9 Comparison of nearfield, earthquake ground motion
duration (from Chang and Krinitzsky, 1977)
DURATION, seconds
(acceleration > 0.05 g)
70 
EPICENTRAL DISTANCE, km
Figure 3.10 Duration versus epicentral distance and
magnitude for rock sites (from Chang and
Krinitzsky, 1977)
DURATION, seconds (acceleration > 0.05 g)
82
Figure 3.11
Duration versus epicentral distance and magnitude
for soil sites (from Chang and Krinitzsky, 1977)
83
Table 5 Related Coefficients of Equations (3.3) and (3.4)
Acceleration a b c a A B E
Vertical 6.29 2.90 0.172 10.89 4.21 0.0672 7.25
Horizontal 4.88 2.33 0.149 10.67 2.92 0.0830 7.17
D as + bM + cR + a (3.3)
where D is strong motion duration in seconds; s = 0, 1, and 2 for
alluvium deposits, intermediate type rocks and basement rock,
respectively; M is earthquake magnitude; R is epicentral distance in
km; and a is one standard deviation.
Assuming the standard deviation increases with increasing
epicentral distance due to the disturbance happened in seismic wave
propagation, then a can be estimated by using Equation (3.4).
a A + BR + E (3.4)
where A, B and E values are shown in Table 5. In Trifunac and
Brady's study (1975), a limitation should be pointed out. Over one
half of the data from San Fernando, California earthquake of 1971
were used in formulating the above correlations. The formulation
may be biased toward the earthquakes near San Fernando area, and may
not be applicable to other areas.
3.4.3 Strong Motion Duration Database in EOGEN
In this research, in order to assist a user in choosing an
appropriate strong motion duration input, a strong motion duration
database is implemented in EQGEN computer program. The strong motion
duration database is divided into five groups as described below.
84
The first one is specially for the case when the Trifunac (1976)
strong motion database is chosen to retrieve a design spectrum; and
the others, where rock sites and stiff soil sites are classified as
hard sites; deep cohesionless soil sites and soft to medium clay and
sand sites are classified as soft sites, are for general purposes.
(1) For applying Trifunac*s model (1976'). Equation (3.3),
proposed by Trifunac and Brady (1975), is recommended to
estimate strong motion duration. Using this database,
an upper bound, a lower bound and a mean value of dura
tion are provided. It is suggested that several
accelerogram simulations are performed by assigning
different duration values that fall within the upper and
lower bounds.
(2) For nearfield (R < 10 km), hard sites. The curve
shown in Figure 3.12 is applied to estimate the strong
motion duration. This curve is basically modified from
the curve proposed by Donavan (1972), and it is close to
an upper bound of the proposed correlations shown in
Figures 3.8 and 3.9.
(3) For farfield (R 10 km). hard sites. The Equation
(3.5) shown below is applied to estimate the strong
motion duration:
Duration (1 + r )D (3.5)
where D is the time interval determined by using Figure
3.10 proposed by Chang and Krinitzsky (1977), and the
rvalue, as shown in Table 6, is determined according
85
Figure 3.12 Curve for estimating the earthquake ground motion
duration under nearfield (R<10 km), hard site
condition
86
to the userspecified earthquake magnitude. The deter
mination of the rvalue in Table 6 is discussed in
Section 5.2.1 (see Equation 5.1).
Table 6 Coefficient r of Equation (3.5)
Earthquake Magnitude, Ml r
4.0  5.0 0.938
5.1  6.0 0.930
6.1  7.0 0.716
7.1  8.0 .0.455
8.1  8.5 0.251
(4) For nearfield (R < 10 km). soft sites. As shown in
Figure 3.9, the solid curve B, proposed by Chang and
Krinitzsky (1977), is used to estimate the strong motion
duration.
(5) For farfield (R 10 km), soft sites. The curves
proposed by Chang and Krinitzsky (1977), as shown in
Figure 3.11, are implemented in the database.
It is noted that the implementation of the strong motion
duration database in EQGEN is related to the development of the
revised envelope function presented in Chapter 5.
3.5 Peak Ground Acceleration (For Chang's model only)
When using Chang's model (1981) in the strong motion data
base of EQGEN, it is required to input a peak ground acceleration
