TY - JOUR
AU - Wang,, Shangxu
AB - Abstract Attenuating random noise while preserving edges and structures of seismic signals is of great significance for seismic interpretation and inversion. An edge-preserving and signal-preserving noise reduction method is presented in this paper. The method regards noise reduction problem as an inverse problem based on a Bayesian framework. A priori information that the difference data of seismic signals obey Cauchy distribution, a kind of non-Gaussian distribution, is used as a constraint to control the noise reduction. The performance of the method is mainly dependent on a trade-off parameter and a scale parameter, which make it easier to obtain a relatively perfect signal preservation and noise reduction result. With directional constraints, the method recovers linear events and quasi-linear events without dip steering. In particular, the method can adaptively protect edges of discontinuous events. A synthetic data example and two field post-stack data sets illustrate the effectiveness of the method. 93.85.Bc, 93.85.Rt noise reduction, Bayesian inversion, non-Gaussian distribution, edge-preserving, inverse problem 1. Introduction An objective of exploration geophysics is to reveal the subsurface structural features and physical properties by using seismic signals. However, the effective seismic signals are always smeared by random noise. Therefore, reducing random noise while preserving useful seismic events, especially for edges of events which often correspond to important geological features such as faults, pinch-outs, channels, lens or fractures, is a critical task of seismic data processing (Luo et al2002, Bonar and Sacchi 2012). In the past decades, a large number of random noise reduction methods have been proposed. For example, f–x deconvolution (Canales 1984), based on the assumption that signals are predictable in the f–x domain while random noise is not, is widely used in real seismic data processing. However, this method predicts noise as source noise rather than as additive white noise, and it cannot protect edges of discontinuous events well. Some improved methods, such as f–x projection (Soubaras 1994) and f–x ARMA filter (Sacchi and Kuehl 2001), overcome inconsistency between source noise and additive noise, but they are still not good at the preservation of signal edges. f–x Cadzow filtering (Trickett 2008) is essentially a singular value decomposition (SVD) filtering method in the f–x domain. Compared with local SVD filtering methods in the t–x domain (Lu 2006, Bekara and van der Baan 2007), f–x Cadzow filtering can perfectly recover conflicting events without dip steering (Yuan and Wang 2011). However it still fails to preserve edges of discontinuous events, because it is based on the same assumption as f–x deconvolution. Recently, multichannel singular spectrum analysis was presented to simultaneously reduce random noise and reconstruct missing seismic traces (Oropeza and Sacchi 2011). It is essentially an equivalent algorithm to f–x Cadzow filtering. The edge-preserving smoothing (EPS) method (Luo et al2002) focuses on preserving edges, but with an implicit assumption that all seismic events are horizontal (Liu et al2003). The dipping events will be strongly distorted if they are directly processed without dip steering (Yuan et al2012). In addition, the EPS method pays more attention to edge preservation than signal preservation (Lu and Lu 2009), thus it is usually used as a preprocessing stage for seismic data interpretation. In the earlier work, Yuan et al (2012) introduced an edge-preserving noise reduction method using Bayesian inversion. For this method, a priori information that total variation (TV) data of seismic data satisfy the Laplace distribution, a form of non-Gaussian distribution, is adopted to control the noise reduction. In this paper, we propose an adaptively edge-preserving noise reduction method with directional difference constraints. The method can directly retrieve dipping events and adaptively preserve edges of discontinuous events. With appropriate parameters, the method will achieve a good trade-off between signal fidelity and noise reduction. Compared with our earlier work, the method we proposed in this paper has three main improvements. Firstly, we use the Cauchy distribution, another kind of non-Gaussian distribution, to control noise reduction, since using TV function as a constraint readily makes peaks and troughs of events blocky or flat. And we change matrix expression and computation into vector expression and computation, since vector form is more straightforward and easier to extend to other kinds of non-Gaussian distribution about directional difference data. Secondly, we extend two directional difference constraints in earlier works to more directional difference constraints. This improvement ensures our method is more flexible and better at protecting dipping events. Lastly, we try to theoretically clarify the motivation behind our method. We start this paper with theory (section 2), including noise reduction based on Bayesian inversion framework (section 2.1), prior information (section 2.2), edge-preserving noise reduction (section 2.3) and numerical scheme (section 2.4). Then we use a synthetic data example and two field post-stack data sets (section 3) to test the performance of our method. Finally, we give some discussions (section 4) and conclusions (section 5). 2. Theory 2.1. Noise reduction based on Bayesian inversion framework As we know, the seismic data is often regarded as a superposition of coherent signals and random noise. So the 2D seismic data with m time samples and k traces can be expressed as 1 where matrices D, S and N denote observed seismic data, effective seismic signals and random noise, respectively. Because vector form is more straightforward and readily linked to directional difference operator, we change matrix operation into vector operation. By rearranging the elements of 2D matrices D, S and N, we rewrite equation (1) with 1D column vectors d, s and n as 2 where d = [D11, …, Dm1, D12, …, Dm2, …, D1k, …, Dmk]T, s = [S11, …, Sm1, S12, …, Sm2, …, S1k, …, Smk]T, n = [N11, …, Nm1, N12, …, Nm2, …, N1k, …, Nmk]T and the superscript T denotes transpose. According to Bayes' rule, the posterior probability of signals is derived by combining the likelihood of observed data p(d|s) and the prior of signals p(s) with the normalizing term p(d): 3 Assuming that n is approximated as uncorrelated zero-mean Gaussian noise with variance σ2, the likelihood function in equation (3) can be expressed as 4 where l = m × k is the total sample number of seismic data. By maximizing the above likelihood function (equation (4)), we get the optimum inverted result ⁠. In this case, the noise fully remains in the inverted result. Therefore, for the purpose of reducing noise, a priori information about signals must be added. In this section, it is assumed that the prior of signals is written as the following exponent function form 5 where c is a normalized constant, and ϕ(s) is a function of s. This prior form is usually utilized in seismic inverse problems (e.g. Sacchi and Ulrych 1995), and also presented as the framework of a Gibbs distribution (e.g. Geman and Geman 1984). By substituting equations (4) and (5) into equation (3), we have 6 After taking the negative logarithm for both sides of equation (6), we get 7 Since the logarithm function is monotonic, the maximization of p(s|d) is identical to the minimization of -log p(s|d). Thus, we wish to minimize the following objective function that establishes a trade-off between signal fidelity and noise reduction 8 with a positive parameter λ = 2σ2. The solution of minimizing equation (8) is the noise reduction result. Obviously, the performance of this method depends on how to set parameter λ and select function ϕ(s). 2.2. Prior information As seen in equation (8), a priori information about signals is necessary for noise reduction. However, how to obtain reasonable prior knowledge is difficult. In this section, we present a kind of non-Gaussian prior information about difference data of seismic signals. Suppose that seismic signals are composed of smooth signals, separated by edges. The difference data of smooth signals distribute near zero value, while those of edges distribute away from zero. If smooth signals are dominant in a seismic profile, the difference data of seismic signals approximately obey a non-Gaussian distribution. In our earlier work, TV has been successfully applied to noise reduction as a constraint (Yuan et al2012), where the prior knowledge we adopted is that TV data of seismic signals obey the Laplace distribution, a kind of non-Gaussian distribution. In this paper, we introduce another kind of non-Gaussian distribution about directional difference data of seismic signals, the Cauchy distribution, as a priori information. For convenience, we consider difference operation along directions 1, 2, 3 and 4 (figure 1), where 1, 2, 3 and 4 represent horizontal space direction, 45° space direction, time direction and 135° space direction, respectively. Along these four directions, we define four difference operators C1, C2, C3 and C4, respectively. Then we combine them into a large difference operator C = [α1C1, α2C2, α3C3, α4C4]T, where α1, α2, α3 and α4 are weights balancing direction constraints. In this paper, we only consider that αu (u = 1,…,4) is 0 or 1. If there are strong coherent signals along or near spatial directions 1, 2 or 4, we set the weight αu of this direction as 1. Otherwise, set αu as 0 and force the corresponding Cu to be an empty matrix. Moreover, we always set α3 = 1 in order to enhance the signals along the time direction and the coherent signals between the time direction and the chosen spatial direction(s). Figure 1. Open in new tabDownload slide Four difference directions, where 1, 2, 3 and 4 represent horizontal space direction, 45° space direction, time direction and 135° space direction, respectively. Figure 1. Open in new tabDownload slide Four difference directions, where 1, 2, 3 and 4 represent horizontal space direction, 45° space direction, time direction and 135° space direction, respectively. We use Cs to represent difference data of seismic signals. Assuming Cs can be described by the Cauchy distribution and every element in Cs is uncorrelated with every other, we have 9 where μ is the location parameter, specifying the location of the peak of the distribution, and γ is the scale parameter. Considering μ = 0, we can simplify equation (9) as 10 In fact, equation (10) is a specific form of equation (5) with a constant c = 1/(πγ)K and the Cauchy function ⁠. 2.3. Edge-preserving noise reduction In this section, we try to theoretically clarify motivation for edge-preserving noise reduction by using the Cauchy function as a priori information. By substituting the Cauchy function into equation (8), we obtain the objective function 11 By letting the first-order derivative of the function with respect to the unknown s be zero, we obtain a necessary condition for solving equation (11) as follows: 12 where A(s) = I + λB(s), is a sparse matrix, is a diagonal matrix with diagonal entries composed of vector and I is an identity matrix. We can minimize equation (11) and obtain the noise-reduction data by solving equation (12). However, equation (12) is a nonlinear problem, as matrix A contains unknown s. Therefore we transform this nonlinear problem into a series of linear inversion problems. For each linear problem, matrix A is required to be known. That means vector s in matrix A is required to be known. Assuming that the known matrix A is ⁠, we describe the linear problem as 13 where ⁠, and ⁠. The identity matrix I plays a key role in stably solving equation (13), as matrix B is a singular matrix. For convenience, we replace element in vector form with element in matrix form, where i ∈ {1, m}, j ∈ {1, k}. From equation (13), it can be derived that the sample point Di,j (i ∈ {1, m}, j ∈ {1, k}) is the linearly weighted sum of Si-1,j-1, Si-1,j, Si-1,j+1, Si,j-1, Si,j, Si, j+1, Si+1,j-1, Si+1,j and Si+1,j+1: 14 where and ⁠. For simplicity, we take αu = 0 (u = 2, 3 and 4) and only consider C1. Then in equation (13) can be replaced by ⁠. In this case, the sample point Di,j is a linearly weighted sum of points Si,j-1, Si,j and Si,j+1. The weights are the elements A(j-1)*m+i,(j-2)*m+i (or ⁠), A(j-1)*m+i,(j-1)*m+i (or ⁠) and A(j-1)*m+i,(j+1)*m+i (or 1 - μ) in matrix A, respectively. These weights depend on coefficients B(j-1)*m+i,(j-2)*m+i (or ⁠), B(j-1)*m+i,(j-1)*m+i (or ⁠) and B(j-1)*m+i,(j+1)*m+i (or -μ) in matrix B, which are related to the features of seismic data. When there is large amplitude variation at sample point ⁠, and are large, thus coefficients B(j-1)*m+i,(j-2)*m+i, B(j-1)*m+i,(j-1)*m+i and B(j-1)*m+i,(j+1)*m+i are small. Supposing that and are infinity, coefficients B(j-1)*m+i,(j-2)*m+i, B(j-1)*m+i,(j-1)*m+i and B(j-1)*m+i,(j+1)*m+i are close to 0, thus weights A(j-1)*m+i,(j-2)*m+i, A(j-1)*m+i,(j-1)*m+i and A(j-1)*m+i,(j+1)*m+i are close to 0, 1 and 0, respectively. Therefore, Si,j approximately equals Di,j. The large amplitude variation is naturally preserved. When there is a discontinuity between sample point and sample point ⁠, weight μF1 is small. When is infinity, weight μF1 vanishes and there is no smoothing in this direction, thus preserving discontinuity in this direction. The same conclusion can be obtained when all directional differences are considered. In this way, discontinuities or large amplitude variations, which often correspond to important geology features, such as faults, pinch-outs, channels, lens or fractures, can be adaptively preserved. 2.4. The numerical scheme In section 2.3, we point out that the solution of equation (12) can be obtained by solving a series of linear problems. In this section, an iteratively reweighted least-squares algorithm (Scales and Gersztenkorn 1988) is adopted to solve the equation. Firstly, we set parameters αu (u = 1,…,4), λ as well as γ, and choose an starting model s0 = d. Then we solve equation (13) by using s0 as and obtain solution s1. After that, using s1 as to solve equation (13) again, we obtain solution s2. Like this, we iteratively solve equation (13) until convergence condition or maximum number of iterations is achieved. The final iterative solution is the noise-reduction result. The whole algorithm flow is summarized as follows: Fix weights αu (u = 1, …, 4), tolerance error ε and maximum number of iterations itmax, and set the trade-off parameter λ and the scale parameter γ. Initialize s0 = d. For the gth iteration, solve sg using equation (13). Calculate ‖sg - sg - 1‖2. If ‖sg - sg - 1‖2 > ε and u < itmax, go to step 3; otherwise, go to step 5. Output noise-reduction result. The performance of the method is affected by parameters αu (u = 1, …, 4), ε, itmax, λ and γ. We have introduced how to choose parameter αu in section 2.2. In general, ε and itmax are readily determined. In this paper, we set ε = 10-6 × ||d||2, and itmax = 20. As for λ and γ, they play an important role in balancing noise reduction and signal preservation. The larger λ or the smaller γ, the better the noise attenuation, but the worse the signal fidelity and vice versa. Therefore, we should choose large λ and large γ, or small λ and small γ, in order to achieve a good trade-off between noise attenuation and signal fidelity. 3. Examples The method is applied to a synthetic data set and two field post-stack seismic data sets to test its performance. We also illustrate the influence of trade-off parameter λ and scale parameter γ, as well as different directional constraints, on the noise-reduction result. 3.1. Synthetic data We make a noise-free synthetic seismic data set (figure 2(a)), consisting of 50 traces with a sample interval of 1 ms. This profile contains two curved events with different curvatures (A), a quasi-linear discontinuous event (B), an isolated event (C), a broken event (D), a horizontal discontinuous event (E), two conflicting events (F), a event with peak frequency of wavelet changing (G), a event with amplitude variations (H), and a event with phase of wavelet changing (I). In general, the curved events and conflicting events are mainly observed in pre-stack data set, and the other events we designed are always observed in post-stack data set. The main aim for designing the data set is to illustrate the effectiveness of the method for enhancing events with different features. Figure 2. Open in new tabDownload slide The synthetic data. (a) Noise-free synthetic data and (b) noisy data with 50% Gaussian noise. Figure 2. Open in new tabDownload slide The synthetic data. (a) Noise-free synthetic data and (b) noisy data with 50% Gaussian noise. By adding 50% (the noise-to-signal ratio N/S = 0.5, defined as ||n||22/||s||22) white Gaussian noise to the noise-free profile, a noisy profile (figure 2(b)) is generated. To quantitatively evaluate the noise-reduction result, we define a relative error as an evaluation criterion. By taking directions 1 and 3 constraints to control the noise attenuation, we try a large number of values from 0 to 1 with step 0.01 for parameter pair (λ, γ), in order to test their influence on our method and find the optimum result. Figure 3 is the relative error Q versus λ and γ. When λ = 0.06 and γ = 0.17, we can get the optimum noise-reduction result with minimum Q 8.3%. The noise-reduction result and difference profile are shown in figures 4(a) and (b), respectively. As the two subfigures show, the noise is reduced adequately while the signals are preserved well, especially for the discontinuities and amplitude variations, such as B, C, D, E, G and H. For comparison, we also show the results (figures 4(c)–(f)) using two other groups of parameters (0.8, 0.5) and (0.1, 0.7). When (λ, γ) is (0.8, 0.5), we obtain a very clean profile (figure 4(c)), but the signals are strongly affected (figure 4(d)), especially for the two curved events (A) and the dipping event (F). When (λ, γ) is (0.1, 0.7), we obtain a good signal-preserving profile (figure 4(e)), but less noise is reduced (figure 4(f)). Figure 3. Open in new tabDownload slide The relative error Q variations versus λ and γ when directions 1 and 3 are used as constraints. The red asterisk (λ = 0.06, γ = 0.17) corresponds to the minimum Q (8.3%). Q denoted by the two blue circles (0.8,0.5) and (0.1,0.7) are 19.6% and 21.5%, respectively. The big oval indicates the good trade-off region between signal fidelity and noise reduction. Figure 3. Open in new tabDownload slide The relative error Q variations versus λ and γ when directions 1 and 3 are used as constraints. The red asterisk (λ = 0.06, γ = 0.17) corresponds to the minimum Q (8.3%). Q denoted by the two blue circles (0.8,0.5) and (0.1,0.7) are 19.6% and 21.5%, respectively. The big oval indicates the good trade-off region between signal fidelity and noise reduction. Figure 4. Open in new tabDownload slide The noise-reduction results and difference profiles by using directions 1 and 3 as constraints with different parameter pairs (λ, γ). (a) Noise-reduction result when (λ, γ) takes (0.06, 0.17), (b) difference profile between figures 2(b) and (a), (c) noise-reduction result when (λ, γ) takes (0.8, 0.5), (d) difference profile between figures 2(b) and (c), (e) noise-reduction result when (λ, γ) takes (0.1, 0.7), and (f) difference profile between figures 2(b) and (e). Figure 4. Open in new tabDownload slide The noise-reduction results and difference profiles by using directions 1 and 3 as constraints with different parameter pairs (λ, γ). (a) Noise-reduction result when (λ, γ) takes (0.06, 0.17), (b) difference profile between figures 2(b) and (a), (c) noise-reduction result when (λ, γ) takes (0.8, 0.5), (d) difference profile between figures 2(b) and (c), (e) noise-reduction result when (λ, γ) takes (0.1, 0.7), and (f) difference profile between figures 2(b) and (e). In addition, we compare our method (figures 4(a) and (b)) with f–x deconvolution (figures 5(a) and (b)) and conventional EPS filtering (figures 5(c) and (d)). For f–x deconvolution, the noise reduction result is the average of forward and backward prediction filtering with the filter operator length 4. For EPS filtering, we use a 5 × 5 window with nine panels: four hexagonal, four pentagonal and one square. As figures 4(a), (b) and 5 show, our method is better at protecting edges than f–x deconvolution, and better at protecting signals than EPS filtering with lower Q. Although none of the methods accurately retrieve curved events or quasi-linear events, our method has the least signal loss. Figure 5. Open in new tabDownload slide The results obtained by using other noise reduction methods. (a) Noise-reduction result of f–x deconvolution, (b) difference profile between figures 2(b) and (a), (c) noise-reduction result of conventional EPS filtering, and (d) difference profile between figures 2(b) and (c). The relative errors for f–x deconvolution and EPS are 16.0% and 25.1%, respectively. Figure 5. Open in new tabDownload slide The results obtained by using other noise reduction methods. (a) Noise-reduction result of f–x deconvolution, (b) difference profile between figures 2(b) and (a), (c) noise-reduction result of conventional EPS filtering, and (d) difference profile between figures 2(b) and (c). The relative errors for f–x deconvolution and EPS are 16.0% and 25.1%, respectively. 3.2. Field data To test the performance of our method for retrieving real seismic signals, we choose a ‘clean’ profile (figure 6(a)), consisting of 160 traces with a sample interval of 1 ms, which has been processed using f–x deconvolution. By adding 50% white Gaussian noise, we get a noisy profile (figure 6(b)). We use different parameters and direction constraints to conduct noise attenuation. Since there is no obvious visual difference among the noise-reduction profiles, we only show the difference profiles. Figures 7(a)–(c) are the difference profiles by using single direction constraint (directions 1, 2 and 3) respectively, with their optimum parameter pairs (λ, γ). Figure 7(d) is the difference profile by simultaneously using directions 1, 2 and 3 as constraints with its optimum parameter pair. As figures 7(a)–(d) show, the noise-reduction result via direction 3 constraint has the maximum signal loss. The result via direction 1 constraint preserves horizontal or nearly horizontal events well, while that via direction 2 constraint preserves dipping events along or near direction 2 well. Compared with the results via single directional constraint, the noise-reduction result via three directional constraints has the best signal preservation as well as noise reduction. However, slight indications of coherent events also remain in the difference profile, mainly because the noise-reduction result is just the optimum trade-off between noise attenuation and signal fidelity. Figure 7(e) is the difference profile by using conventional EPS filtering. From this figure, it can be observed that conventional EPS preserves edges well, but with the lowest signal preservation. Figure 6. Open in new tabDownload slide Real post-stack data. (a) Profile without numerical noise, and (b) profile with 50% Gaussian noise. Figure 6. Open in new tabDownload slide Real post-stack data. (a) Profile without numerical noise, and (b) profile with 50% Gaussian noise. Figure 7. Open in new tabDownload slide The difference profiles using different schemes with their optimum parameter pairs. (a) Result by using direction 1 constraint, (b) result by using direction 2 constraint, (c) result by using direction 3 constraint, (d) result by using directions 1, 2 and 3 constraints, and (e) result of conventional EPS. The corresponding relative errors are 17.7%, 17.3%, 18.5%, 12.3% and 34.6%, respectively. Figure 7. Open in new tabDownload slide The difference profiles using different schemes with their optimum parameter pairs. (a) Result by using direction 1 constraint, (b) result by using direction 2 constraint, (c) result by using direction 3 constraint, (d) result by using directions 1, 2 and 3 constraints, and (e) result of conventional EPS. The corresponding relative errors are 17.7%, 17.3%, 18.5%, 12.3% and 34.6%, respectively. Figure 8 is the relative error Q versus λ and γ when directions 1, 2 and 3 are simultaneously used as constraints. Similar to figure 3, it also shows a good trade-off region (denoted by the big oval) between signal fidelity and noise reduction. For example, we can choose small λ and small γ, or large λ and large γ. Figure 8. Open in new tabDownload slide The relative error Q variations versus λ and γ when directions 1, 2 and 3 are used as constraints. The red asterisk (λ = 0.06, γ = 0.34) corresponds to the minimum Q (12.3%). The big oval indicates the good trade-off region between signal fidelity and noise reduction. Figure 8. Open in new tabDownload slide The relative error Q variations versus λ and γ when directions 1, 2 and 3 are used as constraints. The red asterisk (λ = 0.06, γ = 0.34) corresponds to the minimum Q (12.3%). The big oval indicates the good trade-off region between signal fidelity and noise reduction. Another field post-stack seismic data set (figure 9(a)), consisting of 200 traces with a sample interval of 2 ms, is adopted to test the application potential of our method. The data set mainly includes discontinuous events, dipping events, events with amplitude and phase variations, quasi-linear events and background noise. It is obvious that the signals, especially for the discontinuities and lateral variations of the events, are blurred by random noise. Figure 9. Open in new tabDownload slide Field post-stack seismic data example. (a) Original seismic data, (b) noise-reduction result, and (c) difference profile between figures 9(a) and (b). Figure 9. Open in new tabDownload slide Field post-stack seismic data example. (a) Original seismic data, (b) noise-reduction result, and (c) difference profile between figures 9(a) and (b). We choose directions 1, 3 and 4 as constraints to implement inversion, as the predominant directions of seismic signals are near directions 1 and 4. Based on the experience from the above two examples (figures 3 and 8), we choose the parameters as λ = 0.4, γ = 0.7 by several trials. Figure 9(b) is the noise-reduction result and figure 9(c) is the difference profile. As figure 9(b) shows, almost all signals are enhanced. The edges of discontinuous events, as well as the amplitude and phase variations of events, become clearer. Moreover, dipping events are perfectly retrieved. In the difference profile (figure 9(c)), there are no obvious indications of coherent events removed. 4. Discussions By using Cauchy distribution of directional difference data of seismic signals as a priori information, we propose a noise reduction method based on Bayesian inversion framework. The idea is well used for seismic noise reduction, even though non-smooth regularizations or constraints via non-Gaussian distribution have been widely presented for acoustic impedance inversion (e.g. Oldenburg et al1983), sparse deconvolution (e.g. Sacchi 1997), seismic data interpolation (e.g. Herrmann and Hennenfent 2008), seismic imaging (e.g. Herrmann and Li 2012) and full-waveform inversion (e.g. Guitton 2012). Besides the Cauchy distribution we utilized in this paper, we do not deny that other non-Gaussian distributions with heavy tails will be appropriate. Based on our constraint or a priori information (equation (11)), every sample point of seismic signals is implicitly only related to its first-order adjacent points, but the relationship can be passed to any other sample points (equation (13)). It is demonstrated by the fact that when a certain spatial difference constraint is used, the signals along the direction can be perfectly preserved. The performance of the method mainly depends on parameters λ and γ, as well as weights of directional constraints. We define a relative error as the criterion to get some experience for choosing appropriate parameters λ and γ by a synthetic data example and the first field post-stack data example. In a word, by choosing small λ with small γ, or large λ with large γ, we can obtain a good trade-off result between signal fidelity and noise reduction. If we visually observe that seismic signals are seriously harmed, we guess parameter λ is set too large and parameter γ is set too small. If we visually observe that too much noise still remains in the noise reduction profile, we guess parameter λ is set too small and parameter γ is set too large. As for the weights of directional constraints, they are simply determined as 0 or 1 mainly according to the energy distribution of coherent signals along spatial directions 1, 2 and 4. In general, if there are strong coherent signals along or near a certain spatial direction, the weight αu of this direction is set as 1. Moreover, α3 is fixed as 1 in order to enhance the signals along the time direction and the coherent signals between the time direction and the chosen spatial direction(s). Although a series of seismic dataset examples in this paper are used to illustrate the effectiveness of multiple direction constraints with weights 0 or 1, further work about the choice of weights and more directional constraints are worth studying in the future. The method can preserve high-amplitude spike noise, but it does not break down the preservation of signals and does not lead to spurious events. The method can also be directly used to simultaneously reduce spike noise by replacing the first term in equation (11) with functions such as L1 norm (Koko and Jehan-Besson 2010) or Huber norm (Guitton and Symes 2003). Furthermore, it is easy to be extended to 3D case by changing the difference operator C in equation (11), as an additive 3D synthetic data example (figure 10) shows. The potential 3D method strongly enhances effective signals, especially for perfectly recovering edges of 3D discontinuous events, with the relative error Q 1.2%. Figure 10. Open in new tabDownload slide The 3D synthetic data example with a sample interval of 1 ms. (a) Noise-free synthetic data, (b) noisy data with 30% Gaussian noise, and (c) 3D noise-reduction result by using inline direction, crossline direction and time direction as constraints with parameter pairs (0.04, 0.14). Figure 10. Open in new tabDownload slide The 3D synthetic data example with a sample interval of 1 ms. (a) Noise-free synthetic data, (b) noisy data with 30% Gaussian noise, and (c) 3D noise-reduction result by using inline direction, crossline direction and time direction as constraints with parameter pairs (0.04, 0.14). 5. Conclusions Taking features of seismic data into account, we present a novel random noise reduction method with directional difference constraints. The method can adaptively preserve edges of discontinuous events. It also allows us to achieve a good trade-off between signal fidelity and noise reduction by properly selecting directional difference operators as well as parameters λ and γ. Moreover, the method can directly preserve dipping events via directional difference constraints, without a dip steering step. Acknowledgments We are grateful to two anonymous reviewers for their constructive comments on this paper. We also would like to thank Felix J Herrmann for providing some good suggestions and thank the Seismic Laboratory for Imaging and Modelling for providing us parallel computation platform. This research was supported by National 973 program (2013CB228600), Program for Changjiang Scholars, Innovative Research Team (PCSIRT) in University, and China Postdoctoral Science Foundation Funded Project (2012M510771). References Bekara M , van der Baan M . , 2007 Local singular value decomposition for signal enhancement of seismic data , Geophysics , vol. 72 (pg. V59 - V65 ) 10.1190/1.2435967 Google Scholar Crossref Search ADS WorldCat Crossref Bonar D , Sacchi M D . , 2012 Denoising seismic data using the nonlocal means algorithm , Geophysics , vol. 77 (pg. A5 - A8 ) 10.1190/geo2011-0235.1 Google Scholar Crossref Search ADS WorldCat Crossref Canales L L . , 1984 Random noise reduction 54th Annu. Int. Meeting, SEG, Expanded Abstracts (pg. 525 - 527 ) Geman S , Geman D . , 1984 Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images , IEEE Trans. Pattern Anal. Mach. Intell. , vol. 6 (pg. 721 - 741 ) 10.1109/TPAMI.1984.4767596 Google Scholar Crossref Search ADS PubMed WorldCat Crossref Guitton A . , 2012 Blocky regularization schemes for full-waveform inversion , Geophys. Prospect. , vol. 60 (pg. 1 - 15 ) 10.1111/j.1365-2478.2011.00957.x Google Scholar Crossref Search ADS WorldCat Crossref Guitton A , Symes W W . , 2003 Robust inversion of seismic data using the Huber norm , Geophysics , vol. 68 (pg. 1310 - 1319 ) 10.1190/1.1598124 Google Scholar Crossref Search ADS WorldCat Crossref Herrmann F J , Hennenfent G . , 2008 Non-parametric seismic data recovery with curvelet frames , Geophys. J. Int. , vol. 173 (pg. 233 - 248 ) 10.1111/j.1365-246X.2007.03698.x Google Scholar Crossref Search ADS WorldCat Crossref Herrmann F J , Li X . , 2012 Efficient least-squares imaging with sparsity promotion and compressive sensing , Geophys. Prospect. , vol. 60 (pg. 697 - 712 ) 10.1111/j.1365-2478.2011.01041.x Google Scholar Crossref Search ADS WorldCat Crossref Koko J , Jehan-Besson S . , 2010 An augmented Lagrangian method for TVg+L1-norm minimization , J. Math. Imaging Vis. , vol. 38 (pg. 182 - 196 ) 10.1007/s10851-010-0219-1 Google Scholar Crossref Search ADS WorldCat Crossref Liu Y K , Luo Y , Deng Z Y . , 2003 Edge preserving smoothing for oblique images , Geophys. Res. Lett. , vol. 30 pg. 1699 10.1029/2003GL017399 OpenURL Placeholder Text WorldCat Crossref Lu W K . , 2006 Adaptive noise attenuation of seismic images based on singular value decomposition and texture direction detection , J. Geophys. Eng. , vol. 3 (pg. 28 - 34 ) 10.1088/1742-2132/3/1/004 Google Scholar Crossref Search ADS WorldCat Crossref Lu Y , Lu W K . , 2009 Edge-preserving polynomial fitting method to suppress random seismic noise , Geophysics , vol. 74 (pg. V69 - V73 ) 10.1190/1.3129907 Google Scholar Crossref Search ADS WorldCat Crossref Luo Y , Marhoon M , Al-Dossary S , Alfaraj M . , 2002 Edge-preserving smoothing and applications , Leading Edge , vol. 21 (pg. 136 - 158 ) 10.1190/1.1452603 Google Scholar Crossref Search ADS WorldCat Crossref Oldenburg D W , Scheuer T , Levy S . , 1983 Recovery of the acoustic impedance from reflection seismograms , Geophysics , vol. 48 (pg. 1318 - 1337 ) 10.1190/1.1441413 Google Scholar Crossref Search ADS WorldCat Crossref Oropeza V , Sacchi M . , 2011 Simultaneous seismic data denoising and reconstruction via multichannel singular spectrum analysis , Geophysics , vol. 76 (pg. V25 - V32 ) 10.1190/1.3552706 Google Scholar Crossref Search ADS WorldCat Crossref Sacchi M D . , 1997 Reweighting strategies in seismic deconvolution , Geophys. J. Int. , vol. 129 (pg. 651 - 656 ) 10.1111/j.1365-246X.1997.tb04500.x Google Scholar Crossref Search ADS WorldCat Crossref Sacchi M D , Kuehl H . , 2001 ARMA formulation of FX prediction error filters and projection filters , J. Seismic Explor. , vol. 9 (pg. 185 - 197 ) OpenURL Placeholder Text WorldCat Sacchi M D , Ulrych T J . , 1995 High-resolution velocity gathers and offset space reconstruction , Geophysics , vol. 60 (pg. 1169 - 1177 ) 10.1190/1.1443845 Google Scholar Crossref Search ADS WorldCat Crossref Scales J A , Gersztenkorn A . , 1988 Robust methods in inverse theory , Inverse Problems , vol. 4 (pg. 1071 - 1091 ) 10.1088/0266-5611/4/4/010 Google Scholar Crossref Search ADS WorldCat Crossref Soubaras R . , 1994 Signal-preserving random noise attenuation by the f–x projection 64th Annu. Int. Meeting, SEG, Expanded Abstracts (pg. 1576 - 1579 ) Trickett S R . , 2008 F–xy Cadzow noise suppression CSPG CSEG CWLS Convention (pg. 303 - 306 ) Yuan S Y , Wang S X . , 2011 A local f–x Cadzow method for noise reduction of seismic data obtained in complex formations , Petrol. Sci. , vol. 8 (pg. 269 - 277 ) 10.1007/s12182-011-0144-y Google Scholar Crossref Search ADS WorldCat Crossref Yuan S Y , Wang S X , Li G F . , 2012 Random noise reduction using Bayesian inversion , J. Geophys. Eng. , vol. 9 (pg. 60 - 68 ) 10.1088/1742-2132/9/1/007 Google Scholar Crossref Search ADS WorldCat Crossref © 2013 Sinopec Geophysical Research Institute
TI - Edge-preserving noise reduction based on Bayesian inversion with directional difference constraints
JO - Journal of Geophysics and Engineering
DO - 10.1088/1742-2132/10/2/025001
DA - 2013-04-01
UR - https://www.deepdyve.com/lp/oxford-university-press/edge-preserving-noise-reduction-based-on-bayesian-inversion-with-Ev6D6WbXv0
VL - 10
IS - 2
DP - DeepDyve
ER -