Polarization Filtering Method for Suppressing Surface Wave in Time-Frequency Domain

Show more

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] .

$WT\left(b,a\right)=\left({g}_{b,a},x\right)={\displaystyle {\int}_{-\infty}^{+\infty}\frac{1}{a}{g}^{\ast}\left(\frac{t-b}{a}\right)x\left(t\right)\text{d}t}$ (1)

Let

${g}_{b,a}\left(t\right)=\frac{1}{a}g\left(\frac{t-b}{a}\right)$ (2)

where $WT\left(b,a\right)$ 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:

$g\left(t\right)={\text{e}}^{2\text{\pi}i{f}_{0}t}{\text{e}}^{-\frac{{r}^{2}}{2{a}^{2}}}$ (3)

where e is a natural constant; i is an imaginary number; f_{0} is the center frequency.

The data $WT\left(b,a\right)$ can be recovered by the formula (4).

$x\left(t\right)=\frac{1}{{c}_{g}}{\displaystyle {\int}_{-\infty}^{+\infty}\frac{1}{{a}^{2}}g\left(\frac{t-b}{a}\right)WT\left(b,a\right)\text{d}b\text{d}a}$ (4)

where ${c}_{g}={\displaystyle {\int}_{0}^{+\infty}\frac{{\left|\stackrel{\xaf}{g\left(f\right)}\right|}^{2}}{f}\text{d}f}$ , $\stackrel{\xaf}{g\left(f\right)}$ is the Fourier transform of wavelet g(t).

The scale factor a is related to the center frequency f_{0} of g(t) and can be related to the physical frequency by a = f_{0}/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 u_{i}(t) (i = x, y, z), its wavelet transform in the time-frequency domain is
$W{T}_{i}\left(b,a\right)$ . When the time is b, its local approximation can be expressed as
${W}_{i}\left(b+\tau ,a\right)$ , this can be expressed as formula (5) by the modulus and phase of the wavelet transform:

${W}_{i}\left(b+\tau ,a\right)\approx \left|W{T}_{i}\left(b,a\right)\right|\cdot \mathrm{cos}\left[{\Omega}_{i}\left(b,a\right)\tau +\mathrm{arg}W{T}_{k}\left(b,a\right)\right]$ (5)

where ${\Omega}_{i}\left(b,a\right)$ is the instantaneous frequency function of i component at scale a and time b, it is defined as:

${\Omega}_{i}\left(b,a\right)=\frac{\partial}{\partial b}\mathrm{arg}W{T}_{i}\left(b,a\right)$ (6)

In order to avoid numerical derivation, the wavelet g(t) and its derivative is used to obtain the exact equation of ${\Omega}_{i}\left(b,a\right)$ .

${\Omega}_{i}\left(b,a\right)=-\mathrm{Im}\left[\frac{W{T}_{i}^{{g}^{\prime}}\left(b,a\right)}{W{T}_{i}^{g}\left(b,a\right)}\right]$ (7)

where Im[.] represents the complex imaginary part.

Using the wavelet transform $W{T}_{i}\left(b,a\right)$ of three-component seismic data, we can construct a covariance matrix in the time-frequency domain.

$M\left(b,a\right)=\left[\begin{array}{l}{I}_{xx}\left(b,a\right){I}_{xy}\left(b,a\right){I}_{xz}\left(b,a\right)\\ {I}_{xy}\left(b,a\right){I}_{yy}\left(b,a\right){I}_{yz}\left(b,a\right)\\ {I}_{xz}\left(b,a\right){I}_{yz}\left(b,a\right){I}_{zz}\left(b,a\right)\end{array}\right]$ (8)

Here

$\begin{array}{l}{I}_{k,m}\left(b,a\right)=\left|W{T}_{k}\left(b,a\right)\right|\left|W{T}_{m}\left(b,a\right)\right|\times \{\text{sinc}\left[\frac{{\Omega}_{k}\left(b,a\right)-{\Omega}_{m}\left(b,a\right)}{2}{T}_{k,m}\left(b,a\right)\right]\\ \text{}\times \mathrm{cos}\left[\mathrm{arg}W{T}_{k}\left(b,a\right)-\mathrm{arg}W{T}_{m}\left(b,a\right)\right]\\ \text{}+\text{sinc}\left[\frac{{\Omega}_{k}\left(b,a\right)-{\Omega}_{m}\left(b,a\right)}{2}{T}_{k,m}\left(b,a\right)\right]\\ \text{}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\begin{array}{c}\text{\hspace{0.05em}}\\ \text{\hspace{0.05em}}\end{array}\times \mathrm{cos}\left[\mathrm{arg}W{T}_{k}\left(b,a\right)+\mathrm{arg}W{T}_{m}\left(b,a\right)\right]\}-{\mu}_{k}{\mu}_{m}\end{array}$ (9)

where k, m = (x, y, z), μ is the mean of the k component in the time window:

$\begin{array}{l}{\mu}_{k}\left(b,a\right)=\frac{1}{{T}_{k,m}\left(b,a\right)}{\displaystyle {\int}_{-T\left(b,a\right)/2}^{T\left(b,a\right)/2}{W}_{k}\left(b+\tau ,a\right)\text{d}\tau}\\ \text{}=\mathrm{Re}\left[W{T}_{k}\left(b,a\right)\right]\text{sinc}\left[T\left(b,a\right){\Omega}_{k}\left(b,a\right)/2\right]\end{array}$ (10)

where Re(∙) represents a complex imaginary part, Sinc(∙) represents a singular function. $T\left(b,a\right)$ 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 $T\left(b,a\right)$ by the instantaneous frequencies of multiple components.

$T\left(b,a\right)=\frac{6\text{\pi}N}{{\Omega}_{x}\left(b,a\right)+{\Omega}_{y}\left(b,a\right)+{\Omega}_{z}\left(b,a\right)}=\frac{2\text{\pi}N}{{\Omega}_{av}^{xyz}\left(b,a\right)}$ (11)

where N is a positive integer, N = 1 or 2; ${\Omega}_{av}^{xyz}$ is the average of the instantaneous frequencies of the three components.

The covariance matrix
$M\left(b,a\right)$ 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 V_{i}. 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 V_{1}(V_{1,x}, V_{1,y}, V_{1,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

$R\left(b,a\right)=\sqrt{{\lambda}_{1}\left(b,a\right)}\frac{{V}_{1}\left(b,a\right)}{\Vert {V}_{1}\left(b,a\right)\Vert}$ (12)

Instantaneous intermediate polarization vector

${r}_{s}\left(b,a\right)=\sqrt{{\lambda}_{2}\left(b,a\right)}\frac{{V}_{2}\left(b,a\right)}{\Vert {V}_{2}\left(b,a\right)\Vert}$ (13)

Instantaneous sub-polarization vector

$r\left(b,a\right)=\sqrt{{\lambda}_{3}\left(b,a\right)}\frac{{V}_{3}\left(b,a\right)}{\Vert {V}_{3}\left(b,a\right)\Vert}$ (14)

2) Ellipticity

Primary ellipticity

$\rho \left(b,a\right)=\Vert {r}_{s}\left(b,a\right)\Vert /\Vert R\left(b,a\right)\Vert $ (15)

Secondary ellipticity

${\rho}_{\text{1}}\left(b,a\right)=\Vert r\left(b,a\right)\Vert /\Vert {r}_{s}\left(b,a\right)\Vert $ (16)

3) Dip angle

$\delta \left(b,a\right)=\mathrm{arctan}\left(\sqrt{{V}_{1,x}{\left(b,a\right)}^{2}+{V}_{1,y}{\left(b,a\right)}^{2}}/{V}_{1,z}\left(b,a\right)\right)$ (17)

4) Azimuth

$\alpha \left(b,a\right)=\mathrm{arctan}\left({V}_{1,y}\left(b,a\right)/{V}_{1,x}\left(b,a\right)\right)$ (18)

2.3. Filter Function

After analysis and testing, the ellipticity $\rho \left(b,a\right)$ and dip angle $\theta \left(b,a\right)=\text{\pi}/2-\delta \left(b,a\right)$ can separate the body wave and the surface wave more effectively. This paper combines these two parameters to construct a filter function.

$WT\left(b,a\right)={F}_{e}\left(b,a\right)W{T}_{i}\left(b,a\right),\text{}a>0$

${F}_{e}\left(b,a\right)=\{\begin{array}{l}1\text{}\rho \left(b,a\right)\in {P}_{\rho}\u4e14\theta \left(b,a\right)\in {P}_{\theta}\\ 0\text{others}\end{array}$ (19)

where F_{e} 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:

$t\left(x\right)={t}_{0}+\frac{x}{v}$ (20)

where t(x) is the start time of the surface wave suppression; x is the offset; t_{0} 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.

References

[1] Shimshoni, M. and Smith, S.W. (1964) Seismic Signal Enhancement with Three-Component Detectors. Geophysics, 29, 664-671.

https://doi.org/10.1190/1.1439402

[2] Flinn, E.A. (1965) Signal Analysis Using Rectilinearity and Direction of Particle Motion. Proceedings of the IEEE, 53, 1874-1876.

https://doi.org/10.1109/PROC.1965.4462

[3] Samson, J.C. (1977) Matrix and Stokes Vector Representation of Detectors for Polarized Wave-Forms. Geophysical Journal of the Royal Astronomical Society, 51, 583-603.

https://doi.org/10.1111/j.1365-246X.1977.tb04208.x

[4] Benhama, A., Cliet, C. and Dubesset, M. (1988) Study and Application of Spatial Directional Filtering in Three-Component Recordings. Geophysical Prospecting, 36, 591-613.

https://doi.org/10.1111/j.1365-2478.1988.tb02182.x

[5] Jurkevics, A. (1988) Polarization Analysis of Three-Component Array Data. Bulletin of the Seismological Society of America, 78, 1725-1743.

[6] Shieh, C. and Herrmann, R.B. (1990) Ground Roll: Rejection Using Polarization Filters. Geophysics, 55, 1216-1222.

https://doi.org/10.1190/1.1442937

[7] Huang, Z.Y., Gao, L., Xu, Y.M., et al. (1996) The Analysis of Polarization in Three-Component Seismic Data and Its Application. Geophysical Prospecting for Petroleum, 35, 9-16.

[8] Chen, Y., Zhang, Z.J. and Tian, X.B. (2005) Complex Polarization Analysis Based on Windowed Hilbert Transform and Its Application. Geophysics, 48, 889-895.

https://doi.org/10.1002/cjg2.736

[9] Liu, J.H., Liu, F.T. and Xu, Y. (2006) Polarization Analysis of Three-Component Seismic Data. Progress in Geophysics, 21, 6-10.

[10] Christoffersson, A., Husebye, E.S. and Ingate, S.F. (1985) A New Technique for 3-Component Seismogram Analysis. Semiannual Technical Summary, 2, 85-86.

[11] Jackson, G.M., Mason, I.M. and Greenhalgh, S.A. (1991) Principal Component Transforms of Triaxial Recordings by Singular Value Decomposition. Geophysics, 56, 528-533.

https://doi.org/10.1190/1.1443068

[12] Vidale, T. (1986) Complex Polarization Analysis of Particle Motion. Bulletin of the Seismological Society of America, 76, 1393-1405.

[13] Rene, R.M., Fitter, J.L., Forsyth, P.M., et al. (1986) Multi-Component Seismic Studies Using Complex Trace Analysis. Geophysics, 51, 1235-1251.

https://doi.org/10.1190/1.1442177

[14] Stewart, R.R. (1990) Ground-Roll Filtering Using Local Instantaneous Polarization. Crewes, 97, 28-33.

[15] Morozov, L.B. and Smithson, S.B. (1996) Instantaneous Polarization Attributes and Directional Filtering. Geophysics, 61, 872-881.

https://doi.org/10.1190/1.1444012

[16] Zhang, Y. and Zhou, J. (1999) Analysis of Factors Affecting Spatial Direction Filtering Effect. Oil and Gas Geology, 20, 212-215.

[17] Tatham, R.H. and McCormack, M.D. (1991) Multicomponent Seismology in Petroleum Exploration. Society of Exploration Geophysicists.

[18] Ma, J. and Li, Q. (2015) Suppression of Seismic Surface Waves by Time-Frequency Domain Polarization Filtering. Oil Geophysical Prospecting, 50, 1089-1097.