Polarization Filtering Method for Suppressing Surface Wave in Time-Frequency Domain ()
1. Introduction
In seismic exploration, surface waves are typical regular disturbances, and surface waves are usually suppressed based on the apparent velocity and frequency difference between the surface wave and the effective wave. In the three-component seismic exploration, the converted wave energy is weak. In the shallow layer, the converted wave apparent velocity overlaps with the surface wave. In the deep layer, the converted wave frequency overlaps with the surface wave. Therefore, the converted wave is greatly disturbed by surface waves. According to the different polarization characteristics between body wave and surface wave, a polarization filtering method is studied to suppress the surface wave.
Initially, Shimshoni and Smith separated the signals in natural and nuclear explosions according to the polarization characteristics of seismic waves [1] . Flinn calculated the three-component seismic data covariance matrix in a certain time window, constructed a series of polarization parameters that can describe the polarization characteristics of seismic waves, and separated the wave fields with different polarization characteristics. The disadvantage of this method is that it depends on the selection of window length, and the polarization parameters are not time-varying [2] . Samson proposed the concept of polarization coefficient, which provides a basis for designing various polarization filter factors. With the application of three-component technique in practice, geophysicists have developed a variety of polarization analysis techniques [3] . Benhama, Jurkevics, Shieh C, Huang, Chen, Liu, etc. perform polarization analysis in the time domain and frequency domain [4] - [9] . Christoffersson proposed the maximum likelihood method [10] . Jackson et al. proposed the singular value decomposition method [11] . Vidale, Rene, Stewart, Morozov, etc. developed instantaneous polarization analysis, etc. [12] [13] [14] [15] . Huang Zhongyu and Gao Lin studied the three-component polarization analysis and its practical application [7] . Zhang Yufen and Zhou Jianxin analyzed the factors affecting filtering effect in the spatial direction [16] . Tatham et al. pointed out that due to the complexity of underground media, the use of the elliptic rate alone does not effectively separate the surface waves [17] . The above mentioned methods all use fixed time windows, so the effect of suppressing surface waves is not ideal.
In order to solve the problem of time windows, this paper proposes an adaptive window function, the length of which is adaptive to the instantaneous frequency of the 3D seismic record, so polarization characteristic parameters can be obtained at each time sampling point of the three-component seismic record. In this time-varying time window, the wavelet transform is used to construct the time-frequency domain instantaneous frequency equation and the time-frequency domain covariance matrix. The combination of the elliptic rate and the dip polarization parameter is used as the constraint condition to construct the filter function to separate the surface wave. At last, the separated surface wave was transformed into the time domain. Then, it is adaptively subtracted from the raw data to suppress the surface wave. The actual data processing result shows that this method protects the effective signal while suppressing the surface wave.
2. Method Principle
2.1. Wavelet Transform
The wavelet transform decomposes signals x(t) by telescopic translation in multi-scales, and finally achieves the result that the time subdivision at high frequency and the frequency subdivision at low frequency [18] .
(1)
Let
(2)
where
is the wavelet transformed signal; wavelet g(t) produces a “time-frequency” window that changes with frequency by stretching a and translation b, a is the frequency scale factor, b is the time factor, the symbol * indicates the conjugate complex number. The wavelet transform of seismic signals generally uses Morlet wavelets:
(3)
where e is a natural constant; i is an imaginary number; f0 is the center frequency.
The data
can be recovered by the formula (4).
(4)
where
,
is the Fourier transform of wavelet g(t).
The scale factor a is related to the center frequency f0 of g(t) and can be related to the physical frequency by a = f0/f. The scale factor a can stretch the basic wavelet, the wavelet width is proportional to a, and the amplitude is proportional to 1/a, but the shape (area) of the wavelet remains unchanged. The larger the scale, the longer the wavelet function is in time. Then the wavelet function processes the longer part of the signal, and the low frequency characteristic of the signal is obtained. On the other hand, when the scale factor is smaller, the wavelet function only processes the shorter part of the signal, so the high frequency characteristic of the signal is obtained.
2.2. Parameters Calculation
For the three-component seismic signal ui(t) (i = x, y, z), its wavelet transform in the time-frequency domain is
. When the time is b, its local approximation can be expressed as
, this can be expressed as formula (5) by the modulus and phase of the wavelet transform:
(5)
where
is the instantaneous frequency function of i component at scale a and time b, it is defined as:
(6)
In order to avoid numerical derivation, the wavelet g(t) and its derivative is used to obtain the exact equation of
.
(7)
where Im[.] represents the complex imaginary part.
Using the wavelet transform
of three-component seismic data, we can construct a covariance matrix in the time-frequency domain.
(8)
Here
(9)
where k, m = (x, y, z), μ is the mean of the k component in the time window:
(10)
where Re(∙) represents a complex imaginary part, Sinc(∙) represents a singular function.
is the adaptive time window length determined by the instantaneous frequency of time b and scale a. In fact, different components of the signal often have different instantaneous frequencies, so it is necessary to calculate an average time window
by the instantaneous frequencies of multiple components.
(11)
where N is a positive integer, N = 1 or 2;
is the average of the instantaneous frequencies of the three components.
The covariance matrix
is defined at each time-frequency point in the time-frequency domain, and its physical meaning represents the instantaneous energy at each time-frequency point. Calculate its eigenvalues λi(i = 1, 2, 3; λ1 > λ2 > λ3) and the corresponding eigenvector Vi. The instantaneous polarization parameters in the time-frequency domain are derived from the results, and the energy magnitude and direction of the particle in the time window are described. For linearly polarized waves, the covariance matrix has only one non-zero eigenvalue, and the polarization direction of the seismic wave can be judged by the corresponding eigenvector. For seismic waves vibrating in one plane, the covariance matrix has two non-zero eigenvalues, which determine the polarization plane. In the actual data, the polarization direction exists in the entire 3 space, so the three eigenvalues are not zero. The largest eigenvalue λ1(b, a) corresponds to the eigenvector V1(V1,x, V1,y, V1,z) as the main eigenvector. Among the three eigenvalues, the seismic signal energy is focused on the largest eigenvalue.
The polarization parameters calculated from the eigenvalues and eigenvectors of the time-frequency domain covariance matrix are mainly as follow:
1) Instantaneous polarization vector
Instantaneous main polarization vector
(12)
Instantaneous intermediate polarization vector
(13)
Instantaneous sub-polarization vector
(14)
2) Ellipticity
Primary ellipticity
(15)
Secondary ellipticity
(16)
3) Dip angle
(17)
4) Azimuth
(18)
2.3. Filter Function
After analysis and testing, the ellipticity
and dip angle
can separate the body wave and the surface wave more effectively. This paper combines these two parameters to construct a filter function.
(19)
where Fe is the polarization filter function applied to each of the original components, Pρ and Pθ define the range of variation of ρ and θ. Set the appropriate Pρ and Pθ to obtain the desired output. In the actual data processing, the parameters were tested multiple times to ensure that the surface wave was suppressed properly.
3. Filtering Steps
According to the polarization filtering principle and the actual data, the processing flow of the polarization filtering suppression surface wave is summarized in Figure 1.
The specific steps of polarization filtering are as follows:
1) Define the surface wave suppression window
Since the propagation characteristics of surface waves are only related to the geological features of the surface, the surface wave velocity of the same area are almost consistent. On a single shot record, the surface waves are fan-shaped. According to this feature of the surface wave, the window area for suppressing surface waves can be defined as follow:
(20)
where t(x) is the start time of the surface wave suppression; x is the offset; t0 is the start time of the time window; v is the apparent velocity of the surface wave.
2) Wavelet transform, convert the seismic data to the time-frequency domain.
3) Generate the covariance matrix of multi-component seismic trace.
4) Calculate the polarization characteristic parameters in the time-frequency domain.
5) Use polarization filtering to separate surface wave. The body wave propagating in a homogeneous isotropic medium is a linearly polarized wave, and the surface wave is an elliptical polarized wave. However, in the case of anisotropic medium or interference, the particle vibration of the wave becomes very complicated. The wave often appears as an irregular ellipse with a large flattening rate. Therefore, in the polarization filtering, the polarization filtering parameters should be reasonably set according to the characteristics of the actual data.
6) Inverse wavelet transform, convert the surface wave energy to the time domain.
7) Subtract the surface wave from the original record adaptively.
Directly subtracting the surface wave in the time-frequency domain will be too stiff and may damage the effective signal. So we separate the surface wave
Figure 1. The flowchart of polarization filtering.
from raw converted wave in the time-frequency domain, then the surface wave was subtracted from the raw record adaptively in the time domain.
4. Examples
Figure 2(a) is a shot record of converted wave acquired in a certain area of Saudi Arabia. The surface wave has a wide spatial distribution and strong energy in the original single shot, which seriously pollutes the effective signal. Figure 2 and Figure 3 show the single shot and stack profiles before and after applying the polarization filtering proposed in this paper. Before the noise is suppressed, the event on the stack profile is covered by the surface wave energy. It is difficult to identify a continuous event, which has a great impact on subsequent processing. After the polarization filtering, as shown as Figure 2(b) and Figure 3(b), both the shot and stack profiles have clear events, and the signal-to-noise ratio is greatly improved. As we can see from Figure 2(c) and Figure 3(c), the noise suppressed, the polarization filtering does not damage the effective reflection energy.
Figure 4 shows the time-frequency spectrum before and after the surface wave attenuation. Before the surface wave is suppressed, there is noise energy between 0 - 15 Hz in the low frequency band. After the surface wave is attenuated, the noise energy of the frequency band is greatly reduced.
5. Conclusions
Compared with the conventional surface wave suppression method, polarization filtering can effectively separate converted waves and surface waves. Due to
Figure 2. The single shot comparison before and after polarization filtering (a) shot gather before polarization filtering; (b) shot gather after polarization filtering; (c) and the surface wave removed.
Figure 3. The stack section comparison before and after polarization filtering (a) stack section before polarization filtering; (b) stack section after polarization filtering; (c) the surface wave removed.
factors such as underground geological conditions, the polarization characteristics of body waves or surface waves are not ideal linear polarization or elliptical polarization. It is necessary to test the polarization parameters based to construct an appropriate filter function.
Figure 4. Comparison of spectrum between before and after polarization filtering: (a) the spectrum before surface wave attenuation; (b) the spectrum after surface wave attenuation.
By transforming the separated surface wave from time-frequency domain to the time domain, and then subtracting the surface wave from the raw record in time domain adaptively, this processing make the denoising effect more natural.