Seismic Periodic Noise Attenuation Based on Sparse Representation Using a Noise Dictionary

Abstract

Periodic noise is a well-known problem in seismic exploration, caused by power lines, pump jacks, engine operation, or other interferences. It contaminates seismic data and affects subsequent processing and interpretation. The conventional methods to attenuate periodic noise are notch filtering and some model-based methods. However, these methods either simultaneously attenuate noise and seismic events around the same frequencies, or need expensive computation time. In this work, a new method is proposed to attenuate periodic noise based on sparse representation. We use a noise dictionary to sparsely represent periodic noise. The noise dictionary is constructed based on ambient noise. An advantage of our method is that it can automatically suppress monochromatic periodic noise, multitoned periodic noise and even periodic noise with complex waveforms without pre-known noise frequencies. In addition, the method does not result in any notches in the spectrum. Synthetic and field examples demonstrate that our method can effectively subtract periodic noise from raw seismic data without damaging the useful seismic signal.

1. Introduction

Noise decreases the signal-to-noise ratio (SNR) of seismic data and affects the quality of subsequent processes [1,2]. There are six types of noise recorded using geophones [3], in which periodic noise is a kind of noise caused by power lines, pump jacks [4], engine operation [5], or other interferences, shown as monochromatic noise or multitoned noise. Sometimes, periodic noise is so strong that seismic records are severely contaminated. However, it is not easy to attenuate periodic noise, since it overlaps with seismic waves in the time domain and the frequency domain.

For periodically monochromatic noise, the conventional method of attenuation is notch filtering, which requires exact knowledge of the frequency of monochromatic noise [6]. Obviously, the notch filtering method can attenuate seismic waves at the cutoff frequency. Model-based approaches [7,8,9,10] such as sinusoidal approximation are used to remove power line noise, though this method requires accurate estimation of noise frequency and needs significant computation time [11]. Henley [12] presented a spectral clipping method to detect monochromatic noise automatically; however, it is not applicable to weak periodic noise. Karsli and Dondurur [13] used an improved mean filtering method to attenuate power line harmonic noise without noise frequency estimation; however, this requires knowledge of the rough frequency band.

In recent years, following the application of wide-azimuth, broadband, high-density seismic acquisition technology—including micro-seismic observation [14] and time-lapse monitoring based on fiber sensing—the size of 3D seismic data is increasing, and the frequency band of reflections is becoming wider. Besides monochromatic noises suppression, recognition of multitoned noises and their automatic suppression are urgently needed. Methods based on sparsity representation, such as S-transform, singular spectrum analysis and empirical mode decomposition, extend realization of automatic denoising [15,16,17]. Xu et al. [4] proposed a method based on morphological diversity of monochromatic noise and seismic waves, which assumed that the raw data are only composed of two kinds of signals, the monochromatic noise and seismic waves. It is not applicable to seismic data with strong white Gaussian noise or multitoned noise.

In this paper, a new method is proposed to attenuate periodic noise. The proposed method is based on sparse representation using a noise dictionary. The novelty is the construction of a noise dictionary, which can represent periodic noise sparsely. First, the noise dictionary is constructed using ambient noise. Next, the noise dictionary is used to sparsely represent periodic noise. Then, the de-noised data are obtained by subtracting the periodic noise from the raw seismic data. The method is applied to synthetic and field seismic data. The effectiveness of the proposed method is that it can subtract periodic noise from raw seismic data without any notches in the spectrum compared with the notch filtering method.

2. Method

Periodic noise is recurring in raw seismic data, and its amplitude is constant within the recording time [6]. Based on these characteristics, we can differentiate and estimate periodic noise from ambient noise. Because ambient noise, such as periodic noise, white Gaussian noise and other non-stationary [1,18] random noise, is ubiquitous in seismic records before the first arrivals, we can extract noise features and construct the noise dictionary from it. The algorithm is outlined in the flow skeleton (Figure 1).

2.1. Noise Period Scanning

The waveform of periodic noise is similar in a single trace. Based on the similarity, we scan the noise period from ambient noise, which exists in the records before the first arrivals. First, given the scanned period interval $\left. [T_{\min },T_{\max }\right]$, the noise period $T$ is changed from the initial scanned period $T_{\min }$ to the final scanned period $T_{\max }$. The data $S_j$ on the single trace $h_j(j=1,2,\dots ,m)$ are split into many time windows $s_i(i=1,2,\dots ,n)$, denoted as

$$S_j=[s_1,s_2,\dots ,s_n]$$

(1)

where the length of time window $s_i$ is equal to the scanned period $T$. The size of $S_j$ is $N\times 1$ ($N\ge nT$). Next, the correlation coefficients for the adjacent two time windows $s_i$ and $s_{i+1}$ are calculated:

$$\mathrm{Corr}(s_i,s_{i+1})=\frac{\mathrm{Cov}(s_i,s_{i+1})}{\sqrt{\mathrm{Var}(s_i)}\sqrt{\mathrm{Var}(s_{i+1})}}$$

(2)

where $\mathrm{Cov}\left(s_i,s_{i+1}\right)$ is the covariance of the two time windows and $\mathrm{Var}(s)$ is the variance of $s$. For accuracy, we average the correlation coefficients to obtain a coefficient to evaluate the similarity of all time windows.

$$C(T)=\frac{1}{n-1}{\displaystyle \sum _{i=1}^{n-1}\mathrm{Corr}(s_i,s_{i+1})}$$

(3)

Then, the scanned period $T$ is increased and the former procedures are repeated until all periods in the interval $\left. [T_{\min },T_{\max }\right]$ are scanned. The period which matches the maximum value of the correlation coefficients is our estimated period:

$$\tilde{T}=\underset{T}{\arg \max }\left\{C\left(T\right)\right\}$$

(4)

where $\tilde{T}$ is the period of periodic noise on the trace $h_j$.

2.2. Waveform Estimation by Stacking

The time window, whose length is $\tilde{T}$, is chosen and is represented as ${\tilde{s}}_i$. The chosen time window ${\tilde{s}}_i$ is approximated as the noise waveform; however, it is affected by white Gaussian noise. A stacking method is widely used to improve the SNR of the seismic profile [19], because it weakens white Gaussian noise and emphasizes seismic waves. We stack the approximate waveforms for an accurate noise waveform:

$$w_j={\displaystyle \sum _{i=1}^n{\tilde{s}}_i}$$

(5)

where $w_j$ is the noise waveform on the trace $h_j$ and its length $\tilde{T}$ is equal to that of ${\tilde{s}}_i$. The waveforms of different traces are similar to the near traces $h_1,h_2,\dots ,h_m$ when the traces are interfered with by the same noise source. Therefore, we can stack the similar waveforms $w_1,w_2,\dots ,w_m$ along different traces:

$$w={\displaystyle \sum _{j=1}^m{w}_j}$$

(6)

where $w$ is the noise waveform estimated from the traces $h_1,h_2,\dots ,h_m$. We note that waveforms $w_1,w_2,\dots ,w_m$ obtained by Equation (5) are not in phase. Before stacking, waveforms $w_j(j\ne 1)$ need to be cyclically shifted to the same phase as $w_1$ by scanning the shift length in an interval $\left. [0,1,\dots ,\tilde{T}-1\right]$. The suitable shift length corresponds to the maximum correlation coefficient between $w_j(j\ne 1)$ and $w_1$.

2.3. Periodic Noise Representation

Based on the stationarity of periodic noise, the basis function ${\phi }_0$ of periodic noise is obtained by extending the waveform $w$ iteratively to the length of seismic data on a single trace and then energy-normalized:

$$\phi =[w,w,\dots ,w]$$

(7)

$${\phi }_0=\frac{\phi }{{‖\phi ‖}_2}$$

(8)

where ${‖•‖}_2$ is the $l_2$ norm. The size of $\phi$ is $N\times 1$, which is equal to that of $S_j$. The noise dictionary is constructed by basis functions of different phases:

$$D=\left. [\begin{array}{ccccc}{\phi }_0& {\phi }_1& {\phi }_2& \cdots & {\phi }_{T_j-1}\end{array}\right]$$

(9)

where ${\phi }_1,{\phi }_2,\dots ,{\phi }_{\tilde{T}-1}$ are obtained by cyclically shifting $1,2,\dots ,\tilde{T}-1$ time samples. The size of $D$ is $N\times \tilde{T}$. The dictionary provides a sparse representation of periodic noise. In sparse representation theory [20], the signals can be efficiently explained as linear combinations of prespecified basis functions, where the linear coefficients are sparse. Based on sparse representation theory, the mathematical model of our noise representation is

$$S_j=Dx$$

(10)

where $S_j(j=1,2,\dots ,m)$ is the actual periodic noise and $x$ is its coefficient represented by the dictionary $D$. The size of $x$ is $\tilde{T}\times 1$. An approach to solve Equation (10) is using the sparsity constraint via an $L_0$ regularization term. Then, the periodic noise of multi-trace seismic data is obtained by solving the following optimization problem:

$${\tilde{S}}_j=\underset{S_j}{\arg \min }{‖S_j-Dx‖}_2^2 s.t.{‖x‖}_0\le 1$$

(11)

where ${\tilde{S}}_j(j=1,2,\dots ,m)$ is the approximate periodic noise and $x$ is its coefficient represented by the dictionary $D$. In the optimization problem, the sparsity of this representation is 1. Because only one basis function corresponds to a single trace, the coefficient $x$ has only one non-zero component. This is the reason for the condition on the right in Equation (11) ${‖x‖}_0\le 1$. Equation (11) is solved by the matching pursuit algorithm, which entails computing the inner products between the residual and the dictionary elements, updating the coefficient, and updating the residual iteratively [20,21]. Finally, the de-noised data are obtained by subtracting the periodic noise from the raw seismic data.

3. Synthetic Example

A synthetic dataset consisting of 21 traces is shown in Figure 2 with a 1 ms sampling rate and a 10 m geophone interval. Multitoned noises of 40 Hz and 50 Hz are incorporated to the noise-free data. In addition, a small amount of white Gaussian noise is added. Figure 3 shows the contaminated seismic dataset and its frequency–wavenumber (f-k) spectrum. The reflections are hardly identifiable owing to the strong periodic noise. The contaminations in the seismic data appear as two horizontal bands at 40 and 50 Hz in the f-k domain. The signal on the 11th trace (Figure 4a) shows the influence of multitoned noise clearly. Its amplitude spectrum (Figure 4b) also shows high spikes at frequencies of 40 Hz and 50 Hz.

The first 0.4 s of data is dominated by periodic noise marked by a red rectangle in Figure 3a. In addition, it does not consist of seismic waves. Therefore, it is chosen to construct the noise dictionary. The waveform estimated from the chosen data is shown in Figure 5. After sparse representation, the result of noise attenuation is shown in Figure 6. The main noise is attenuated and reflections can be clearly seen in Figure 6a. The eliminated noise is periodic noise as shown in Figure 6b. Their spectra show that the multitoned noise is attenuated and no seismic waves are eliminated. The de-noising result on the 11th trace and the corresponding amplitude spectra are shown in Figure 7. The de-noised seismic signal is close to the theoretical signal in both time sequence and amplitude spectrum, except for weak white Gaussian noise.

A synthetic signal is shown in Figure 8a, where the main noise is stationary and the waveform is recurring except for weak white Gaussian noise. The eliminated noise and de-noised signal are shown in Figure 8b,c, respectively. Their amplitude spectra are shown in Figure 9. Both the time sequence and amplitude spectrum of the de-noised signal are consistent with the theoretical signal. Because the amplitude spectrum is complex, it is hard to estimate its amplitudes, frequencies and phases, which are important for model-based approaches [8,9]. Therefore, the proposed method can effectively attenuate the multitoned noise and even periodic noise with a complex waveform but not influence the seismic events.

4. Field Example

The proposed method was tested using field data. The field data are a land seismic shot recorded with a 2 ms sample rate and a 40 m geophone interval, as shown in Figure 10a. The data contain periodic noise which may be caused by power lines, engine operation and other interferences. The noise is so strong that it affects the quality of the subsequent processes, especially for the first 12 traces. As shown in Figure 10c, we cannot see any waves after the time 1.5 s. There are three horizontal bands around the frequencies of 7 Hz, 14 Hz and 50 Hz in the f-k domain (Figure 10b). The periodic noise at the frequency of approximately 14 Hz, highlighted in Figure 10b, is the weakest, and the periodic noise at the frequency of 50 Hz is the strongest.

We use the ambient noise marked by a red rectangle in Figure 10a to construct the noise dictionary. The estimated waveform is shown in Figure 11. The result of noise attenuation by the proposed method is shown in Figure 12. The periodic noise is attenuated to a large degree. The horizontal bands of the f-k spectrum are almost eliminated, except for weak residual noise of about 7 Hz frequency. In addition, there are no spectral notches in the f-k spectrum (Figure 12b).

The notch filtering method is applied to the field data for a comparison. The narrow stop band of the notch filter is the frequency range $[f_0-\Delta f,f_0+\Delta f]$, where $f_0$ is the noise frequency and $2\times \Delta f$ is the noise bandwidth. The values of $f_0$ for the three noise bands are set to 7 Hz, 14 Hz and 50 Hz, respectively. The noise bandwidth is estimated to be 2 Hz and $\Delta f=1\mathrm{Hz}$ is set. The f-k spectra of the filtering result are shown in Figure 13. The horizontal bands around the frequencies of 7 Hz, 14 Hz and 50 Hz are separated from the seismic data. However, the seismic waves are also eliminated at those frequencies. This causes spectral notches of seismic waves (Figure 13a). To further show the effectiveness of our method, the spectra of the single-trace de-noised signals by the two methods are compared in Figure 14. For the notch filtering method, the spectral notch causes amplitude loss of seismic events around the frequencies of 7 Hz, 14 Hz and 50 Hz. However, amplitude loss does not occur using our proposed method. The eliminated noise (Figure 15a) is stationary. This is consistent with the characteristics of periodic noise caused by the power lines or engine operation. For comparison, the noise eliminated by the notch filtering method is not stationary because it contains seismic signals.

5. Discussions

The proposed de-noising method is based on sparse representation of periodic noise. The key to our method is the construction of the noise dictionary. Because ambient noise contains no seismic waves but dominant periodic noise and other random noise, periodic noise can be estimated without the influence of seismic reflections. Therefore, our method is useful regardless of whether the seismic waves are strong or weak. A scanning method is used to estimate the noise period. The accuracy of the noise period largely influences our de-noising result. To obtain an accurate noise waveform, the waveforms in the time domain and the space domain are stacked shown as the Equations (5) and (6), respectively. It must be emphasized that our de-noising method is only applicable to stationary noise with a constant period, waveform and amplitude. It can be used to attenuate power line harmonic noise, pump jack noise and engine operation noise in land or oceanic seismic exploration.

Based on the proposed method, ambient noise detection is urged for noise estimation. However, ambient noise has not specifically been detected in oil exploration. Therefore, we suggest that ambient noise should be acquired for one second before source excitation. Long-term ambient noise will be helpful for building a perfect noise dictionary.

Wind- or water-induced noise is not strictly stationary, and the stationarity becomes weak with increasing recording time [22,23]. Our method cannot be used to attenuate this kind of noise. To broaden the method for non-stationary noise, higher-order statistics [24] need to be considered further. In addition, methods of noise feature extraction from ambient noise based on machine learning [25] will draw increased attention in the future.

6. Conclusions

A new method is proposed to attenuate periodic noise based on sparse representation. The novelty is the construction of a noise dictionary based on ambient noise. Our method can attenuate monochromatic or multitoned periodic noise automatically without pre-known noise frequencies. The noise is assumed to be stationary noise with a constant period, waveform and amplitude. Synthetic and field tests show the effectiveness of the proposed method. Compared with the conventional notch filtering method, the proposed method can obtain de-noised data with no distortion in the time and frequency domains. Therefore, our method can attenuate periodic noise without damaging the seismic events.