Bulletin of the Seismological Society of America, Vol. 95, No. 4, pp. 1447–1456, August 2005, doi: 10.1785/0120030183

Estimation of Site Response in Time Domain Using the Meyer–Yamada
Wavelet Analysis
by Gülüm Birgören and Kojiro Irikura

Abstract In most applications of the stochastic Green’s function simulation, amplitude characteristics of seismic responses at target sites are taken into consideration,
but phase information of them is often disregarded or just assumed to be random.
However, when seismic waves from the source pass through different geologicalstrata, they take multiple paths and, as result, the arrival times of different seismicphases to the receiver vary. Especially in basin structures, seismic waves trapped in
unconsolidated deposits or secondarily induced at the basin edges elongate the duration of ground motion. In the present study, a simple empirical method using theMeyer–Yamada wavelet procedure is proposed for the phase-dependent site response
estimation in time domain. The method reveals the amplitude and phase effects ofextended duration of ground motion on site response by utilizing the wavelet transform. Because the wavelet transform yields a localization of a signal in both frequency and time domains, the phase information of the signal stays intact. Application of this method has been performed using the aftershock records of the 1999
Düzce, Turkey Earthquake (Mw 7.1). Amplitude spectra of absolute site responses
were also calculated by a spectral calculation method at each station for comparison
purpose. Good agreement between the amplitude spectra calculated with the wavelet
method and those of the spectral calculation method supports the efficiency of the
proposed method.
Introduction
Conventional methods for estimating the site response
are mainly based on the calculation of the amplitude spectrum of the ground motion and do not effectively reflect the
information about the duration (Borcherdt, 1970; Andrews,
1986; Iwata and Irikura, 1988; Fletcher and Boatwright,
1991). However, in sedimentary basins the impact of wave
propagation effects on ground motion can be more significant than that due to the surfacial layers. The long-period
ground motions recorded at the sites underlain by deep basins may be associated not with the direct body waves butwith secondarily generated surface waves, which elongates
the ground motion. These long-period surface waves cancause large amount of damage in flexible structures. Therefore, extended duration of the ground motion is an inevitableconsequence of the wave generation in sedimentary-filledbasin structures and has to be considered together with site
amplification for realistic simulation of the ground motion
for assessment of a future earthquake.
Arrival time information of seismic waves to a site is
kept in the phase spectrum of a signal. Derivation of theunwrapped phase spectrum gives the group delay time versus frequency behavior of the site. Sawada (1998) has suggested that the group delay time calculation is an efficientmethod to understand the phase property of the seismic motion. Following the method of Sawada (1998), Beavual etal. (2003) showed the variety of the lengthening of theground motion at basin stations using observed and simulated data. They attempted to estimate synthetic signals at
stations located inside a basin by modifying the amplitude
and phase spectra of the signal at a reference site. However,this method demands very close distance between the reference site and the sites of interest such that the source andpath effects cancel out.
The method presented here focuses on the calculationof the site response at basin stations, taking into consideration the elongation of the ground motion due to basin effects,
and proposes an empirical method to estimate the phasedependent site response directly in the time domain utilizing
the Meyer–Yamada wavelet analysis. Because the method
estimates the site response at a station directly from its time
domain data by removing source and path effects at eachwavelet level, the phase information is automatically kept in
the analysis; therefore, the method necessitates neither anear-distant reference site nor additional analyses to estimate
the elongation of ground motion.
The 1999 Turkey earthquakes (17 August, Mw 7.4 and

1447

1448

G. Birgören and K. Irikura

12 November, Mw 7.1) caused extensive damage, in particular, in Düzce City and several towns located in the Düzce
Basin. To perform a realistic ground-motion simulation tounderstand the nonuniform damage distribution in the basin,site amplification and extended duration of the ground motion due to basin-generated waves should be computed. To
achieve this objective the proposed method was applied tothe aftershock data set of the 1999 Düzce, Turkey, earthquake.

the user to modify the signal in each frequency band. Inverse
wavelet transform reassembles the signal from modified coefficients.
For a discrete time series f(t) with sampling rate Dt andthe number of data points N (N  2n where n is an integer)
inverse discrete wavelet transform is defined as (Yamada
and Ohkitani, 1991);
N1 2j1

f(t) 

Time Domain Representation of Site Response
Discrete Wavelet Transform
Because it has been used for several decades, a nonstationary signal can be treated as the combination of cosine
and sine waves, which comprises the basis of the Fourier
analysis. The signal can be more accurately characterized by
the usage of basis functions, in other words, wavelets (Daubechies, 1990; Mallat, 1989; Meyer 1989). The wavelet transform decomposes a signal into its high- and low-frequency
components by using shifted and scaled versions of themother wavelet. Orthogonality is the most desired property
in any signal analysis. Besides, localization or compactnessof the wavelet is required to avoid the redundancy of thesignal. Hence, in most of wavelet transform applications inseismology, the discrete wavelet analysis is preferred because the wavelet function of the discrete wavelet analysis
satisfies the orthonormality condition (Meyer, 1989). A timeseries f(t) can be represented by an inverse discrete wavelet
transform asf(t) 

兺j 兺k ␣j,k Wj,k (t)(1)

Wj,k (t) represents the analyzing wavelet (basis function), t
is the time, j and k are the integer values denoting the spatial
scale and position of the wavelet, respectively (Yamada andOhkitani, 1991). Wavelet coefficients (␣j,k) are calculated byconvolving the time series by the analyzing wavelet as␣j,k 兺 兺 ␣j,k Wj,k (t) .

Among several analyzing wavelet definitions, the Meyer–
Yamada procedure (Meyer, 1989; Yamada and Ohkitani,1991) has been employed in the present study (Fig. 1). The
priority of the Meyer–Yamada wavelet is that amplitudes of
the wavelet spectrum approximately correspond to the amplitudes of the power spectrum of a signal at the jth scale.
Ej 兺k |␣j,k|2

⬁

Wj,k (t) * f(t)dt

2j
2j2
ⱕfⱕ
3Td
3Td

(2)

(3)

Because wavelet transform is a linear technique wavelet coefficients are modified in the time domain even though they
are associated with frequency levels (wavelet levels). In
other words, data operations can be performed with wavelet
coefficients of the corresponding frequency band instead of
the signal itself. Hence, this property of the transform allows(6)

Td is the signal length in seconds (Td  N Dt). Note that
the Meyer–Yamada Wavelet is nonzero only between thefrequency ranges given at (6). Wavelet coefficients havehigh resolution in time and low resolution in frequency as
the wavelet level increases and vice versa. The wavelet decompositions of a strong-motion accelerogram and its power
and wavelet amplitude spectra are illustrated in Figure 2 as
an example.

Asterisk denotes the complex conjugate.
The analyzing wavelet is derived by translation and dilatation from a mother wavelet function (W) as following;
Wj,k (t)  2j/2 W(2j t  k)(5)Frequency range of the jth scale is defined as

⬁

冮

(4)

j0 k0

Figure 1.

Meyer mother wavelet in the time domain.

Estimation of Site Response in Time Domain Using the Meyer–Yamada Wavelet Analysis

1449

Figure 2.

(left) Wavelet decompositions of a strong-motion accelerogram. (right)
Power spectrum of the same accelerogram is shown in light gray. Amplitudes of eachwavelet level (each frequency band) are drawn at geometric mean frequencies (2j1/
3Td) with filled circles. Thick black line represents 20% logarithmic smoothing of thepower spectrum.

Site Response Calculation with Wavelet Method
The Fourier amplitude spectrum of an observed seismic
signal is represented as the multiplication of source, propagation path, and site response (Iwata and Irikura, 1988).
Oij ( f )  S( f )i P( f )ij G( f )j(7)

where S(f)i is the seismic source of the ith event, P(f)ij isthe propagation path effect from the ith event to the jth site,
and G(f)j is the site response at the jth site.
The simplest empirical way to calculate local site responses is to remove the source and propagation path effects
from the observed record by taking the spectral ratio between the site of interest and a reference site for the same

event, as proposed by Borcherdt (1970). This method is effective if the distance between the stations is much smaller
than the hypocentral distance.
The empirical estimation of the site response in time
domain is performed essentially in the same way as the site
amplification is calculated from the spectral ratio technique.
This method utilizes the discrete wavelet transform to localize a seismic signal in both time and frequency domains
simultaneously and enables retention of amplitude and phase
properties of the observed signal. The method is based on
the modification of wavelet coefficients of seismic signals at
each discrete wavelet level to yield wavelet coefficients ofsite response. In the method each horizontal component of
the seismic signal, starting from S-wave onset to the endof the data, is decomposed into its wavelet levels, wavelet

1450

G. Birgören and K. Irikura

coefficients at each level are divided by the corresponding
source and propagation factors as follows:␣j,Gk, i ␣Oj, k, i| S(x)j, i  P(x)j, i |

(8)

O
(␣G
j, k, i) and ␣j, k, i are the wavelet coefficients of the site response and the wavelet coefficients of the signal at the jthscale of the ith event The denominators S(w)j,i and P(w)j,i in(8) are the source and propagation path scalar terms at thegeometric mean frequency (2j1/3Td) at the jth scale of theith event, respectively. In the calculation, the wavelet scales
whose frequency range is smaller than the noise levels of
signals are suppressed by forcing the related coefficients tozero (i.e., ␣j,k,i  0). Schematic illustration of the calculation
method is given in Figure 3.
The source spectrum required in (7) and consequentlyin (8) is calculated from the estimation of the low-frequency
flat level (X) and corner frequency (fc) of the displacement
spectra with w2 point source assumption at reference stations(Aki, 1967; Brune, 1970).

S( f )i 

Xi /2
f 2
1
fci

冢冣

(9)

where f is the frequency. Flat level is divided by 2 to accountfor the free surface effect. Regarding propagation path ef-

Figure 3.

fects, radiation of S waves from the source are taken into
account. Propagation path effect has a general equation asP( f )ij 

1 pRij f /Q( f )VS
e
Rij

(10)

where R is the hypocentral distance between the ith eventand the jth station, Q(f) is the quality factor, and VS is the
S-wave velocity of the medium.
Equation (8) is valid for the small earthquakes in shortdistance with similar source durations. If this condition is
fulfilled, phase delay due to the source effect and the propagation path effects can be assumed negligibly small compared with those due to the site effects. Finally, the averagecoefficients of the site response at a station are calculated asn␣Gj, k, ini1

␣j,Gk  兺

(11)

where ␣ji ⬆ 0, n  {number of events whose ␣j,k,i ⬆ 0}.
The modified average wavelet coefficients calculated
from the preceding operations together with wavelet function yields the site response in the time domain through theinverse wavelet transform defined in (4). According to themean value theorem of statistics, the noise should disappear
as the number of records used in the analysis increases. Thus,
the stacking technique in (11) enhances common site effects(coherent part) and reduces noises (incoherent part).

Schematic representation of the site response in the time domain by the wavelet method.

Estimation of Site Response in Time Domain Using the Meyer–Yamada Wavelet Analysis

Probably one of the most important priorities of the
wavelet analysis over the Fourier analysis is that it can yieldtime and frequency information simultaneously in accordance with the uncertainty principle, where the Fourier analysis can yield only the frequency information. Current anaysis can also be done in the complex domain, because it
does not explicitly use the time-frequency calculation. However, Fourier coefficients are not stable for each event recorded at a station and the smoothing functions used for the
analysis do not always guarantee the orthonormality condition. In addition, wavelet analysis makes it possible to explain the amplitudes with a limited number of positions
rather than using all data points.
Spectral Calculation Method (SCM)with Reference Event
Absolute site-amplification factors are estimated by aspectral method to evaluate the validity of the site-amplification calculation with the wavelet method.
The site-amplification factors are calculated by removing the source and propagation path effects from the observed record. A diversion from the technique of Borcherdt(1970) arises from the fact that, the source spectra, estimatedfor the events recorded at the reference site, and propagation
path spectra of all events are calculated empirically according to (9) and (10) as long as the source spectra are explainedwith the w2 model. Therefore, the absolute site effects canbe estimated at each station with respect to the theoreticalS(f) and P(f) spectra.
The definition of average site response at the jth station
calculated for N events is expressed as:N

Gj( f ) 

O (f)

/N
兺 j
i1 Si( f )Pij( f )(12)

In Figure 4 the calculation method is expressed in a schematic way. The calculation method will be addressed as

1451“the spectral calculation method (SCM) with reference
event,” hereafter.
A Case of Düzce Basin
During the 1999 Kocaeli and Düzce earthquakes, urban
areas located in the Düzce Basin, in particular, Düzce andKaynaşlı, suffered from heavy structural damage. The damage is mainly associated with the malfunction of the structural system as well as the effect of the soil characteristics.
Borings drilled up to 15 m and seismic refraction measurements performed in Düzce City and Kaynaşlı Town show
similar characteristics. The areas underlain by the Quaternary deposit consists of silt, clay, and sandy gravels. The
groundwater levels estimated from borehole logs vary between 2 and 7 m (Kayabalı et al., 2001; Özdemir et al.,2002). A microtremor study carried out by Kudo et al.
(2002) in Düzce City indicates that Düzce is a sedimentarybasin with 1.5 km depth, consisting of thick and soft sediments at surface (Vs ⬃ 260 m/sec) and intermediate depths(Vs ⬃ 460–510 m/sec).
The site responses at the observation stations inside andoutside of the Düzce Basin were analyzed by using thestrong-motion accelerograms of the 1999 Düzce, Turkey,
earthquake aftershocks, recorded by temporary networks deployed by Kandilli Observatory and Earthquake ResearchInstitute of Boğaziçi University (EQE) and LDEO of Columbia University (CU). The strong-motion stations are listed inTable 1. Fourteen shallow earthquakes of duration magnitude, 3.2 ⱕ Md ⱕ 5.0 with hypocentral distance of up to
40 km were used. For the estimation of site-amplification
factors, the station CU1058 was chosen as the reference site
based on the topographic and geological data.
Figure 5 shows the location of the stations and aftershocks together with the bandpass-filtered (0.5–2.0 Hz) velocity waveforms of north–south component of the aftershock no. 12, starting from P-wave arrival. Even though therecording time is shorter in the basin stations due to misadjustment of instruments, the ground-motion differences

Figure 4. Schematic representation of the calculation method of site amplificationby the spectral calculation method with reference event (SCM).

1452

G. Birgören and K. Irikura

Table 1
Strong-Motion Recording Stations and Their Specifications
Station Name

Location

CU1058CU1059CU1060CU1061
CU1062
CU0362
CU0375
CU0531EQESAZ

40.75N–31.06E
40.75N–30.87E
40.78N–30.63E
40.72N–30.79E
40.72N- 30.82E
40.67N–30.67E
40.74N–30.87E
40.70N–30.85E
40.84N–31.05E

EQEBALEQEGON
EQEFCM

40.78N–31.10E
40.82N–31.21E
40.83N–31.19E

Instrument Specifications

Reftek-72A-08 24 bit
Terra T 320(full scale  2 g)
S.R: 100 Hz*

Kinemetrics SSA 12(full scale  0.5 g) S.R: 200 Hz
GeoSys GSR-12
Terra T 320(full scale  2.0 g) S.R: 200 Hz*S.R., sampling rate.

between basin stations and the reference station are quite
visible. The waveforms get complicated inside the basin andstrong-motion durations get longer. The extended durationof shaking at basin stations may be the indicator of the basininduced surface waves.
Noise-amplitude spectra of the records were calculated
by using the 5-sec time window before P- wave arrival. Records having signal-to-noise ratios lower than a factor of 2
were removed from the analysis. At most CU stations thesignal-to-noise ratio is above the factor of 2 between 0.4 and
10 Hz. Most of the records at the EQEGON and EQEFCM
stations either have very short duration before P-wave arrival or started after P-wave onset. Therefore, the noise levels of those earthquakes are biased at lower-frequency
ranges, (particularly at f ⱕ 0.5 Hz) Acceleration spectra were
calculated for a 5-sec window of the records, starting 1 sec
before the S-wave arrival, as the horizontal vector summation of north–south and east–west component’s spectra andtransformed in displacement with double integration. Amplitudes were corrected for 1/R geometric spreading, R is the
hypocentral distance. The low-frequency flat level (X) andcorner frequency (fc) of the displacement spectra of 14 earthquakes recorded at the reference site was estimated by a
Genetic Algorithm (Holland, 1975) assuming w2 source
model. The seismic-moment values (M0) were derived fromM0  4pqVS3 X/Rh␾

(13)

S-wave velocity of the medium (Vs), density (q), andradiation-pattern coefficient for S waves (Rh␾) are taken as
3.2 km/sec, 2.7 g/cm3 (Mindevalli and Mitchell, 1989), and
冪2/5 (Andrews, 1986), respectively. The source parametersof the 14 earthquakes were given in Table 2.
Moya (personal comm., 2003) estimated the quality factor of the Düzce area as Q  154f 0.15 through a spectral
inversion method by using a larger aftershock data set recorded at the same stations. Propagation path spectra were
calculated with equation (10), based on the Q value mentioned above.

Figure 6 shows the average site response results together
with site response results of all individual earthquakes usedin the calculation for the frequency range of 0.5–10 Hz. The
data preparation for the calculation of phase-dependent site
responses were done in the following order: All the records
were corrected to zero baseline, cosine tapered starting from
1 sec before the S-wave onset until the end of data, and zero
padded to equalize the number of wavelet levels of each
earthquake at a station. Decomposition levels were chosen
according to the noise level of records. In this aspect, thedata selection differs from the SCM. For instance, the data
having noise contamination beyond 0.7 Hz were not included in the estimation with the SCM; however, it canbe included in the wavelet method suppressing the lowfrequency level coefficients and including higher-frequency
level coefficients in the estimation. The largest event (Md 5)
was not used in the wavelet calculation to keep the source
duration of events similar. Time domain representations of
the site responses, bandpass filtered between 0.5 and 5 Hz,
are drawn next to the stations in Figure 7. Low-frequency
level is due to the reliability of the observed data. Averagesite response spectra of the stations calculated from bothmethods were compared in Figure 8. Note that the site response estimation with the SCM comprises mainly the Swave portion of the data where the wavelet method covers
all available waves arriving after S-wave onset.

Results
Site response estimations using spectral calculation andwavelet methods are presented. Some of the remarkable
points can be summarized as follows:• Flat site amplification in low-frequency range and fundamental frequency observed in higher-frequency range at
the station CU1058 imply harder soil site. This result supports the appropriateness of the reference site selection.
• Two close basin stations EQEFCM and EQEGON show
similar response shape and amplification levels. The stations EQEBAL and EQESAZ located in the edge and center of the basin have amplification levels as high as 10 at
about 2 Hz and 1 Hz, respectively. This is the largest site
amplification level estimated for basin stations at midfrequency range (0.5–2 Hz)• Time domain representation of site response shown in Figure 7 reflects the difference in soil characteristics inside
the basin and outside the basin. The stations located on
harder soil (CU1058, CU1060) show a site response almost like an impulse function; however, stations inside thebasin, in particular, EQEBAL and EQESAZ, exhibit
longer and more complex response functions and are richin high-frequency contents.
• At stations CU1059 and CU0375 the difference in amplification levels is significant even though they are locatedapproximately 500 m from each other. The station
CU0375 has the highest amplification levels with shortEstimation of Site Response in Time Domain Using the Meyer–Yamada Wavelet Analysis

1453

Figure 5.

Location of the strong-motion stations (filled triangles) and aftershocks(filled circles). North–South component strong-motion velocity records of the earthquake
no. 12 (99/12/20/03:27 Md 4.2, depth  10.0 km), bandpass filtered between 0.5 and2 Hz are drawn next to the stations to show the variation of ground-motion records.

Table 2
Source Parameters of Reference Events
No.

Event Date(yy/mm/dd/h:min)

Depth
(km)

Moment
(N m)

fc
(Hz)

Magnitude
(Md)

1
2
3
4
5
6
7
8
9
10
11
12
13
14

99/11/19/14:01
99/11/19/19:59
99/11/21/04:31
99/11/21/04:33
99/12/03/12:40
99/12/03/12:46
99/12/13/19:13
99/12/14/17:38
99/12/15/00:12
99/12/17/05:29
99/12/17/12:12
99/12/20/03:27
99/12/20/03:41
99/12/23/14:40

10.0
9.0
10.0
7.0
5.0
3.0
7.0
8.0
8.0
8.0
3.5
10.0
7.0
5.0

3.9E15
2.6E166.0E15
4.8E14
3.5E15
4.4E14
5.7E155.2E14
3.6E142.4E14
4.5E14
4.3E15
6.1E14
5.0E14

1.7
0.9
1.6
2.1
2.1
2.3
1.2
1.9
2.4
2.5
2.2
2.3
4.1
4.5

4.1
5.0
4.2
3.5
3.7
3.3
4.3
3.6
3.3
3.2
3.5
4.2
3.4
3.2

duration. This reflects a strong velocity contrast between
thin uppermost sedimentary layer and engineering basement layer under the station.
• The average site amplification spectra estimated by the
wavelet analysis method show good agreement with theestimates obtained by the SCM.

Conclusion and Discussion
An empirical method has been proposed to estimate thesite response in time domain. The main asset of the method
is that it allows the calculation of site response without ig-

noring the phase properties of a seismic signal. The method
was applied to the aftershock data set of the 1999 Düzce
Earthquake. The fairly good fitting of both amplitude levels
and shape of the site response calculated with the SCM andthe wavelet method supports the validity of the proposed
method. The efficiency of the method should be tested withthe stochastic Green’s function simulations and comparedwith the damage distribution observed at the Düzce Basin.
The success of the method highly depends on the existence of the reference site to define source effects with highaccuracy. S-wave onset picking is also a major factor affecting results, because misalignment of records highly increases the error in estimation of phase information. Moreover, the method requires a high signal-to-noise ratio tocalculate low-frequency site response realistically.
Probably the most important limitation of this approachis related to the calculation of source and propagation path
spectra. The method can not be performed successfully without accurate estimation of the seismic moment, corner frequency of reference events, and inelastic attenuation factorof the medium. Azimuthal dependence of the site response
was not considered in the present study because of limited
number of the earthquakes. For further estimation this factor
should be taken into consideration, in particular, for the basin edge stations.
Beside these, many extensions of this research still deserve further consideration; the main assumption followedin this research is that site response mainly consists of thecoherent part of the signal (as in equation 11). Spectra fittings shown in the application confirm this assumption. A
full development of the site response considering both co-

1454

G. Birgören and K. Irikura

Figure 6.

Average site response results versus frequency calculated with the spectral
calculation method with reference event (SCM). Gray lines show the events used for
calculation.

Figure 7. Time domain representations of the site response bandpass filtered between 0.5 and 5 Hz. Upper and lower waveforms are the north–south and east–west
components of the site response at each station, respectively.

1455

Estimation of Site Response in Time Domain Using the Meyer–Yamada Wavelet Analysis

Figure 8.

Comparison of amplitude spectra of site response calculated with spectral
calculation (in black) and wavelet methods (in gray).

herent and incoherent parts of the signal is beyond the scopeof the present study. For further direction of the research,
the incoherent part of site effects should be calculated basedon amplitude spectra, assuming random phases to prevent
underestimation of the site response.
Another issue raised by the current results is the overestimation of the high-frequency ground motion inside thebasin, possibly because of the nonlinear effect of the soil.
The response of a specific site may change according to different deterministic scenarios under consideration. Special
attention should be given to the nonlinearity of the soil inbasin structures. For more realistic prediction of the groundmotion, synthesized waveforms should be reanalyzed, considering nonlinear effects and assuming a shallow velocity
structure model of the region.

Acknowledgments
We thank Aaron Moya for the calculation of the frequency-dependent
quality factor and for his help using the Generic Algorithm method. We
also thank the scientists and technicians at the Boğaziçi University andColumbia University, who are in charge of the strong-motion stations
whose recordings are used in the present study. This work was supportedby the Special Project for Earthquake Disaster Mitigation in Urban Areas
from the Ministry of Education, Culture, Sports, Science, and Technologyof Japan. Some of the graphs were prepared by using GMT v3.1 (Wessel
and Smith, 1998).

References
Aki, K. (1967). Scaling law of seismic spectrum, J. Geophys. Res. 72,
1217–1231.
Andrews, D. J. (1986). Objective determination of source parameters andsimilarity of earthquakes of different size, in Earthquake Source Mechanics, D. Shamita, J. Boatwright, and C. H. Scholtz (Editors), Vol.
6, American Geophysical Union, Washington, D.C., 259–267.
Beavual, C. Y., P.-Y. Bard, P. Moczo, and J. Kristek, (2003). Quantificationof frequency-deprendent lengthening of seismic ground-motion duration due to local geology: applications to the Volvi area (Greece),
Bull. Seism. Soc. Am. 93, 371–385.
Borcherdt, R. D. (1970). Effects of local geology on ground motions near
San Francisco Bay, Bull. Seism. Soc. Am. 60, 29–61.
Brune, J. N. (1970). Tectonic stress and the spectra of seismic shear waves
from earthquakes, J. Geophys. Res. 75, 4997–5009.
Daubechies, I. (1990), The wavelet transform, time-frequency localizationand signal analysis, IEEE Trans. Inf. Theory 36, no. 5, 961–1005.
Fletcher, J., and J. Boatwright (1991). Source parameters of Loma Prieta
aftershocks and wave propagation characteristics along the San Francisco peninsula from a joint inversion of digital seismograms, Bull.
Seism. Soc. Am. 81, 1783–1812.
Holland, J. H. (1975). Adaptation in Nature and Artificial Systems, University of Michigan Press, Ann Arbor, Michigan.
Iwata, T., and K. Irikura, 1988. Source parameters of the 1983 Japan Sea
earthquake sequence, J. Phys. Earth. 36, 155–184.
Kayabali, K., E. Ulugergerli, and E. Tokgöz. (2001) Düzce İli Merkezinde
246 Hektarlik Alanin İmar Planina Esas Jeolojik-Jeofizik-Jeoteknik
Etüd Raporu, Ankara (in Turkish).
Kudo, K., T. Kanno, H. Okada, O. Özel, M. Erdik, T. Sasatani, S. Higashi,
M. Takahashi, and K. Yoshida (2002). Site-specific issues for strong

1456
ground motions during the Kocaeli, Turkey, earthquake of 17 August
1999, as inferred from array observations of microtremors and aftershocks, Bull. Seism. Soc. Am. 92, 448–465.
Mallat, S. G. (1989). A theory for multiresolution signal decomposition:
the wavelet representation, IEEE Trans. Pattern Anal. Machine Intell.
11, no. 7, 674–693.
Meyer, Y. (1989). Orthonormal wavelets, in Wavelets, J. M. Combes, A.
Grossman, and Ph. Tchamitchian (Editors). Springer, Berlin.
Mindevalli, O. Y., and B. J. Mitchell (1989). Crustal structure and possible
anisotropy in Turkey from seismic surface wave dispersion, Geophys.
J. Int. 98, 93–106.
Özdemir, P., B. Taşkin, C. Vatansever, A. Sezen, A. İlki, and H. Boduroğlu(2002). November 12, 1999 Düzce earthquake: a reconnaissance report for the town of Kaynasli, in 12th Proc. European Conference ofEarthquake Engineering, London, United Kingdom.
Sawada, S. (1998). Phase characteristics on the site amplification of layered
ground with irregular interface, in The Effect of Surface Geology on
Seismic Motion, K. Irikura, K. Kudo, H. Okada, and T. Sasatani (Editors), Balkema, Rotterdam, 1009–1013.

G. Birgören and K. Irikura
Yamada, M., and K. Ohkitani (1991). Orthonormal wavelet analysis of
turbulence, Fluid Dyn. Res. 8, 101–115.
Wessel, P., and W. H. F. Smith (1998). New, improved version of Generic
Mapping Tools released, EOS Trans. AGU 79, 579.
K.O.E.R.I, Bogazici University
Cengelkoy, Istanbul
34684 Turkey
birgore@boun.edu.tr(G.B.)
D.P.R.I., Kyoto University
Gokasho, Uji
Kyoto, 611-0011, Japan
irikura@egmdpri01.dpri.kyoto-u.ac.jp(K.I.)

Manuscript received 4 September 2003.