Dispersion analysis of MASW

Three common algorithms are applied for extracting the dispersion relationship from multi-stations seismic records. They are (1) frequency-wavenumber, f-k, transformation (Capon, 1969; Yilmaz, 1987; Alleyne et al, 1990; Forchap et al, 1998; Lin et al, 2003;Lu et al, 2004) (2) frequency- slowness, f-p, transformation (McMechan et al, 1981) (3) Phase shift method (Park et al, 1998b). These methods are reviewed as follows.

2.3.3.1 f-k transformation

Due to the frequency-dependant characteristic, geophysicists always need to map data from time domain to frequency domain. The Discrete Fourier Transform (DFT) is the common one used for any sequential series. For a discrete 2-D wavefield collected in the field, u(tm,xn) with m samples in time domain and n samples in space domain, the algorithm starts with a DFT and its DFT spectra at multiple stations is,

( )

where u is the velocity or acceleration measured by the receiver, U is the DFT of u, j=√-1, tm

= mΔt, fi = iΔf = i/[(M-1) Δt], and xn = nΔx. The i, n, and m in (2-56) are integer indices to represent respectively discrete points in the frequency, space, and time domain.

For each frequency component, the wavefield is a harmonic function of space. The wavenumber k (i.e. spatial frequency) can be determined from the wavenumber analysis (spectral analysis in space). The wavenumber analysis of the multi-station signals can be performed using the discrete-space Fourier Transform (Lin et al, 2003):

( )

where Y(fi,k) represents the wavefield in the frequency-wavenumber domain, as shown in Fig.

2-20.

Obviously the number of stations is much less than the number of samples in the time domain. Hence the aliasing problem in space domain is the first great concern when f-k transformation applied. Besides the adjustment of spatial testing parameters, focusing on the aspect of signal processing, the common way for enhancing the resolution in space domain is to add zero-value traces after the raw data. However the algorithm presented here (Lin et al, 2003) applies the discrete-space Fourier Transform in space domain.

The discrete-space Fourier Transform is different from the discrete Fourier Transform in that the wavenumber remains continuous but the fast algorithm (FFT) cannot be used. Due to deficiency of space sampling, the discrete-space Fourier Transform rather than discrete Fourier Transform is used in the space domain, such that the resolution in the wavenumber domain can be arbitrarily chosen. The wavenumber (k) of the surface wave can be identified at the peaks of the amplitude spectrum of Y(fi,k). The phase velocity is then determined by the definition v = 2πf/k. An example of f-k transformation is shown in Fig. 2-20.

Fig. 2-20 An example of f-k transformation

2.3.3.2 p-f transformation

The f-p transform is equivalent to a plane-wave decomposition of the wavefield, where p is the horizontal slowness (i.e. the inverse of velocity) and the intercept τ is a transformed time (linearly moved out). The transform is often called Slant-Stack because considering a wavefield the basic operation is that of stacking all values along each inclined (slant) line (Robinson, 1982). The sum of all the values along the line is then associated to a point in the

new domain having as coordinates the slope p0 and intercept τ0 of the line as shown on Fig.

2-21 (Foti, 2000).

Fig. 2-21 Exemplification of the Slant-Stack transform concept (Foti, 2000)

Fig. 2-22 An example of p-f transformation (Beaty et al, 2003)

For a recorded wavefield, u(xi,t), the slant-stack transform, as shown in Fig. 2-22, is given in the following expression: (Beaty, 2000; Forbriger, 2003)

) , (

) , (

∑

=

+

=

= ^N

i

i i x px t

u p

U τ τ (2-59)

where τ is an intercept time, p=1/v is slowness or the inverse of velocity, xi is the offset value

for the channel i and N is the number of the seismic channels used in the τ-p domain.

If the seismic signal u(t, xi) is written as a Fourier integral:

+∞

∫

The slant stack becomes:

∑ ∫

From (2-61) the Fourier coefficients of the slant stack:

∑

To gain full resolution it is best to normalize the signal energy to the offset dependence of plane wave in elastic media which is (Forbriger, 2003):

+∞

∫

=

0

2(t,x)dt constant

u (2-63)

Then (2-64) becomes:

∑

fi are appropriate factor to scale the seismogram to match (2-63). Substantial components of the seismic wavefield that travel with phase slowness p at an angular frequency ω will produce an amplitude maximum of G^SLS at (ω, p). The dispersion relation p(ω) of the surface wave will become apparent from these maxima. An example of p-f transformation is shown in Fig. 2-22.

2.3.3.3 Phase shift method

A phase shift method, as shown in Fig. 2-23, was presented by Park et al (1998). A wavefield in time-space domain is represented as u(x,t). Applying the Fourier Transform with

respect to time gives,

∫

= u x t e dt x

U( ,ω) ( , ) ^iωt (2-65)

U(x,ω) can be expressed as the multiplication of two separate terms:

Fig. 2-23 An example of Phase shift method (Park et al, 1998)

where P(x,ω) and A(x,ω) are phase and amplitude spectrum respectively. P(x,ω) contains all the information about dispersion properties and A(x,ω) contains the information about all other properties such as attenuation and spherical divergence. Therefore, U(x,ω) can be expressed as follows:

)

Applying the following integral transformation to U(x,ω) in (2-67) we obtain V(ω,φ):

( ) dx

The integral transformation in (2-68) can be thought of as the summing over offset of

wavefields of a frequency after applying offset-dependent phase shift determined for an assumed phase velocity cω (=ω/φ) to the wavefields in (2-67). This process is identical to applying a slant stack to the equivalent time-domain expression of U(x,ω)/∣U(x,ω)∣ for a single frequency. To insure equal weighting during analysis of wavefields from different offsets U(x,ω) is normalized with respect to offset compensating for the effects of attenuation and spherical divergence. Therefore, for a givenω, V(ω,φ) will have a maximum if

ω ω

φ = /C

=

Φ (2-69)

Because A(x,ω) is both real and positive. For a value of φ where a peak of V(ω,φ) occurs, the phase velocity cω can be determined. If higher modes get appreciable amount of energy, there will be more than one peak.

Dispersion curves result from transforming of V(ω,φ) to obtain I(ω,cω) through changing the variables such that cω=ω/φ . In the I(ω,cω) wavefields, there will be peaks along the cω-axis that satisfy (2-69) for a given ω. The locus along these peaks over different values of ω that permits the images of dispersion curves to be constructed (Park et al, 1998). An example of f-v transformation is shown in Fig. 2-23.