Phase Velocity Continuous Wavelet Transform Array

Summary

In this paper, we propose a method of surface waves characterization based on the deformation of the wavelet transform of the analysed signal. An estimate of the phase velocity (the group velocity) and the attenuation coefficient is carried out using a model-based approach to determine the propagation operator in the wavelet domain, which depends nonlinearly on a set of unknown parameters. These parameters explicitly define the phase velocity, the group velocity and the attenuation. Under the assumption that the difference between waveforms observed at a couple of stations is solely due to the dispersion characteristics and the intrinsic attenuation of the medium, we then seek to find the set of unknown parameters of this model. Finding the model parameters turns out to be that of an optimization problem, which is solved through the minimization of an appropriately defined cost function. We show that, unlike time–frequency methods that exploit only the square modulus of the transform, we can achieve a complete characterization of surface waves in a dispersive and attenuating medium. Using both synthetic examples and experimental data, we also show that it is in principle possible to separate different modes in both the time domain and the frequency domain.

1 Introduction

Surface wave propagation in heterogeneous media can provide a valuable source of information about the subsurface structure if properly processed and interpreted (Campillo 1996; Shapiro 1997, 2000). For example, surface waves can be used to invert for the subsurface rigidity through inversion of the shear-wave velocity. This parameter is essential for engineering purposes and seismic hazard assessment.

In a recent paper (Scherbaum 2003), assuming that the fundamental-mode Rayleigh wave is the dominant wave type in observed ambient vibration records in the vicinity of the city of Cologne, the authors were able to estimate the shear-wave velocity profile of the superficial layers by matching the extracted Rayleigh wave ellipticity and phase-velocity dispersion curves to those obtained from theoretical earth models.

Because of the non-uniqueness of earth models that can be fitted to a given dispersion curve, the inversion for the average shear velocity profile is usually treated as an optimization problem where one tries to minimize the misfit between experimental and theoretical dispersion curves computed for a given earth model that is assumed to best represent the subsurface under investigation. With this inversion approach, one has to deal with two potential sources of errors. The first source of errors is related to the measurement of correct dispersion curves, be it for the phase or group velocity. Errors at this stage would lead to an erroneous inverted shear velocity profile regardless of the goodness of the input earth model. The second type of errors is connected to the sensitivity of the optimization algorithm with respect to the parametrization of the earth model (number of layers, shear and P-wave velocities, density, layer thickness, etc.), as well as to the initial starting model employed. Both parametrization and choice of the starting model will considerably affect the inversion outcome.

Regarding the estimation of velocity dispersion, which is the main focus of the present research, several methods exist. In general, estimation of the phase velocity from single station records is possible but requires the knowledge of the source mechanism and the epicentral distance (Udías 1999). Group-velocity estimates, however, can be derived from single station records without prior information on the source mechanism if the epicentral distance is known (Herrmann 1973; Keilis-Borok 1989). Both group and phase velocities can be obtained from two closely spaced stations along the great-circle path from event and stations (Aki & Richards 2002).

Time–frequency representations of non-stationary signals such as the Gabor transform, the continuous wavelet transform (CWT) or bilinear transforms like the Wigner–Ville or smoothed Wigner–Ville transform are commonly used for the analysis of non-stationary signals (Levshin 1972; Prosser & Seale 1999; Abbate & Das 1995; Niethammer & Jacobs 2001; Pedersen 2003b).

In the analysis of dispersive surface waves, most of the methods commonly used, are essentially based on the study of the square modulus of the representation, which can be interpreted as an energy density in the time–frequency plane for linear transforms (Delprat 1992). Restricting the analysis to the square modulus of the representation provides only partial information that cannot be used to fully characterize the different propagation modes. In the context of Rayleigh waves, the ability to distinguish between the different modes is very important especially for analysis of polarization or the interpretation of H/V spectral ratios from ambient noise recordings (Fäh 2001; Scherbaum 2003).

Having pointed out some major sources of errors in the estimate of dispersion characteristics of seismic surface waves, which magnify the trade-off between inverted and actual earth structure on the one hand and the difficulty related to mode separation when only the square modulus of the time frequency is considered on the other, we propose a new approach here which overcomes some of the above mentioned shortcomings. The method is based on the CWT time–frequency representation. Let us mention at this point that the same kind of analysis may be carried over to Gabor transforms in a straightforward way. The reader may safely replace wavelet transform in this paper by Gabor transform. We, however, stick to wavelets for definiteness reasons.

The paper is organized as follows. After the introduction and a short overview of the wavelet transforms method, we present our concept of wave propagation modelling in the wavelet domain. To this end we use the wavelet propagator (Kulesh 2005) to link the wavelet transforms of signals at different stations with the assumption that the propagating signal is affected only by the dispersion and attenuation (intrinsic) characteristics of the medium. The operator explicitly incorporates the phase, and the attenuation coefficient. We limit ourselves to a first-order approximation, which is a pointwise operator. Higher-order approximations would involve integral kernels, which are numerically unfavourable.

In we develop our inversion algorithm for the characterization of surface waves. and are devoted to the application of our method to both numerical and experimental data. Using seismograms from two or more stations, we show how to take advantage of the wavelet operator to separate different wave groups in both the time and the frequency domain. We also show how to extract the dispersion curves and attenuation coefficient for each coherent wave group separately, using a parametrization of the phase and attenuation curve in conjunction with an optimization algorithm that seeks to minimize a properly defined merit function. Here again, wavelet analysis is useful to define suitable merit functions that are localized in wavelet phase space. To assess the validity of our inversion approach, our results are compared to those from independent methods. is devoted to discussions and concluding remarks.

2 Continuous Wavelet Transform

The CWT of a signal with respect to a mother wavelet g(t) is defined as follows,

where 〈·, ·〉 is the L ²–scalar product and is generated from g through dilation (a > 0) and translations. The symbol * denotes the usual complex conjugate. The wavelet g is assumed to be a function well localized in the time and frequency domains and obeys the oscillation condition

The wavelet transform can be expressed in terms of the Fourier transform of s,

From this we see that the inverse 1/a of the scale may be associated with a frequency measured in units of the central frequency of g. If the central frequency of the wavelet is assumed to be f ₀, then each scale a can be related to the physical frequency f by a=f ₀/f. Therefore, if we select a wavelet with a unit central frequency, we can directly obtain the physical frequency by taking the inverse of the scale. For the sake of clarity, we will consider a wavelet with unit central frequency and proceed with the physical frequency instead of the scale in the next sections. The signal s(t) can be recovered from its wavelet transform as follows,

where the constant C_g is defined below.

In this paper we limit ourselves to real-valued signals and admissible, progressive analysing wavelets, that is, we assume for f≤ 0, and that the constant is finite.

3 Modelling Seismic Wave Propagation Using Cwt

3.1 Estimating The Wavelet Propagator From Recorded Signals At Two Stations

Let us assume that s ₁(t) and s ₂(t) represent two signals observed at two stations, a distance d apart. We also assume a single source and a geometrical configuration where the source and the receiver stations are aligned. If dispersive and dissipative characteristics of the medium are represented by the frequency-dependent wavenumber k(f) and attenuation coefficient α(f), the relation between the Fourier transforms of the two signals reads

1

As previously mentioned, obtaining the phase velocity and attenuation from the equation above is only possible when the observed signals consist of single non overlapping wavelet arrivals. In addition, a well-behaved phase-velocity estimate requires the knowledge of the integral number of full wavelength cycles between the two stations or the consideration of more than two stations in the analysis.

In order to offer an alternative approach more suitable for the analysis of actual seismic surface waves that may incorporate multimode arrivals superimposed by body wave arrivals, we use the wavelet propagation operator to link the wavelet transforms of signals observed at two different stations. With the assumption that the attenuation function α(f) and the phase function k(f) are slowly varying with respect to the effective size of the spectrum of the wavelet transform, we can use the methodology developed by Kulesh (2005) to express the spectral relation of eq. (1) in terms of wavelet transforms of the source signal, and that of the propagated signal, as,

2

The equation above shows that remapping of the initial wavelet image into through the deformation is characterized by the wavelet phase alteration and energy dissipation characterized by the factor and a change in location indicated by the shift from the position t to

Note that the shift k′(f) is constant for a given scale (frequency) and can be identified as the group time that can be used to compute the group velocity from the ridge of the transforms. Deriving the group velocity in this way only requires the use of the modulus of the transform and as such would not provide any additional advantage over other methods such as the classical use of Gabor transform moduli and the Wigner–Ville distribution.

Note also that in the special case, with the assumption that the analysing wavelet has linear phase (with time derivative approximately equal to 2π, as it is the case for the Morlet wavelet), eq. (2) can be written in terms of the phase and group velocities as

3

The advantage of using eqs (2) and (3) is that the dispersion and attenuation characteristics are expressed explicitly and could, therefore, be extracted. In addition, the accuracy of the estimated parameters can be controlled by comparing the actual signal with the signal resulting from the inverse wavelet transform of eq. (2).

The accuracy of our formulae can be improved by adjusting the width of the wavelet. Indeed, for an arbitrary a priori fixed frequency range, in the limit of Δf/f going to zero, our formulae become exact. Clearly there is a trade-off between accuracy of these simple formulae and the time resolution of our analysis.

3.2 Estimating The Wavelet Propagator From The Cross-Correlations

We consider a set of seismic data consisting of M traces s_i (t) from M aligned stations. We also assume that the change of the waveform's shape and amplitude through the different stations results only from the dispersion and attenuation properties of the medium described by k(f) and α(f). With these assumptions, it is possible to express the cross-correlation of any chosen pair of traces in terms of the autocorrelation of a selected reference trace. For instance, taking a trace s_r (t) at a distance d_r as reference and two other traces s_i (t) at a distance d_i and s_j (t) at a distance d_j (distance with respect to the position of the source), we can write the cross-correlation in the Fourier domain as

4

where. By analogy to the derivation of eq. (2) from eq. (1), we can express the wavelet transform of T_ij in terms of the wavelet transform of T_rr as follows,

5

eqs (1)–(5) are equivalent in the sense that they can all be used to model the propagation of a source wavelet in the medium when the attenuation and dispersion characteristics are known. The advantage of using a wavelet propagator derived from the cross-correlation becomes more perceptible when dealing with the inverse problem, namely, the attempt to estimate the dispersion and attenuation characteristics from recorded signals in an actual seismic survey. However, we should note that even with the wavelet propagator (eqs 2 and 5), if only two stations are considered in the analysis, the resolution of the phase velocity will still be bound by the ambiguity related to the unknown integral of 2π-cycle skips between the stations, which needs to be determined otherwise (Aki & Richards 2002).

Note that we could have alternatively taken the correlations of the voices of the wavelet transforms instead of the wavelet transforms of the correlations. This is due to the fact that

and that correlation (g, g) is again a wavelet. Numerically, however, the first approach is less demanding.

4 Inversion For The Dispersion And Attenuation Characteristics Of The Medium

Given the recorded signals at two or more stations we now intend to develop an approach to simultaneously estimate k(f) (and hence k′(f)) and the attenuation coefficient α(f). With the assumptions made at the beginning of, the knowledge of these functions completely determine the dispersion and dissipation characteristics of the medium under investigation. To avoid the need of searching for the integral number of 2π-cycle skips for the estimate of the phase velocity, several stations will be included in the analysis. As it is not possible to guess the correct functions straight away, we proceed with a model-based approach for these functions. Once these models have been selected we need to fix the set of parameters to be varied in order for the model to match the observation. Finding an acceptable set of parameters can be thought of as an optimization problem that seeks to minimize a cost function χ² that involve eqs. (1), (2), (4) or (5) depending on the nature of the signal to be analysed and can be formulated as follows,

where P is the number of parameters used to model the frequency-dependent attenuation α(f) and M the number of parameters used to model the frequency-dependent wavenumber k(f). p and m represent the vector of parameters describing the frequency-dependent wavenumber and attenuation, respectively. In the remainder of the contribution, the sensitivity of all quantities to the parameters p and m is implied through the dependence on α (i.e. α(p)) and k (i.e. k(m)) on these parameters.

In the present contributions we parametrize α and k using polynomials or B-spline functions. For the polynomial parametrization we write

For the parametrization with B-spline functions of order four we define α and k as

where [k1, k2] is the frequency interval in which we want to estimate α and k. Constants Δ_α and Δ _k are defined as

In the following, we intend to discuss the different steps involved in our inversion algorithm as depicted in the chart of Fig. 1.

1

Flowchart showing the sequence of optimizations to be performed in order to extract the dispersive and dissipative characteristics of individual modes (that can be identified with the continuous wavelet transform) from multimode surface wave records.

4.1 Step 1

The first step will consist of seeking a good initial condition by performing an image matching using the modulus of the wavelet transforms of a pair of traces. In this case the cost function χ² is defined as

6

where L is the number of scales considered in the wavelet transform and N is the number of time samples. The optimization is carried out over the whole frequency range of the signal. This step is equivalent to an image matching problem and finding the global minimum is a feasible task in most situations. At the end we get the frequency-dependent derivative of the wavenumber k′(f) and attenuation α(f). At this stage we need to distinguish between the case where the analysed signal consists only of one coherent arrival from the case where it is made of several coherent arrivals. In the former, the derived functions are meaningful and characterize those of the analysed event. However, in the latter, these functions cannot be easily interpreted since on one hand they were obtained from an optimization that only takes the modulus of the transform into account, which means that the phase information is lost and on the other hand, the signals involved consist of many overlapping arrivals.

4.2 Step 2

If only one single phase is observed in all the traces s_i (t), it will be enough to minimize a cost function that involves some selected seismic traces in order to estimate the attenuation and phase velocity. A possible definition of this cost function may be

7

where N is the number of samples in each signal and indicates the inverse Fourier transform, where the relation between s_i and s_j is described by eq. (1).

In order to reduce the effect of uncorrelated noise on our estimates, it is preferable to use a cost function based on the cross-correlations,

8

The sets I and J indicate the index of traces used for the cross-correlations.

Minimizing the cost function based on the cross-correlations is more advantageous for two reasons. On one hand, the effect of random noise is cancelled, on the other, with geophones laid out symmetrically around the source in a 2-D seismic survey, cross-correlations of traces from seismic waves propagating in opposite directions can be combined in the optimization. For single phase arrival in all traces, the output at the end of step 2 yields the desired result.

4.3 Step 3

In the case where the observed signals consist of a mixture of different wave types and modes, a cascade of optimizations in the Fourier domain and the wavelet domain will be necessary in order to fully determine the dispersion and attenuation characteristics specific to each coherent arrival.

In the present situation, using the frequency-dependent wavenumber k(f, p₀ ) and attenuation α(f, m₀ ) obtained from the optimization with the cost functions (7) or (8), we can approximate the cross-correlation of any pair of traces as follows:

9

Despite the good phase matching that can be achieved between the actual cross-correlations and their approximations using eq. (9), we are not yet able to extract the dispersion and the attenuation curves of specific arrivals. This limitation is due to the fact that the attenuation and the wavenumber estimated using eq. (7) or (8) represent some sort of averaged attributes from the different modes present in the signal.

In order to overcome this limitation, we carry out the optimization using the wavelet transforms of the cross-correlations in a frequency and time window around the coherent arrival of interest.

Firstly, we perform the optimization on the modulus of the transforms in which case the attenuation for the specified event is derived:

10

Next, we perform an optimization involving the argument of the wavelet transforms, which will finally provide the phase- and group-velocity curves of the analysed coherent arrival:

11

In the optimization with the cost functions (10) and (11) we use k(f, p₀ ) and α(f, m₀ ) estimated using eq. (7) or (8) as starting models. This optimization can be repeated to characterize each coherent arrival separately.

Note that if the different wave groups of interest live in completely different frequency bands and time windows, an appropriate pre-filtering followed by an optimization using the cost functions (7) or (8) will be sufficient to characterize each individual event. However, when the delimitation between the different wave groups are not along lines that are parallel to either the frequency or the time axes, one will need to take into account the variation of the frequency band with time. In such a situation, one will need to determine the width of the frequency band that delimits the range of f_q in eqs (10)–(11) for each time t_l .

An alternative approach would be to determine the wavelet ridge associated to each event. The wavelet ridge, is defined (Delprat 1992) as the position (t, f) in the time–frequency domain such that

At this position, the magnitude of the wavelet transform reaches its local maximum. Once this local curve, f _ridge(t) is determined for the event of interest, the summation over the frequency in eqs (10) and (11) can be restricted over the range of f_q such that |f _ridge (t_l ) −f_q | ≤Δf (Daubechies & Maes 1996). Δf is the wavelet frequency bandwidth at the position (t, f _ridge(t)).

Since the dependence of the cost functions (6)–(10) on the parameters p and m is highly non-linear, each function may have several local minima. To obtain the global minimum that corresponds to the true parameters, a non-linear least-squares minimization method that proceeds iteratively from a reasonable set of initial parameters is required. In the present contribution, we use the Levenberg–Marquardt algorithm (Press 1992).

5 Numerical Examples

In order to assess the performance of our method in characterizing waves with mixed coherent arrivals, we simulate seismic traces from two interfering Ricker (Ricker 1953) pulses with different dispersion and attenuation characteristics as shown in Fig. 2. Given the phase velocity, from which the wavenumber is computed and the frequency-dependent attenuation of each pulse, propagation modelling was performed using eq. (1); the synthetic traces are then formed by adding the pulses obtained from the inverse Fourier transform of eq. (1). In this manner, five seismic traces were generated to simulate the observation at five successive stations. From Fig. 2(d) the manifestation of 'effective' dispersion and attenuation are readily perceptible from the simple fact of decreasing amplitude and changing wave train shape as it can be observed from s ₀(t) to s ₄(t). Here we use the term 'effective' to express the fact that the derived parameters characterize the wave train which may consist of different surface wave modes (Lai & Rix 1998). Note that since the initial step of the method involves an optimization based on the modulus of the transform, it is basically the high-energy events (i.e. surface wave arrivals) that will determine the goodness of the fit and not those arrivals with relatively low energy content such as the body waves.

2

Modelling of a seismic arrivals that consist of two interfering Ricker wavelets of different attenuation and dispersion characteristics. (a) The phase velocities, (b) the group velocities, (c) the corresponding frequency-dependent attenuation for each pulse, (d) overlapping wave trains propagated at different stations and (e) the Fourier spectra of the different traces in (d).

From the synthesized traces at the different stations, we seek to recover the phase velocity, the group velocity and the attenuation function of each pulse by successive minimization of the cost functions as outlined in the previous section. For the parameter vectors we use polynomials of order 5 at most. The reason is that the parameters associated with higher-order terms are negligible in comparison to those from the lower terms.

For the given example we selected s ₁(t) and s ₂(t) from Fig. 2 for this purpose. The output of an optimization based on minimizing the misfit between the modulus of the wavelet transform of s ₁(t) (Fig. 3b) and the modulus of the wavelet transform of s ₂(t) (Fig. 3d) using eq. (6) below provides the effective attenuation function α _e (f) and the derivative k′ _e (f) of the effective wavenumber function.

3

Application of the wavelet optimization based on the modulus on s ₁(t) and s ₂(t) to obtain the initial values. (a) Plot of s ₁(t), (b) the modulus of the wavelet transform of s ₁(t). (c) the signal s ₂(t) (dashed line) is compared to the one obtained from the inverse wavelet transform of the estimated wavelet coefficient (solid line) using derived parameters. Because the optimization is based only on the modulus, the phase of the signal cannot be predicted correctly. (d) The modulus of the wavelet transform of s ₂(t) and (e) the modulus of the wavelet transform of the estimated signal in (c).

Fig. 3(c) shows a comparison between the actual s ₂(t) and the estimated s ₂(t) from the derived parameters. We observe that while the group arrivals of the two wave trains are in agreement, their phases show a significant mismatch. This observation is consistent with the fact that the cost function considered for the optimization is only sensitive to the modulus of the transform and, therefore, does not incorporate any phase information. Comparing the moduli of the wavelet transforms of the actual and predicted s ₂(t) reveals a good match of the wavelet images (Figs 3d and e).

To correct for the observed phase misfit, we consider the minimization of a new cost function described by eq. (7) in which the effective parameters estimated previously are taken as initial values. The actual cost function (eq. 7) outperforms the previous one (eq. 6) in the sense that it takes into account the phase information and furthermore uses all the traces for the minimization process. Thus, minimization of eq. (7) provides effective parameters (wavenumber k_e (f) and attenuation α _e (f)) that are representative of the medium under investigation. Fig. 4 depicts a comparison between the actual traces at each receiver with its predicted counterpart from the estimated model parameters. Even though a very good fit can be observed for all traces, the estimated parameters are not yet indicative of the characteristics of the individual pulses hidden in composite wave trains, as can be seen from the comparison of the inverted effective wavenumber with those of the individual pulses. This curve seems to fit more or less that of the pulse with higher phase velocity but tends to disagree with the latter at high frequency. The difficulty of associating the effective parameters with those of the individual pulse is more striking if one attempts to compare the effective attenuation with those of the individual pulses as show in Fig. 4(c). Here we observe that the effective attenuation does not fit any of the individual pulses except for some very narrow frequency band. In order to circumvent the difficulty of extracting the characteristics of each event, we need to design a filter that would restrict the optimization to a specific frequency band. The design of such a cost function based on the cross-correlations of selected traces and its advantages has been discussed at the end of the preceding section.

4

Comparison between the original signals (dashed lines) with those produced at the end of the optimization indicated by step 2 on the flowchart in Fig. 1 (solid lines). (b) Comparison between the derived effective wavenumber (solid line) with those of the individual pulses (dashed line). (c) Comparison between the estimated frequency-dependent attenuation (solid line) with those of the individual pulses (dashed line). Even though a good fit of the waveforms has been achieved with this optimization, the inverted parameters do not reveal explicitly the characteristics of the individual pulses.

We note for instance that in Fig. 3(d) two coherent arrivals that correspond to the two interfering pulses can be identified. This gives an indication about the frequency band to be used for the filter. For the present example, two frequency bands are considered in order to separate the interfering pulses and determine their dispersion characteristics. The first one ranging from 5 to 40 Hz and the second one from 45 to 80 Hz. For each frequency band, the parameters resulting from the previous optimization are used as initial conditions with the cost functions (10) and (11), respectively, to first derive the attenuation and the frequency-dependent wavenumber from which the group and phase velocities are then computed. The estimated dispersion and attenuation results for each pulse in the composite traces are summarized in Figs 5(e)–(g) and 6(e)–(g). The phase velocity, the group velocity and the attenuation are plotted along side the estimated ones, showing an excellent agreement.

5

Characterization of the low-frequency pulse. (a) The modulus of the wavelet transform for the cross-correlation between s ₀(t) and s ₃(t), (b) the corresponding phase image. (c) The estimated modulus of the wavelet transform for the cross-correlation between s ₀(t) and s ₃(t) and (d) the corresponding estimated phase image. (e) Original phase velocity (dashed line) and the estimated phase velocity (solid line), (f) original group velocity (dashed line) and the estimated group velocity (solid line) and (g) original attenuation (dashed line) and the estimated attenuation (solid line). Indeed the estimated phase velocity, group velocity and the attenuation fit the initial model curves very well.

6

Characterization of the high-frequency pulse. (a) The modulus of the wavelet transform for the cross-correlation between s ₀(t) and s ₃(t) and (b) the corresponding phase image. (c) The estimated modulus of the wavelet transform for the cross-correlation between s ₀(t) and s ₃(t) and (d) the corresponding estimated phase image. (e) Original phase velocity (dashed line) and the estimated phase velocity (solid line), (f) original group velocity (solid line) and the estimated group velocity (dashed line) and (g) original attenuation (dashed line) and the estimated attenuation (solid line). As for the low-frequency pulse, here again we observe that the estimated phase velocity, group velocity and attenuation fit the initial model curves very well.

Now that we have the dispersive and dissipative characteristics of each event, we may assess the goodness of the estimation by comparing the original cross-correlations with those obtained from the estimated parameters in the time domain. To do so, we reconstruct the cross-correlation from each pulse at all stations. The results at each station are added together and compared to the original cross-correlation. The results are shown in Fig. 7. The agreement between the original and reconstructed cross-correlation is very good.

7

Comparison between the modelled cross-correlations (dashed lines) with those obtained adding the cross-correlations derived form the inverse wavelet transform (solid lines). The agreement of the time-series is very good both in terms of phase and amplitude.

6 Application To Experimental Data

Field data are used to verify the ability of the proposed method to estimate phase velocities and attenuation from seismic data.

The experimental data (Fig. 8) consist of a 2-D shallow seismic survey (stations along a line) at Kerpen, a particular site in the Lower Rhine embayment where the buried scarp of a historically active fault is presumed. Several profiles of 48 channels with 2-m interreceiver spacing were collected using hammer blows as seismic source. From the analysis of the first arrival data, a 2-D tomographic model of the investigated site was derived, however, no shear-wave velocity information is available for the shallow structure. We selected a seismogram profile (Fig. 8) with prominent low-frequency, high-amplitude arrivals that correspond to the surface wave arrivals we intend to characterize.

8

Observed seismograms obtained from a shallow seismic experiment using a sledgehammer as source. The distance between consecutive stations is 2 m. Traces in subsections A and B are delimited by the extent of the bold lines.

6.1 Characterization Of The Surface Wave With The Cwt-Based Method

We selected two subsections from the seismograms for our analysis. These subsections are labelled 'subsection A' and 'subsection B' in Fig. 8. In subsection A, the surface wave arrival consists of one coherent energy arrival while in subsection B two distinct coherent energy arrivals can be observed. The occurrence of these two different energy patches in the wavelet image are not indicative of the typical modes of the observed Rayleigh wave as in such a case a certain frequency overlap between the different modes is to be expected. A very good example from experimental data that illustrates such a behaviour is shown in (Beaty & Schmitt 2003; Forbriger 2003a,b). In the presented example the two coherent energy arrivals do not show any overlap in frequency content. Therefore, we speculate that existence of these two coherent arrivals in the frequency band between 5 and 40 Hz in 'subsection B' is most likely due to the strong attenuation characteristic of the medium under investigation for a certain frequency range (around 30 Hz), which happens to lie within the frequency band of the observed surface wave. Therefore, the attempt to separate these two coherent arrivals in the present section is not aimed at the extraction of phase-velocity curves of the fundamental and higher modes of the Rayleigh arrival but simply as an illustration of the performance of the proposed method in estimating the dispersive and dissipative characteristics of coherent energy arrivals that can be distinguished clearly in the time–frequency plane.

For the analysis of the data in each subsection, we followed the different steps depicted in the chart in Fig. 1 and illustrated with the synthetic example in. In the present and next example we will skip the comment on the intermediate steps and give only the final results. The same line of reasoning is followed as in the numerical examples.

Fig. 9(a) shows selected cross-correlations for subsection A and their respective predicted counterparts. Note the excellent match between the observed and estimated cross-correlation. This emphasizes the advantage of using the cost function based on the cross-correlations of the traces instead of the one based on the traces themselves as the contribution of events that are not coherent on all traces would not significantly contribute to the minimization. The estimated phase velocity, group velocity and attenuation for this coherent arrival are depicted in Figs 9(b) and (c).

9

Characterization of the observed surface wave as identified from the traces in subsection A. (a) Comparison between the actual cross-correlations and the estimated ones for selected pairs of traces. (b) The estimated phase velocity (solid line) and group velocity (dashed line) and (c) the estimated attenuation.

Figs 10–12 summarize the results of the analysis of those traces in the subsection 'B'. The surface wave arrival in this subsection is made of two coherent arrivals that can be identified from the modulus of the wavelet transforms of the traces at stations 23 and 27 shown in Figs 10(b) and (d). Aside from the filtering needed to separate the different events, the processing steps are similar to those of the previous example.

10

Details of the wavelet transform for signals s ₂₃(t) and s ₂₇(t) from subsection B of the observed seismograms. (a) Plot of s ₂₃(t), (b) the modulus of the wavelet transform of s ₂₃(t), (c) the signal s ₂₇(t) and (d) the modulus of the wavelet transform of s ₂₇(t). Note that two coherent arrivals that correspond to the surface wave arrivals can be observed in the frequency between 20 and 40 Hz.

11

Characterization of the low-frequency (a, b) and the high-frequency (c, d) surface wave events identified in the seismograms from subsection B. (a, c) The estimated phase velocity (solid line) and group velocity (dashed line) and (b, d) the estimated attenuation.

12

(a) Comparison between the reconstructed signals (solid line) and the observed seismic traces (dashed lines). (b) Comparison between the reconstructed cross-correlations (solid line) and those from the observed seismograms (dashed lines). In both cases a good match between the time-series can be observed. Note that the cross-correlations compare much better than the signals (see text for explanation).

The final results for the high- and low-frequency coherent arrival are summarized in Figs 11 and 12, respectively. As in the numerical example, we reconstruct the actual observed traces (computed cross-correlation) by adding the signals (cross-correlations) that can be obtained for each coherent arrival, at each station by using the estimated dispersion and attenuation characteristics. Fig. 12(a) shows the comparison for the signals and Fig. 12(b) the comparison for the cross-correlations in subsection B. In both cases, the actual and estimated time-series compare very well. Note that the estimated signals do not replicate the high-frequency component present in the actual traces. To avoid any synergetic effect of such misfits into the parameters' estimation, it is preferable to use the cost function based on the cross-correlation for the minimization process.

6.2 How Reliable Are The Dispersion Curves Obtained From The Cwt-Based Method?

An independent way to check the quality of the phase-velocity estimates from the CWT-based method would be to perform a multimode forward modelling of the Rayleigh wave. This would require accurate information about the shear-wave velocity distributions, which is unfortunately not available. However, we still have the possibility of cross-checking the CWT-based phase velocity by comparing it with the results obtained from commonly employed array techniques applied to the linear configuration of receivers. We evaluated the phase velocities using Capon's high-resolution minimum variance f−k (Capon 1969) and the MUltiple Signal Classification (MUSIC) method introduced by Schmidt (Schmidt 1986). Both of these methods exploit the information content of the cross-spectral matrix :

formed from the pairwise cross-spectra of the signals s(r _i ,) recorded by the N sensors at position vectors r_i .

An estimate of the phase-velocity v_p (ω) of a plane wave front crossing the sensor array by Capon's algorithm is obtained by maximizing the expression

12

where a (k, ω) are the steering vectors accounting for the travel time delays of the plane front with components denotes the horizontal wavenumber vector with |k| =ω/v_p (ω) and the direction of arrival is given by . For the linear receiver configuration of the hammer shot data set we perform the maximization of eq. (12) in a grid search manner along the 1-D wavenumber grid projected along the source–receiver direction.

The MUSIC algorithm (Schmidt 1981, 1986) belongs to the class of subspace-methods. Those methods rely on the orthogonality between signal and noise subspaces spanned by the eigenvectors of the cross-spectral matrix to deduce signal propagation characteristics for multiple signal contributions. The phase velocities of multiple sources (arrivals) can be derived by searching for the roots of the signal steering vectors projected onto the noise subspace or equivalently by maximizing the function

13

E _l (ω), l= 1, … , N are the normalized eigenvectors of denotes the number of acting sources and a(k, ω) are the steering vectors as above. The maximization of eq. (13) is performed again by a grid search approach and the number of sources q is estimated by an information theoretical approach (Wax & Kailath 1985) based on Akaike's criterion (Akaike 1973).

We selected the complete waveform windows of subarrays corresponding to subsections A and B along the shot profile to compute P(ω) in eqs (12) and (13) for horizontal wavenumbers k sampled equidistantly in one dimension and for a set of 200 discrete frequencies spaced equidistantly on a logarithmic frequency scale from 10 to 50 Hz. In Fig. 13 we show the results of this analysis for the two subarrays. Each subfigure pictures the results of the MUSIC approach as greyscale background image. For better visibility we normalized the maxima to one for each individual analysis frequency. The superimposed contour lines give the distribution of P(ω) from Capon's analysis method. In order to compare these results with the phase-velocity estimates from the CWT-based method, we finally superimposed the corresponding estimates as white solid curves.

13

Comparison of the slowness estimates obtained for subsection A (a) and subsection B (b) from the CWT method (solid curve), Capon's high-resolution f−k method (contour lines) and the MUSIC algorithm (greyscaled background image).

For subsection A all three methods show a consistent estimate of the phase velocities from the seismic records. For subsection B, we recognize less consistency between the different estimates of the phase velocities. It is interesting to note that the contour plot of Capon's analysis results do not reach the frequency range where the surface wave is highly attenuated. Also the greyscale background depicting the result from the MUSIC approach shows a decreased resolution of the velocity in this frequency range, indicated by the peak broadening along the vertical scale. The region where the phase-velocity estimates from all methods come quite close correspond to those where the energy of the surface arrival are significantly high (around 25 and 35 Hz). The relatively large misfit between the velocity estimates from the presented method on one hand and those from MUSIC and Capon's analysis on the other is probably due to the sensitivity of the methods to the reduced signal to noise ratio in this frequency range.

7 Discussions And Conclusions

A model-based method for velocity dispersion curves estimate is developed. The method allows simultaneous computations of both phase and group velocity and, therefore, the risk of error amplification by computing the group velocity from the phase velocity and vice versa is eliminated. The method owes its robustness to the fact that the minimization process involves not only the modulus but also the phase of the wavelet transform thus making it possible, in principle, to reconstruct the dispersed signal from the manipulated wavelet coefficients. Based on the accuracy of inverted group and phase velocities and on the reconstructed signal we show that the optimization based on a cost function using either the argument or the modulus of the wavelet transform where FWT and IWT is performed at each iteration step, provides the best results. The velocity dispersion curves estimate is tested on experimental data and compared with the results from independent methods such as MUSIC and Capon's spectral analysis. The results from the different methods show a consistent estimate of the phase velocities.

For multimode signals where wave groups are not easily separated in the CWT representation, a pre-processing step is required for the application of the proposed method as illustrated in Fig. 1. Such a step will consist on the application of the reassignment technique (Auger & Flandrin 1995; Chassande-Mottin 1998; Pedersen 2003a) to the CWT or the Gabor transform of the signal, which allows a much better reading of the group velocity associated to each mode. For each mode where the corresponding ridge can be at least partially identified in the reassigned representation, only step 1 and step 2 from the chart in Fig. 1 will be needed to obtain an estimate for the phase velocity and attenuation curve associated to that mode. The approach depicted above is the subject of a forthcoming contribution.

Acknowledgments

This project is supported by the Deutsche Forschungsgemeinschaft (DFG) within the framework of the priority program SPP 1114, 'Mathematical methods for time series analysis and digital image processing'. We thank Gabi Laske and other anonymous referee for the constructive review of the manuscript.

References

1

,

1995

.

Wavelet transform signal processing for dispersion analysis of ultrasonic signals

,

IEEE Ultrasonics Symposium

, pp.

751

–

755

.

2

,

1973

.

Information theory and an extension of the maximum likelihood principle

,

Proc. 2nd Int. Symp. Inform. Theory

, pp.

267

–

281

.

3

,

2002

.

Quantitative seismology

, 2nd Ed., University Science Books, Sausalito, California.

4

,

1995

.

Improving the readability of time-frequency time-scale representation by the reassignment method

,

IEEE Transaction on Signal Processing

,

5

,

1068

–

1089

.

5

,

2003

.

Repeatability of multimode Rayleigh-wave dispersion studies

,

Geophysics

,

68

(

3

),

782

–

790

.

6

,

1996

.

Crustal structure south of the Mexican Volcanic Belt, based on the group velocity dispersion

,

Geophys. Int.

,

152

,

361

–

370

.

7

,

1969

.

High-resolution frequency-wavenumber spectrum analysis

,

Proc. IEEE

,

57

,

1408

–

1418

.

8

,

1998

.

Méthode de réallocation dans le plan temps-fréquence pour l'analyse et le traitement des signaux non stationnaires

,

PhD thesis

, University of Cergy-Pontoise, France.

9

,

1996

.

A nonlinear squeezing of the continuous wavelet transform based on auditory nerve models

, in

Wavelets in Medecine and Biology

, eds ,

CRC Press

, Boca Raton.

10

,

1992

.

Asymptotic wavelet and Gabor analysis: Extraction of instantaneous frequencies

,

IEEE Trans. Info. Theor.

,

38

(

2

),

644

–

664

.

11

,

2001

.

A theoretical investigation of the average H/V ratios

,

Geophy. J. Int.

,

145

,

535

–

549

.

12

,

2003

.

Inversion of shallow-seismic wavefields: I. wavefield transformation

,

Geophys. J. Int

,

153

,

719

–

734

.

13

,

2003

.

Inversion of shallow-seismic wavefields: II. inferring subsurface properties from wavefield transforms

,

Geophys. J. Int

,

153

,

735

–

752

.

14

,

1973

.

Some aspects of band-pass filtering of surface waves

,

Bull. seism. Soc. Am

,

63

,

663

–

671

.

15

,

1989

.

Seismic Surface Waves in Laterally Inhomogeneous Earth

,

Kluwer Academic Publishers

, London.

16

,

2005

.

Modeling of wave dispersion using continuous wavelet transforms

,

Pure and Applied Geophysics

,

162

,

843

–

855

.

17

,

1998

.

Simultaneous Inversion of Rayleigh Phase Velocity and Attenuation for Near-Surface Site Characterization, Report No. GIT-CEE/GEO-98-2

, Georgia Institute of Technology, School of Civil Engineering, Georgia.

18

,

1972

.

On a frequency-time analysis of oscillations

,

Ann. Geophys.

,

28

,

211

–

218

.

19

,

2001

.

Time-frequency representations of Lamb waves

,

J. acost. Soc. Am.

,

109

(

5

),

1841

–

1847

.

20

,

2003

.

Improving group velocity measurements by energy reassignment

,

Geophysics

,

5

,

677

–

684

.

21

,

2003

.

Improving surface-wave group velocity measurement by energy reassignment

,

Geophysics

,

68

(

2

),

677

–

684

.

22

,

1992

.

Numerical Recipe in C: The art of Scientific Computing

,

Cambridge University Press, Cambridge

.

23

,

1999

.

Time-frequency analysis of the dipsersion of lamb-modes

,

J. acoust. Soc. Am.

,

105

(

5

),

2669

–

2676

.

24

,

1953

.

The form and laws of propagation of seismic wavelets

,

Geophysics

,

18

,

10

–

40

.

25

,

2003

.

Determination of shallow shear wave velocity profiles in the Cologne, Germany area using ambient vibration

,

Geophys. J. Int.

,

152

,

597

–

612

.

26

,

1981

.

A signal subspace approach to multiple emitter location and spectral estimation

,

PhD thesis

, Stanford University, California.

27

,

1986

.

Multiple emitter location and signal parameter estimation

,

IEEE Trans. Ant. Prop.

,

34

(

3

),

276

–

280

.

28

,

1997

.

Surface waves propagation across the Mexican Volcanic Belt and origin of the long-period seismic-wave amplification in the Valley of Mexico

,

Geophys. J. Int.

,

128

,

151

–

166

.

29

,

2000

.

Average shear-wave velocity structure of the Kamchatka Peninsula from the dispersion of surface waves

,

Earth Planets Space

,

52

,

1049

–

1056

.

30

,

1999

.

Principles of Seismology

,

Cambridge University Press, Cambridge

.

31

,

1985

.

Detection of signals by information theoretic criteria

,

IEEE Transactions on ASSP

,

33

(

2

),

387

–

392

.

© 2005 RAS

jacksonthinty1951.blogspot.com

Source: https://academic.oup.com/gji/article/163/2/463/2056781

Share this post