On Friday afternoon of 28th, May 2004, at 12:38:46 Coordinated Universal Time, UTC, (17:08:46 local time scale) a strong earthquake occurred at the Western part of Mazandaran Province of Iran. The studies show that the location of the earthquake was in the South Western part of Alborz range, somewhere between Baladeh and Marzan-Abad villages. This earthquake is also known as Kojour-Firoozabad earthquake. The Alborz Mountains constitute a broad arch of parallel anticlines and Synclines are located at the Southern part of Caspian Sea (Fig. 1).
||The locations of the mainshock and aftershocks
The study of historical seismicity of Alborz (4th century B.C. to 1900 Data) shows that almost all the active areas in the historical period have also been seismically active in the recent century. These studies also indicate that many areas of this fold like Rasht, Rudbar, Lahijan, Fasham, Jajrud, Damavand, Amol, Babol, Sari, Behshahr and Gorgan have been experienced destructive earthquakes.
This study deals with the study of seismological parameters of Kojour-Firoozabad main event and using them for synthesizing process. Obtaining seismological parameters is one of the main steps in estimating more reliable strong motion of expected future events and therefore it can essentially affect the design of seismic resistant structures in the regions with lack of reliable data like the region under study. Minimizing the uncertainties that inherently exist in these parameters such as hypocenter location, focal mechanism (strike, dip and rake), rupture length and width is one of the main objectives of this study.
The three components of acceleration time histories of this earthquake are synthesized at one station far from causative source and compared with those of recorded data through their elastic response spectra. A Genetic Algorithm (GA) approach is developed for generating a number of seismological parameter set. Each set of these parameters is applied to simulate acceleration time histories, after which a fitness function is used for evaluation each parameter set by comparing the elastic response spectra of synthesized and recorded data. This comparison is done via three components, two horizontal and one vertical. This process is repeated through the GA several times to find the best set of seismological parameters regarding the fitness value.
There are a number of methods for calculating strong ground motions in literature;
theoretically, the first calculation of strong motions was made by Aki
(1968) and Haskell (1969). They used kinematic source
model propagating dislocation over a fault plane in an infinite homogeneous
medium which had five parameters: the fault lengths, fault width, rupture velocity,
final offset of dislocation and rise time (Irikura, 1983).
Point source method introduced by Brune (1970) and extended
by Boore (1983) is widely used in ground motion modeling
for engineering purposes. Despite the popularity of this approach, the point
source approximation is unable to characterize key features of large earthquakes
such as their long duration (Beresnev and Atkinson, 2002).
The finite fault method that is developed by Irikura (1983)
and Beresnev and Atkinson (1997) can solve some of the
problems of point source model such as the above aforementioned one; while it
is believed that both methods suffer from modeling the dependence of amplitudes
on the azimuth to the observation point, i.e., directivity effect and fling
step (Beresnev and Atkinson, 2002).
In many regions, the subsoil medium is not sufficiently known to simulate the
waves propagation in a suitable frequency band for earthquake-engineering
purposes (between 0.1 and 20 Hz) (Kohrs et al., 2005).
To overcome this problem, Hartzell (1978) proposed the
use of recordings of small earthquakes that are summed with time delays so as
to reproduce the rupture propagation effect. These small events represent all
the propagation effects between the source, path and the receiver and known
as empirical greens function (Hartzell, 1978).
So, they can properly apply modifications to earthquake ground motions, especially
in the high frequency range depending on the shear wave velocity and depth of
the soil (Lam et al., 2000). Another noticeable
point is that, although simple teleseismic body waves (P and SH direct waves)
can be precisely and deterministically estimated, most parts of the seismograms
can not and that is the reason of being helpful to use a smaller event, EGF,
to take the propagation between the source and the receiver into account (Vallee,
2004). It is worth noting that, because the synthesis approach is based
on Valleesite response from very small ground motion, it is only valid for linear
The EGF approach can be applied both in direct and inverse problems. In the
inverse problems, EGF helps to find the source information such as slip distribution
on fault (Plicka, 2003). Some investigators such as Mueller
(1985), Frankel et al. (1986), Mori
and Frankel (1990), Velasco et al. (1994),
Courboulex et al. (1997) and Plicka
and Zahradnik (2002) have used this application of EGF. In direct problems,
this technique is widely used in strong motion simulation of expected future
events (Plicka, 2003). Some seismologists like Irikura
(1983), Joyner and Boore (1986), Boatwright
(1988), Kanamori et al. (1993), Zeng
et al. (1994), Tumarkin et al. (1994)
and Hutchings et al. (2007) have utilized this
approach to synthesize strong motions. Both of these methods are adopted in
this study; the GA technique is used for estimating suitable source-path parameters
that are applied in the EGF approach to simulate time histories for the region
under study in a direct approach.
Among a number of researchers that have used EGF method, Wennerberg
(1990) used Joyner and Boore (1986) stochastic summation
approach to describe the source with the empirical Green's functions. Whereas,
probability distributions do not provide an extensive description of the phase
effects of the finite source and do not generate similar time histories respect
to the observed records. Hutchings (1991) used empirical
Green's functions to constrain propagation path and site response information
and proposed a range of simple kinematic rupture models to describe the source
in predicting strong ground motion for the full time history. This simple earthquake
model was tested by Jarpe and Kasameyer (1996) for the
1989 Loma Prieta earthquake with source parameters determined from different
studies of the earthquake. They figured out that for periods 0.05 to 0.4 sec,
the standard error between observed and predicted response spectra is less than
or equal to the errors from other methods. For periods between 0.5 and 2 sec,
this standard error was significantly less than regression methods based on
measured strong motion data. Hutchings (1994) reported
that Wu (1978) and Hutchings and
Wu (1990) examined the theoretical relation between rupture parameters and
synthesized seismograms using an exact solution to the representation relation,
in the frequency and time domain that uses empirical Green's functions. This
solution uses empirical Green's functions to obtain the Green's functions of
the representation relation, allowing for completely synthetic rupture models.
The computer Code EMPSYN originally written by Hutchings
(1988) was used in this study. It calculates synthetic seismograms by numerically
computing the discredited representation relation with Empirical Greens
Functions (EGF). The formula it uses has the following form (Hutchings
et al., 2007):
where, X and t are position and time in space relative to the hypocenter and
the origin time of the synthesized earthquake. N is the number of elements and
i denotes the values at an element. Ai is the area of an element
of the fault such that sum of all A is equals the total rupture area. μi
is the rigidity at an element. S(t´)i is the desired slip function
at an element analytically deconvolved with the step function. en(X,
t´)i is the recording of a small earthquake with effectively a step
source time function and interpolated to have a source and origin time at the
location of the ith element. t´ is relative to the origin time of
the synthesized earthquake. tr is the time that rupture reaches from
the hypocenter to the element. It is the integral of radial distance from the
hypocenter of the synthesized earthquake divided by the rupture velocity, which
can be a function of position on the fault. Me0i is the
scalar seismic moment of the source event and * stands for the convolution operator.
U has the same units as en (Hutchings et al.,
As mentioned in, this study presents a new method to study the seismological parameters of an earthquake and using them for synthesizing process. In fact, uncertainties inherently exist in seismological parameters are reduced. Sources of uncertainty in engineering safety problems are classified into two major groups known as aleatory uncertainty (due to inherent randomness) and epistemic uncertainty (due to limited knowledge of mathematical form of the model used). The aleatory uncertainty is the result of the estimated seismological source parameters and the factors affecting the wave propagation and site amplification and the epistemic uncertainty is arisen from the mathematical formulation of the model.
GA is one of the most effective techniques for optimizing/reducing the above
mentioned two families of uncertainties (Naeim et al.,
2004). In this study, GA is used to extract seismological parameters. Here,
the parameters are bounded in specific ranges. Furthermore, different sub-faults
dimensions and Brune correction (Hutchings et al.,
2007) are used. Incorporating the aforementioned parameters within predefined
ranges is an example of optimizing/reducing the first type of uncertainty. Utilizing
a number of predefined different sub-faults dimensions and Brune correction
as well, are samples of confronting the 2nd of uncertainty.
MATERIALS AND METHODS
This research was conducted on June 2007-May 2009 at Iran University of Science and Technology, Tehran, Iran.
Strong Ground Motion Data
The Main Shock Data
The main shock of the earthquake has been recorded by some of the accelerographs
of Iran Strong Motion Network (ISMN). The nearest station to the epicenter of
the earthquake was Poul (Kojour) Station, located at geographical coordinates
of 36.38 N and 51.72 E. The maximum acceleration recorded on the L (Longitudinal)
component equals to 0.305 g, on V (Vertical) component is 0.272 g and on T (transverse)
component is 0.170 g (These values are uncorrected). The distance of this station
to hypocenter is about 20 km.
Investigatin of different stations shows the effect of geological and topographical conditions. As an example, the Boumehen station which is located on alluvial, has recorded the earthquake with the PGA of 0.021 g and at a distance of 4 km away, the stone bedding Roudehen station has not recorded the earthquake. The Kondour station located on the heights (1890 m from the sea level and near to Karaj Dam), on the geographical coordinates of 51.11 E and 35.84 N, has not been able to record this shock due to its geological and topographic conditions.
Data of the Aftershocks
The Poul station acclerograph has recorded 18 aftershocks, after recording the main shock of the earthquake. The first aftershock has been recorded about 2 min after the main shock with the PGA of 0.012 g. The largest aftershock reported by IGTU is (Mb=5) occurred at 13:53:47 local time (9:23:47 international time scale) on 29th May 2004. This aftershock has been also recorded by Poul station accelerograph with the PGA of 0.083 g. The focal distance of aftershocks is 16-29 km from Poul (Kojour) station.
Site Soil Conditions
The strong motions were synthesized based on linear-soil-response ground
motion. The site soil condition of two stations, Poul and Baladeh, are not well
documented. The shear wave velocity is about 165 m sec-1 for upper
30 m depth of soil at Noshahr station (BHRC) (Table 1). The
selected stations are sufficiently far away from the causative fault, so their
recorded time histories were not greatly influenced by near field problems such
as directivity and fling step. As the small events have already included the
properties of the propagation path and local site effects, the EGF method can
be used for synthesizing earthquake at these stations. It is worth mentioning
that the geology beneath some of these sites may be such that it may respond
non-linearly during strong shaking; because the synthesis approach is based
on site response from very small ground motions, it is only valid for linear
material response. Therefore, soil nonlinearity is not taken into account in
this study. This is the same assumption made also by Hutchings
et al. (1997).
The Proposed Algorithm
One of the best methods to achieve a nearly optimal solution in problems
dealing with uncertainty and a lot of imprecise variables is known to be Genetic
Algorithm (GA). A genetic algorithm is a computer simulation of natural evolutionary
processes to solve search and optimization problems. Genetic Algorithms was
initially developed by John Holland over the course of 1960s and 1970s and was
extended by his student, Goldberg (1989a, b).
The applicability of genetic algorithm in engineering problems stems in the
efforts done by Goldberg. Genetic Algorithm as an optimization procedure has
also been successfully applied to structural design (Camp
et al., 1998). The GA technique has also previously been used in
seismology. Wu et al. (2008) have used GA to
determine the focal mechanism of some earthquakes in Taiwan during 1991 to 2005.
||The specifications of stations considered in this study
||The proposed GA based model
In this study, this approach is used to find the best parameters for simulating
strong ground motions in the region of interest. The 2004 Kojour Earthquake
M6.2 at the Poul station is simulated, using the EGF approach. The GA technique
is used for reducing uncertainties inherently exist in seismological source
parameters by comparing the 5% elastic response spectra of the synthesized and
recorded data at this station. Thereafter, to know the accuracy of the estimated
parameters, using empirical greens function, these parameters were used
to synthesize the three components of strong motion recorded at Noshahr station.
The main steps shown in the flowchart of Fig. 2 were taken
to precisely establish the proposed GA-based algorithm.
The initial population is generated randomly. Starting with this population, the main GA loop is initiated. Each chromosome in the population contains the seismological parameters. These parameters are objectives of the GA to optimize. In each step, the population individuals are evaluated with the fitness function. The fitness calculation process uses the parameters in the chromosome to simulate the strong ground motion. The parameters are fed into the EMPSYN program, which operates based on the EGF approach. The output of this stage is the acceleration time history. The time histories are analyzed to generate the elastic response spectra. The fitness value will be inversely relative to the error between the spectra of the synthesized and recorded acceleration time histories. The use of EGF method and GA technique is not new in seismology; but this study uses these two approaches together with the proposed method which is novel.
In general, 12 parameters of the earthquake are optimized in this process,
which are listed as follows: 3 parameters for hypocenter location (Table
2), 4 parameters for distances of edges of the fault from hypocenter (Table
3), 1w parameter for the sub-fault dimension, 1 parameters for the magnitude
of the earthquake (Table 2) and 3 parameters for the focal
mechanism (Strike, Dip and Rake) of the ground motion (Table 2).
||Comparison of hypocenter location and focal faulting type
estimated through this study and those of other references for Kojour Earthquake
||Estimated Hypocenter location and causative fault dimensions
||Binary representation of a single parameter cooperated into
the GA procedure
For all these parameters upper and lower bounds are assumed. The bounds for
values of the parameters are determined according to those of other references.
These bounds are used for chromosome representation illustrated in the following
A sample binary representation of a bit string that is used as a part of
the proposed chromosome design is demonstrated in Fig. 3.
Each chromosome consists of several bit string parameters cooperated into the
GA procedure using this representation.
In Fig. 3 g1j˜gjnj represents the values of jth parameter. nj is the number of bits representing the jth parameter. The source parameters such as longitude, latitude and depth of hypocenter, fault plain dimensions and orientations (Strike, Dip and Rake) surveyed in this study.
The binary strings are converted to real numbers using Eq. 2.
||Phenotype of bit string
||Maximum value of the parameter
||Minimum value of the parameter
||Decimal value of bit string
||Length of bit string
The mathematical expressions to evaluate the fitness of individuals in a
generation used for the two horizontal components L, T and vertical component
V of Poul station are as follows:
||L component error which is the sum of the squares of the differences
between synthesized and observed acceleration response spectra at any period
||T component error which is the sum of the squares of the differences between
synthesized and observed acceleration response spectra at any period
||Z component error which is the sum of the squares of the differences between
synthesized and observed acceleration response spectra at any period
The three-point crossover is used as the crossover operator. After the two
parent chromosomes are selected, three random numbers, less than or equal to
chromosome length, are generated and the bits associated with these numbers
are used as the boundaries for dividing chromosomes into four sections. Thereafter,
the second and forth sections of the two selected chromosomes are exchanged.
Figure 4 shows the crossover operation process.
A simple mutation operator is used, in which a random number less than or equal to the chromosome length is generated and the value of associated bit is flipped.
The selection operator is used to generate the next generation. The bigger
the fitness value of the chromosome, the more chance to be selected. The elitism
is used to save the best results of the previous generation. This means that
the best chromosome of previous generation (the chromosome with the biggest
fitness value between all of the population) is randomly replaced into one of
the chromosome of population in the new generation.
Termination Criterions and Parameter Settings
The genetic algorithm will be terminated if the number of generations that
population is no longer improved meets a specific number (ng).
The proposed genetic algorithm was applied on the simulation process of EMPSYN computer software to find adoptive seismological parameters of the Kojour earthquake. Each chromosome of the population was considered as a series of seismological parameters that used for synthesizing earthquake. Population size was set to 200, crossover rate to 0.8 and mutation rate to 0.08 according to the experimental values that is often proposed for these rates. Three-point crossover, as mentioned above, together with a roulette wheel selection was taken into use. The number of bits for each chromosome was set to nj = 8. The maximum number of generations used as the termination criterion was ng = 50 and finally, as discussed before, 12 parameters are cooperated into the algorithm to be optimized.
RESULTS AND DISCUSSION
The three components of Kojour-Firoozabad earthquake, Mw6.2, at Poul (Kojour)
station were simulated using the widely used EGF approach. This station is relatively
20 km far away from the causative fault. The GA was used as a tool for inversion
solution of the problem. The seismological source parameters such as hypocenter
location, focal mechanism (strike, dip and rake), rupture length and width,
explained before, were estimated through a Genetic Algorithm (GA). To do this,
the proposed genetic algorithm was used to optimize the seismological parameters
for minimizing the differences between elastic response spectra of the synthesized
and recorded data at Poul station. Some of these parameters have been previously
reported by some other references such as USGS and HRVD. The estimated seismological
parameters are compared with those of other references, previously studied the
2004 Kojour-Firoozabad earthquake (Table 2). The basic theory
of simulation procedure used in this study was that of Hutchings
(1994). A computer FORTRAN program, EMPSYN, originally written by Hutchings
(1988) was used. The Kostrov slip function with variable rise time and constant
stress drop has been incorporated in the model. Figure 5a-c
show the observed and simulated acceleration time histories at Poul station.
Furthermore, the elastic response spectra of simulated strong motions at this
station were compared with those of recorded data as shown in Fig.
6a-c. Also the sub-faults dimension parameter was determined
to be 0.5 km, after the GA process.
Validation and Prediction Process
In order to validate the estimated source parameters and find the accuracy
of the proposed method, the main event was synthesized at Noshahr station located
at about 22 km far away from the epicenter. The source parameters acquired from
GA analysis at Poul station were used to synthesize the strong ground motion
at Noshahr station. The inverse problem was not solved at this station and a
forward simulation was made incorporating the estimated source parameters and
using this site specific EGFs. The satisfactory match of elastic response spectra
and the observed data at these two stations (Poul and Noshahr stations) confirms
the reliability of the technique in estimating the seismological parameters.
Figure 7a-c shows the elastic response spectra
of synthesized strong motions at Noshahr station along with those of recorded
Also, using the source parameters resulted from GA analysis, the main shock
at Baladeh station was predicted at which because of technical destroys, there
was not any recorded data during the main shock and only the aftershocks was
recorded at this station. Figure 8a-f shows
the three component acceleration time histories and response spectra at Baladeh
||(a-c) Recorded and synthesized strong motions at Poul st
As already mentioned, the values of previously estimated parameters were used as criteria for the upper-lower bounds of parameters used in GA. Table 2 shows the estimated source parameters such as Longitude, Latitude and Depth of hypocenter, fault plain orientation (Strike, Dip, Rake) and those have already surveyed by other references. Table 3 shows the obtained Hypocenter location and Fault dimensions.
The error value for each chromosome is calculated to find minimum error in each generation. Figure 9 shows the minimum error value, introduced before, in each generation which contains 200 chromosomes. As it can be seen, the differences between the estimated and observed data are reduced as the generation number is increased, confirming the applicability and effectiveness of the technique in estimating/modifying the aforementioned seismological parameters. As is shown in Fig. 9, the curve becomes asymptotic to a horizontal line expressing the adequacy of synthesizing/modifying process.
In this study, an approach is proposed by which it is possible to predict the
three components of strong motion time histories at selected station as a result
of inversion solution at other stations. The EGF method is used to synthesize
acceleration time histories and a GA approach is developed throughout the paper
to optimize the seismological parameters such as hypocenter location, focal
mechanism (strike, dip and rake), rupture length and width.
||(a-c) Comparison between elastic response spectra of synthesized
and recorded data at poul station
The optimization is objected to minimize the differences between the elastic
response spectra of the synthesized and the recorded data. The applicability
of the technique was shown through estimating the aforementioned parameters
at Poul station and validating the estimated parameters by simulating the three
components of acceleration time histories at Noshahr station located at about
22 km far away from the epicenter. Also, the acceleration time histories of
this earthquake at another station, Baladeh, were predicted at which the strong
motion was not recorded during the main earthquake. All the selected stations
are sufficiently far away from the causative fault, so as the results are not
greatly influenced by the near source problems such as directivity and fling
||Comparing the predicted three components of strong motion
and the recorded data at Noshahr St., in the form of elastic response spectra
It is notable that, in spite of the fact that the estimated elastic response
spectra corresponding to six components of synthesized strong motions are comparatively
in a good agreement with those of observed data, it is not claimed that the
proposed technique ends with the results free from uncertainty, but it is an
adaptable approach for the regions with lack of reliable earthquake data, like
the one used in this study. The GA approach and validation process in this research
work show the applicability of the proposed approach for synthesizing the happened
earthquake in the region under study, which is the aim of this study. As the
GA and validation procedures have been done on the elastic response spectra,
the resulted time histories may be different from the real one to some extent.
In fact, the proposed approach can be used to find the acceleration response
spectra for regions with lack of recorded time histories during earthquakes.
||(a-c) Proposed strong motion and (d-f) acceleration response
spectra at Baladeh station (L., T. and V. Components)
||Minimum Error value in each generation
So the structural behavior can be estimated during happened earthquakes. Also,
this technique can be utilized to extract time histories that are compatible
with the synthesized response spectra for time history analysis of structures
during earthquakes in the regions under study. The estimated values of parameters
are in fact the adopted values of these parameters for synthesizing process,
where their suitability was checked during validation process. Also the use
of EGF and GA methods is not new in seismology; the proposed model of combination
of these two approaches is new and has not been presented before.
Many thanks to Professor Hutchings for kindly offering the valuable simulation codes and Dr. Yasin Fahjan for giving us the EMPSYN FORTRAN Codes.