Research and application of Rayleigh wave 3-D phase velocity imaging based on wavelet transform

Rayleigh wave detection is a recently developed method for shallow seismic exploration. Current Rayleigh wave data processing and interpretation methods can only provide the transverse average wave velocity of rock-soil bodies under the geophone array range, resulting in a low lateral resolution of wave velocity. To solve this problem, this paper presents a Rayleigh wave data processing method based on wavelet transform. First, the Hankel matrix is constructed from the intercepted Rayleigh wave, and the effective singular value is preserved by singular value decomposition to filter the Rayleigh wave. Then, the appropriate center frequency is selected and the corresponding relationship between the time and frequency of the Rayleigh wave is obtained via wavelet transform. The waveform of each frequency component can be extracted and the complete time difference of each frequency component between two geophones will be obtained and used to calculate the phase velocity-depth profile of the Rayleigh wave in a rock-soil body. This method is applied to examine unfavorable geological bodies that are underground in a yard. By combining the phase velocity-depth profiles of several survey lines, the 3-D image of phase velocity of Rayleigh wave underground can be obtained. This method can provide the phase velocity distribution of the formation below the survey line by only one measurement, which greatly improves upon the work efficiency and lateral resolution of the traditional Rayleigh wave data processing method.


Introduction
Rayleigh wave exploration is a recently developed method of shallow seismic exploration. A Rayleigh wave is a type of surface wave that propagates along a free surface [Rayleigh, 1885]. It is a secondary wave composed of two types of body waves, namely longitudinal and shear waves (SV wave), whose energy mainly concentrates near the surface [Chen, 2010]. Rayleigh waves have many vibration modes, and different vibration modes correspond to different frequencies, wavelengths, and velocities. That is, Rayleigh waves propagate in media with dispersion characteristics. Simultaneously, the effective penetration depth of a Rayleigh wave with different frequencies is proportional to its wavelength [Haskell, 1953]. The relationship between the phase velocity and detection depth of a Rayleigh wave can be established from the correlation between wave velocity, frequency, wavelength, and effective penetration depth, thus yielding the wave velocity of rock-soil bodies at different depths.
Scientific research on Rayleigh wave exploration began in the 1970s [Chang, 1973]. In 1982, the VIC Corporation of Japan developed the GR-810 Sato automatic underground prospecting machine based on steady state Rayleigh wave prospecting to solve shallow engineering geological problems; this was the first practical application of Rayleigh wave exploration [Sato, 1986]. In 1996, Liu Yunzhen and others constructed their own multi-channel seismic data acquisition and processing system. It can increase the effective depth of Rayleigh wave exploration to 30 m and even reach 50 m under good conditions, which basically met the requirements for geotechnical engineering investigation. Rayleigh wave exploration is typically divided into two types, namely the transient method and the steady state method. Because the operational efficiency of the transient method is much higher than that of steady state method, transient Rayleigh wave exploration is widely used in engineering practice [Liu et al., 1996;Wang et al., 2000;Zhang et al., 2001;Chen et al., 2006;Pan et al., 2011Pan et al., , 2015Cui et al., 2012]. The Rayleigh wave excited by transient impulse contains a variety of frequency components, which overcomes the inefficiency of the steady state method by which one excitation can only generate a single-frequency signal.
However, because of this, there are many difficulties in processing and interpreting transient Rayleigh wave data.
Presently, multi-channel transient Rayleigh wave exploration can only calculate the lateral average velocity of rocksoil body under the geophone arrays. Furthermore, problems such as low efficiency and insufficient lateral resolution of the rock-soil body remain [Qi, 2001]. To address the weaknesses of the existing multi-channel transient Rayleigh wave exploration, this paper presents a time difference method based on wavelet transform. From the overall time difference of Rayleigh wave waveforms of different frequencies received by two geophones, the phase velocity of rock and soil mass between the two geophones is calculated. By combining several survey lines, 3-D wave velocity imaging can be completed.

Rayleigh wave data processing based on phase difference method
Presently, the phase difference method is the most commonly used in Rayleigh wave data processing. f 1 (t) and f 2 (t) are the time-domain records of two different observation points. By applying Fourier transform to f 1 (t) and f 2 (t), the spectrum ₁( ) and ₂( ) can be obtained. Thus, the cross power spectra of ₁( ) and ₂( ) can be expressed as follows: In the formula, ₁ * ( ) and ₂ * ( ) are the complex conjugations of ₁( ) and ₂( ), respectively. From Formula (1), the phase difference can be obtained as follows: (2) Under the premise of satisfying the spatial sampling theorem (group interval should be less than or equal to half of the wavelength), the phase velocity of the Rayleigh wave can be calculated as follows: (3)

Introduction of wavelet transform
For signal f(t), the expression of wavelet transform can be written as follows: where a is the scale parameter and a > 0; b is the translation parameter. ( ) is a function called mother wavelet or basic wavelet. Defining wavelet basis function , ( ), Then, the wavelet transform of f(t) can also be written as follows: Let the Fourier transform of ( ) be ( ) and the Fourier transform of ( ) be ( ); then, the frequency domain expression of the wavelet transform is (7) In the wavelet transform, the scale parameter a is used to produce a dilating transform for ( ), and the translation parameter b is used to determine the time position of the analysis. Thus, combining a and b helps determine the central position and time width of the wavelet transform for ( ). Therefore, the wavelet transform can be understood as a method for analyzing signals with a family of wavelet basis functions, whose central position and time width are constantly varying. This feature solely satisfies the basic requirement that different components of the signal require different resolution analyses.
From the definition of the abovementioned wavelet transform, the mother wavelet ( ) directly affects the result of the transform. It is generally constructed artificially according to certain requirements and exhibits diversity. At present, several mother wavelets, such as the following, are commonly used: Haar wavelet, Morlet wavelet, Mexican Hat wavelet, Daubechies wavelet, and Meyer wavelet.
The instantaneous frequency of the signal can be obtained using the wavelet transform; in other words, the corresponding time of each frequency component of each geophone can be obtained. Because the group interval is then known, the velocity distribution of rock-soil body can be calculated using a time difference method.

Rayleigh wave 3-D processing method based on wavelet transform
The traditional Rayleigh wave data processing calculates the dispersion curve using the phase difference method. In general, only one dispersion curve (Figure 1) can be obtained by one survey line, which causes the low lateral resolution.
In this paper, a new Rayleigh wave processing method based on the wavelet transform is proposed. Using this method, velocity-depth profiles can be obtained by one survey line; further, if several surveying lines are combined, underground 3-D velocity distribution can also be obtained, which helps solve the problem of low lateral resolution of traditional methods.

Interception of Rayleigh Wave
To extract the dispersion curve of Rayleigh wave more accurately, it is necessary to separate Rayleigh waves from seismic waves. This step can be accomplished directly using the seismic waveform record. Compared with the velocities of P-wave and S-wave, that of the Rayleigh wave is slower and its amplitude is larger. Based on this, the Rayleigh wave can be intercepted ( Figure 2).
After interception, the first point of the first waveform should be considered as the initial time, and the last point of the last waveform as the end time. The missing part of each waveform should be compensated by 0 to ensure that the length of each waveform is the same, to facilitate the subsequent calculation.

Rayleigh wave filtering based on improved singular value decomposition
Shallow seismic exploration is susceptible to the influence of external environments, causing noise in the seismic waveform. Simultaneously, the velocity of Rayleigh wave is similar to that of shear wave, and the Rayleigh wave obtained is only intercepted in the time domain; thus, it is difficult to avoid mixing noise with the shear wave.
Therefore, it is necessary to weaken the interference. Presently, the most commonly used method is to filter the intercepted Rayleigh wave in the frequency domain. In general, it is believed that the random noise is of high frequency; thus, a low pass filter is constructed to remove the high frequency components. However, there are two problems in frequency domain filtering. First, the cut-off frequency is difficult to judge: the extremely high cut-off frequency will mix the effective signal with noise; the extremely low cut-off frequency will render the effective signal incomplete. Second, the effective signal also has high-frequency components, and the noise also has lowfrequency components. The frequency domain filtering judges whether the signal components are effective or not only by the cut-off frequency, which cannot directly distinguish noise from an effective signal. In addition, singular value decomposition is widely used in noise reduction, filtering, and other related processes [Regadíoa et al., 2019;Yue et al., 2019]. In this study, an improved singular value decomposition is introduced to desiccate the intercepted Rayleigh waves. First, assuming that the Rayleigh wave intercepted in 3.1 is (n), the sampling points are (1), (2), (3)…, (N). A Hankel matrix can be constructed according to the sampling points as follows: The singular value decomposition of is performed.
where U is a unitary matrix of order N-L, is a unitary matrix of order L+1 and is a positive semi-definite matrix of order (N-L)·(L+1) whose diagonal line comprises singular values of and other elements are 0.
For a signal, the bigger the singular value, the more it can reflect the characteristics of the signal. Therefore, the core idea of denoising via singular value decomposition is to delete smaller singular values and retain larger singular values. In general, a singular value is validated by comparing it with the maximum singular value. Let the singular value be σ w , w = 1, 2, 3…, and the maximum singular value is σ max .
Setting threshold θ 0 , if θ>θ0, the corresponding singular value represents the effective component of the signal that should be retained, and if θ<θ0, the corresponding singular value represents the noise of the signal that should be discarded. This method uses the maximum singular value as the reference standard to validate the other singular values. However, for different seismic signals from different detectors, the maximum singular value varies, rendering the criteria inconsistent for choosing the effective singular value. Simultaneously, θ0 still exists. In addition, θ0 are also difficult to determine and unite. In the existing research, θ0 = 10 -2 or 10 -3 ; however, neither can be applied to all situations. Different seismic signals are often required to be analyzed one by one.
In general, in the process of mathematical operation, the singular values will be arranged from large to small, and eventually tend to zero. In other words, for the Rayleigh wave from a detector, if the singular values are marked in rectangular coordinates, and then connected into a smooth curve, the curve will approach a straight line. Thus, starting from a point, the first derivative of the curve will be approximately equal to 0. These points are the invalid singular values. The invalid singular values are judged using the following equation: σ′ = 0 (σ′ should be round up and round down numbers) where I = c, c+1, c+2…, σ c is the first invalid singular value.
For the Rayleigh wave from the first seismic trace in Figure 3, the first derivative of the singular value change curve is drawn in Figure 4. that the figure shows that, starting from the eleventh point, i.e., c = 11, the first derivative of the singular value change is always approximately 0. These singular values can be considered invalid.
The singular values are validated using the above two methods, as shown in Figure 5. The point marked with red denotes the result of the common method: when θ0 = 10 -2 , c = 7; when θ0 = 10 -3 , c = 59. The point marked with green denotes the result of the method proposed in this paper: c = 11. If c = 59, the singular values near σ 59 are generally small and no significant difference is observed among them. Thus, it is impossible to remove the noise accurately. If c = 7, there are still large singular values behind the point, which may cause the loss of effective signals. In contrast, c = 11 is a more accurate boundary between larger and smaller singular values. Therefore, the method proposed in this paper to validate singular values is more effective for the Rayleigh wave.
After determining the effective singular value, only the first c-1 row and c-1 column of matrix are retained, which is denoted as ′, a c-1 order diagonal matrix; the first c-1 column of matrix U are retained, which is denoted as ′, a (N-L)·(c-1) order matrix; and the first c-1 row of matrix V are retained, which is denoted as ′, a (c-1)·(L+1) order matrix. The matrix ′ is calculated as follows: ′ is a matrix of order (N-L)·(L+1). The denoised Rayleigh wave y'(n) can be reconstructed by combining the first row of ′ and the last column of ′ that the first row is deleted. As shown in Figure 6, by singular value decomposition filtering, the noise is weakened, and the waveform is smoothed.
Xiang Min et al.

Rayleigh wave data processing based on wavelet transform
In previously reported studies, the instantaneous frequency of the signal was calculated by wavelet transforms and the dispersion curve was subsequently deduced with the instantaneous frequency.
Usually, the Morlet wavelet is used to analyze Rayleigh waves. The general expression of the Morlet wavelet is as follows: (13) where c = -¼ .
When the center frequency of the Morlet wavelet is 1 Hz, the time-frequency distribution of each Rayleigh wave can be obtained, as shown in Figure 7.
The instantaneous frequency of each component is determined by peak detection, that is, by detecting the maximum amplitudes of each of the peaks, whose corresponding time is transformed to the frequency component.
However, in practice, this causes problems in accurate estimation. The Heisenberg uncertainty principle is satisfied between time and frequency, that is, time and frequency are a pair of contradictions, and absolutely accurate time and frequency cannot be obtained simultaneously. Thus, it is impossible to determine the accurate time of a certain frequency component, and only the correspondence between the time periods and frequency bands of signal components can be obtained. Therefore, the instantaneous frequency in the absolute sense cannot be calculated.
Taking the first Rayleigh wave in Figure 3 as an example, Morlet transforms with the center frequencies of 1, 3, 5, and 7 Hz are used to compute the time-frequency distribution.
As shown in Figure  the energy of the signal cannot be concentrated in the time and frequency domains at the same time. The shorter the frequency band, the longer is the corresponding time period, and the time resolution decreases. Therefore, accurate instantaneous frequency cannot be obtained by peak detection; hence, the dispersion curve cannot be calculated accurately as well.
Considering the Heisenberg uncertainty principle, this work extends the peak detection method for deriving the dispersion curve from instantaneous frequency to time-period frequency, i.e., replacing the time variation of the maximum amplitude with the complete time variation of a frequency component.
First, it is necessary to determine the center frequency of the Morlet wavelet. The center frequency needs to ensure that it has the highest frequency resolution under the premise that the components can be resolved in the direction of the time axis, that is, when Rayleigh waves of different seismic traces can be distinguished in the direction of the time axis, the maximum central frequency should be selected to obtain the optimal frequency resolution. In practice, as shown in Figure 8, when the center frequency of the Morlet wavelet is 5 Hz, the component near 40 Hz is spread throughout the figure, which will render the time characteristics of this component unobvious, and it is difficult to extract the time differences between each seismic trace. In comparison, the center frequency of 5 Hz is appropriate.
The Raleigh waves in Figure 3 are processed by the Morlet wavelet with a center frequency of 5 Hz to obtain Figure 9.
where N is the length of time series ₁(n) and ₂(n).
The similarity r is calculated once every translation of ₁(n). When after translating m 0 times and r is minimum, ₁(n) and ₂(n) are most similar. Let the sampling interval be t 0 . The time differences of the components with the frequency of f 1 in the Rayleigh wave can be considered as m 0 q 0 t 0 between the first and second seismic traces.
After obtaining the time difference of each frequency component between every seismic trace, the phase velocity can be calculated by (15).
where represents the group interval. Then, the relationship between phase velocities, frequencies, and wavelengths is obtained.
Because the wavelength of the Rayleigh wave is proportional to the detection depth, the ratio coefficient called depth coefficient is recorded as .
where ℎ represents the detection depth, represents the wavelength, and is a constant related to the Poisson's ratio of the medium (see Table 1).  It should be noted that the space sampling theorem requires that the group interval should be less than half the wavelength.
where Δ represents group interval and represents the wavelength.
The effective detection depth of the Rayleigh wave can thus be calculated by (18).
According to the above steps, the phase velocities of each frequency component between each seismic trace can be calculated, and the phase-velocity-depth profile of the Rayleigh wave can be obtained (as shown in Figure 11).
Compared with the traditional Rayleigh wave processing (as shown in Figure 1), this method presented here has better lateral resolution and can directly obtain the phase velocity distribution of the formation below the arrangement through one shot excitation, which greatly improves the work efficiency.

Engineering Application
The method introduced in this paper is applied to the detection of abnormal bodies in the pipe jacking construction of a yard. Ten parallel survey lines are arranged, and the distance between each pair of lines is 2 m.
Rayleigh wave detection uses an 18-pound hammer to excite the waves. Each line is equipped with 10 detectors with a group interval of 2 m.  The 3-D phase velocity distribution of the formation below the yard can be obtained by processing the Rayleigh waves of 10 surveying lines with the method described above (as shown in Figure 12). There is a corresponding relationship between the type of rock-soil body and the phase velocity of the Rayleigh wave. When the wave velocity is lower than 800 m/s, the geological body is usually soil or soft rock. When the wave velocity is higher than 800 m/s, the geological body is hard and complete rock, which is not conducive to pipe jacking. Accordingly, adjusting the angle of the velocity distribution figure, the position suitable for pipe jacking can be selected.

Conclusion
In this paper, a Rayleigh wave data processing method based on wavelet transform is proposed and verified by an engineering example. The following conclusions are drawn: 1) the filtering method based on improved singular value decomposition proposed in this paper can distinguish effective signals from noises; 2) in this paper, the shortcomings of spectrum peak detection are explained by the Heisenberg uncertainty theory.
Meanwhile, it is proposed to replace the time variation of the peak with the time variation of the whole waveform of each frequency component, and then to calculate the variation of phase velocity with depth between the adjacent geophones; 3) the phase velocity-depth profiles of the formation below the seismic array can be obtained by this method.
Compared with the traditional Rayleigh wave processing, only one dispersion curve can be gotten by a seismic array, which is more intuitive and efficient, and has higher lateral resolution; 4) introducing this method to solve an engineering problem, combing the phase velocity-depth profiles of several surveying lines, the 3-D phase velocity underground can be imaged, which can judge the geological body that is not conducive to construction, and verify the feasibility of this method.
Xiang Min et al.