Abstract: This study addresses wavelet transform for multichannel surface wave method to overcome the limitations of conventional seismic signal analysis. Seismic surface wave method is familiar nondestructive seismic test to reveal the soil profile in geotechnical engineering. The spectral analysis surface wave method using two receivers is constraint due to interferences of other signals and various types of noises. To reduce these problems, the multichannel analysis of surface wave method is developed as more reliable surface wave technique with the performance of mode discrimination. But this method is tedious and slower for phase velocity extraction using Fourier transform. However, because of characteristic of fourier transform used in the conventional phase extraction method, dispersion curve is sensitive to background noise and body waves in the low frequency range. Furthermore, under some field conditions such as pavement site, the phase extraction method can lead to erroneous dispersion curve. To overcome these problems, a new method of determining the group velocities and phase velocities using time frequency analysis based on Complex Morlet wavelet transform is highlighted with the algorithm in this study. To estimate the applicability of the proposed method, the dispersion curve at soil layer profile is revealed and the dispersion curve by proposed method is more reliable than those by the usual phase extraction method. The significance of this study is to obtain a robust and consistent performance with noise reduction for surface wave analysis regarding multichannel and Morlet wavelet transform.
INTRODUCTION
Seismic surface wave profiling is gaining popularity in engineering practice for determining shear-wave velocity profile. The main advantage of surface wave testing is essentially related to its non-destructive and non-invasive nature that allows the characterization of hard-to-sample soils without the need for boreholes that makes the subsurface seismic methods (such as down-hole and cross-hole methods) expensive and time consuming. One of the familiar seismic method for determining the shear wave velocity is the Spectral Analysis of Surface Waves (SASW) method, based on dispersive characteristics of surface waves. The SASW method is widely established as a subsurface investigative tool and is implemented in a wide variety of geotechnical environments, including pavements, solid waste landfills and sea beds (Roesset, 1998).
A unique method of collecting and calculating data to produce phase plots was utilized for the SASW testing based on two receivers. Phase plots were calculated and shown in the field by the dynamic signal analyzer, using information from the receivers as well as the source. Using the source function verifies that the phase shift between receivers is correlated to the source signal. The testing procedure of SASW method is revealed in Fig. 1.
The wrapped phase angle and unwrapped phase angle corresponding to SASW method are plotted in Fig. 2.
| Fig. 1: | General SASW field testing source-receiver geometry |
| Fig. 2: | (a) Wrapped phase angle in deg, (b) Unwrapped phase angle in degrees |
The surface wave phase velocity is calculated from phase angle as:
| (1) |
where, D is the distance between two receivers and φ12 is the phase angle between two receivers. These calculations yield an experimental dispersion curve from SASW method. The evaluation method of the dispersive phase and group velocities for SASW technique was also proposed using the harmonic wavelet transformation (Park and Kim, 2001; Kim and Park, 2002).
During the field testing various kinds of noises are in existence. The conventional phase unwrapping method with filter criteria usually leads to the correct dispersion curve. However, under some irregular profiles where the more modes dominate the surface motion, the conventional phase unwrapping method can lead to incorrect dispersion curve. For the phase unwrapping method, it is difficult to evaluate correct number of 3600 cycles in the phase spectrum with noises. Furthermore, due to the short distance coverage of SASW method, the spatial resolution is limited.
Park et al. (1999a, b, 2000) has shown that Multichannel Analysis of Surface Waves (MASW) method represents an improvement over SASW, overcoming the few but significant weaknesses of the SASW method. MASW is also a fast method of evaluating near-surface vs profile because the entire range of investigation depth is covered by one or a few generation of ground roll without changing receiver configuration. The main advantage of the MASW method is its ability to take a full account of the complicated nature of seismic waves that always contain higher modes (Strobbia and Foti, 2006; Penumadu and Park, 2005) of surface waves, body waves, scattered waves, traffic waves. Incorporating multi-station receivers and 2-D wavefield transformation (Lin and Chang, 2004; Park et al., 1999a) improves inherent difficulties in evaluating signal from noise with only a pair of receivers. The MASW technique has proven to mute the interfering seismic waves in the shot records and filter noisy surface-wave modes and thus significantly improve the range and resolution of multimodal dispersion curves in the phase-velocity-frequency domain (Ivanov et al., 2005; Penumadu and Park, 2005). Observation and detection of higher mode surface wave is the pioneering success of the MASW method (Park et al., 2000, 2001a, b). MASW method is also convenient to apply in a shallow marine environment (Kaufmann et al., 2005), to evaluate stiffness of water-bottom sediments (Park et al., 2005), to determine a sinkhole impact area as well as collapse and subsidence feature (Xia et al., 2004), to represent near surface anomalies (Park et al., 1999b) and to map bedrock surface (Miller et al., 2000, 2001a, b).
The MASW method is arranged as data acquisition, data processing for dispersion, finally inversion analysis shown in Fig. 3. This method simplifies the testing procedure and automates the data analysis. The determination of the dispersion spectrum S (w, v) from multichannel surface wave data is based on the Eq. 2 in which A (x, ω) is the normalized energy spectrum for each receiver, v is the assumed phase velocity and V is the phase velocity for a given frequency.
| (2) |
| Fig. 3: | Scheme of the three steps of the MASW method |
S(ω,v) is maximized when the two velocities are equal.
However, the usual MASW method is constraint to adjust irregularities and presence of noises for phase velocity extraction using Fourier transform. However, because of characteristic of Fourier transform used in the conventional phase extraction method, dispersion curve is sensitive to background noise and body waves in the low frequency range. Furthermore, under some field conditions such as pavement site, the phase extraction method can lead to erroneous dispersion curve. To overcome these problems, a new method of determining the group velocities and phase velocities using time frequency analysis based on Complex Morlet wavelet transform is highlighted with the algorithm.
In this study, the new evaluation method of phase and group velocities is proposed based on complex Morlet wavelet transform. Wavelet decompositions and windowed Fourier decompositions have emerged in the last few years as tools of great importance for the analysis of signals in which local frequencies can be extracted (Torresani, 1993). Wavelet theory can be viewed as a modern improvement and extension of the Fourier theory and hence allow a flexible alternative to the Fourier method in non-stationary signal analysis (Cho and Chon, 2006). Wavelet transform has recently been through considerable development adding new instruments and special functions and considered more robust with noise deduction performance. The aim of this study is to develop the consistent surface wave analysis technique with noise lessening performance using complex Morlet wavelet transform.
MATERIALS AND METHODS
Location of study: The study of the proposed method is conducted at University Kebangsaan Malaysia since January, 2008 to March, 2009. The seismic signal for soil profile collected by MASW testing at University Kebangsaan Malaysia (UKM) in Bangi, Selangor, Malaysia is used in this study. The schematics arrangement of the test is shown in Fig. 4. The propagation of the waves were detected using multiple receiving geophones transferred digitally a notebook computer. The recorded time domain signals were saved in notebook computer and were then converted to frequency domain where noise contaminations were analyzed by using MATLAB. The data acquisition and analysis are performed since April, 2008 to January, 2009 at University Kebangsaan Malaysia (UKM) in Bangi, Selangor, Malaysia.
Proposed method using wavelet transform: Wavelet Transformation (WT) is equivalent to the convolution of the wavelet function and the signal under investigation. The WT is performed by projecting a signal s(t) onto a family of zero mean functions deduced from an elementary function ψ by translations and dilations (Yves et al., 1992).
| (3) |
The variables a and b control the scale and position of the wavelet, respectively. For Seismic test, we are required to work with complex signals. The Fourier Transformation (FT) coefficients of such signals are no longer symmetric.
| Fig. 4: | Data acquisition using MASW method |
Likewise it is convenient to work with analytic wavelets which can separate amplitude and phase components and allow the measurement of the time evolution of frequency transients. The WT for processing of surface waves can be implemented in such a way that only the coefficients resulting from the forward flow components are obtained when the scale is positive and only the coefficients resulting from the reverse flow components are obtained when the scale is negative. Ignoring the translational parameter b, the scale dependent complex Morlet wavelet is given by equation:
| (4) |
where, ω0 is the nondimensional frequency and usually assumed to 5 to 6 satisfy the admissibility condition. The FT is given by equation:
| (5) |
where, H stands for the Heaviside step function. From Eq. 5, one can observe that a frequency spectrum of an upper analytic signal is obtained for a>0 and a frequency spectrum of a lower analytic signal is obtained for a<0.
The WT of a signal s(t) with the Morlet wavelet is given by Eq. 6 as:
| (6) |
If the number of scales is J, a complete set of directional wavelet coefficients can be mapped over the scales from a = -J to b = J, excluding a = 0. In order to see how negative scales have been utilized to obtain directional wavelet coefficients, let us evaluate Eq. 6 for one positive and one negative scale. Assuming b = 0 for simplicity, the WT of the signal for a = 1 and b = -1 are given, respectively, by equation:
| (7) |
| Fig. 5: | Schematic diagram for the wavelet analysis in MASW method |
| (8) |
where, the constant C = π-1/4|a|G1/2. Taking the FT of Eq. 7 and 8 yield, respectively
| (9) |
| (10) |
where, S (ω) is the FT of the signal and C is a constant. The inverse FT of Eq. 9 and 10 yields complex wavelet coefficients formed only by the forward and the reverse flow components respectively. The wavelet coefficients are used to determine the phase delay and group delay for MASW method. The crucial attention of the proposed method is dispersion curve analysis for the representation of phase velocity versus frequency curve. The design of the new method is comparable to the phase velocity extraction using Fourier transform for MASW method (Park et al., 2000, 2001a, b). The implementation of wavelet transform to determine phase velocity and group velocity for multi-channels surface wave method is revealed in the Flow-chart of wavelet analysis in MASW method shown in Fig. 5.
Determination of phase and group velocity: The wavelet transform is used to decompose a signal into wavelets, small oscillations that are highly localized in time. Whereas the Fourier transform decomposes a signal into infinite length of sines and cosines and effectively losing all time-localization information, to which the WT's basis functions are scaled and hence shifted versions of the time-localized mother wavelet. The WT is used to construct a time-frequency representation of a signal that offers very good time and frequency localization.
The WT is an excellent tool for mapping the changing properties of non-stationary signals. The WT is also an ideal tool for determining whether or not a signal is stationary in a global sense. When a signal is judged non-stationary, the WT can be used to identify stationary sections of the data stream.
The definition of the continuous WT is given by Eq. 11 and 12:
| (11) |
| (12) |
| (13) |
| (14) |
The WT is a convolution of the data sequence with a scaled and translated version of the mother wavelet, the psi function. This convolution can be accomplished directly, as in the first equation, or via the Fast Fourier Transformation (FFT)-based fast convolution in the second equation. Note that the WT is a continuous function except for the discrete data series x and its discrete Fourier transform. In these equations the *symbolizes complex conjugation, N is the data series length, s is the wavelet scale, dt is the sampling interval, n is the localized time index and omega is the angular frequency. Each of the equations contains a normalization so that the wavelet function contains unit energy at every scale.
In the WT, for each value of the scale used, the correlation between the scaled wavelet and successive segments of the data stream is computed. Unless reconstruction is needed, there are no restrictions in the WT as to the number of scales are to be used, nor of the spacing between the scales. A WT spectrum can use linear or logarithmic scales of any density desired. If needed, a high resolution spectrum can be generated for a narrow range of frequencies. The convolutions can be done up to N times at each scale and must be done all N times if the FFT is used. The WT consists of N spectral values for each scale used, each of these requiring an inverse FFT. The computational load of the WT and its memory requirements are thus considerable. The benefit from this high measure of redundancy in the WT is an accurate time-frequency spectrum. The most commonly used WT wavelet is the Gaussian wavelet is defined as:
| (15) |
by taking the pth derivative of f. The integer p is the parameter of this formula, Cp is the such that:
| (16) |
where f (p) is the pth derivative of f.
The Morlet wavelet depicted as Gaussian-windowed complex sinusoid is also commonly used WT wavelet that is defined as following in the time and frequency domains:
| (17) |
| (18) |
In the Eq. 17 and 18, ψ, eta is a non-dimensional time parameter, m is the wavenumber and H is the Heaviside function. In the time domain plot that follows, the complex Morlet wavelet is shown in Fig.6 with an adjustable parameter m (wavenumber) of 8. In addition, the complex Gaussian wavelet is defined as:
| (19) |
| Fig. 6: | Representation of the complex morlet wavelet |
So, the complex Morlet wavelet is defined by Eq. 20.
| (20) |
where, fb is a bandwidth parameter.
fc is the wavelet centre frequency.
The complex Gaussian wavelet and complex Morlet wavelet is expedient for the decomposition of the complex wave like as the seismic wave.
We can use Euler’s formula to define the analytical signal as,
| (21) |
where, i is the imaginary unit to give a more concise formula.
Now,
| (22) |
The Fourier coefficients are then given by equation:
| (23) |
The Fourier coefficients an, bn, cn are related via
| (24) |
And
| (25) |
However, after decomposition of the seismic signal basis on the complex Morlet wavelet, the wavelet coefficient can be represented as the decomposed signal by Eq. 26:
| (26) |
| (27) |
Here
| (28) |
where, H represents the Hilbert transform; the magnitude of ψ (t) is
| (29) |
and the phase of ψ (t) is:
| (30) |
The Eq. 29 and 30 show the magnitude and phase extraction through the complex Morlet wavelet transform. The phase spectrum using complex Morlet wavelet transform is revealed in the Fig. 7. This phase extraction using complex Morlet wavelet is the significant step to determine phase velocity.
The determination of group velocity and phase velocity is correlated. The group velocity is obtained from the group delay which is the guiding principle to obtain phase velocity. Figure 8 shows the group delay corresponding to the time domain signal between two receivers.
| Fig. 7: | Determination of phase spectrum through complex morlet wavelet |
| Fig. 8: | Determination of group delay tg1 and tg2 at receiver No. 1 and 2 |
| Fig. 9: | Determination of Phase delay tph1 and tph2 at receiver No. 1 and 2 |
| Fig. 10: | Multi-channels surface wave signal |
The phase delay determination correlated with group delay is the key step to determine phase velocity. The phase delay is also shown by the Fig. 9 which is used to extract phase velocity. Extraction of phase velocity through MASW method is the key concentration of this study. The concept of MASW method is developed on using multiple receivers for acquisition (Park et al., 1999a, b). The acquired multiple receiver corresponding signal is processed to determine phase velocity as well as group velocity with robust performance. The concept of the extraction of phase velocity and group velocity is implemented corresponding to multiple-channels surface wave signal shown in Fig. 10. Every adjacent pair is regarded to attain phase velocity corresponding to multiple receivers. The phase velocity versus frequency curve known as dispersion curve is an attempt to deliver the information of soil characteristics through inversion analysis. The multiple dispersion curve-phase velocity versus frequency curve is compared to reduce the presence of noise and erroneous outcome.
The procedure to determine the group and phase velocities with frequency is as follows:
| (1) | Compute Complex Morlet wavelet transform of signals obtained at receiver 1 and 2. |
| (2) | Determine phase and group delays in frequency at fc frequency which is the centre frequency of complex Morlet wavelet. |
| • | The group delays at receivers 1 and 2 which are tg1 and tg2 are obtained. The group delay is a time corresponding to the maximum of magnitude of ψ1 and ψ2 |
| • | From phase information of, ψ1 (t) is taken as a phase corresponding to tg1 |
| • | The tL and tR are obtained from phase information of ψ2. The tL is the time corresponding to θ1 which is the most close to tg2 on the left side of tg2 and tR is the time corresponding to θ1 which is the most close to tg2 on the right side of tg2 |
| • | tph1 is defined as tg1 and tg2 is either tL or tR depending upon which is closer to tg2 |
| (3) | To determine the phase and group delays at whole frequencies, the procedure (2) has to be repeated for all complex Morlet wavelet coefficients. |
| (4) | If the distance between receivers 1 and 2 is D, then the group velocity Vgr and phase velocity Vph at each frequency are obtained as follow: |
| (31) |
| (32) |
| (5) | The procedure is continued until all receivers corresponding signal is regarded iteratively for processing. |
RESULTS
The new method is more robust and persistent than conventional phase extraction method because the reduction of noises during analysis as well as easier and faster analysis. The dispersion curve extracted using proposed method is shown by Fig. 11 to estimate the performance. At frequencies larger than about 150 Hz, it can be noticed that the dispersion curve obtained by the conventional MASW method in Fig. 12a coincides well with the dispersion curve obtained by the proposed method in Fig. 11. With comparing the result of new method with previous finding relates through the high frequency corresponding low phase velocity. The good agreement of both finding is also supportive by the dispersive characteristic of soil because high frequency penetrates near surface layer consisting low shear wave velocity.
| Fig. 11: | Phase velocity versus frequency curve depicted from phase delay |
However, at low frequencies below 150 Hz, both dispersion curves do not match. This disagreement may be owing to the near field effect which comes from the coupling between surface and body waves as noises. The disagreement at very low frequency like below 10 Hz is not so concentric for the measurement of near surface layer profile. Among the frequencies about 50 to 100 Hz the irregularities is found regarding the dispersive trend in the conventional extracted phase velocity in Fig. 12a and this portion of previous findings contradicts with the phase velocity of those ranges in new findings of Fig. 11. According to conventional result in Fig. 12a, low phase velocity 200 m sec-1 is shown at low frequency 50Hz and phase velocity about 500 m sec-1 is pointed at 100 Hz. This outcome means deep layer consisting soft soil surrounded by hard soil which is not supportive by dispersive characteristics of soil. According to this usual MASW method, the irregularities at frequency range 50-100 Hz in Fig. 12a is amended by curve fitting in Fig. 12b. The adjusted phase velocity by curve fitting criteria in Fig. 12b is well suited with the extracted phase velocity by proposed method in Fig. 11 except the very low frequencies below 10 Hz. For the near surface layer profile, the phase velocity at very low frequencies is regardless for the half-space constrains.
The significant outcome for the new method is representation of consistent dispersion curve avoiding irregularities and noise impressions. The reliable dispersion curve is the key success of this study because this stage is ineffort to characterize soil profile through inversion method. The dispersive outcome of the proposed method is sustained with outcome of MASW method regarding noise filter criteria. The presence of noises during signal processing is kept away from the extracted phase velocity based on new method.
| Fig. 12: | (a) Extracted phase velocity using conventional surface wave method, (b) adjusted dispersion curve with erroneous |
Furthermore, the analysis of the proposed method is robust and persistent for using the wavelet transformation. The outcome of this new method is obtained easily and rapidly avoiding Fourier transformation. The main benefit from this high measure of redundancy in the WT is an accurate time-frequency spectrum.
DISCUSSION
The performance of new method is estimated by comparing the dispersive output with that of conventional MASW method. For the accuracy and pioneering performance, the new MASW method possesses several advantages for different Geotechnical application. Anomalies that include fracture zones within bedrock, dissolution/potential subsidence features, voids associated with old mine works and erosional channel eched into the bedrock surface have been effectively identified by proposed MASW method. Advantages of using wavelet transform in MASW method are to detect, delineate, or map anomalous subsurface materials including insensitivity of cultural noise, ease of generating and propagating surface wave energy and its sensitivity to change in velocity.
The dispersion curve represents phase velocity corresponding to frequencies in geotechnical characteristics for surface wave analysis. Multi-station recording permits a single survey of a broad depth range, high levels of redundancy with a single field configuration and the a bility to adjust the offset (Lin and Chang, 2004). Park et al. (1999a, b, 2000) extracted the dispersion curve with modal discrimination using frequency domain analysis for MASW method. This development is also important for representing the mode separation in dispersive outcome. But, the noise reduction performance is not improved through this conventional MASW method. Miller et al. (2000) has compared shear-wave velocity profiles from usual MASW method with Borehole Measurements and resulted in differences between inverted S-wave velocities between the MASW method and borehole measurements to be as low as 18%, with potential improvement as low as 9%. The overall difference between Shear (S) wave velocities derived from the MASW (multi-channel analysis of surface wave) technique and borehole measurements shown by Xia et al. (2002a, b) is about 15%.
Due to the analysis corresponding to the new method, the extracted phase velocity data is achievable to compare and adjust if any irregularities occur. After the comparing and irregularities adjustment, the final dispersion curve shown in Fig. 11 is obtained through proposed method regarding complex Morlet wavelet transform. Usually, the noise reduction technique or irregularities adjustment is not maintained for the final phase velocity representation (Park et al., 2000).
However, in the conventional MASW method, the phase velocity is extracted using Fourier transformation and cross spectral density corresponding to all received signal which is tedious for analysis. Some noise and interferences added during Fourier transformation shown in Fig. 13 as well as wrapped and unwrapped phase angle extraction shown in Fig. 14. Furthermore, the Fourier transformation of each receiver corresponding signal, cross spectral density between each pair of receivers and additional filtering technique of noises make the conventional method slower.
On the contrary, the proposed method is robust with faster, easier and noise reduction performance during analysis for using WT avoiding Fourier transformation and conventional phase velocity extraction.
| Fig. 13: | Representation of noises during Fourier transformation of seismic signal |
| Fig. 14: | Presence of noises in (a) wrapped phase angle, (b) unwrapped phase angle |
For the utilization of the wavelet transform in multichannel analysis, the performance is superior because wavelet theory is one of the most successful tools to analyze, visualize and manipulate complex non-stationary data, for which the traditional Fourier methods cannot be applied to directly (Yves et al., 1992). In each of these fields, the wavelets are applied for data compression, noise removal, feature extraction, classification and regressions face (Saito, 2004).
For multichannel analysis in the surface wave method with Morlet wavelet transform, the phase velocity outcomes for each pair of receivers are compared. In the WT, for each value of the scale used, the correlation between the scaled wavelet and successive segments of the data stream is computed. Unless reconstruction is needed, there are no restrictions in the WT as to how many scales are used, nor of the spacing between the scales. If needed, a high resolution spectrum can be generated for a narrow range of frequencies. The unique success of the proposed method is estimated by the rapidness, easiness and noise reduction performance during analysis rather than conventional surface wave method.
CONCLUSION
In this study, a wavelet decomposition technique is implemented on MASW method to obtain an evaluation of dispersion curve. A more stable and consistent outcome is obtained with reference to the de-noising performance of Morlet wavelet analysis. In the evaluation method, the phase and group velocity are determined in correspondence to multiple receivers at each frequency component using the information obtained around the time at which the signal energy is concentrated. Therefore, the new method is less affected by noise, effect of body wave or higher mode surface wave. Furthermore, the imaging of soil layer profile in geotechnical subsurface exploration can be developed inaugurating wavelet analysis in MASW method.
ACKNOWLEDGMENTS
This research is sponsored by Research Project of Science Fund No. 01-01-02-SF0338 from Ministry of Science, Technology and Innovation of Malaysia.