## 1 Introduction

The reflected events contained in the seismic data play an important role in the study of geological structures. However, the noise contamination corrupts the quality of the reflected events and disturbs the identification of geological information. Therefore, the noise attenuation is crucial to seismic data analysis (Klemperer and Brown, 1985; Zhang and Klemperer, 2005; Zhang and Ulrych, 2003). Usually, random noise can be generated during data acquisition by various sources and this noise is unpredictable in space and time. So, it is difficult to remove random noise from the seismic data. In recent years, many efforts have been made for seismic random noise attenuation (e.g., Abma and Claerbout, 1995; Cao and Chen, 2005; Jones and Levy, 2006; Wang, 1999, 2002). However, most of the existing denoising methods do not work well in the case of low signal-to-noise ratios (SNRs). It performs poorly in signal preservation, which leads to a decline of the fidelity in seismic data processing. Thus, one of the tasks is to design a denoising method to balance the noise attenuation and signal preservation, especially when the SNR is low.

TFPF is a one-dimensional (1-D) time-frequency filtering algorithm. It can give an unbiased estimation where the signal is linear in time and embedded in white Gaussian noise (Barkat and Boashash, 1999; Boashash and Mesbah, 2004; Zahir and Hussain, 2002). It has been successfully applied to the seismic data processing field (Li et al., 2009; Lin et al., 2007, 2013; Zhang et al., 2013). This method enhances non-stationary signals from random noise through frequency modulation and instantaneous frequency estimation by taking the peak in the time-frequency distribution. In general, we use pseudo Wigner–Ville distribution (PWVD) to realize local linearity. The unbiased estimation condition of the TFPF can be better satisfied if the signal has a high degree of linearity. The nonlinearity of signals will result in amplitude loss in the TFPF. So how to increase the linearity is crucial to the TFPF. One of the effective solutions is to filter along the reflected events. It also takes the coherence between the adjacent channels into consideration, which has a positive effect on the continuity of reflected events. With this in mind, an improved TFPF algorithm has been developed and filtering along a radial trace with an angle to the channel direction (Wu et al., 2011). However, a fixed radial filtering trace cannot fit with the curve reflected events effectively, which leads to some limitations in actual application. A similar trace to the reflected events for filtering is important (Tian and Li, 2014), but the field seismic data is too complicated to find a suitable trace for all the reflected events.

This paper introduces a Multiple Directional Time-Frequency Peak Filtering technique. To do the TFPF along the reflected events, we decompose the 2-D matrix of seismic data into multiple direction components. The direction decomposition is equivalent to the use of some linear segments to approximate the curve reflected events. The distributions of the events in each direction component show as the straight line and have the same slope; the slope direction is the direction of decomposition. So, TFPF could be done along a group of parallel radial filtering traces with the same slope as the decomposed reflected events in each direction component. The sum of the filtering results from all the given directions as the final processed result. The direction decomposition is realized using the directional filter bank. Our method makes the filtering along the reflected events possible by a novel idea without finding the accurate traces of all the reflected events.

## 2 Time-frequency peak filtering

Let the noisy signal $s\left(t\right)$ be modeled by the equation:

$$s\left(t\right)=x\left(t\right)+n\left(t\right)\text{,}$$ | (1) |

First, we encode the noisy signal $s\left(t\right)$ as the instantaneous frequency of the analytic signal ${z}_{\text{s}}\left(t\right)$ via frequency modulation; ${z}_{\text{s}}\left(t\right)$ can be expressed as:

$${z}_{\text{s}}\left(t\right)={\text{e}}^{\text{j}2\pi \mu {\int}_{0}^{t}s\left(\lambda \right)\text{d}\lambda}\text{,}$$ | (2) |

Then we calculate the Wigner–Ville distribution (WVD) of the analytic signal ${z}_{\text{s}}\left(t\right)$ through

$$WV{D}_{{z}_{\text{s}}}\left(t,f\right)={\int}_{-\infty}^{\infty}{z}_{s}\left(t+\tau /2\right){{z}_{s}}^{\ast}\left(t-\tau /2\right){\text{e}}^{\u2013\text{j}2\pi f\tau}\text{d}\tau \text{,}$$ | (3) |

Finally, we estimate the peak of $WV{D}_{{z}_{s}}\left(t,f\right)$ to obtain the instantaneous frequency of z_{s}(t). Since we encode the noisy signal s(t) as the instantaneous frequency of z_{s}(t) by Eq. (2), the estimated instantaneous frequency is the estimation of the signal x(t), which is expressed as:

$$\stackrel{\u02c6}{x}\left(t\right)=\underset{f}{\mathrm{arg}\text{\hspace{0.28em}}\text{max}}\left[WV{D}_{{z}_{s}}\left(t,f\right)\right]/\mu \text{.}$$ | (4) |

The TFPF method is unbiased for a linear signal $x\left(t\right)=\alpha t+C$ ($\alpha $ and C are constants) from white Gaussian noise. The bias B(t) is defined as $B\left(t\right)=E\left[\stackrel{\u02c6}{x}\left(t\right)-x\left(t\right)\right]$, E denotes the mathematical expectation. If a signal is embedded in a stationary white Gaussian noise background, the final derivation result of the bias $B\left(t\right)$ of the TFPF algorithm is written as (Boashash and Mesbah, 2004):

$$B\left(t\right)=\underset{f}{\mathrm{arg}\text{\hspace{0.28em}}max}(WV{D}_{{z}_{x}}\left(t,f\right)\underset{f}{\times}{\left(4{\pi}^{2}{k}_{{\text{n}}_{2}}{\mu}^{2}/((2{\pi}^{2}{k}_{{\text{n}}_{2}}{\mu}^{2}\right)}^{2}+{\left(2\pi f\right)}^{2})))-x\left(t\right)\text{,}$$ | (5) |

$$\begin{array}{l}B\left(t\right)=\underset{f}{\mathrm{arg}\text{\hspace{0.28em}}max}\left(\delta \left(f-x\left(t\right)\right)\underset{f}{*}\left(4{\pi}^{2}{k}_{{n}_{2}}{\mu}^{2}/\left({\left(2{\pi}^{2}{k}_{{n}_{2}}{\mu}^{2}\right)}^{2}+{\left(2\pi f\right)}^{2}\right)\right)\right)-x\left(t\right)\\ \text{\hspace{0.28em}}\text{\hspace{0.28em}}\text{\hspace{0.28em}}\text{\hspace{1em}}=\underset{f}{\mathrm{arg}\text{\hspace{0.28em}}max}\left(4{\pi}^{2}{k}_{{n}_{2}}{\mu}^{2}/\left({\left(2{\pi}^{2}{k}_{{n}_{2}}{\mu}^{2}\right)}^{2}+{\left(2\pi \left(f-x\left(t\right)\right)\right)}^{2}\right)\right)-x\left(t\right)=0.\end{array}$$ | (6) |

However, the reflected seismic signals are nonlinear. So, there must be a bias in the filtering processing by the TFPF. To reduce the bias, TFPF adopts the PWVD, which is the windowed version of the WVD to make the signal approximately linear within the window. The PWVD is defined as (Barbarossa, 1997):

$$PWV{D}_{z}\left(t,f\right)={\int}_{-\infty}^{\infty}h\left(\tau \right)z\left(t+\tau /2\right){z}^{\ast}\left(t-\tau /2\right){\text{e}}^{-\text{j}2\pi f\tau}\text{d}\tau \text{,}$$ | (7) |

The PWVDs of an analytic Ricker wavelet contaminated by an additive white Gaussian noise (SNR = 0 dB) are shown on Fig. 1. On this illustration, we can see that the peak of the PWVD is mainly located at the instantaneous frequency of the analytic signal. So, we can use the peak of the PWVD to recover the reflected signal.

## 3 Principle of multiple directional TFPF

To obtain an unbiased result, one effective way is to improve the degree of linearity of the desired signal, which is an unbiased estimation condition of TFPF. The solution can be implemented by filtering along the reflected events. Fig. 2b shows five representative radial traces corresponding to the numbered trajectories on Fig. 2a (Wu et al., 2011). We can see that a common-origin seismic reflected event can be mapped into a comparatively higher degree of linearity when the slope of the trajectory is close to that of the event. So, the radial traces such as direction 1 or 2 are suitable for the unbiased TFPF because the signals extracted from these radial traces have a much higher degree of linearity. Thus, constructing the reasonable filtering traces of events is the key of TFPF method to reduce the bias. However, the traces suitable for all the reflected events in the field data are too difficult to construct because of complexity and noise interference. However, we may change our idea to implement trace construction for all the reflected events by decomposing the seismic reflected events into some linear segments. In other words, we use some linear segments to approximate the events. Then, TFPF can be done along multiple linear segments that are the local direction of events reflected by a group of parallel radial filtering traces. Considering that the linear segmentation is difficult to realize in the time domain, we can perform the decomposition in the frequency domain. The spectrum of the noisy synthetic seismic record (Fig. 3a) in the frequency-wavenumber domain is shown on Fig. 3b. We divide the region with energy distribution into some directions. The direction division in the frequency domain is equivalent to approximating the events by some linear segments in the time domain. The more directions we divide in frequency domain, the higher the linearity of local reflected events in each direction component is. In other words, if the curvature of reflected events is larger, more linear segments are used to approximate the events so that we can do filtering along the local reflected events by a group of parallel radial filtering traces.

These direction components are obtained by a 2-D directional filter bank. This directional filter bank can decompose a 2-D matrix of seismic data according to directions while achieving perfect reconstruction (Do and Vetterli, 2005; Feilner et al., 2005). The directional filter bank is efficiently implemented via an l-level binary tree decomposition that leads to ${2}^{l}$ sub-bands with wedge-shaped frequency partitioning. The decomposition of a noisy seismic data $s(n)$ of time-offset domain $\left(n=\left(d,t\right),d\text{\hspace{0.28em}}\text{is}\text{\hspace{0.28em}}\text{offset}\right)$ by the directional filter bank can be written as:

$$s\left(n\right)=\sum _{k}{s}_{k}\left(n\right)=\sum _{k=0}^{{2}^{l}-1}\sum _{m\in {Z}^{2}}{y}_{k}\left(m\right){g}_{k}\left(n-{S}_{k}m\right)\text{,}$$ | (8) |

$${y}_{k}\left(m\right)\le s\left(n\right),{h}_{k}\left({S}_{k}m-n\right)>\text{,}$$ | (9) |

$${S}_{k}=\left\{\begin{array}{l}diag\left({2}^{l-1},2\right)for0\le k<{2}^{l-1},\hfill \\ diag\left(2,{2}^{l-1}\right)for{2}^{l-1}\le k<{2}^{l},\hfill \end{array}\right.$$ | (10) |

${s}_{k}\left(n\right)$ is a decomposed direction component of $s\left(n\right)$ and $k$ the direction index. ${h}_{k}\left(n\right)$ is the impulse response of analysis filter and ${g}_{k}\left(n\right)$ is the impulse response of synthesis filters. $m$ is used to denote location of noisy seismic data.

Next we give the whole filter process of our multiple directional TFPF. First, we decompose a noisy seismic data $s\left(d,t\right)$ into some direction components by Eq. (8). Each decomposed component ${s}_{k}\left(d,t\right)$ has its own orientation depending on the index $k$. We take the noisy synthetic seismic record (Fig. 3a) as an example. The decomposed direction components of Fig. 3a are shown on Fig. 4.

Then the data in each component need to be resampled along this known decomposed direction. The processing can be realized by transforming a 2-D matrix of seismic data into a new data matrix along some parallel traces in this decomposed direction. The corresponding multiple directional filtering traces are shown on Fig. 5. The equation of the parallel trace that we used is expressed as

$${t}_{i}={a}_{k}{d}_{i}+b\text{.}$$ | (11) |

_{i}according to Eq. (11). After resampling, the new data matrix in each direction component is filtered by the TFPF method.

The parameter ${a}_{k}$ denotes the direction of filtering method. Our method applies the TFPF along different ${a}_{k}-\text{directions,}$ which are also the decomposition directions of the reflected events. The final result is obtained by putting the filtering results from all the directions of reflected events together. The effect of the whole processing is equivalent to the filtering along seismic reflected events and avoids finding the directions of all the events.

## 4 Implementation and applications

### 4.1 Synthetic data example

To verify the effectiveness of the proposed MD–TFPF, we gave it to process synthetic seismic data. The synthetic example is a seismic model with three hyperbolic reflected events as shown on Fig. 6a. We add a white Gaussian noise to the synthetic model. Fig. 6b is a noisy version (SNR = −5 dB) of the noise-free Fig. 6a. In this example, MD–TFPF is compared with conventional TFPF. Fig. 6c and d show the denoising results by the two methods. As can be seen from Fig. 6c and d, we can achieve the visual qualitative performance. Both two methods can suppress some random noise and the three events can be identified. While the MD–TFPF achieves clean background and the reflected events are continuous and recovered, especially in the rectangular region.

Detailed effects of the two processing consequents can be obtained by analyzing the single trace and its amplitude spectrum. We extract the 36th trace of the synthetic records and redisplay it as an example. The amplification waveforms in time and frequency domains are shown on Fig. 6e and f. Fig. 6e exhibits the denoising result of a seismic wavelet from the 36th trace of the synthetic record in the time domain. Notice that the filtering result of MD–TFPF is close to the original signal and the amplitude value of noise on both sides of the signal is significantly lower than conventional TFPF does. Through the calculation, the correlation coefficient between the original signal and the denoised signal by using MD–TFPF is up to 0.83, but using a conventional TFPF method, the correlation coefficient is only 0.54. Furthermore, we also compare the filtering performance in the frequency domain, which is shown on Fig. 6f. It can also be seen that the MD–TFPF provides the best-matching amplitude spectrum to the noise-free spectrum for the signal components in the 0–70 Hz band. The spectrum of the convention TFPF has larger high-frequency components after 70 Hz than our method. So the random noise is effectively suppressed.

In the following, we ran extensive tests on data with white Gaussian noise at different levels to show the effectiveness of the MD–TFPF on SNR (in dB) and mean square error (MSE).

Table 1 gives the resulting SNR and MSE for synthetic seismic data of Fig. 6a. The better ﬁltering effect has larger SNR value and lower MSE value. As can be seen in Table 1, the SNRs after filtering by the MD–TFPF are larger and the MSEs are smaller than when applying conventional TFPF in any SNR. This advantage of the MD–TFPF is obvious when the SNR is below −8 dB. The result implies that MD–TFPF provides accurate estimation and improvement to noisy signals.

**Table 1**

SNR and MSE after applying conventional and MD–TFPF with different noise levels.

SNR (in dB) computed from the original signal | Processed signal by conventional TFPF | Processed signal by MD–TFPF | ||

SNR (dB) | MSE | SNR (in dB) | MSE | |

−20 | −13.21 | 0.8652 | −10.56 | 0.4283 |

−16 | −9.83 | 0.6831 | −5.33 | 0.2946 |

−12 | −6.1 | 0.4314 | −1.46 | 0.1887 |

−8 | −1.16 | 0.3427 | 2.78 | 0.1153 |

−4 | 3.24 | 0.2135 | 5.92 | 0.0832 |

0 | 7.86 | 0.1682 | 9.32 | 0.0667 |

4 | 9.17 | 0.1026 | 12.29 | 0.0465 |

8 | 12.93 | 0.0868 | 14.45 | 0.0329 |

Our method has another advantage, which is that it can recover signals from non-stationary random noise with very low SNR. We add a non-stationary random noise with high noise levels to the synthetic model (Fig. 6a). The noisy record with the SNR −15 dB is shown on Fig. 7a, where the three hyperbolic reflection events are almost submerged completely in the strong noise. In this example, the MD–TFPF is compared with conventional TFPF and the curvelet denoising method. Fig. 7b, c and d show the denoising results by curvelet denoising method, conventional TFPF and our method, respectively. Notice that the events on Fig. 7b and c can hardly be identified while they are well recovered on Fig. 7d. Thus, the MD–TFPF is also very effective in the low SNR.

### 4.2 Field data application

A common shot point seismic data that we test is shown on Fig. 8a. The receiver interval is 30 m and the sampling frequency is 1000 Hz. The distance between the source and the receivers changes from 600 to 3600 m. The continuities of seismic reflected events are badly damaged and the events are obscured because of the interference from random noise. To show the effectiveness of our method, we compare the results of MD–TFPF with conventional TFPF and F–X deconvolution methods. It can be found that F–X deconvolution (Fig. 8b) and conventional TFPF methods (Fig. 8c) cannot reduce random noise as effectively as MD–TFPF, whose results are displayed on Fig. 8d, particularly in the ellipse area. To demonstrate the visual quality of the events, we show an enlarged view between traces 60–110 and double-time 500–1200 ms (Fig. 9) in the rectangle on Fig. 8. In contrast, the reflected events processed by MD–TFPF (Fig. 9d) are clearly recovered and the events continuity is enhanced (e.g., in the rectangles). According to the experimental results, the proposed MD–TFPF method has a great advantage in eliminating random noise and efficiently protecting the reflected events compared with the other two methods.

We also apply our method to a real seismic section to test the applicability. Fig. 10a shows a section with high background random noise. This background random noise overwhelms weak reflected events and disrupts the continuity of events over the entire section. We compare the MD–TFPF with conventional TFPF and the curvelet denoising method. As seen on Fig. 10b and c, some random noise still exists after filtering by curvelet denoising method and conventional TFPF method, especially in the rectangular region. However, the application of MD–TFPF method shows that most of the background noise is suppressed without the edges of the reflected events being blurred and the method does not make any significant changes to the shapes of reflected events. Furthermore, the weak reflected events continuity is also enhanced. This method greatly improved SNR, which is clearly identified on the seismic section (Fig. 10d). The result on Fig. 10d shows that even with strong random noise, the proposed method is still effective in suppressing noise.

## 5 Conclusions

In this paper, a multiple directional time-frequency peak filtering has been proposed to improve the filtering result in the conventional time-frequency peak filtering. It filters the seismic data along multiple decomposed directions of the reflected events instead of the channel direction. The coherence between adjacent channels is fully considered and utilized, which recovers the reflected events more clearly and continuously. The direction decomposition is realized based on a directional filter bank. In this novel method, the directions of reflected events are indirectly obtained by a directional filter bank without accurately tracing all the events, which makes filtering along the reflected events in the complex field seismic data possible. Therefore, our method has a high practicability in seismic data denoising. Texts on both synthetic and field seismic data have demonstrated that the proposed method can preserve signals and effectively enhance the continuity of reflected events. Furthermore, the proposed method makes it possible to recover the desired seismic signal in the case of SNR down to −15 dB. It has great practical value, especially when the SNR of seismic data is very low.

## Acknowledgments

This research is financially supported by the National Natural Science Foundation of China (under grants 41130421 and 41274118); Science and Technology Development program of Jilin Province of China (under grant 20110362) and the Development and Reform Commission of Jilin Province (grants No. 2011006-2, JF2012C013-2).