Blind deconvolution methodology for site-response evaluation exclusively from ground-surface seismic recordings

Blind deconvolution methodology for site-response evaluation exclusively from ground-surface seismic recordings

Soil Dynamics and Earthquake Engineering 18 (1999) 47–57 Blind deconvolution methodology for site-response evaluation exclusively from ground-surface...

367KB Sizes 0 Downloads 6 Views

Soil Dynamics and Earthquake Engineering 18 (1999) 47–57

Blind deconvolution methodology for site-response evaluation exclusively from ground-surface seismic recordings A. Zerva a,*, A.P. Petropulu b, P.-Y. Bard c a

Department of Civil and Architectural Engineering, Drexel University, Philadelphia, PA 19104, USA Department of Electrical and Computer Engineering, Drexel University, Philadelphia, PA 19104, USA c Observatoire de Grenoble, Joseph Fourier University, LGIT/IRIGM BP 53 X, 38041 Grenoble Cedex, France b

Received 8 January 1998; accepted 16 January 1998

Abstract A novel blind deconvolution methodology for identification of the local site characteristics based on two seismograms recorded on the free surface of a sediment site is presented. The approach does not require recordings at depth nor at a nearby rock outcrop, and eliminates the need for any prior parameterization of source and site characteristics. It considers that the surface recordings are the result of the convolution of the ‘input motion at depth’ with transfer functions (channels) representing the characteristics of the transmission path of the waves from the input location to each recording station. The input motion at depth is considered to be the common component in the seismograms (same input in a statistical sense). The channel characteristics are considered to be the part in the seismograms that is non-common, since the travel path of the waves from the input motion location at depth to each recording station is different, due to spatially variable site effects. By means of blind deconvolution, the algorithm eliminates what is common in the seismograms, namely the input motion at depth, and retains what is different, namely the transfer functions of the site from the input location to each recording station. It estimates the site response in both frequency and time domains, and identifies the duration of the site’s transfer functions. The methodology is applied herein to synthetic data at realistic sites for performance validation. The blindly estimated results are in almost perfect agreement with the actual site characteristics, indicating that the approach is a promising new tool for seismic site-response identification from recorded data. 䉷 1998 Elsevier Science Ltd. All rights reserved. Keywords: Seismic recordings; Site response analysis; Blind deconvolution; Time and frequency domain; Duration; Non-minimum phase channels

1. Introduction Local site effects can significantly amplify seismic motions and increase the potential for damage and collapse of engineering structures. Site amplification has been observed during many earthquakes (for example, the 1985 Michoacan, Mexico, and the 1989 Loma Prieta, California, earthquakes) and can reach values as high as 40 [1]. A lot of attention has been focused recently on establishing site response characteristics from the analysis of recorded data (for example, Refs [1–12]); Bard [2], Field and Jacob [5] and Safak [13] presented extensive reviews on the subject. Although these ‘site-specific’ analyses provide information only for the site under investigation, they improve our understanding of the effects of site amplification. Generally, the available methods for establishing the sitespecific response using recorded data can be classified as * Corresponding author.

‘reference’ and ‘non-reference’ site techniques. The most common reference site techniques are based on the evaluation of spectral ratios (spectral, cross-spectral and response spectral ratios), that require recorded data at two locations: the first recording (input motion), ideally, at the base of the sediment site and the second (output motion) at its free surface. However, because the true input motion cannot be recorded, it is approximated by downhole recordings, which may be contaminated by free surface effects, or by recordings at the nearby rock outcrop, for which it is assumed that the source effects, azimuth and incident angles of waves are the same as those at the base of the site [5]. The non-reference site techniques involve a generalized inversion scheme that estimates source, path and site effects [14], but may, in cases, require prior parameterization of the source or path characteristics. As to the single-station techniques, based on the horizontal-to-vertical spectral ratio of earthquake recordings (‘receiver function’ technique [15]) or microtremor recordings [16], it is generally accepted that

0267-7261/99/$ - see front matter 䉷 1998 Elsevier Science Ltd. All rights reserved. PII: S0267-726 1(98)00005-0

48

A. Zerva et al. / Soil Dynamics and Earthquake Engineering 18 (1999) 47–57

they do provide a reliable estimate of the fundamental frequency of soft soils, but their ability to provide information on the amplification level is still questionable [17]. Blind deconvolution techniques have also been used in site response analyses in order to identify both the seismic input and the characteristics of the travel path of the seismic waves based on a recorded seismic motion. In order to guarantee identifiability of the input signal and the travel path characteristics of the waves (channel), existing methods impose certain conditions (restrictions), which, however, may not always be realistic. A commonly imposed condition requires that the seismic input and the channel are separable in the cepstrum domain [18]. A blind deconvolution scheme, which is based on higher-order statistics, has been applied to seismic data recorded at the SMART-1 array [19]. This paper employs a recently proposed blind deconvolution methodology [20] for reconstruction of the seismic input at depth and the estimation of the characteristics of the site through transfer functions between the seismic input location and two surface stations. The required data are seismograms recorded at two different stations on the free surface of the sediment site. In this way, the approach eliminates the need for either a recording at a nearby rock site or a downhole recording, and does not require any prior parameterization of source and path effects. Furthermore, the blind deconvolution methodology identifies not only the frequency domain characteristics of the site response, but, also its true time domain characteristics, including the actual duration of the site transfer function. The approach, which is described in the following two sections, is applied herein to synthetic data for performance validation.

2. Problem set-up Consider two time series, that were obtained in the following manner. A random excitation of unknown properties was passed through two different filters, again with unknown characteristics, to yield the two time series. Only the resulting two time series are available for analyses. From this information, the blind deconvolution methodology identifies what is ‘common’ in the time series, namely the excitation, and what is ‘non-common’, namely the characteristics of the two filters by means of their transfer functions (channels in the signal processing terminology); the filter characteristics are identified in both time and frequency domains. In the site response analyses scenario, the two time series are the two recordings on the free surface of the sediment site. It is reasonable to consider that these recordings resulted from a ‘common excitation at some depth’, that was transmitted through different paths to the two recording stations on the ground surface. The application of the blind deconvolution methodology to the recorded data on the ground surface can identify ‘blindly’ the common excitation

(‘input motion at depth’) and the characteristics of the travel paths of the waves from the input motion location to each recordings station. It is noted, that the terms ‘input motion’ and ‘at depth’ are vague: if the characteristics of the travel paths of the waves from the source to each recording station were different, then the input motion for the two seismograms would be their common excitation at the source, and the identified transfer functions would represent the transmission path of the waves from the source to each recording station. However, for recordings on the free surface of a sediment site, it is more reasonable to expect that the ‘common excitation’ will be the incident motion at the bedrock–site layers’ interface, and the ‘location’ of this input motion will be that interface. In this case, the transfer functions will represent the (different) characteristics of the travel paths of the waves from the interface to each recording station, since these characteristics, except for the ideal case of horizontal layers, will be different due to spatially variable site effects. Still, the exact location of the input motion cannot be known a priori. It is recalled, that the algorithm evaluates the common and non-common characteristics in the seismograms blindly. If, for example, the separation distance between the recording stations is short, then the ‘location’ of the common excitation could be at some intermediate location within the sediments, and the input motion would then incorporate the effects of, possibly, the lower layers within the sediment site. Mathematically, the set-up of the approach is as follows. In linear site response analyses, it is generally assumed that the ground surface seismic motions are the result of the convolution of the ‘input motion at depth’ and the site transfer function. A similar modeling is used herein in the blind deconvolution approach: let x1(n) and x2(n) be two seismograms recorded at two free surface stations; the seismograms follow the model: xi …n† ˆ hi …n† ⴱ s…n† ⫹ hi …n†;

i ˆ 1; 2

…1†

i.e. it is considered that the surface recordings x1(n) and x2(n) are the convolution of the motion at depth (input), s(n), with the transfer functions (channels) h1(n) and h2(n), representing the transmission of the waves from the input location to each recording station; n in Eq. (1) indicates discrete time in samples, i.e. n ˆ t/Dt, t being time and Dt the time step; ‘ ⴱ ’ in Eq. (1) indicates convolution; and h1(n) and h2(n) indicate noise at the two locations. The transfer functions, h1(n) and h2(n) in Eq. (1), are considered to be the part in the seismograms that is non-common, since, as indicated before, the travel path of the waves from the input motion location to each recording station is different, due to spatially variable site effects. The input motion s(n), which is considered to be the common component in the seismograms (same input in a statistical sense), is described by a stationary, generally non-white, zero-mean random process: s…n† ˆ e…n† ⴱ r…n†

…2†

A. Zerva et al. / Soil Dynamics and Earthquake Engineering 18 (1999) 47–57

with e(n) indicating a zero-mean white process, and r(n) the color of the process, i.e. the filter that ‘colors’ the white input to produce the random process. Although real seismic signals are not stationary, it is common practice in earthquake engineering to view short segments of these signals as stationary in a limited or weak sense [21]. The blind deconvolution algorithm eliminates what is common in the seismograms, namely the input motion at depth, s(n), and retains what is different, namely the transfer functions of the site from the input location to each recording station. The approach considers that the transfer functions (channels) are ‘non-minimum phase’. If the channel is viewed as the superposition of impulses, then its convolution with the seismic excitation will produce a superposition of delayed and attenuated versions of the excitation. Assuming that the channel is minimum phase implies that the seismic excitation versions that will arrive faster at the surface will be stronger. This, however, is not always true: for example, P-waves, which travel faster through the media, are not stronger than S-waves, which arrive later. The analysis based on non-minimum phase channels is more general and includes the minimum phase channels as a special case. For both cases, the amplitude spectra of the sequences will be identical, but their phase spectra will be different. Waiving the minimum phase restriction on the channels allows, for the first time, the evaluation of the true time domain characteristics of the site response. Once the minimum, ii(n), and maximum, oi(n), phase components of the channels are identified, as will be described in the following section, the channels, within a time shift and a scalar constant, are expressed in terms of their minimum and maximum phase components as [22]: h1 …n† ˆ i1 …n† ⴱ o1 …n†;

h2 …n† ˆ i2 …n† ⴱ o2 …n†

…3†

Finally, by means of deconvolution of each seismogram with the corresponding estimated channel, the input motion can also be reconstructed.

3. Blind deconvolution methodology The blind deconvolution methodology [20] is described briefly in the following. The approach assumes that: 1. the channels hi(n), i ˆ 1,2, have no common zeros, and no zeros on the unit circle; 2. there are no zero-pole cancellations between the channels, hi(n), i ˆ 1,2, and convolutional components of s(n); 3. there are no common zeros between convolutional components of s(n) and each of the channels; and 4. hi(n), i ˆ 1,2, are noise processes, not necessarily Gaussian, uncorrelated to each other and to s(n). Under the above conditions, the channels h1(n) and h2(n) are identifiable within a constant and a delay (time shift). From Eqs. (1) and (2), the recorded sequences are

49

rewritten as: xi …n† ˆ e…n† ⴱ gi …n† ⫹ hi …n†;

i ˆ 1; 2

…4†

where gi(n) is the convolution of the color of the input, r(n), with the respective channel, hi(n), i ˆ 1,2. Let g˜i(n) denote ˜ i(v) be its the minimum phase equivalent [22] of gi(n), and G Fourier transform; v, throughout the approach, indicates normalized frequency. The minimum phase equivalent can be obtained from the power cepstrum of the recorded data as: g~i …n† ˆ F ⫺1 {exp‰Ci …v†Š};

ci …n† ˆ F ⫺1 {ln‰Sxi …v†Š}u…n† …5†

where Sxi …v† is the power spectrum of xi(n), Ci(v) is the Fourier transform of ci(n), u(n) equals one for n ⬎ 0 and equals zero otherwise, and F ⫺1 indicates inverse Fourier transform. Let Smin …v† W Sx1 x2 …v†G~ 1 …v†G~ ⴱ2 …v†

…6†

and Smax …v† W Sx1 x2 …v†G~ ⴱ1 …v†G~ 2 …v†

…7†

where Sx1 x2 …v† is the cross-spectrum of x1(n) and x2(n), and superscript ⴱ indicates conjugation. Based on the phase properties of g˜1(n) and g˜2(n) in relation to those of the convolutional components of x1(n) and x2(n) (Eq. 4), it can be shown [20] that: arg{Si1 i2 …v†} ˆ

1 1 arg{Smin …v†} ⫹ kv 2 2

…8†

and arg{So1 o2 …v†} ˆ

1 1 arg{Smax …v†} ⫹ kv 2 2

…9†

where k is an integer, Si1 i2 …v† is the cross-spectrum of the minimum phase parts of h1(n) and h2(n), and So1 o2 …v† is the cross-spectrum of the maximum phase parts of h1(n) and h2(n); the contribution of the minimum, ii(n), and maximum, oi(n), phase parts to the respective channels, hi(n), has been described in Eq. (3). The noises, due to the above stated assumptions, do not appear in the cross-spectrum of the observations. They do, however, contribute to the minimum phase equivalent of the observations. The effect of the noise contribution and ways of controlling it have been studied by Pozidis and Petropulu [20]. Let hmin …n† W i1 …n† ⴱ i2 …⫺n†

…10†

be the time equivalent of Si1 i2 …v† (Eq. 8), and hmax …n† W o1 …n† ⴱ o2 …⫺n†

…11†

be the time equivalent of So1 o2 …v† (Eq. 9). It can be shown [20], that the sequences hmin(n) and hmax(n) can be reconstructed from the phase of Si1 i2 …v† and So1 o2 …v†, respectively, within a scalar. When the phase and the duration (length) of any sequence are known exactly, the sequence can be

50

A. Zerva et al. / Soil Dynamics and Earthquake Engineering 18 (1999) 47–57

Fig. 1. Time histories and Fourier transform amplitudes of the data at sites 1 and 2 of Example 1. Input is a zero-mean white process and channels are nonminimum phase; time is given in time samples and frequency is normalized.

obtained as the least-squares solution of a system of equations, which involves known phase samples, i.e. NX ⫺1

x…n† sin…vn ⫺ ux …v†† ˆ sin…ux …v††

nˆ1

{v ˆ

2p k; k ˆ 0; …; L ⫺ 1} L

…12†

where x(n), in this case, plays the role of hmin(n) or hmax(n), ux(v) and N are the corresponding Fourier phase and sequence length, respectively, L( ⬎ N ⫺ 1) is the number of discrete frequencies used in the least-squares solution, and it was assumed that x(0) ˆ 1. However, from Eqs. (8) and (9), the phases of Si1 i2 …v† and So1 o2 …v†, and thus of the Fourier phases of hmin(n) and hmax(n), can be computed up to a linear phase. Also, the lengths of hmin(n) and hmax(n) are not known. Therefore, in order to identify hmin(n) and hmax(n), the reconstruction from phase methodology has to be applied inside a double loop as explained in the sequel. Let Lmin be the true length of hmin(n), which is unknown, and let l be an estimate for it. For l ˆ 2,3,… For m ˆ ⫺ l:0.5:l (i.e. m increases between ⫺ l and l with step 0.5) Rewrite Eq. (12) based on the phase:

um …v† ˆ

1 arg{Smin …v†} ⫺ mv 2

…13†

and length l. Let xm,l be the corresponding least-squares solution, and uxm;l …v† be its Fourier phase. If emin …v† W 1 2 arg{Smin …v†} ⫺ uxm;l …v† is a straight line in v, then xm,l is the right solution, i.e. xm,l(n) ˆ hmin(n), and the iteration stops here. If not, then the loop in m proceeds for the next value of m. end end To determine whether emin(v) is a straight line in v, we perform a least-squares fit of emin(v) to the equation of a line. Let Emin(l,m) be the least-squares (LS) error of the fit corresponding to length l and phase um(v), and let E(l) ˆ minm{Emin(l,m)}. For each value of l, the reconstructed sequence corresponding to the location of the minimum in Emin(l,m) is a potential solution. However, for l ⬍ Lmin, no solution exists. For l ⱖ Lmin, the LS error exhibits a profound minimum, that stays almost constant as l increases, and which is several orders of magnitude smaller than the minima corresponding to l ⬍ Lmin. We can monitor E(l) as l increases, and, when it stabilizes to a low value, i.e. E(l0) ⬇ E(l0 ⫹ 1) ⬇ E(l0 ⫹ 2)…, choose as solution the one corresponding to length l0. It has been shown [20], that the smallest length solution, indicated by an associated low least-squares error, is the correct solution. The same procedure is followed for evaluation of the length of hmax(n). The reconstructed sequence hmin(n) (Eq. 10) is decomposed into the minimum phase parts of the channels, i1(n)

A. Zerva et al. / Soil Dynamics and Earthquake Engineering 18 (1999) 47–57

51

Fig. 2. Time domain sequences hmin(n) and hmax(n) identified from the reconstruction from phase methodology for Example 1.

and i2(n). This decomposition is achieved through polynomial rooting, so that the minimum phase zeros yield i1(n) and the maximum phase zeros i2(⫺n). Similarly, hmax(n) (Eq. 11) is decomposed into the maximum phase parts of the channels, o1(n) and o2(n). Finally, from the information on their minimum and maximum phase parts, the channels are computed either in the time domain from Eq. (3), or in the frequency

domain as H1 …v† ˆ I1 …v†O1 …v†

H2 …v† ˆ I2 …v†O2 …v†

…14†

where Ii(v) and Oi(v) are the Fourier transforms of ii(n) and oi(n), i ˆ 1,2, respectively. It is noted that, because the true length of the sequences hmin(n) and hmax(n) is identified in the process, the time domain characteristics of the channels (Eq. 3) reflect their true duration.

Fig. 3. Comparison of actual and blindly estimated response of sites 1 and 2 in time and frequency domains for Example 1. — — —, blindly estimated transfer functions; - - -, actual transfer functions; time is given in time samples and frequency is normalized.

52

A. Zerva et al. / Soil Dynamics and Earthquake Engineering 18 (1999) 47–57

Fig. 4. Time histories and Fourier transform amplitudes of the data at sites 1 and 2 of Example 2. Input is a colored zero-mean white process and channels are non-minimum phase; time is given in time samples and frequency is normalized.

4. Application to synthetic seismic data The blind deconvolution methodology described in the previous section is applied herein to examples based on synthetic data for performance validation; the degree of complexity in the examples increases, and the last example analyses a realistic site response. 4.1. Example 1 The recordings at two sites, 1 and 2, are evaluated in this example by means of convolution of two different transfer functions with a white excitation. The ‘recording’ at site 1 is the result of the convolution of a realization of the white process with the following transfer function: h1 …n† ˆ ‰⫺0:306; 0:2446; 1:0; ⫺0:237; ⫺0:231Š

…15†

The ‘recording’ at site 2 is obtained from the convolution of the same realization of the white process used for the first site with the transfer function: h2 …n† ˆ ‰⫺0:357; ⫺0:250; 1:0; 0:327; ⫺0:360Š

…16†

Physically, the white excitation may be viewed as the incident motion from the bedrock at the bedrock–surface layer interface of two sites with different layer characteristics. The two transfer functions are selected so as to incorporate both minimum and maximum phase components: it can be seen from Eqs. (15) and (16), that the strongest part of the signal will not arrive first at the surface of both sites. The time series (over 200 samples) ‘recorded’ at the two sites,

and the amplitudes of their corresponding Fourier transforms are presented in Fig. 1; Fig. 1 suggests that the white excitation dominates the shape of the ‘recordings’. From the information provided in Fig. 1 (time series and their Fourier transforms), the blind deconvolution algorithm is requested to identify the non-common characteristics in the seismograms, namely the functions of Eqs. (15) and (16). Fig. 2 presents the identified hmin(n) and hmax(n) (Eqs. (10) and (11)) from the reconstruction from phase methodology; as expected, the sequences hmin(n) and hmax(n) have finite (non-zero) length, indicating that both minimum and maximum phase components are present in the channels. These minimum and maximum phase components of the channels are obtained from polynomial rooting of hmin(n) and hmax(n) (Fig. 2), and the channels are reconstructed in the time (Eq. 3) and frequency (Eq. 14) domains. The blindly predicted transfer functions are shown in Fig. 3 in the time and frequency domains, together with the actual transfer functions (Eqs. (15) and (16)), and their Fourier transforms. The solid lines indicate the blindly estimated results and the dashed ones the actual transfer functions. All quantities in the figures have been normalized with respect to their maximum amplitude, as the blind deconvolution algorithm identifies the functions within a scalar; for this reason, the time domain blindly estimated transfer functions may appear as the negative image of the actual ones (Fig. 3 and, for the next example, Fig. 5). Furthermore, since the algorithm identifies the time domain characteristics within a time shift, the ‘starting points’ of the estimated results are shifted compared to those of the actual ones.

A. Zerva et al. / Soil Dynamics and Earthquake Engineering 18 (1999) 47–57

53

Fig. 5. Comparison of actual and blindly estimated response of sites 1 and 2 in time and frequency domains for Example 2. — — —, blindly estimated transfer functions; - - -, actual transfer functions; time is given in time samples and frequency is normalized.

Fig. 3 suggests that the blindly estimated results are in almost perfect agreement with the actual ones in terms of frequency content, shape—in both time and frequency domains—and actual duration in the time domain. It is noted again that the only information given to the algorithm were the time series in Fig. 1. 4.2. Example 2 This example is similar to the previous one, except that in this case the input is colored, i.e. the white process is convolved with a transfer function (color), common at the two sites, before it is convolved with the transfer functions of Eqs. (15) and (16). Physically, this example may represent the case of two sites, both consisting of two layers overlaying the same bedrock. The bottom layer is the same at both sites, and its transfer function is given, in the example, by the well-known impulse response function of a single-degree-of-freedom oscillator (see, for example, Ref. [23]): q …17† r…n† ˆ exp…⫺zvnDt† sin{v …1 ⫺ z2 †nDt} with v and z indicating the frequency and damping coefficient, respectively, of the oscillator; the numerical values selected for the parameters were v ˆ 5p/4 (rad s ⫺1), z ˆ 0.65, with a time step of Dt ˆ 0.5 s, which yield: r…n† ˆ ‰0:0; 1:0; 0:04; ⫺0:08; ⫺0:01Š

…18†

It is noted that this impulse response function is minimum phase, i.e. it transfers the strongest part of the signal first.

The characteristics of the top layers of the sites, however, are different, and given by Eqs. (15) and (16); i.e. the ‘recordings’ at the two sites are evaluated from: xi …n† ˆ hi …n† ⴱ s…n† ˆ hi …n† ⴱ {r…n† ⴱ e…n†} i ˆ 1; 2

…19†

where hi(n), i ˆ 1,2, are different for the two sites, but r(n) is the same (Eq. 18), and e(n) is the same realization of the zero-mean white process. The outcome of the convolution (Eq. 19) in the time and frequency domains is presented in Fig. 4. Eq. (17) has been used extensively in earthquake engineering (e.g. the Kanai–Tajimi spectrum [24,25]) to model the transfer function of a single layer overlaying a rigid bedrock; furthermore, passing a white excitation through two filters (transfer functions) has been common practice in modeling seismic excitations (e.g. the Clough– Penzien spectrum [23]). It is noted that, although the process of convolutions in Eq. (19) does not conform with the actual transmission of waves through a layered medium, it may be viewed as a reasonable approximation for the first application of the algorithm in site response analysis. From the information presented in Fig. 4, the algorithm is requested to blindly estimate the non-common characteristics of the time histories. The outcome of the process is presented in Fig. 5. Again, the solid lines represent the blindly estimated (non-common) transfer functions and the dashed ones the actual ones (Eqs. (15) and (16)). The agreement is again almost perfect in terms of frequency content, shape in both time and frequency domains, and actual duration in the time domain. hmin(n) and hmax(n) (Eqs. (10) and (11)) for this example, identified from the

54

A. Zerva et al. / Soil Dynamics and Earthquake Engineering 18 (1999) 47–57

Table 1 Properties of site 1 (Ashigara Valley, station KS2) Layer

Depth (m)

Density (g cm ⫺3)

P-wave velocity (m s ⫺1)

S-wave velocity (m s ⫺1)

1 2 3 4 5 6

0.0–4.0 4.0–11.0 11.0–78.0 78.0–95.0 95.0–280.0 280.0

1.40 1.50 1.70 1.90 2.00 2.50

1500.0 1500.0 1800.0 2000.0 2200.0 3000.0

70.0 150.0 400.0 700.0 800.0 1500.0

reconstruction from phase methodology, are not presented, as they are identical to those identified in the first example (Fig. 2). 4.3. Example 3 In this example, the blind deconvolution methodology is applied to synthetic data that bear close resemblance to actual recorded time histories. Two sites with similar properties (Tables 1 and 2) are used in the evaluation of the synthetic seismograms. The sites consist of horizontal layers overlaying the bedrock, which has the same properties in both cases. Site 1 (Table 1) consists of five horizontal layers overlaying the bedrock; the characteristics of this site (layers and bedrock) are essentially those reported for Ashigara Valley, station KS2 [26]. Site 2 consists of three layers overlaying the bedrock; the total depth of the layers is the same as that of site 1, and the layer properties of site 2 are rough average values of those of site 1. The incident motion at the bedrock–layers interface is a segment (30 sample points, from 70 to 100) of the time history recorded during the strong motion shear wave window at station I03 of the SMART-1 array (Fig. 6); in Fig. 6, the amplitudes of the motions are given as dimensionless numbers, and 2048 units in the figure correspond to 2g. The synthetic (normalized) time and frequency domain results at the free surface of the two sites are presented in Fig. 7, and are considered to be the ‘two recorded seismograms on the free surface of a sediment site’. From the information in Fig. 7, the algorithm is asked to blindly estimate the characteristics of the sites in both time and frequency domains. Fig. 8 presents the reconstructed sequences hmin(n) and hmax(n), that were identified from the smallest length solution associated with low leastsquares error in the reconstruction from phase methodology. It is noted, from the comparison of hmin(n) and hmax(n) in the

figure, that the maximum phase component of the channels is significantly shorter and ‘simpler’ than the minimum phase one. hmin(n) and hmax(n) of Fig. 8 are used for the identification of the minimum and maximum phase components of the channels (Eqs. (10) and (11)), from which the sites’ transfer functions in the time (Eq. 3) and frequency domains (Eq. 14) are reconstructed. Fig. 9 presents the comparison of the actual impulse response function of sites 1 and 2 and their Fourier transforms with the channel characteristics reconstructed by means of blind deconvolution; again the solid lines represent the blindly estimated results, whereas the dashed ones the actual transfer functions of the sites. It is noted that the blindly estimated results are compared with the impulse response of each site, because, from the set-up of this example, it follows that the ‘input motion’ (common component) is the incident motion at the bedrock–layers interface, and its location is the interface. The agreement of the actual and reconstructed site characteristics in Fig. 9 is almost perfect in the time and frequency domains. In particular, it is noted that the duration of the actual impulse response functions is fully matched by that of the reconstructed ones. The frequency content of the site response characteristics is also fully captured by that of the blindly identified ones. Only small differences are observed: the blindly estimated results are slightly richer in frequency content than the actual site characteristics in the higher frequency range. It is also noted that the time step used in the evaluation of the ‘recordings’ in this example was Dt ˆ 0.02 s, and hence, the actual cut-off frequency in Fig. 9 corresponds to 25 Hz.

5. Conclusions A novel blind deconvolution approach for the identification of site response characteristics based on two

Table 2 Properties of site 2 Layer

Depth (m)

Density (g cm ⫺3)

P-wave velocity (m s ⫺1)

S-wave velocity (m s ⫺1)

1 2 3 4

0.0–20.0 20.0–120.0 120.0–280.0 280.0

1.40 1.70 2.00 2.50

1500.0 1800.0 2200.0 3000.0

150.0 400.0 800.0 1500.0

A. Zerva et al. / Soil Dynamics and Earthquake Engineering 18 (1999) 47–57

55

Fig. 6. Recorded time history at station I03 of the SMART-1 array in the N–S direction of Event 5. Time is given in time samples; the figure presents the (windowed) strong motion shear-wave window of actual duration 5.12 s.

seismograms recorded on the free surface of the sediment site is presented. The approach considers that the seismic surface motions result from the convolution of a common input excitation at depth with transfer functions, representing the characteristics of the travel path of the waves from the input motion location to each recording station. It is based on the reconstruction of signals from phase information only, and allows the site response to be non-minimum

phase. Unique features of the approach include the following: • It does not require seismic records at depth nor at a nearby rock site, as is commonly the case in site response analyses, and eliminates the need of any prior parameterization of source or site characteristics. • It provides a two-dimensional picture of the site, since it

Fig. 7. Time histories and Fourier transform amplitudes of the data at sites 1 and 2 of Example 3. Input is a segment of an actual recorded time history, and the transfer functions represent the characteristics of actual sites (Ashigara Valley); time is given in time samples and frequency is normalized.

56

A. Zerva et al. / Soil Dynamics and Earthquake Engineering 18 (1999) 47–57

Fig. 8. Time domain sequences hmin(n) and hmax(n) identified from the reconstruction from phase methodology for Example 3.

identifies channel characteristics from the ‘input motion location’ to each recording station; current approaches estimate a one-dimensional site response. • It fully identifies the true time domain characteristics of the site transfer functions, including their actual shape and duration.

The blind deconvolution methodology was applied herein to synthetic data that used simple site transfer functions as well as transfer functions of realistic sites, and incident motions at the bedrock–layer interface described by white processes, colored white processes and recorded seismic time histories. An essentially perfect match of the actual

Fig. 9. Comparison of actual and blindly estimated response of sites 1 and 2 in time and frequency domains for Example 3. — — —, blindly estimated transfer functions; - - -, actual transfer functions; time is given in time samples and frequency is normalized.

A. Zerva et al. / Soil Dynamics and Earthquake Engineering 18 (1999) 47–57

and the blindly estimated site response characteristics in both time and frequency domains resulted from the application of the method in all examples. As with any new approach, the proposed algorithm requires further testing before it is fully utilized in actual recorded seismograms on the free surface of sediment sites; currently, it is being applied to additional synthetic data. The present results strongly indicate that it will yield a promising new tool in seismic site response identification from recorded data.

[11]

[12]

[13] [14]

Acknowledgements APP would like to acknowledge the support of NSF under grant MIP-9553227.

[15]

[16]

References [1] Archuleta RJ, Seale SH, Sangas PV, Baker LM, Swain ST. Garner Valley downhole array of accelerometers: instrumentation and preliminary data analysis. Bulletin of the Seismological Society of America 1992;82:1592–1621. [2] Bard P-Y. Effects of surface geology on ground motion: recent results and remaining issues. In Proceedings of the Tenth European Conference on Earthquake Engineering, Vienna, Austria, 1994. [3] Blakeslee S, Malin P. High-frequency site effects at two Parkfield downhole and surface stations. Bulletin of the Seismological Society of America 1991;81:332–345. [4] Cramer CH. Weak-motion observations and modeling for the Turkey Flat, U.S., site effects test area near Parkfield, California. Bulletin of the Seismological Society of America 1995;85:440–451. [5] Field EH, Jacob KH. A comparison and test of various site response estimation techniques, including three that are not reference-site dependent. Bulletin of the Seismological Society of America 1995;85:1127–1143. [6] Field EH, Jacob KH, Hough SE. Earthquake site response estimation: a weak motion case study. Bulletin of the Seismological Society of America 1992;82:2283–2307. [7] Gagnepain-Beyneix J, Lepine JC, Nercessian A, Hirn A. Experimental study of site effects in the Fort-de-France area (Martinique Island). Bulletin of the Seismological Society of America 1995;85:478–495. [8] Kato K, Aki K, Takemura M. Site amplification from coda waves: validation and application to S-wave site response. Bulletin of the Seismological Society of America 1995;85:467–477. [9] Lermo J, Chavez-Garcia FJ. Site effect evaluation using spectral ratios with only one station. Bulletin of the Seismological Society of America 1993;83:1574–1594. [10] Liu H-P, Warrick RE, Westerlund RE, Sembera ED, Wennerberg L.

[17]

[18] [19]

[20]

[21]

[22] [23] [24]

[25]

[26]

57

Observation of local site effects at a downhole-and-surface station in the Marina district of San Francisco. Bulletin of the Seismological Society of America 1992;82:1563–1591. Seale SH, Archuleta RJ. Site amplification and attenuation of strong ground motion. Bulletin of the Seismological Society of America 1989;79:1673–1696. Su F, Aki K. Site amplification factors in central and southern California determined from coda waves. Bulletin of the Seismological Society of America 1995;85:452–466. Safak E. Models and methods to characterize site amplification from a pair of records. Earthquake Spectra 1997;13:97–129. Boatwright J, Fletcher JB, Fumal TE. A general inversion scheme for source, site and propagation characteristics using multiply recorded sets of moderate-seized earthquakes. Bulletin of the Seismological Society of America 1991;81:1754–1782. Langston CA. Structure under Mount Renier, Washington, inferred from teleseismic body waves. Journal of Geophysical Research 1979;84:4749–4762. Nakamura Y. A method for dynamic characteristics estimation of subsurface using microtremor on the ground surface. Quarterly of Railway Technology Research Institute, Japan 1998; 30(1). Lachet C, Bard P-Y. Numerical and theoretical investigations on the possibilities and limitations of the Nakamura’s technique. Journal of Physics of the Earth 1994;42:377–397. Arya VK, Aggarwal JK, editors. Deconvolution of seismic data. Hutchinson Ross Publishing Company, 1982. Zerva A, Petropulu AP, Abeyratne U. Blind deconvolution of seismic signals in site response analysis. In: Proceedings of the Fifth International Conference on Seismic Zonation, Nice, France, 1995. Pozidis H, Petropulu AP. Cross-spectrum based blind channel identification. IEEE Transactions on Signal Processing 1997;45:2977– 2993. Zerva A. Spatial variability of seismic motions recorded over extended ground surface areas. In: Kausel E, Manolis GD, editors. Wave problems in earthquake engineering. Southampton: CMP, in press. Oppenheim AV, Schafer RW. Discrete-time signal processing. Englewood Cliffs, NJ: Prentice Hall, 1989. Clough RW, Penzien J. Dynamics of structures. New York: McGrawHill, Inc., 1975. Kanai K. Semi-empirical formula for the seismic characteristics of the ground. Bulletin of the Earthquake Research Institute, University of Tokyo, Japan 1957;35:309–325. Tajimi H. A statistical method of determining the maximum response of a building structure during an earthquake. In: Proceedings of the Second World Conference on Earthquake Engineering, Tokyo and Kyoto, Japan, 1960. Sawada Y. (The Subcommittee for the Geotechnical Survey of the Ashigara Valley Blind Prediction Test). Geotechnical data. In: Proceedings of the International Symposium on the Effects of Surface Geology on Seismic Motion, Odawara, Japan, 1992.