TY - JOUR
AU - Zhang,, Dong
AB - SUMMARY Amplitude-preserving data processing is an important and challenging topic in many scientific fields. The amplitude-variation details in seismic data are especially important because the amplitude variation is directly related with the subsurface wave impedance and fluid characteristics. We propose a novel seismic noise attenuation approach that is based on local plane-wave assumption of seismic events and the amplitude preserving capability of the orthogonal polynomial transform (OPT). The OPT is a way for representing spatially correlative seismic data as a superposition of polynomial basis functions, by which the random noise is distinguished from the useful energy by the high orthogonal polynomial coefficients. The seismic energy is the most correlative along the structural direction and thus the OPT is optimally performed in a flattened gather. We introduce in detail the flattening operator for creating the flattened dimension, where the OPT can be applied subsequently. The flattening operator is created by deriving a plane-wave trace continuation relation following the plane-wave equation. We demonstrate that both plane-wave trace continuation and OPT can well preserve the strong amplitude variation existing in seismic data. In order to obtain a robust slope estimation performance in the presence of noise, a robust slope estimation approach is introduced to substitute the traditional method. A group of synthetic, pre-stack and post-stack field seismic data are used to demonstrate the potential of the proposed framework in realistic applications. Image Processing, Seismic Noise 1 INTRODUCTION Seismic noise attenuation is one of the most significant steps in the whole seismic data processing and imaging workflow. It has great influence to many subsequent processing tasks, such as amplitude-variation-offset inversion, reverse time migration, full waveform inversion and automatic interpretation for oil and gas detection (Wang et al., 2015; Gan et al., 2016; Gao et al., 2016; Huang et al., 2016; Qu et al., 2016; Xue et al., 2016b; Asgedom et al., 2017; Zeng et al., 2017; Zhang et al., 2017; Chen, 2018). In the past several decades, a large number of algorithms have been developed for seismic noise attenuation. Stacking the seismic data along the spatial directions, for example, the offset direction, can enhance the energy of spatially coherent useful waveform signals as well as mitigate the spatially incoherent random noise (Liu et al., 2009a; Xie et al., 2016; Wu & Bai, 2018c). One of the commonly used state-of-the-art algorithms is the prediction-based method, including t–x predictive filtering (Abma & Claerbout, 1995), f–x deconvolution (Canales, 1984), the polynomial fitting based approach (Liu et al., 2011) and non-stationary predictive filtering (Liu et al., 2012; Liu & Chen, 2013). This type of methods utilize the predictive property of useful signals along spatial direction to create a regression-like model for distinguishing between signal and noise. Another type of commonly used methods are based on data decomposition. This type of methods assume that noisy seismic data can be decomposed into different components where signal and noise are separated based on their frequency difference or morphological difference (Huang et al., 2017). Empirical mode decomposition (EMD; Huang et al.1998; Chen 2016) and its improved version, for example, ensemble empirical mode decomposition (EEMD; Wu & Huang 2009), complete ensemble empirical mode decomposition (CEEMD; Colominas et al.2012;), have been used intensively for reducing the noise in seismic data (Chen et al., 2016a) . Variational mode decomposition was proposed by Dragomiretskiy & Zosso (2014) for substituting EMD because of its explicit control on the decomposition performance. It has been utilized for noise attenuation in Liu et al. (2017) and for time-frequency analysis by Liu et al. (2016a). Regularized non-stationary decomposition (Yang et al., 2014; Wu et al. , 2016) is another decomposition method which is also based on a solid mathematical model. Singular value decomposition (SVD) can also be used to extract the most spatially coherent components from the multidimensional seismic data. Sparse transform based approaches assume that multidimensional seismic data can be compressed in a sparse transformed domain, where the signal is represented by high-amplitude coefficients and the noise is represented by small-amplitude coefficients (Gholami, 2013; Bai & Wu, 2018). Hence, by transforming data to the sparse domain, a soft thresholding can be applied to reject those small-amplitude coefficients that correspond to noise, which is followed by an inverse transform from the thresholded coefficients to time–space domain (Chen, 2017). This type of methods are closely connected with the compressive sensing paradigm (Lorenzi et al., 2016). Widely used sparse transforms are Fourier transform, curvelet transform (Candès et al., 2006; Herrmann et al., 2007; Herrmann & Hennenfent, 2008); Wang et al., 2011; Zu et al., 2017), seislet transform (Fomel & Liu, 2010; Chen & Fomel, 2015a; Gan et al., 2015a; Gan et al., 2015b), shearlet transform (Kong et al., 2016), Radon transform (Foster & Mosher, 1992) and a variety of sparse wavelet transforms (Mousavi & Langston, 2016b; Anvari et al., 2017), for example, synchrosqueezing (Daubechies et al., 2011; Mousavi et al., 2016; Mousavi & Langston, 2016a, 2017) or empirical wavelet transforms (Liu et al., 2016b), etc. Recently, the adaptive dictionary learning has gained a lot of attention in the seismic data processing field (Chen, 2017; Wu & Bai, 2018b). The dictionary learning based sparse representation differs from the traditional sparse transforms in that the basis functions for the sparse transform are adaptively learned from the data itself, instead of being fixed in the traditional transforms. Rank reduction methods are one of the most effective methods in the seismic data processing community, which includes the Cadzow filtering (Trickett, 2008), singular spectrum analysis (SSA) (Vautard et al., 1992; Wu & Bai, 2018a; Zhou et al., 2018), damped SSA (Chen et al., 2016b; Zhang et al., 2016a,b), and multistep SSA (Zhang et al. , 2016c). There are two least-squares projection step in the damped SSA method. The first step can be considered as a rank reduction method while the second step can be interpreted as a compensation step for the non-optimal performance of the rank-reduction method, that is, the approximated signalsubspace in the traditional rank-reduction framework is a mixture of both signal and noise subspaces. From a different aspect, Xue et al. (2016a) proposed a rank-increasing method for iteratively estimating the spike-like noise instead of estimating signals in deblending of simultaneous-source data (Zu et al. , 2016a,b; Bai & Wu, 2017; Zhou & Han, 2018; Wu & Bai, 2018d; Bai et al., 2018a,b). Mean and median filters utilize the statistical difference between signal and noise to reject the Gaussian white noise or impulsive noise (Liuet al., 2009b; Liu, 2013). In addition to these classic noise attenuation methods, some advanced denoising methods have been proposed in recent years. Time-frequency peak filtering (Kahoo & Siahkoohi , 2009; Lin et al., 2013, 2015) based approaches utilize the high-resolution property of time-frequency transform to distinguish between useful signals and random noise. Instead of developing a standalone denoising strategy, Chen & Fomel (2015b) proposed a two-step denoising approach that tries to solve along-existing problem in almost all denoising approaches: the signal leakage problem. By initiating a new concept called local orthogonalization, Chen & Fomel (2015b) successfully retrieved the coherent signals from the removed noise section to guarantee no signal leakage in any denoising algorithms. For all the aforementioned state-of-the-art noise attenuation algorithms, none of them are specifically designed for preserving the strong amplitude-variation details in seismic data. As we know, the amplitude variations in seismic data greatly affect the subsurface oil and gas exploration and production. Hence, the amplitude preserving capability is one of the backbone features we need to keep in mind when designing a new denoising algorithm. In this paper, we are solving the serious problem that is often neglected in traditional seismic data processing by proposing the plane-wave orthogonal polynomial transform (OPT) method. Here we want to clarify that the amplitude variation we mention here refers to strong amplitude variation, not simply the edge details or weak signals that are often mentioned in the literature. We first introduce the basic knowledge of the OPT, which is the key component that brings us the amplitude-preserving capability in the proposed framework. We then introduce the theory of plane-wave trace continuation that is used forflattening the seismic events without damaging the amplitude information. We show that both plane-wave trace continuation and OPT can well preserve the amplitude variation details in the seismic data, which accounts for the superb performance of preserving the amplitude details in the real data applications. Considering the strong influence of the slope estimation to the plane-wave flattening, we introduce a robust slope estimation method that can substitute the traditional plane-wave destruction (PWD; Fomel 2002) based methods in the presence of strong noise. A group of synthetic, pre-stack and post-stack field seismic data are used for demonstrating the performance of the proposed framework. 2 THEORY 2.1 Orthogonal polynomial transform In a seismic profile, the amplitude of time t and space x can be expressed as: \begin{equation*} A(t,x) = \sum _{j=0}^{M-1} C_j(t) P_j(x), \end{equation*} (1) where {Pj(x), j = 0, 1, 2, ⋅⋅⋅, M − 1} is a set of orthogonal polynomials and M is the number of basis functions and {Cj(t), j = 0, 1, 2, ⋅⋅⋅, M − 1} is a set of coefficients. The Pj(x) is a unit basis function that satisfies the condition: \begin{equation*} P_j(x)P_i(x)=\delta _{ij}, \end{equation*} (2) where δij denotes the Kronecker delta. The spectrum defined by Cj(t) denotes the energy distribution of the t–x domain data in the orthogonal polynomials transform domain. Besides, the low-order coefficients represent the effective energy and the high-order coefficients represent the random noise energy. We provide a detailed introduction about how we construct the orthogonal polynomial basis function in Appendix A. In a matrix-multiplication form, eq. (1) can be expressed as the following equation: \begin{equation*} \mathbf {A} = \mathbf {C}\mathbf {P}, \end{equation*} (3) where |$\mathbf {A}$| is constructed from A(t, x), |$\mathbf {C}$| is constructed from Cj(t) and |$\mathbf {P}$| is constructed from Pj(x). |$\mathbf {A}$| is known and |$\mathbf {P}$| can be constructed using the way introduced in Appendix A. The unknown is |$\mathbf {C}$|⁠. |$\mathbf {C}$| can be obtained by inverting the eq. (3), \begin{equation*} \mathbf {C}=\mathbf {A}\mathbf {P}^H(\mathbf {P}\mathbf {P}^H)^{-1}, \end{equation*} (4) where [ · ]H denotes matrix transpose. In this paper, we choose M = 20, which indicates that we select 20 orthogonal polynomial basis function to represent the seismic data. Hence, inverting equation |$\mathbf {P}\mathbf {P}^H$| is simply inverting a 20 × 20 matrix and is computationally efficient. In the OPT method, we need to define the order of coefficients we want to preserve, the process of which corresponds to applying a mask operator to the orthogonal polynomial coefficients. Mask operator can be chosen to preserve low-order coefficients and reject high-order coefficients. It takes the following form: \begin{equation*} \mathcal {M}_{\tau }(C_j(t)) = \left\lbrace \begin{array}{ll}C_j(t) & \text{for}\quad j \le \tau \\ 0 & \text{for}\quad j \gt \tau \end{array}\right., \end{equation*} (5) where |$\mathcal {M}$| denotes the mask operator, Cj(t) denotes the orthogonal polynomial coefficients at time t and order j. τ denotes a threshold for coefficients rejection. In this paper, we chose τ = 2. The coefficients after applying the mask operator 5 become \begin{equation*} \hat{\mathbf {C}}=\mathcal {M}_{\tau }(\mathbf {C}). \end{equation*} (6) The useful signals can be reconstructed by \begin{equation*} \hat{\mathbf {A}} = \hat{\mathbf {C}}\mathbf {P}, \end{equation*} (7) where |$\hat{\mathbf {A}}$| denotes the denoised data. 2.2 Plane-wave trace continuation In this section, we will derive a plane-wave flattening operator so that the seismic data can be flattened locally and the OPT can then be applied to the flattened gather. The key question here is how to map the curved events into flattened events. We do the data mapping by a recursively predicting strategy. Each trace in a seismic gather can be predicted using neighbour traces. Given a reference trace, each trace in the gather can predict the reference trace in some ways, for example, by recursive trace continuation. Ideally, arranging the predicted reference traces (from all other traces) into a gather constructs a flattened gather. Next, we will introduce the theory of how we predict traces following the plane-wave equation, which we call plane-wave trace continuation. For simplicity, we always treat the first trace in the gather as the reference trace. To flatten the gather, we need to predict the first trace from all other traces and arrange them together. In a brief mathematical way, predicting the first trace from the jth (j ≠ 1) trace can be expressed \begin{equation*} \mathbf {d}_1=\mathbf {P}_{2,1}\mathbf {P}_{3,2}\cdots \mathbf {P}_{j,j-1}\mathbf {d}_{j} , \end{equation*} (8) where |$\mathbf {P}_{p,q}$| denotes a prediction operator to predict trace |$\mathbf {d}_q$| from trace |$\mathbf {d}_p$|⁠. Specifically, |$\mathbf {P}_{p,p-1}$| denotes the prediction between two traces from right to left. In the inverse process, the jth trace can be predicted in a similar recursive formula from left to right: \begin{equation*} \mathbf {d}_j=\mathbf {P}_{j-1,j}\cdots \mathbf {P}_{2,3}\mathbf {P}_{1,2}\mathbf {d}_{1}. \end{equation*} (9) Predicting from the jth trace to (j + 1)th trace (or from the (j + 1)th trace to jth trace) requires solving the plane-wave equation: \begin{equation*} \frac{\partial u}{\partial x} + \sigma \frac{\partial u}{\partial t} = 0, \end{equation*} (10) where u(t, x) is the seismic record and σ is local slope. In the case of the constant local slope, eq. (10) has the following solution: \begin{equation*} u(t,x) = f(t-\sigma ), \end{equation*} (11) where f is the waveform function. In the variable-slope case, we can solve eq. (10) by discretizing it. Let |$u_p^v$| denote u(vΔt, pΔx), and then we obtain: \begin{equation*} \begin{split} &\frac{u_{p+1}^{v+1} - u_p^{v+1}}{2\Delta x} + \frac{u_{p+1}^{v} - u_p^{v}}{2\Delta x}\\ &+\sigma ^v_p \left(\frac{u_{p+1}^{v+1}-u_{p+1}^v}{2\Delta t}+ \frac{u_{p}^{v+1}-u_{p}^v}{2\Delta t}\right) = 0. \end{split} \end{equation*} (12) Rearranging the terms in eq. (12), we get \begin{equation*} \begin{split} &\left(\frac{1}{\Delta x}+\frac{\sigma _p^v}{\Delta t}\right) u_{p+1}^{v+1} + \left(-\frac{1}{\Delta x}+\frac{\sigma _p^v}{\Delta t}\right)u_{p}^{v+1} +\\ &\left(\frac{1}{\Delta x}-\frac{\sigma _p^v}{\Delta t}\right)u_{p+1}^{v}+ \left(-\frac{1}{\Delta x}-\frac{\sigma _p^v}{\Delta t}\right)u_{p}^{v}=0. \end{split} \end{equation*} (13) Then we have the following point-to-point recursion from |$u_{p}^{v}$| to |$u_{p+1}^{v+1}$|⁠: \begin{equation*} \begin{split} &u_{p+1}^{v+1}= \left(\frac{1}{\Delta x}+\frac{\sigma _p^v}{\Delta t}\right)^{-1}\left[\left(\frac{1}{\Delta x}-\frac{\sigma _p^v}{\Delta t}\right)u_{p}^{v+1} +\right.\\ &\left.\left(-\frac{1}{\Delta x}+\frac{\sigma _p^v}{\Delta t}\right)u_{p+1}^{v}+ \left(\frac{1}{\Delta x}+\frac{\sigma _p^v}{\Delta t}\right)u_{p}^{v}\right]. \end{split} \end{equation*} (14) Eq. (14) can be used for the continuation from the first trace to all other traces, in a similar way, all other traces can predict the first trace by an inverse trace continuation process. Since the trace continuation process is achieved using the local plane-wave assumption of seismic data, we call the relation in eq. (14) the plane-wave trace continuation relation. Fig. 1 shows an example of the trace prediction process. We start from the first trace in a gather, and predict each trace in the gather from the first trace, following a given local slope field. Fig. 1(a) shows the initial status of the trace prediction process, where only the first trace is shown. Following the slope field shown in Fig. 1(b), we can predict a complete gather from the first trace using the plane-wave equation. The complete gather is shown in Fig. 1(c). It is clear that the morphology of the predicted gather is consistent with the local slope shown in Fig. 1(b). Fig. 1(d) shows a flattened gather from the curved events shown in Fig. 1(c). We flatten the events by predicting the first trace from each trace shown in Fig. 1(c). For a clear view of the first trace, we plot it in Fig. 1(e). Figure 1. Open in new tabDownload slide Trace prediction for clean data. (a) Reference trace. (b) Slope field. (c) Predicted gather from the reference trace. (d) Flattened gather by predicting the reference trace from each trace in Fig. (c). (e) First trace in Fig. (a). Figure 1. Open in new tabDownload slide Trace prediction for clean data. (a) Reference trace. (b) Slope field. (c) Predicted gather from the reference trace. (d) Flattened gather by predicting the reference trace from each trace in Fig. (c). (e) First trace in Fig. (a). We then show an example in the presence of random noise. The noisy data is simulated with |$\mathbf {d}=\mathbf {s}+\mathbf {n}$|⁠, where |$\mathbf {s}$| is signal, that is, the solution to a wave equation in some random medium. The signal has a certain mean and variance, and a certain spatial and temporal correlation structure. |$\mathbf {n}$| is the noise, for which, presumably the expectation 〈n〉 is zero. The overall spatial variance of the noise is a certain number, and its covariance with the signal is zero. The noise is distributed in the 2-D plane following a Gaussian rule and has no spatial correlation. Fig. 2(a) shows the noisy reference trace. The details of the noisy trace are shown in Fig. 2(e). Following a given slope field shown in 2(b), we predict a complete gather from the first trace, and show the gather in Fig. 2(c). It can be seen that the noise is preserved during the prediction process. Fig. 2(d) shows the flattened gather from Fig. 2(c). We can conclude from this test that trace prediction can preserve any component in the starting trace. Figure 2. Open in new tabDownload slide Trace prediction for noisy data. (a) Reference trace. (b) Slope field. (c) Predicted gather from the reference trace. (d) Flattened gather by predicting the reference trace from each trace in Fig. (c). (e) First trace in Fig. (a). Note that during trace prediction, the noise is preserved as coherent signal. Figure 2. Open in new tabDownload slide Trace prediction for noisy data. (a) Reference trace. (b) Slope field. (c) Predicted gather from the reference trace. (d) Flattened gather by predicting the reference trace from each trace in Fig. (c). (e) First trace in Fig. (a). Note that during trace prediction, the noise is preserved as coherent signal. We then use another example to show the amplitude preserving feature of the flattening operator. Fig. 3(a) shows a complete gather containing three curved events and amplitude variation along the spatial direction. Fig. 3(b) shows the flattened gather from the curved events, where we can see that the amplitude variation is well preserved during the trace prediction process. Figure 3. Open in new tabDownload slide (a) Curved events. (b) Flattened events by predicting the first trace from each trace in (a). Note that during trace prediction, the amplitude is well preserved. Figure 3. Open in new tabDownload slide (a) Curved events. (b) Flattened events by predicting the first trace from each trace in (a). Note that during trace prediction, the amplitude is well preserved. 2.3 Robust slope estimation Another important factor in plane-wave OPT is the local slope calculation. The accuracy of the slope estimation affects performance of the flattening operation and the following OPT. In this part, we will introduce a robust slope estimation method that is based on the Hilbert transform (Liu et al., 2015). Rearranging eq. (10) we get \begin{equation*} \sigma = -\frac{\frac{\partial u}{\partial x}}{\frac{\partial u}{\partial t}}. \end{equation*} (15) Eq. (15) can be further derived such that \begin{equation*} \begin{split} \sigma = -\frac{\frac{\partial u}{\partial x}}{\frac{\partial u}{\partial z}}=-\frac{F_x^{-1} [ H_{\mathrm{D}\mathrm{X}}[ F_x u ] ]}{F_t^{-1} [H_{\mathrm{D}\mathrm{T} }[ F_t u ] ]}, \end{split} \end{equation*} (16) where HDX is frequency response function of the partial derivative in the x direction, and HDT is frequency response function of the partial derivative in the t direction. Fx and Ft denote the Fourier transform along the x and t directions, respectively. It can be straightforwardly derived that \begin{equation*} \sigma =-\frac{\mathcal {H}_x(u)}{\mathcal {H}_t(u)}, \end{equation*} (17) where |$\mathcal {H}_x(u)$| denotes the Hilbert transform of u along x direction and |$\mathcal {H}_t(u)$| denotes the Hilbert transform of u along t direction. Fig. 4 shows a slope calculation test. We calculate the slope from the noisy data using the traditional PWD method and the robust slope calculation method, respectively. As a comparison, an accurate slope estimation from the clean data using the PWD algorithm is used to evaluate the robustness of different slope estimation approaches in the case of noise. Fig. 4(a) shows the clean data, and Fig. 4(c) shows the slope estimated from the clean data using the PWD algorithm, which is deemed to be the accurate slope. Fig. 4(b) shows the noisy data by adding some Gaussian white noise. Fig. 4(d) shows the slope calculated using the robust slope estimation. It is salient that the slope estimated from the noisy data is fairly close to the accurate slope field. However, using the traditional PWD algorithm, it is difficult to obtain an acceptable slope estimation from the noisy data, as can be seen from the result shown in Fig. 4(e). From this test, we conclude that the robust slope estimation can be used to obtain robust slope estimation performance even in the presence of strong random noise. Figure 4. Open in new tabDownload slide Slope calculation test. (a) Clean data. (b) Noisy data. (c) Slope calculated from the clean data using PWD algorithm. (d) Slope calculated from the noisy data using the robust slope calculation algorithm. (e) Slope calculated from the noisy data using the PWD algorithm. Figure 4. Open in new tabDownload slide Slope calculation test. (a) Clean data. (b) Noisy data. (c) Slope calculated from the clean data using PWD algorithm. (d) Slope calculated from the noisy data using the robust slope calculation algorithm. (e) Slope calculated from the noisy data using the PWD algorithm. It is worth mentioning that, by eqs (10) and (15), we do not consider the spatial gradient of amplitude. In the case of smooth spatial amplitude change (e.g. small spatial gradient), the slope estimation method also works, since the calculation is done locally and the small spatial gradient almost has no influence. However, in the case of sharp spatial amplitude change (e.g. large spatial gradient), the method cannot be adopted. This drawback can be hopefully overcome in the future work. In addition, the problem of spatial gradients of amplitude and the implications for non-plane wave solutions was mentioned in Wielandt (1993). 2.4 Plane-wave orthogonal polynomial transform We have introduced in detail the theory of plane-wave trace continuation, that is, how we predict an arbitrary trace in a seismic gather from a random starting trace. We have shown that by discretizing the plane-wave equation, we can derive the spatial trace continuation relation, which can be used for trace prediction. We have presented that during trace continuation, the amplitude of seismic waveforms can be well preserved. Regarding the slope estimation, which is an important factor in the plane-wave trace continuation operator, we introduce the robust slope estimation approach. We also show that the robust slope estimation approach can obtain robust slope estimation in the presence of strong noise. Considering the amplitude-preserving capability of the OPT in a flattened dimension, we can cascade the plane-wave trace continuation operator and the OPT together to obtain a twofolds amplitude-preserving performance during a complete workflow. Thus, we name the cascaded framework as the plane-wave OPT. The complete framework for noise attenuation using the plane-wave OPT is shown in Algorithm 1. Open in new tabDownload slide Open in new tabDownload slide The forward OPT corresponds to inverting |$\mathbf {P}^H(\mathbf {P}\mathbf {P}^H)^{-1}$|⁠. The inverse OPT corresponds to multiplying the orthogonal polynomial coefficients by |$\mathbf {P}$|⁠. In Algorithm 1, the detailed implementations of the forward plane-wave flattening operator and the inverse plane-wave flattening operator are shown in Algorithms 2 and 3, respectively. Open in new tabDownload slide Open in new tabDownload slide Open in new tabDownload slide Open in new tabDownload slide In Algorithms 2 and 3, note that N denotes the number of spatial traces. 1 and −1 in the operator PWTC() denote predicting from a trace to the first trace and predicting the first trace to another trace, respectively. |$\mathbf {D}(i)$| and |$\hat{\mathbf {D}}(i)$| denote the ith column (or trace) in the matrix |$\mathbf {D}$| and |$\hat{\mathbf {D}}$|⁠. 3 EXAMPLES The first example is a synthetic example, as shown in Fig. 5. We apply the Karhunen–Loève (KL) filtering method (Jones & Levy, 1987) and the proposed method to a flattened data set with strong amplitude variation. Fig. 5(a) shows the clean data and Fig. 5(d) shows the noisy data. Figs 5(b) and (c) show the denoised data using the KL filtering method and the proposed method, respectively. Figs 5(e) and (f) show the removed random noise using two approaches. We can observe clearly from Figs 5(b) and (c) that the KL filtering causes significant damages to the events, while the proposed method preserves the amplitude-variation details successfully. Figure 5. Open in new tabDownload slide Synthetic example. (a) Clean data. (b) Denoised data using KL filtering. (c) Denoised data using the proposed method. (d) Noisy data. (e) Noise section corresponding to (b). (f) Noise section corresponding to (c). Figure 5. Open in new tabDownload slide Synthetic example. (a) Clean data. (b) Denoised data using KL filtering. (c) Denoised data using the proposed method. (d) Noisy data. (e) Noise section corresponding to (b). (f) Noise section corresponding to (c). In order to compare the amplitude between different seismic profiles in detail, we compare the amplitude for a single trace from each section shown in Fig. 5. The trace is chosen as the 20th trace in each section of Fig. 5. The comparison is presented in Fig. 6. A zoom-in comparison is shown in Fig. 6(b). The black line is from the clean data. The red line is from the noisy data. The blue line corresponds to the KL method. The green line corresponds to the proposed method. It is apparent that the green line is very close to the black line while the blue line deviates from the black line too much in most areas. This trace amplitude comparison further confirms the superior performance of the proposed algorithm. Figure 6. Open in new tabDownload slide Comparison of the 20th trace amplitude of each seismic gather in Fig. 5. The black line is from the clean data. The red line is from the noisy data. The blue line corresponds to the KL method. The green line corresponds to the proposed method. (a) Comparison of the whole trace. (b) Zoom-in comparison. Note that the black and green lines are very close to each other, thus the reconstruction error using the proposed approach is much less than the traditional method for most parts. Figure 6. Open in new tabDownload slide Comparison of the 20th trace amplitude of each seismic gather in Fig. 5. The black line is from the clean data. The red line is from the noisy data. The blue line corresponds to the KL method. The green line corresponds to the proposed method. (a) Comparison of the whole trace. (b) Zoom-in comparison. Note that the black and green lines are very close to each other, thus the reconstruction error using the proposed approach is much less than the traditional method for most parts. In order to numerically compare the denoising performance, we use the commonly used signal-to-noise ratio (SNR) defined as follows to quantitatively measure the performance (Chen & Fomel, 2015b): \begin{equation*} \mathrm{SNR}=10\log _{10}\frac{\Vert \mathbf {s} \Vert _2^2}{\Vert \mathbf {s} -\hat{\mathbf {s}}\Vert _2^2}. \end{equation*} (18) where |$\mathbf {s}$| denotes the noise-free data and |$\hat{\mathbf {s}}$| denotes the denoised data. In addition, to quantitatively measure the noise removal in the case of no discernable signal damage, we define the metric as the root-mean-square (RMS), \begin{equation*} \mathrm{RMS}=\Vert \mathbf {n}\Vert _2, \end{equation*} (19) where |$\mathbf {n}$| denotes the removed noise. Although the amplitude range varies a lot for different data sets, the RMS metric provides us a quantitative way to evaluate the noise removal performance for one specific data set among different denoising methods, for example, how much better method A performs than method B. In order to compare the performance of two methods in different noise level, we increase the variance of noise from 0.1 to 1.0, and calculate the SNRs of denoised data of both methods and show them in Table 1. To see the varied SNRs more vividly, we plot the data from Table 1 in Fig. 7. The black line shows the SNRs varying with input noise variances. The red line shows the SNRs corresponding to the KL method. The blue line shows the SNRs of the OPT method. It is obvious that both methods obtain large SNR improvement for all noise levels and the SNRs of the OPT method are always higher than the KL method. We can also observe clearly that the difference between the proposed OPT method and the KL method increases as noise variance becomes larger, which indicates that the proposed method outperforms the KL method more when the seismic data becomes noisier. Figure 7. Open in new tabDownload slide SNR diagrams of synthetic example. Figure 7. Open in new tabDownload slide SNR diagrams of synthetic example. Table 1. Comparison of SNRs in dB for different input noise level. The diagram corresponding to this table is shown in Fig. 7. Noise variance Input data (dB) KL (dB) OPT (dB) 0.1 2.60 13.08 17.58 0.2 −3.42 6.92 11.56 0.3 −6.94 3.05 8.04 0.4 −9.44 0.19 5.54 0.5 −11.38 −1.89 3.60 0.6 −12.97 −3.54 2.02 0.7 −14.30 −4.94 0.68 0.8 −15.46- −6.17 −0.48 0.9 −16.49 −7.26 −1.50 1.0 −17.40 −8.25 −2.42 Noise variance Input data (dB) KL (dB) OPT (dB) 0.1 2.60 13.08 17.58 0.2 −3.42 6.92 11.56 0.3 −6.94 3.05 8.04 0.4 −9.44 0.19 5.54 0.5 −11.38 −1.89 3.60 0.6 −12.97 −3.54 2.02 0.7 −14.30 −4.94 0.68 0.8 −15.46- −6.17 −0.48 0.9 −16.49 −7.26 −1.50 1.0 −17.40 −8.25 −2.42 Open in new tab Table 1. Comparison of SNRs in dB for different input noise level. The diagram corresponding to this table is shown in Fig. 7. Noise variance Input data (dB) KL (dB) OPT (dB) 0.1 2.60 13.08 17.58 0.2 −3.42 6.92 11.56 0.3 −6.94 3.05 8.04 0.4 −9.44 0.19 5.54 0.5 −11.38 −1.89 3.60 0.6 −12.97 −3.54 2.02 0.7 −14.30 −4.94 0.68 0.8 −15.46- −6.17 −0.48 0.9 −16.49 −7.26 −1.50 1.0 −17.40 −8.25 −2.42 Noise variance Input data (dB) KL (dB) OPT (dB) 0.1 2.60 13.08 17.58 0.2 −3.42 6.92 11.56 0.3 −6.94 3.05 8.04 0.4 −9.44 0.19 5.54 0.5 −11.38 −1.89 3.60 0.6 −12.97 −3.54 2.02 0.7 −14.30 −4.94 0.68 0.8 −15.46- −6.17 −0.48 0.9 −16.49 −7.26 −1.50 1.0 −17.40 −8.25 −2.42 Open in new tab For computational cost comparison, the KL method takes 0.62 s for processing the data shown in Fig. 5(d) while the proposed algorithm takes 0.01 s. The data contains 151 samples and 61 traces. The computation is done on a PC station equipped with an Intel Core i7 CPU clocked at 3.1 GHz and 16 GB of RAM. Note that both KL and OPT methods require the events to be flattened in order to obtain the best performance, thus we only compare the cost difference in the filtering stage. The second example is a pre-stack field data example. Fig. 8(a) shows the original data. Figs 8(b)–(d) show the denoised data using EMD method, KL method and the proposed method, respectively. Fig. 9 shows the removed noise sections using three approaches. Fig. 9(a) shows that some low-frequency energy is damaged while Fig. 9(c) shows that the removed noise is stronger. In this example, the calculated RMSs for Figs 9(a)–(c) are 389.49, 449.22 and 506.26, respectively. Thus, the proposed method removes 12.7 per cent more noise than the KL method and 30.0 per cent more noise than the EMD method. Figure 8. Open in new tabDownload slide Denoising comparison. (a) Raw noisy data. (b) Filtered using EMD method. (c) Filtered using KL method. (d) Filtered using the proposed method. Figure 8. Open in new tabDownload slide Denoising comparison. (a) Raw noisy data. (b) Filtered using EMD method. (c) Filtered using KL method. (d) Filtered using the proposed method. Figure 9. Open in new tabDownload slide Noise comparison. (b) Removed noise using EMD method. (b) Removed noise using KL method. (c) Removed noise using the proposed method. In this case, the calculated RMSs for (a)–(c) are 389.49, 449.22 and 506.26, respectively. Thus, the proposed method removes 12.7 per cent more noise than the KL method and 30.0 per cent more noise than the EMD method. Figure 9. Open in new tabDownload slide Noise comparison. (b) Removed noise using EMD method. (b) Removed noise using KL method. (c) Removed noise using the proposed method. In this case, the calculated RMSs for (a)–(c) are 389.49, 449.22 and 506.26, respectively. Thus, the proposed method removes 12.7 per cent more noise than the KL method and 30.0 per cent more noise than the EMD method. In order to comprehensively compare three different approaches, we zoomed four frame boxes (A, B, C and D) to show the detailed difference. Fig. 10 shows the comparison from frame box A. It is obvious that the KL approach causes some residual noise while EMD and OPT approaches obtain good results, more careful observation can show that the proposed method can obtain a more coherent image. Fig. 11 shows the comparison for frame box B. It is obvious that OPT method obtains the cleanest result. Fig. 12 shows the comparison for frame box C. It is still obvious that OPT method can make the events more coherent and more importantly, preserve the amplitude-variation-with-offsets details well. Fig. 13 shows the comparison for frame box D. Both Figs 13(b) and (c) show obvious amplitude artefacts while the OPT result in Fig. 13(d) shows nearly zero amplitude in the zoomed section. Figure 10. Open in new tabDownload slide Zoomed frame box A from Fig. 8. (a) Zoomed noisy field data. (b) Zoomed filtered data using EMD method. (c) Zoomed filtered data using KL method. (d) Zoomed filtered data using the proposed method. Figure 10. Open in new tabDownload slide Zoomed frame box A from Fig. 8. (a) Zoomed noisy field data. (b) Zoomed filtered data using EMD method. (c) Zoomed filtered data using KL method. (d) Zoomed filtered data using the proposed method. Figure 11. Open in new tabDownload slide Zoomed frame box B from Fig. 8. (a) Zoomed noisy field data. (b) Zoomed filtered data using EMD method. (c) Zoomed filtered data using KL method. (d) Zoomed filtered data using the proposed method. Figure 11. Open in new tabDownload slide Zoomed frame box B from Fig. 8. (a) Zoomed noisy field data. (b) Zoomed filtered data using EMD method. (c) Zoomed filtered data using KL method. (d) Zoomed filtered data using the proposed method. Figure 12. Open in new tabDownload slide Zoomed frame box C from Fig. 8. (a) Zoomed noisy field data. (b) Zoomed filtered data using EMD method. (c) Zoomed filtered data using KL method. (d) Zoomed filtered data using the proposed method. Figure 12. Open in new tabDownload slide Zoomed frame box C from Fig. 8. (a) Zoomed noisy field data. (b) Zoomed filtered data using EMD method. (c) Zoomed filtered data using KL method. (d) Zoomed filtered data using the proposed method. Figure 13. Open in new tabDownload slide Zoomed frame box D from Fig. 8. (a) Zoomed noisy field data. (b) Zoomed filtered data using EMD method. (c) Zoomed filtered data using KL method. (d) Zoomed filtered data using the proposed method. Figure 13. Open in new tabDownload slide Zoomed frame box D from Fig. 8. (a) Zoomed noisy field data. (b) Zoomed filtered data using EMD method. (c) Zoomed filtered data using KL method. (d) Zoomed filtered data using the proposed method. For this example, we also demonstrate the plane-wave flattening process in Fig. 14. Fig. 14(b) shows the flattened gather from the original data shown in Fig. 14(a) (or Fig.8(a)). It is clear that most events have been flattened well. Fig.14(c) shows the reconstructed data from the inverse plane-wave flattening. The data is almost the same as the original data. Fig. 14(d) shows the difference between the reconstructed data and the original data. The error section is almost zero everywhere, which demonstrates that the plane-wave flattening process does not introduce extra error. Fig. 15 shows the slope estimated from the original data. We also show a detailed comparison between different data in the flattened dimension in Fig. 16. A zoomed comparison among the flattened gathers after filtering is shown in Fig. 17, where we can conclude that the proposed method obtains the smoothest result while best preserving the reflection amplitude. Figure 14. Open in new tabDownload slide Pre-stack field data example. (a) Field data. (b) Flattened field data. (c) Reconstructed field data. (d) Reconstruction error. Figure 14. Open in new tabDownload slide Pre-stack field data example. (a) Field data. (b) Flattened field data. (c) Reconstructed field data. (d) Reconstruction error. Figure 15. Open in new tabDownload slide Local slope estimation of the pre-stack field data. Figure 15. Open in new tabDownload slide Local slope estimation of the pre-stack field data. Figure 16. Open in new tabDownload slide Denoising comparison in the flattened dimension. (a) Raw noisy data. (b) Filtered using EMD method. (b) Filtered using KL method. (c) Filtered using the proposed method. Figure 16. Open in new tabDownload slide Denoising comparison in the flattened dimension. (a) Raw noisy data. (b) Filtered using EMD method. (b) Filtered using KL method. (c) Filtered using the proposed method. Figure 17. Open in new tabDownload slide Zoomed sections from Fig. 16. (a) Zoomed noisy field data. (b) Zoomed filtered data using EMD method. (c) Zoomed filtered data using KL method. (d) Zoomed filtered data using the proposed method. Figure 17. Open in new tabDownload slide Zoomed sections from Fig. 16. (a) Zoomed noisy field data. (b) Zoomed filtered data using EMD method. (c) Zoomed filtered data using KL method. (d) Zoomed filtered data using the proposed method. The next example is a real post-stack seismic image shown in Fig. 18(a). There are 140 spatial traces and 194 temporal samples. The seismic image contains highly curved events and the amplitude along the events is not continuous, which will make the seismic interpretation difficult. After using three approaches, the EMD method, the KL method and the proposed method, the denoised images are shown in Figs 18(b)–(d), respectively. It is obvious that the proposed OPT based filtering approach can obtain a well smoothed seismic image with the continuity and the amplitude of events enhanced greatly. The EMD based approach, however, cannot effectively smooth the seismic events, and still leaves a lot of discontinuity in the image. The KL filtering approach obtains a much better filtering performance compared with the EMD based approach, however, it is not as successful as the performance of the proposed OPT based approach. Figure 18. Open in new tabDownload slide First post-stack field data example. (a) Field data. (b) Filtered data using EMD method. (c) Filtered data using KL method. (d) Filtered data using the proposed method. Figure 18. Open in new tabDownload slide First post-stack field data example. (a) Field data. (b) Filtered data using EMD method. (c) Filtered data using KL method. (d) Filtered data using the proposed method. We can find the mechanism that caused the tremendous difference of denoised images from the comparison in the flattened domain, as shown in Fig. 19. It is even more obvious that the OPT based approach obtains a nearly perfect smoothing along the flattened images (equivalent to along the structure in the original domain). The EMD based approach can achieve some smoothing, but remains less continuous than both KL and OPT based approaches. In this example, the calculated RMSs of removed noise are 0.015 for EMD method, 0.016 for KL method and 0.019 for the proposed method. Thus, the proposed method removes 18.7 per cent more noise than the KL method and 25.0 per cent more noise than the EMD method. Figure 19. Open in new tabDownload slide Comparison in the flattened domain. (a) Field data. (b) Filtered data using EMD method. (c) Filtered data using KL method. (d) Filtered data using the proposed method. Figure 19. Open in new tabDownload slide Comparison in the flattened domain. (a) Field data. (b) Filtered data using EMD method. (c) Filtered data using KL method. (d) Filtered data using the proposed method. The next field data example is shown in Fig. 20(a), which is also a post-stack data and contains weak seismic reflection events. Figs 20(b)–(d) show the denoised results using three different methods. For this example, we further compare the performance of the proposed method with that of the f–x predictive filtering method and the SSA method (Vautard et al., 1992). Fig. 21 shows the corresponding noise sections. For this example, it seems that all three methods obtain much improved results and the performance of different methods is quite similar. In order to compare the performance in detail and more fairly, we plot the F–K spectra of different denoised results. The F–K spectrum of the raw data is shown in Fig. 22(a). The F–K spectra corresponding to different methods are shown in Figs 22(b)–(d). Comparing the F–K spectra of different methods and F–K spectrum of the raw data, it is easy to find that both f–x predictive filtering method and the proposed method preserve the useful signals well, but the f–x predictive filtering method has some residual spectrum energy around the edges (large wavenumber components). SSA method causes significant damages to useful signals. Figure 20. Open in new tabDownload slide Denoising comparison for the second post-stack field data. (a) The second post-stack field data. (b) Filtered data using f–x predictive filtering. (c) Filtered data using SSA. (d) Filtered data using the proposed method. Figure 20. Open in new tabDownload slide Denoising comparison for the second post-stack field data. (a) The second post-stack field data. (b) Filtered data using f–x predictive filtering. (c) Filtered data using SSA. (d) Filtered data using the proposed method. Figure 21. Open in new tabDownload slide Noise comparison for the second post-stack field data. (a) Removed noise using f–x predictive filtering. (b) Removed noise using SSA. (c) Removed noise using the proposed method. Figure 21. Open in new tabDownload slide Noise comparison for the second post-stack field data. (a) Removed noise using f–x predictive filtering. (b) Removed noise using SSA. (c) Removed noise using the proposed method. Figure 22. Open in new tabDownload slide Spectra comparison. (a) Spectrum of the second post-stack seismic data. (b) Spectrum using f–x predictive filtering. (c) Spectrum using SSA. (d) Spectrum using the proposed method. Figure 22. Open in new tabDownload slide Spectra comparison. (a) Spectrum of the second post-stack seismic data. (b) Spectrum using f–x predictive filtering. (c) Spectrum using SSA. (d) Spectrum using the proposed method. In this example, we also calculate the local similarity between the denoised data and removed noise for different methods. The local similarity is an effective way to detect the lost signals in the removed noise. High local similarity indicates that in the noise section, there are significantly similar components as the useful signals, that is, there is lost energy in the noise. The calculation of local similarity is provided in Appendix B. The local similarity maps for different methods are shown in Fig. 23, where we can clearly observe the high similarity anomalies in the f–x predictive filtering and SSA results. Although there are also some similarity anomalies in the result from the proposed method, the similarity value is relatively lower than the other two methods. From this test we conclude that the proposed method causes less damage to useful energy. Figure 23. Open in new tabDownload slide Comparison of local similarity between denoised data and removed noise. (a) Local similarity using f–x predictive filtering. (b) Local similarity using SSA. (c) Local similarity using the proposed method. Figure 23. Open in new tabDownload slide Comparison of local similarity between denoised data and removed noise. (a) Local similarity using f–x predictive filtering. (b) Local similarity using SSA. (c) Local similarity using the proposed method. We also plot a comparison of the average spectrum of all the traces for different data in Fig. 24. The green line corresponds to the proposed approach. The red line corresponds to f–x predictive filtering method. The blue line corresponds to the SSA method. It is quite obvious that the energy preservation of the proposed method in signal frequency band (20–60 Hz) is quite successful. The proposed method mitigates more high-frequency noise than f − x predictive filtering method, which confirms the observation from Fig. 22. We admit that the high-frequency noise of the proposed method is slightly more than the SSA method. However, the proposed method preserves more useful energy than the other two methods in the spectrum. This field data further confirms the superior performance of the presented algorithm. Figure 24. Open in new tabDownload slide Comparisons of the average spectrum of all the traces. The black line denotes the average spectrum of raw data. The green line corresponds to the proposed approach. The red line corresponds to f–x predictive filtering method. The blue line corresponds to the SSA method. Figure 24. Open in new tabDownload slide Comparisons of the average spectrum of all the traces. The black line denotes the average spectrum of raw data. The green line corresponds to the proposed approach. The red line corresponds to f–x predictive filtering method. The blue line corresponds to the SSA method. In this example, to compare the noise removal performance, we need to make sure the removed noise sections do not contain discernable signal energy, as required by the metric defined in eq. (19), and have to adjust the parameters for the f–x and SSA methods. The denoised data and the removed noise sections using the adjusted parameters are shown in Fig. 25. In this case, the calculated RMSs for Figs 25(d)–(f) are 0.059, 0.071 and 0.086, respectively. Thus, the proposed method removes 21.1 per cent more noise than the SSA method and 38.0 per cent more noise than the f–x method. Figure 25. Open in new tabDownload slide Comparison for the second post-stack field data after adjusting the parameters. (a) Filtered data using f–x predictive filtering. (b) Filtered data using SSA. (c) Filtered data using the proposed method. (d) Removed noise using f–x predictive filtering. (e) Removed noise using SSA. (f) Removed noise using the proposed method. In this case, the calculated RMSs for (d)–(f) are 0.059, 0.071 and 0.086, respectively. Thus, the proposed method removes 21.1 per cent more noise than the SSA method and 38.0 per cent more noise than the f − x predictive filtering method. Figure 25. Open in new tabDownload slide Comparison for the second post-stack field data after adjusting the parameters. (a) Filtered data using f–x predictive filtering. (b) Filtered data using SSA. (c) Filtered data using the proposed method. (d) Removed noise using f–x predictive filtering. (e) Removed noise using SSA. (f) Removed noise using the proposed method. In this case, the calculated RMSs for (d)–(f) are 0.059, 0.071 and 0.086, respectively. Thus, the proposed method removes 21.1 per cent more noise than the SSA method and 38.0 per cent more noise than the f − x predictive filtering method. 4 CONCLUSIONS The OPT can be used to effectively separate spatially correlative signals and spatially incoherent noise without losing waveform amplitude. To create a flattened dimension where the OPT can be optimally applied, we derive a plane-wave trace continuation relation for flattening the curved seismic events. Trace prediction in the flattening process and the subsequent OPT are both demonstrated to be amplitude preserving. The robust slope estimation approach can obtain more robust performance than the state-of-the-art PWD method in the presence of strong noise. The proposed framework has been applied to several synthetic, field pre-stack and post-stack seismic data and is shown to better preserve the amplitude variations than other alterative methods. ACKNOWLEDGEMENTS This work was supported by the Starting fund at Zhejiang University, National Natural Science Foundation of China (No. 61401307), Hebei Province Foundation of Returned Overseas Scholars (CL201707), Hebei Province Project of Science and Technology R & D (11213565) and Hebei Province Natural Science Foundation (No. E2016202341). The authors would like to thank Yaru Xue, Shaohuan Zu and Wei Chen for helpful comments and suggestions on the topic of noise attenuation. REFERENCES Abma R. , Claerbout J. , 1995 . Lateral prediction for noise attenuation by t −x and f −x techniques , Geophysics , 60 , 1887 – 1896 .. https://doi.org/10.1190/1.1443920 Google Scholar Crossref Search ADS WorldCat Anvari R. , Siahsar M.A.N. , Gholtashi S. , Kahoo A.R. , Mohammadi M. , 2017 . Seismic random noise attenuation using synchrosqueezed wavelet transform and low-rank signal matrix approximation , IEEE Trans. Geosci. Remote Sens. , 55 ( 11 ), 6574 – 6581 .. https://doi.org/10.1109/TGRS.2017.2730228 Google Scholar Crossref Search ADS WorldCat Asgedom E.G. , Orji O.C. , Soellner W. , 2017 . Rough-sea deghosting of single-sensor seismic data using the knowledge of the sea surface shape , J. Seismic Explor. , 26 , 105 – 123 . WorldCat Bai M. , Wu J. , 2017 . Efficient deblending using median filtering without correct normal moveout - with comparison on migrated images , J. Seismic Explor. , 26 , 455 – 479 . WorldCat Bai M. , Wu J. , 2018 . Seismic deconvolution using iteartive transform-domain sparse inversion , J. Seismic Explor. , 27 ( 2 ), 103 – 116 . WorldCat Bai M. , Wu J. , Xie J. , Zhang D. , 2018a . Least-squares reverse time migration of blended data with low-rank constraint along structural direction , J. Seismic Explor. , 27 ( 1 ), 29 – 48 . WorldCat Bai M. , Wu J. , Zu S. , Chen W. , 2018b . A structural rank reduction operator for removing artifacts in least-squares reverse time migration , Comput. Geosci. , 117 , 9 – 20 .. https://doi.org/10.1016/j.cageo.2018.04.003 Google Scholar Crossref Search ADS WorldCat Canales L. , 1984 . Random noise reduction , 54th Annual International Meeting , SEG, Expanded Abstracts , Atlanta, GA , pp. 525 – 527 . Google Preview WorldCat COPAC Candès E.J. , Demanet L. , Donoho D.L. , Ying L. , 2006 . Fast discrete curvelet transforms , SIAM Multiscale Model. Simul. , 5 , 861 – 899 .. https://doi.org/10.1137/05064182X Google Scholar Crossref Search ADS WorldCat Chen Y. , 2016 . Dip-separated structural filtering using seislet thresholding and adaptive empirical mode decomposition based dip filter , Geophys. J. Int. , 206 , 457 – 469 .. https://doi.org/10.1093/gji/ggw165 Google Scholar Crossref Search ADS WorldCat Chen Y. , 2017 . Fast dictionary learning for noise attenuation of multidimensional seismic data , Geophys. J. Int. , 209 ( 1 ), 21 – 31 . WorldCat Chen Y. , 2018 . Automatic microseismic event picking via unsupervised machine learning , Geophys. J. Int. , 212 ( 1 ), 88 – 102 .. https://doi.org/10.1093/gji/ggx420 Google Scholar Crossref Search ADS WorldCat Chen Y. , Fomel S. , 2015a . EMD-seislet transform , in 85th Annual International Meeting , SEG, Expanded Abstracts , New Orleans, LA , pp. 4775 – 4778 . Google Preview WorldCat COPAC Chen Y. , Fomel S. , 2015b . Random noise attenuation using local signal-and-noise orthogonalization , Geophysics , 80 , WD1 – WD9 .. https://doi.org/10.1190/geo2014-0227.1 Google Scholar Crossref Search ADS WorldCat Chen W. , Chen Y. , Xie J. , Zu S. , Zhang Y. , 2016a . Multiples attenuation using trace randomization and empirical mode decomposition , in 86th Annual International Meeting , SEG, Expanded Abstracts , Dallas, TX , pp. 4498 – 4502 . Google Preview WorldCat COPAC Chen Y. , Zhang D. , Jin Z. , Chen X. , Zu S. , Huang W. , Gan S. , 2016b . Simultaneous denoising and reconstruction of 5-D seismic data via damped rank-reduction method , Geophys. J. Int. , 206 ( 3 ), 1695 – 1717 .. https://doi.org/10.1093/gji/ggw230 Google Scholar Crossref Search ADS WorldCat Colominas M.A. , Schlotthauer G. , Torres M.E. , Flandrin P. , 2012 . Noise-assisted EMD methods in action , Adv. Adapt. Data Anal. , 4 , 1250025 https://doi.org/10.1142/S1793536912500252 . Google Scholar Crossref Search ADS WorldCat Daubechies I. , Lu J. , Wu H.-T. , 2011 . Synchrosqueezed wavelet transforms: an empirical mode decomposition-like tool , Appl. Comput. Harmon. Anal. , 30 , 243 – 261 .. https://doi.org/10.1016/j.acha.2010.08.002 Google Scholar Crossref Search ADS WorldCat Dragomiretskiy K. , Zosso D. , 2014 . Variational mode decomposition , IEEE Trans. Signal Process. , 62 , 531 – 544 .. https://doi.org/10.1109/TSP.2013.2288675 Google Scholar Crossref Search ADS WorldCat Fomel S. , 2002 . Application of plane-wave destruction filters , Geophysics , 67 , 1946 – 1960 .. https://doi.org/10.1190/1.1527095 Google Scholar Crossref Search ADS WorldCat Fomel S. , Liu Y. , 2010 . Seislet transform and seislet frame , Geophysics , 75 , V25 – V38 .. https://doi.org/10.1190/1.3380591 Google Scholar Crossref Search ADS WorldCat Foster D.J. , Mosher C.C. , 1992 . Suppression of multiple reflections using the radon transform , Geophysics , 57 , 386 – 395 .. https://doi.org/10.1190/1.1443253 Google Scholar Crossref Search ADS WorldCat Gan S. , Wang S. , Chen Y. , Chen X. , 2015a . Deblending of distance separated simultaneous-source data using seislet frames in the shot domain , in 85th Annual International Meeting , SEG Expanded Abstracts , New Orleans, LA , pp. 65 – 70 . Google Preview WorldCat COPAC Gan S. , Wang S. , Chen Y. , Chen X. , 2015b . Seismic data reconstruction via fast projection onto convex sets in the seislet transform domain , in 85th Annual international meeting , SEG Expanded Abstracts , New Orleans, LA , pp. 3814 – 3819 . Google Preview WorldCat COPAC Gan S. , Wang S. , Chen Y. , Qu S. , Zu S. , 2016 . Velocity analysis of simultaneous-source data using high-resolution semblance-coping with the strong noise , Geophys. J. Int. , 204 , 768 – 779 .. https://doi.org/10.1093/gji/ggv484 Google Scholar Crossref Search ADS WorldCat Gao Z. , Pan Z. , Gao J. , 2016 . Multimutation differential evolution algorithm and its application to seismic inversion , IEEE Trans. Geosci. Remote Sens. , 54 ( 6 ), 3626 – 3636 .. https://doi.org/10.1109/TGRS.2016.2520978 Google Scholar Crossref Search ADS WorldCat Gholami A. , 2013 . Sparse time-frequency decomposition and some applications , IEEE Trans. Geosci. Remote Sens. , 51 , 3598 – 3604 .. https://doi.org/10.1109/TGRS.2012.2220144 Google Scholar Crossref Search ADS WorldCat Herrmann F.J. , Böniger U. , Verschuur D.J. , 2007 . Non-linear primary-multiple separation with directional curvelet frames , Geophys. J. Int. , 170 , 781 – 799 .. https://doi.org/10.1111/j.1365-246X.2007.03360.x Google Scholar Crossref Search ADS WorldCat Herrmann F.J. , Hennenfent G. , 2008 . Non-parametric seismic data recovery with curvelet frames , Geophys. J. Int. , 173 , 233 – 248 .. https://doi.org/10.1111/j.1365-246X.2007.03698.x Google Scholar Crossref Search ADS WorldCat Huang N.E. et al. , 1998 . The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis , Phil. Trans. R. Soc. Lond., A. , 454 , 903 – 995 . WorldCat Huang W. , Wang R. , Zu S. , Chen Y. , 2017 . Low-frequency noise attenuation in seismic and microseismic data using mathematical morphological filtering , Geophys. J. Int. , 211 , 1318 – 1340 . WorldCat Huang Z. , Zhang J. , Zhao T. , Sun Y. , 2016 . Synchrosqueezing S-transform and its application in seismic spectral decomposition , IEEE Trans. Geosci. Remote Sens. , 54 ( 2 ), 817 – 825 .. https://doi.org/10.1109/TGRS.2015.2466660 Google Scholar Crossref Search ADS WorldCat Jones I.F. , Levy S. , 1987 . Signal-to-noise ratio enhancement in multichannel seismic data via the Karhunen-Loeve transform , Geophys. Prospect. , 35 , 12 – 32 .. https://doi.org/10.1111/j.1365-2478.1987.tb00800.x Google Scholar Crossref Search ADS WorldCat Kahoo A. , Siahkoohi T.R. , 2009 . Random noise suppression from seismic data using time frequency peak filtering , in 71st Annual International Conference and Exhibition , EAGE, Extended Abstracts , doi:10.3997/2214-4609.201400214 . Google Preview WorldCat COPAC Kong D. , Peng Z. , He Y. , Fan H. , 2016 . Seismic random noise attenuation using directional total variation in the shearlet domain , J. Seismic Explor. , 25 ( 4 ), 321 – 338 . WorldCat Lin H. , Li Y. , Yang B. , Ma T. , 2013 . Random denoising and signal nonlinearity approach by time-frequency peak filtering using weighted frequency reassignment , Geophysics , 78 , V229 – V237 .. https://doi.org/10.1190/geo2012-0432.1 Google Scholar Crossref Search ADS WorldCat Lin H. , Li Y. , Ma H. , Yang B. , Dai J. , 2015 . Matching-pursuit-based spatial-trace time-frequency peak filtering for seismic random noise attenuation , IEEE Geosci. Remote Sens. Lett. , 12 , 394 – 398 .. https://doi.org/10.1109/LGRS.2014.2344020 Google Scholar Crossref Search ADS WorldCat Liu C. , Chen C. , Wang D. , Liu Y. , Wang S. , Zhang L. , 2015 . Seismic dip estimation based on the two-dimensional hilbert transform and its application in random noise attenuation , Appl. Geophys. 12 , 55 – 63 .. https://doi.org/10.1007/s11770-014-0474-4 Google Scholar Crossref Search ADS WorldCat Liu G. , Chen X. , 2013 . Noncausal f– x– y regularized nonstationary prediction filtering for random noise attenuation on 3-D seismic data , J. Appl. Geophys. , 93 , 60 – 66 .. https://doi.org/10.1016/j.jappgeo.2013.03.007 Google Scholar Crossref Search ADS WorldCat Liu G. , Fomel S. , Jin L. , Chen X. , 2009a . Stacking seismic data using local correlation , Geophysics , 74 , V43 – V48 .. https://doi.org/10.1190/1.3085643 Google Scholar Crossref Search ADS WorldCat Liu G. , Chen X. , Du J. , Song J. , 2011 . Seismic noise attenuation using nonstationary polynomial fitting , Appl. Geophys. , 8 , 18 – 26 .. https://doi.org/10.1007/s11770-010-0244-2 Google Scholar Crossref Search ADS WorldCat Liu G. , Chen X. , Du J. , Wu K. , 2012 . Random noise attenuation using f-x regularized nonstationary autoregression , Geophysics , 77 , V61 – V69 .. https://doi.org/10.1190/geo2011-0117.1 Google Scholar Crossref Search ADS WorldCat Liu W. , Cao S. , Chen Y. , 2016a . Applications of variational mode decomposition in seismic time-frequency analysis , Geophysics , 81 ( 5 ), V365 – V378 .. https://doi.org/10.1190/geo2015-0489.1 Google Scholar Crossref Search ADS WorldCat Liu W. , Cao S. , Chen Y. , 2016b . Seismic time-frequency analysis via empirical wavelet transform , IEEE Geosci. Remote Sens. Lett. , 13 , 28 – 32 .. https://doi.org/10.1109/LGRS.2015.2493198 Google Scholar Crossref Search ADS WorldCat Liu W. , Cao S. , Wang Z. , 2017 . Application of variational mode decomposition to seismic random noise reduction , J. Geophys. Eng. , 14 , https://doi.org/10.1088/1742-2140/aa6b28 . WorldCat Liu Y. , 2013 . Noise reduction by vector median filtering , Geophysics , 78 , V79 – V87 .. https://doi.org/10.1190/geo2012-0232.1 Google Scholar Crossref Search ADS WorldCat Liu Y. , Liu C. , Wang D. , 2009b . A 1d time-varying median filter for seismic random, spike-like noise elimination , Geophysics , 74 , V17 – V24 .. https://doi.org/10.1190/1.3043446 Google Scholar Crossref Search ADS WorldCat Lorenzi L. , Melgani F. , Mercier G. , 2016 . Missing-area reconstruction in multispectral images under a compressive sensing perspective , IEEE Trans. Geosci. Remote Sens. , 51 ( 7 ), 3998 – 4008 .. https://doi.org/10.1109/TGRS.2012.2227329 Google Scholar Crossref Search ADS WorldCat Mousavi S.M. , Langston C.A. , 2016a . Adaptive noise estimation and suppression for improving microseismic event detection , J. Appl. Geophys. , 132 , 116 – 124 .. https://doi.org/10.1016/j.jappgeo.2016.06.008 Google Scholar Crossref Search ADS WorldCat Mousavi S.M. , Langston C.A. , 2016b . Hybrid seismic denoising using higher-order statistics and improved wavelet block thresholding , Bull. seism. Soc. Am. , 106 ( 4 ), 1380 – 1393 .. https://doi.org/10.1785/0120150345 Google Scholar Crossref Search ADS WorldCat Mousavi S.M. , Langston C.A. , 2017 . Automatic noise-removal/signal-removal based on general cross-validation thresholding in synchrosqueezed domain and its application on earthquake data , Geophysics , 82 ( 4 ), V211 – V227 .. https://doi.org/10.1190/geo2016-0433.1 Google Scholar Crossref Search ADS WorldCat Mousavi S.M. , Langston C.A. , Horton S.P. , 2016 . Automatic microseismic denoising and onset detection using the synchrosqueezed continuous wavelet transform , Geophysics . 81 ( 4 ), V341 – V355 .. https://doi.org/10.1190/geo2015-0598.1 Google Scholar Crossref Search ADS WorldCat Qu S. , Verschuur D. , Chen Y. , 2016 . Full waveform inversion using an automatic directional total variation constraint , in 78th Annual International Conference and Exhibition , EAGE, Extended Abstracts , doi:10.3997/2214–4609.201701340 . Google Preview WorldCat COPAC Trickett S. , 2008 . F-xy cadzow noise suppression , in CSPG CSEG CWLS Convention , Las Vegas, NV , pp. 303 – 306 . Google Preview WorldCat COPAC Vautard R. , Yiou P. , Ghil M. , 1992 . Singular-spectrum analysis: a toolkit for short, noisy chaotic signals , Physica D Nonlinear Phenom. , 58 , 95 – 126 .. https://doi.org/10.1016/0167-2789(92)90103-T Google Scholar Crossref Search ADS WorldCat Wang B. , Wu R. , Chen X. , Li J. , 2015 . Simultaneous seismic data interpolation and denoising with a new adaptive method based on dreamlet transform , Geophys. J. Int. , 201 , 1180 – 1192 . WorldCat Wang Y. , Cao J. , Yang C. , 2011 . Recovery of seismic wavefields based on compressive sensing by an l1-norm constrained trust region method and the piecewise random subsampling , Geophys. J. Int. , 187 , 199 – 213 .. https://doi.org/10.1111/j.1365-246X.2011.05130.x Google Scholar Crossref Search ADS WorldCat Wielandt E. , 1993 . Propagation and structural interpretation of non-plane waves , Geophys. J. Int. , 113 , 45 – 53 .. https://doi.org/10.1111/j.1365-246X.1993.tb02527.x Google Scholar Crossref Search ADS WorldCat Wu G. , Fomel S. , Chen Y. , 2016 . Data-driven time-frequency analysis of seismic data using regularized nonstationary autoregression , in 86th Annual International Meeting , SEG, Expanded Abstracts , Dallas, TX , pp. 1700 – 1705 . Google Preview WorldCat COPAC Wu J. , Bai M. , 2018a . Adaptive rank-reduction method for seismic data reconstruction , J. Geophys. Eng. , 15 , 1688 . https://doi.org/10.1088/1742-2140/aabc74 Google Scholar Crossref Search ADS WorldCat Wu J. , Bai M. , 2018b . Attenuating seismic noise via incoherent dictionary learning , J. Geophys. Eng. , 15 , 1327 . https://doi.org/10.1088/1742-2140/aaaf57 Google Scholar Crossref Search ADS WorldCat Wu J. , Bai M. , 2018c . Fast principal component analysis for stacking seismic data , J. Geophys. Eng. , 15 , 295 – 306 .. https://doi.org/10.1088/1742-2140/aa9f80 Google Scholar Crossref Search ADS WorldCat Wu J. , Bai M. , 2018d . Incoherent dictionary learning for reducing crosstalk noise in least-squares reverse time migration , Comput. Geosci. , 114 , 11 – 21 .. https://doi.org/10.1016/j.cageo.2018.01.010 Google Scholar Crossref Search ADS WorldCat Wu Z. , Huang N.E. , 2009 . Ensemble empirical mode decomposition: a noise-assisted data analysis method , Adv. Adapt. Data Anal. , 1 , 1 – 41 .. https://doi.org/10.1142/S1793536909000047 Google Scholar Crossref Search ADS WorldCat Xie J. , Di B. , Wei J. , Xie Q. , Zu S. , Chen. Y. , 2016 . Stacking using truncated singular value decomposition and local similarity , in 78th Annual International Conference and Exhibition , EAGE, Extended Abstracts , doi:10.3997/2214–4609.201601325 . Google Preview WorldCat COPAC Xue Y. , Chang F. , Zhang D. , Chen Y. , 2016a . Simultaneous sources separation via an iterative rank-increasing method , IEEE Geosci. Remote Sens. Lett. , 13 ( 12 ), 1915 – 1919 .. https://doi.org/10.1109/LGRS.2016.2617338 Google Scholar Crossref Search ADS WorldCat Xue Z. , Alger N. , Fomel S. , 2016b . Full-waveform inversion using smoothing kernels , in 86th Annual International Meeting , SEG Expanded Abstracts , Dallas, TX , pp. 1358 – 1363 . Google Preview WorldCat COPAC Yang W. , Wang R. , Chen Y. , Wu J. , 2014 . Random noise attenuation using a new spectral decomposition method , in 84th Annual International Meeting , SEG Expanded Abstracts , DEnver, CO , pp. 4366 – 4370 . Google Preview WorldCat COPAC Zeng Q. , Guo Y. , Jiang R. , Ba J. , Ma H. , Liu J. , 2017 . Fluid sensitivity of rock physics parameters in reservoirs: quantitative analysis , J. Seismic Explor. , 26 , 125 – 140 . WorldCat Zhang D. , Chen Y. , Gan S. , 2016a . Iterative reconstruction of 3D seismic data via multiple constraints , in 78th Annual International Conference and Exhibition , EAGE, Extended Abstracts , https://doi.org/10.3997/2214–4609.201601241 . Google Preview WorldCat COPAC Zhang D. , Chen Y. , Gan S. , 2016b . Multidimensional seismic data reconstruction with multiple constraints , in 86th Annual International Meeting , SEG, Expanded Abstracts , Dallas, TX , pp. 4801 – 4806 . Google Preview WorldCat COPAC Zhang D. , Chen Y. , Huang W. , Gan S. , 2016c . Multi-step reconstruction of 3D seismic data via an improved MSSA algorithm , in CPS/SEG Beijing 2016 International Geophysical Conference & Exposition , SEG, Expanded Abstracts , Beijing, China , pp. 745 – 749 . Google Preview WorldCat COPAC Zhang P. , Dai Y. , Wang R. , Tan Y. , 2017 . A quantitative evaluation method based on EMD for determining the accuracy of time-varying seismic wavelet extraction , J. Seismic Explor. , 26 , 267 – 292 . WorldCat Zhou Y. , Han W. , 2018 . Multiples attenuation in the presence of blending noise , J. Seismic Explor. , 27 ( 1 ), 69 – 88 . WorldCat Zhou Y. , Li S. , Zhang D. , Chen Y. , 2018 . Seismic noise attenuation using an online subspace tracking algorithm , Geophys. J. Int. , 212 ( 2 ), 1072 – 1097 .. https://doi.org/10.1093/gji/ggx422 Google Scholar Crossref Search ADS WorldCat Zu S. , Zhou H. , Chen Y. , Qu S. , Zou X. , Chen H. , Liu R. , 2016a . A periodically varying code for improving deblending of simultaneous sources in marine acquisition , Geophysics , 81 , V213 – V225 .. https://doi.org/10.1190/geo2015-0447.1 Google Scholar Crossref Search ADS WorldCat Zu S. , Zhou H. , Chen Y. , Pan X. , Gan S. , Zhang D. , Xie C. , 2016b . Recovering the most from big gaps using least-squares inversion , in 86th Annual International Meeting , SEG, Expanded Abstracts , Dallas, TX , pp. 4128 – 4133 . Google Preview WorldCat COPAC Zu S. , Zhou H. , Mao W. , Zhang D. , Li C. , Pan X. , Chen Y. , 2017 . Iterative deblending of simultaneous-source data using a coherency-pass shaping operator , Geophys. J. Int. , 211 ( 1 ), 541 – 557 .. https://doi.org/10.1093/gji/ggx324 Google Scholar Crossref Search ADS WorldCat APPENDIX A: CONSTRUCTION OF POLYNOMIAL TRANSFORMS Let {Pj(x)}, j = 0, 1, ⋅⋅⋅, N denote a set of polynomials, which satisfies the orthogonality condition: \begin{equation*} \sum _{i=0}^{N} P_k(x_i)P_j(x_i) = \delta _{j,k}. \end{equation*} (A1) It is known that as polynomials, Pj(xi) can be expressed \begin{equation*} P_j(x_i) = \sum _{k=0}^j a_{jk}x_i^k, \end{equation*} (A2)ajk denotes polynomial coefficients. It is natural that xj can be expressed based on superposition of different polynomials: \begin{equation*} x^j = \sum _{k=0}^j \beta _{jk} P_k(x). \end{equation*} (A3) Based on eqs (A3) and (A4), jth polynomial can be expressed as lower order polynomials \begin{equation*} P_j(x_i) = \left\lbrace x^j - \sum _{k=0}^{j-1}\beta _{jk}P_k(x_i) \right\rbrace /\beta _{jj}, \end{equation*} (A4) Get squares of eq. (A3) and combine with eq. (A1), we can obtain \begin{equation*} \beta _{jj} = \sqrt{\sum _{i=0}^{N} x_i^{2j} -\sum _{k=0}^{j-1}\beta _{jk}^2} \end{equation*} (A5) and \begin{equation*} \beta _{jk} = \sum _{i=0}^{N} x_i^j P_k(x_i). \end{equation*} (A6) From eqs (A4)–(A6), we can construct the set of polynomials. We first get |$\beta _{00}=\sqrt{N}$| based on eq. (A5), and thus P0 = 1/β00, then compute β10, β11 to construct P1. In the same way, we can construct all polynomials. APPENDIX B: LOCAL SIMILARITY Local similarity between vectors |$\mathbf {a}$| and |$\mathbf {b}$| is defined as \begin{equation*} \mathbf {c}=\sqrt{\mathbf {c}_1\circ \mathbf {c}_2} \end{equation*} (B1) where ○ denotes dot product, |$\mathbf {c}_1$| and |$\mathbf {c}_2$| come from two least-squares minimization problems: \begin{equation*} \mathbf {c}_1 =\arg \min _{\mathbf {c}_1}\Vert \mathbf {a}-\mathbf {B} \mathbf {c}_1 \Vert _2^2 \end{equation*} (B2) \begin{equation*} \mathbf {c}_2 =\arg \min _{\mathbf {c}_2}\Vert \mathbf {b}-\mathbf {A} \mathbf {c}_2 \Vert _2^2 , \end{equation*} (B3) where |$\mathbf {A}$| is a diagonal operator composed of the elements of |$\mathbf {a}$| and |$\mathbf {B}$| is a diagonal operator composed of the elements of |$\mathbf {b}$|⁠. Note that in eqs (B1)–(B3), |$\mathbf {a}$|⁠, |$\mathbf {b}$| and |$\mathbf {c}$| denote vectorized 2-D matrices. Eqs (B2) and (B3) can be solved using shaping regularization with a local-smoothness constraint: \begin{equation*} \mathbf {c}_1 = [\lambda _1^2\mathbf {I} + \mathbf {T}(\mathbf {B}^T\mathbf {B}-\lambda _1^2\mathbf {I})]^{-1}\mathbf {TB}^T\mathbf {b}, \end{equation*} (B4) \begin{equation*} \mathbf {c}_2 = [\lambda _2^2\mathbf {I} + \mathbf {T}(\mathbf {A}^T\mathbf {A}-\lambda _2^2\mathbf {I})]^{-1}\mathbf {TA}^T\mathbf {a}, \end{equation*} (B5) where |$\mathbf {T}$| is a smoothing operator and λ1 and λ2 are two parameters controlling the physical dimensionality and enabling fast convergence when inversion is implemented iteratively. These two parameters can be chosen as |$\lambda _1 = \Vert \mathbf {B}^T\mathbf {B}\Vert _2$| and |$\lambda _2 = \Vert \mathbf {A}^T\mathbf {A}\Vert _2$|⁠. © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices)
TI - Plane-wave orthogonal polynomial transform for amplitude-preserving noise attenuation
JF - Geophysical Journal International
DO - 10.1093/gji/ggy267
DA - 2018-09-01
UR - https://www.deepdyve.com/lp/oxford-university-press/plane-wave-orthogonal-polynomial-transform-for-amplitude-preserving-zF8n9JS7Tb
SP - 2207
VL - 214
IS - 3
DP - DeepDyve
ER -