Application of curvelet denoising to 2D and 3D seismic data — Practical considerations

    Application of curvelet denoising to 2D and 3D seismic data – practical considerations Andrzej G´orszczyk, Anna Adamczyk, Michał Malinowski PII: DOI: Reference:

S0926-9851(14)00078-0 doi: 10.1016/j.jappgeo.2014.03.009 APPGEO 2455

To appear in:

Journal of Applied Geophysics

Received date: Revised date: Accepted date:

12 December 2013 3 March 2014 5 March 2014

Please cite this article as: G´orszczyk, Andrzej, Adamczyk, Anna, Malinowski, Michal, Application of curvelet denoising to 2D and 3D seismic data – practical considerations, Journal of Applied Geophysics (2014), doi: 10.1016/j.jappgeo.2014.03.009

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT Application of curvelet denoising to 2D and 3D seismic data – practical

SC

Andrzej Górszczyk1*, Anna Adamczyk1, Michał Malinowski1

RI

PT

considerations

TE

D

MA

NU

1) Institute of Geophysics, Polish Academy of Sciences, Ks. Janusza 64, 01-452 Warsaw, Poland

Abstract

AC CE P

Contamination of seismic signal with noise of various origins is one of the main challenges encountered during processing and interpretation of seismic data. Several methods exist for eliminating different types of noises like coherent or incoherent noise, multiples etc., but optimal random noise attenuation remains difficult. Here we investigate relatively new technique based on Discrete Curvelet Transform (DCT). Features like multi-resolution, multi-direction and locality of DCT introduce minimal overlapping between coefficients representing signal and noise in curvelet domain which is prime advantage of this method. We present practical application of DCT describing its main features and focusing on useful details, especially more complex thresholding based on analyzing 2D Fourier spectrum and the vector of curvelet coefficients. We demonstrate that better understanding of relations between DCT properties and obtained results in pair with *

Corresponding author: [email protected]

ACCEPTED MANUSCRIPT additional investigation of curvelet domain provides better localization and, in consequence, separation of noise and signal energy. Introduced scale and angle dependent weighting of curvelet

PT

coefficients leads to significant improvements of results with respect to noise attenuation and

RI

signal energy preservation. Effectiveness of our approach is demonstrated both on synthetic 2D

SC

sections with white and colored noise added, as well as on real 2D and 3D post-stack seismic data. Finally, we demonstrate the use of curvelet denoising as the data-preconditioning tool for

NU

frequency-domain full-waveform inversion. Curvelet denoising seems to be much more robust as

MA

compared with traditional filtering (e.g. F-X deconvolution), especially when noise and signal

AC CE P

1. Introduction

TE

D

spectra overlap.

In seismic data processing noise is often understood as the recorded information that is not used to obtain seismic images or models. It interferes with the signal, i.e. the useful part of recorded data, damaging or concealing it yielding difficulties in imaging and leading to misinterpretation. The need to enhance the desired, most informative events in seismic data motivates the researchers to develop increasingly sophisticated denoising techniques. Some of these methods deal with the coherent events, exhibiting phase consistency from trace to trace, some are more often used to remove random noise. Designing a combination of those techniques adequate to our data and the further processing scheme is crucial for successful seismic processing and imaging.

ACCEPTED MANUSCRIPT For direct random noise attenuation, F-X deconvolution (Canales, 1984) based on prediction of linear events in the frequency-space domain is used. Methods developed for coherent noise

PT

removal, such as F-K filtering (Stewart and Schieck, 1989), median filtering (Bednar, 1983), linear

RI

Radon or tau–p transform (Kappus et al., 1990; Turner, 1990; Zhou and Greenhalgh, 1994) can also

SC

be helpful for attenuating random noise by localizing and boosting coherent events. Above listed methods are the most commonly used procedures in seismic processing. A common feature of all

NU

those traditional noise attenuation techniques is that they operate on a single scale, i.e. use signal

MA

or image representation without its segmentation or decomposition into multiresolution elements. However, in recent years multiscale methods have been gaining researchers’ interest

D

with algorithms based on curvelet transform among them. Curvelet transform was first introduced

TE

by Candès and Donoho (2000) as a new multiscale and sparse representation of signals suited for

AC CE P

mapping edge-like objects: smooth along a curve and sharply discontinuous in the normal direction. Since that time curvelets have evolved and found various fields for application. The essence of curvelet-based denoising in the image processing was described by Starck et al. (2002). In the context of seismic imaging, Herrmann et al. (2008) show the application of curvelets to reconstruct data from incomplete measurements, to separate primaries and multiples and to restore migration amplitudes. Neelamani et al. (2008) introduced curvelet-based noise attenuation for 3D seismic data corrupted with random and linear noise. More seismic issues like ground roll attenuation (Yarham and Herrmann, 2008) or enhancing crustal reflection data with sparsity promotion (Kumar et al., 2011) were solved with curvelet-based methods proving their value for seismic processing.

ACCEPTED MANUSCRIPT In this paper we describe our experience with practical noise attenuation by applying curvelet transform to post-stack seismic data in 2D and 3D shots in various acquisition configurations and

PT

processed in a different way. We believe that better understanding of the nature of curvelets and

RI

clarifying the advantages of this noise attenuation method may lead to its popularization and

SC

subsequent development. Experience we gained so far shows that practical application of more complex thresholding leads to better noise attenuation. Such approach is supported by presented

NU

examples. We will start with short description of curvelet transform to show its main features and

MA

introduce the subject, since for many geophysicists this topic may be relatively new. Then we will describe the flow of denoising algorithm in details, supported by demonstrative synthetic

D

examples. Subsequently the results illustrating application of the denoising scheme to synthetic 2D

TE

data as well as 2D and 3D real post-stack data will be presented. We also show the application of

AC CE P

curvelets to condition the data that are subjected to frequency-domain full-waveform inversion. Finally we will conclude with some practical implications and comment on the features of applied method.

2. Curvelet transform

Curvelet transform was first introduced as a mathematical transform capable of representing edges in image processing problems. It was particularly an answer for the difficulties that wavelets revealed when used to represent higher dimension linear singularities. First generation was the ridgelet transform (Candès and Donoho, 1999) applied to blocks of data defined with 2D spatial coordinates and frequency. This implementation, however, had limited application since the geometry of ridgelets is unclear, as they are not true ridge functions in digital images. This led the

ACCEPTED MANUSCRIPT authors to redesign their own construction and make it simpler, faster and less redundant. As a result they substituted block ridgelet transform by direct 2D frequency domain partitioning and

PT

created two digital implementations of, so called, Second Generation Curvelet Transform: Digital Curvelet Transform (DCT) via Unequispaced Fast Fourier Transform,



Digital Curvelet Transform via Wrapping.

SC

RI



each scale and angle (Candès et al., 2006).

NU

Main difference between them comes from the choice of spatial grid used to translate curvelets at

MA

It was just a matter of time when curvelet transform found its application in seismic imaging where recorded signal usually have the form of continuous smooth lines oscillating in the direction

D

in which the wavefront moves. In stacked seismic sections the reflections are represented as linear

TE

objects of similar nature, oscillating in the direction perpendicular to their spatial extension. Even

AC CE P

if subsurface structure is complicated with conflicting dips, intersecting curved events, faults etc., the curvelet transform is still able to decompose the data into linear sum of needle-shaped atoms named curvelets. They are multiscale, multidirectional and local in both spatial and frequency domain as presented in Fig. 1a and Fig. 1b. Moreover, in Fourier domain they occupy strictly localized, directional wedges (Fig. 1c) what is reflected by the fact that they obey the “parabolic scaling” law, which reveal that for each of angular polygons, if its length dimension is proportional to 2j, where j is an integer, then its width dimension is approximately proportional to 2j/2. Those features require increased (in comparison with e.g. Fourier transform) parameterization. In this case each single curvelet is described by characteristic frequency (scale), dip (angle) and two spatial coordinates in the manner shown in Fig. 1d. Such parameterization allows to map particular smoothly continuous events such as wavefronts or reflectors into separate or in practice

ACCEPTED MANUSCRIPT almost separate sets of coefficients, which makes curvelet domain favorable in regard to seismic data processing. Another important feature of the curvelet transform is the sparsity of data

PT

representation. The less components required to represent the signal, the more successful the

RI

transformation. By thresholding the curvelet coefficients, one is able to achieve high level of

SC

sparsity, which might be an advantage in further, easier and faster processing, especially when

NU

dealing with large data sets (i.e. Herrmann et al., 2008).

MA

3. Implementation

D

3.1 General description of the curvelet transform

TE

The issue of noise attenuation stimulates researchers to elaborate new or modify already known

AC CE P

filtrating methods. Soon after curvelets were introduced an approach of image restoration combining curvelets and wavelets was presented (Starck et al., 2001). Kumar et al. (2011) proposed the “cooling method“ which promotes sparsity and involves hard and iterative soft thresholding of curvelet coefficients. We decided to use the simple approach and follow a scheme that is known and widely applied in image processing, namely: 

Forward discrete curvelet transform,



Analysis of obtained curvelet coefficients,



Hard thresholding,



Inverse discrete curvelet transform.

We focus on careful choice of the curvelet coefficients that need to be suppressed in order to enhance the signal. The thresholding strategy includes adapting the parameters according to the

ACCEPTED MANUSCRIPT scale and angle of the data in curvelet domain and requires additional iterations for their adjustment. The whole algorithm is not a closed box and it gives wide independance in choosing

PT

the thresholding parameters, which is an advantage when working with data sets of varying

SC

guidelines we follow are presented later in this article.

RI

quality, data and noise characteristics, etc. Detailed description of this procedure and the

NU

We make use of two assumptions frequently used in noise attenuation elaborations (e.g.

MA

Hennenfent et al. (2011)). The first one is that the recorded seismic data ( ) can be represented as superposition of signal ( ) and noise ( ):

D

(1)

TE

The second corresponds to the fact that in curvelet domain signal and noise maps into

where

AC CE P

approximately disjoint collections of coefficients. Hence after we can write: (2)

denote nonzero vectors of curvelet coefficients.

We would now proceed through some details that might explain and clarify presented scheme. As mentioned in previous chapter there are two implementation of Fast Discrete Curvelet Transform (FDCT). Source codes of both are available from CurveLab project (www.curvelet.org). Since the main focus of this paper is not to present FDCT itself but focus on its practical applications, therefore the reader is referred to Candès et al. (2006) for the omitted details concerning DCT.

ACCEPTED MANUSCRIPT For our needs wrapping-based forward and inverse transform is superior because of its simpler implementation and therefore much faster execution especially when dealing with large portions

PT

of data. As mentioned, curvelet transform has its origins in image processing and thus the input

RI

data is in a form of an image, i.e. seismic traces are treated as equispaced columns of samples (or

SC

pixels). The irregularities in the acquisition have to be dealt with beforehand (resampling of the

four steps:

MA

1. 2D Fast Fourier Transform (FFT),

NU

data, interpolation or filling the empty traces with zeros). In general FDCT architecture consists of

2. Forming a product of scale and angle windows,

D

3. Wrapping this product around origin,

TE

4. Collecting discrete coefficients via inverse FFT 2D.

AC CE P

Before we perform the forward transform we have to choose the number of scales and whether we are going to use wavelets or curvelets for representing the features at finest scale. The appropriate number of scales depends on the characteristics of both signal and noise, e.g. if we deal with high frequency events a higher number of scales might be required in order to correctly map them into the curvelet domain. Random white noise which in curvelet domain is represented by a large number of relatively small coefficients is in most situations easy to remove without much deliberating. However, if the noise is of certain characteristic or coherent and maps into curvelet coefficients of significant values, then densely partitioned frequency spectrum is worth considering, because it leads to better separation and in consequence attenuation of the noise.

ACCEPTED MANUSCRIPT 3.2 Curvelets scale selection

PT

In order to demonstrate influence of scale partitioning on signal and noise separation in curvelet

RI

domain let us consider following situation. Zero-centered white Gaussian noise filtered between

SC

15 and 30 Hz was applied to a synthetic noise-free seismic section. During DCT, 2D frequency spectrum was partitioned in two ways with 4 and 7 scales as demonstrated in Fig. 2a and Fig. 2b.

NU

The following two plots (Fig. 2c and Fig. 2d) present vectors of curvelet coefficients obtained from

MA

noisy section corresponding to each of fragmentation. First vector contains curvelet coefficients (from scales 1-2) received after DCT with sparsely partitioned 2D frequency spectrum. Second

D

vector corresponds to densely partitioned frequency domain (coefficients from scales 1-5 and

TE

partially 6). As it can be clearly seen, increasing number of scales led to extracting coherent energy

AC CE P

in a form of separate, higher coefficients collections at scales 4 and 5, whose energy was primarily mixed with noise and indistinguishable. Usefulness of this observation is explained later in section 3.3.

On the other hand, densely partitioned spectrum may cause certain undesirable effects. If the signal is dominated by incoherent noise it may produce significant artifacts in a form of high frequency coefficients of relatively high value. They result from fitting the finest scale curvelets to the incoherent noise. Furthermore, the computational cost of forward and inverse curvelet transforms increase with the number of scales.

FDCT via wrapping allows to choose between curvelets and wavelets to represent features in the finest scale. To illustrate consequence of employing each of the options we will use another

ACCEPTED MANUSCRIPT synthetic example. In this case zero-centered white Gaussian noise was added to a synthetic seismic section representing part of the Marmousi2 model (Martin et al., 2006) (Fig. 3a – noise

PT

free data, Fig. 3b – data with added noise) and FDCT with wavelets and curvelets at finest scale

RI

were proceeded. Sections presented in Fig. 3c and 3d were recovered using FDCT with wavelets

SC

and curvelets at the finest scale. In both cases global thresholding was applied after estimation of optimal threshold levels which were

for wavelets and

for curvelets.

NU

Selected values allowed 3% and 0.2% of curvelet coefficients to enter inverse DCT respectively.

MA

The results, though both far from ideal, suggest advantage of using curvelets in the finest scale:

D

the image is clearer and the weaker events are easier to follow.

TE

Another aspect of the choice between curvelets and wavelets at the final scale is redundancy of

AC CE P

each approach. Proportion of obtained coefficients varies from about 2.8 in case of wavelet to 7.2 with curvelets at finest scale (Candès et al., 2006). In our works however, despite higher computational requirements, we prefer to use curvelets at finest scales as they provide superior results especially when recovered section contains high-frequency events.

3.3 Scale-adjusted thresholding

All those proprieties of FDCT may significantly help to minimize the overlapping of noise and signal coefficients in curvelet domain. However, correct configuration of the scales is in the background for the most important step of the algorithm which is weighting of the curvelet

ACCEPTED MANUSCRIPT coefficients. During hard thresholding, coefficients of values higher than the selected threshold level are preserved, while smaller coefficients are zeroed.

PT

Therefore picking appropriate threshold values has direct influence on the results’ quality. It

RI

follows from our practice that starting with the threshold equal to the 98 th or 99th percentile of the

SC

curvelet coefficients, depending on the noise level and size of the input section, gives us good first approximation of the desired result. It could happen that such estimation of the threshold level is

NU

sufficient and suitable for larger set of sections (i.e. inlines or crosslines from the same data set).

MA

However, if we deal with a more complicated case of noisy data, further analysis is required. For this purpose it is good practice to plot the vector of coefficients in a manner shown in Fig. 2c and

D

Fig. 2d. Instead of introducing one global threshold level (Fig. 2c), one can set different levels

TE

proportional for each scale (Fig. 2d). In practice all scales might be attenuated when containing

AC CE P

coefficients corresponding to noise of certain characteristic or causing visible artifacts, which is similar to less sophisticated and less costly methods, i.e. band-pass or F-K filtering. For the moment multithreshold approach requires careful inspection of the curvelet coefficients and additional human interaction, however it gives significant improvement of the results when compared to global thresholding.

Fig. 3e presents section restored using scale-adjusted thresholding. Coefficients in curvelet domain were thresholded with 3 different levels: 

– scales 1, 2, 6,



– scales 3, 5,



– scale 4.

ACCEPTED MANUSCRIPT Improvement of recovered section is undeniable comparing to Fig. 3c. Also significant amount of artifacts visible at Fig. 3d was eliminated making image clearer. Most of very weak events,

PT

completely unrecognizable at the noisy input (Fig. 3b) section was recovered and the number of

SC

RI

artifact was minimized.

The next example explains how the analysis in the frequency domain helps to choose the

NU

optimal threshold level. Fig. 4a is the same synthetic section as presented on Fig. 3a, with added

MA

white Gaussian noise band-pass filtered between 1-50 Hz. Additionally 2D Fourier spectrum is placed in the inset. As can be seen only strong reflections are recognizable and the remaining part

D

of the signal is completely dominated by noise. Corresponding inset with 2D Fourier spectrum

TE

proves that coherent energy is mixed with noise of similar frequency characteristics. Fig. 4b . Initially we

AC CE P

presents a section recovered using FDCT with global thresholding

preserved only the highest 0.2% of the coefficients and as a result most of the noise has been removed. However, inspecting the inset with 2D frequency spectrum one can notice that in comparison with the initial spectrum of the signal much of the incoherent energy remained in the central, rectangle-shaped area corresponding to certain scales in the curvelet domain. We can make use of this information and set higher attenuation levels for those scales. For this reason in the second approach we introduced more detailed thresholding with 2 levels: 

– scales 1, 2, 3



– scales 4, 5, 6

Again, a slight intervention in the algorithm provides far better result (Fig. 4c).

ACCEPTED MANUSCRIPT In our consideration about optimal thresholding approach we could even go one step further. Presented example proved superiority of the approach of adjusting the coefficients weighting

PT

according to the scales. Why do not try similar approach with certain angles? Of course

RI

transferring such an approach into angular wedges in practice requires more analysis of the

SC

coefficients in the curvelet domain. However, obtained results are incomparably better than in case of global thresholding. In the presented example, coherent energy of the input section is

NU

concentrated within two wedges symmetrical in relation to the origin of 2D frequency spectrum.

MA

After introducing few threshold levels according to the scales we can still point some residual artifacts. Those are easy to remove by attenuating remaining coefficients localized in certain parts

D

of the 2D frequency domain. Region indicated with a dashed line in the inset in Fig. 4c corresponds

TE

to the specific angular wedges of distinct scales. Multiplying threshold levels for those angular

AC CE P

wedges by 1.5 leads to improvement in result quality (Fig. 4d). Higher attenuation of coefficients localized in detected region results in minimizing the number of artifacts and less damage to coherent energy as compared to the results of thresholding adjusted to scales only (compare Fig. 4c and 4d).

Presented example was highly corrupted by a random noise of frequency characteristics similar to the signal. This resulted in significant overlapping of signal and noise coefficients in curvelet domain and in consequence to large damage to the weak reflections. Nevertheless obtained result proves the effectiveness of our adaptive-thresholding method.

3.4 Signal-to-noise ratio definition

ACCEPTED MANUSCRIPT Until now we have been judging the quality of the restored signal based on the “eyeball norm”. In order to express the efficiency of proposed method we will introduce the signal-to-noise ratio

PT

(SNR) estimation algorithm based on assumption that seismic signal is correlated from trace to

can estimate autocorrelation of signal

– number of traces) we

:

NU

,

and

as the average of cross correlations

MA

between two adjacent traces

(where

SC

Based on this assumption for each trace

RI

trace and the noise is uncorrelated.

(3)

first autocorrelations

and

D

To estimate autocorrelation of the noise for the same trace

,

and

AC CE P

TE

are calculated. Averaging results of the following subtractions:

(4)

(5)

we obtain estimation of the noise autocorrelation ANi for trace Ti as follows:

ANi 





1 A 'Ni  A 'Ni1 . 2

(6)

Since autocorrelation has its maximum value when the lag is zero and because the zero-lag value of the autocorrelation is no more than the sum of squared values of samples, hence calculating:

ACCEPTED MANUSCRIPT SNRi 

max( ASi ) rms of signal  rms of noise max ANi

(7)

 

PT

we estimate signal to noise ratio SNRi at trace Ti . Repeating this flow for all pairs of traces [ Ti . ,

RI

Ti 1 ; Ti 1 , Ti  2 ;; TN 1 , TN ] we receive vector of values ( SNR1 , SNR2 , , SNRN 1 ). Final SNR is

NU

1  SNR1  SNR2  SNRN 1  N 1

(8)

MA

SNR 

SC

obtained by averaging:

The approach described above is based on Dash and Obaidullah (1970) and Ikelle and

TE

D

Amundsen (2005).

AC CE P

4. Application to synthetic and real data

4.1 2D synthetic data

For synthetic data considerations we decided to use time stack section from Marmousi2 data set (available at www.agl.uh.edu). The Marmousi2 model contains wide range of complicated structures and reflectors of very weak to strongly outlined amplitudes and varying dips. This complexity gives us the opportunity to show how curvelet-based noise attenuation copes with heterogeneous environment containing different kinds of noise. For our tests a smaller region of the model was extracted and resampled to the size of 1024x1024 pixels as shown in Fig. 5a. Figure 5b and 5c present amplitude and F-K spectra of selected data. The SNR of the clean synthetic section was estimated as 3,97. We introduce additional parameters to describe and compare the

ACCEPTED MANUSCRIPT quality of the results, the L2 norm N and standard deviation σ0. For the clean section from the Marmousi2 model the norm equals 16.48 and the standard deviation is 0.08.

PT

Testing of the algorithm consisted of processing the input section contaminated with synthetic

N 2 =12.773, N 3 =19.159) were modeled.

SC

(σ1=0.1, σ1=0.2, σ3=0.3) and norms ( N1 =6.386,

RI

noise of varying attributes. For this purpose white Gaussian noises of different standard deviations

Subsequently low- (0-25 Hz) and high-pass (50+ Hz) filters were designed. After applying them to

NU

the mentioned white noises, we obtained another 6 (3 low-frequency and 3 high-frequency) types

MA

of noises. Selected bands have highly different dominant frequency as well as possibly large overlap with the signal spectrum. Since the band-pass filtering is changing standard deviations and

D

norms of the results as compared to initially designed white noises, therefore scaling factors were

TE

applied to low- and high-frequency noises in order to equalize them with white noise in respect of

AC CE P

the norm. Fig. 6 demonstrates frequency characteristics of noises (norm N1 =6.386) comparing with the spectrum of the clean synthetic section. In consequence we constructed 9 types of noises with different frequency and standard deviation characteristics. Applying them to the clean synthetic section we obtained 9 input test sections (parameters summarized in Table 1) presented in Fig. 7. The top row (Figs 7a, b and c) shows sections 1-3 contaminated with white noise, the noise in sections 4-6 (Figs 7d, e and f) is low-frequency and in sections 7-9 (Figs 7g, h and i) – highfrequency. Noise amplitudes increase from left- to right-hand side. The standard deviations of noises (up to the value of 0.3) and initial noise-free section (0.08) show that the signal was highly corrupted. The following two figures (Fig. 8 and 9) present frequency and F-K spectra for the sections shown in Fig. 7.

ACCEPTED MANUSCRIPT The noisy sections were used as input data for the denoising algorithm. We used DCT via wrapping with 6 scales and default number of angles which is 16 at 2nd coarsest and 64 at finest scale.

PT

Initially only the highest 1% of the curvelet coefficients at all scales were selected for inverse FDCT.

RI

Quick analysis of the curvelet coefficients and 2D frequency spectra allowed to adjust threshold

SC

levels individually for each scale. Moreover, for certain angular wedges additional threshold scaling was applied depending on the type of noise. Table 1 summarizes thresholds for all sections

NU

according to the scales. Fig. 10 presents polygons subtracted from curvelet grid containing angular

MA

wedges within which threshold levels were multiplied by constant factor c (see Table 1). Shape of

D

the selected polygons depends on the frequency characteristic of the applied noise.

TE

Figure 11 illustrates the effectiveness of our method. Each of the panels presents the results of the

AC CE P

curvelet denoising of the data from respective panels in Fig. 7. In general, increasing the noise amplitudes makes it more difficult to recover the weakest reflections and causes more artifacts in the final image. In the sections with the smallest values of the standard deviation and the norm of the noise amplitudes, the signal looks similarly well restored, independent of the noise frequency spectrum. However, increasing the noise amplitude causes damage to the weaker events. The most accurate results were obtained for the section polluted with high-frequency noise (Fig. 11g – Fig. 11i) where we observe relatively small signal damage and a minimal number of artifacts. This is also reflected in Fig. 12g – Fig. 12i containing frequency spectra of the recovered sections compared to the noise-free synthetic section. Significant damage to the amplitudes of the dominant signal frequency can be observed as the standard deviation of the noise increases (see Fig. 12d – Fig. 12f). This might be a consequence of the more pronounced overlapping of the low-

ACCEPTED MANUSCRIPT frequency signal and noise coefficients in the curvelet domain. Slightly better results are observed in case of the white noise (see Fig. 11a – Fig. 11c and corresponding frequency spectra in Fig. 12a –

PT

Fig. 12c). Additionally Fig. 13 presents F-K spectra of the corresponding filtered sections. As can be

RI

seen, main areas of signal energy are left well preserved, however some small damages are visible

SC

when the norm of noise increases. Tables 2 and 3 summarize the statistics of the synthetic data before and after denoising. Norm and rms values were calculated in order to quantitatively

NU

describe applied and removed noise (Table 2). Levels of applied and removed noise are very

MA

similar one to another which is a desired effect in our synthetic considerations. Slightly higher rms values in case of removed noise suggest that some part of the coherent energy was damaged

D

during thresholding. In order to determine the level of signal loss, analogous statistics for the clean

TE

basic and recovered sections were calculated (Table 3). As could be expected, increasing level of

AC CE P

the applied noise results in more significant signal loss (especially weaker events) during the attenuation process. The problem is particularly evident in case of the low frequency and white noise. The reason might be the similarity of the noise and signal frequencies and the consequent dominance of the noise over the weaker reflectors.

SNR values of the data before and after denoising were collected and presented in Fig. 14. The robustness of the applied method is reflected by significant increase of the SNR and the fact, that the final SNR of the recovered data is a good approximation of the SNR value estimated for the clean basic section. However, by losing portion of the weakest signals we generate some empty spots within sections in time-space domain, which combined with the strong reflectors, highly correlated from trace to trace, causes the excess of the SNR value above the initial level.

ACCEPTED MANUSCRIPT Moreover, SNR of the input sections contaminated with the low-frequency noise (sections 4-6) is higher than other noisy sections. It seems to be a natural consequence of higher cross-correlation

RI

PT

value between traces containing low frequencies.

SC

We finish this section with the examples demonstrating results of F-X deconvolution applied to sections 3, 6 and 9 (Fig. 15a-c). Comparing these sections with corresponding results of our

NU

curvelet-based method one can easily realize volume of noise which survived during F-X decon

MA

attenuation and simultaneous damage of signal energy especially in case of section corrupted with white noise (Fig. 15a). Results presented on Fig. 15b and 15c (denoised section 6 and 9

D

respectively) exhibit better continuity of signal comparing to Fig. 15a. This is consequence of

TE

better prediction for these signal frequencies which are separable with frequency bands of colored

AC CE P

noises. However there are still visible local loses of reflections and phase shifts combined with significant volume of remaining noise and created artifacts. Fig. 16 demonstrates shape and fitting of the frequency spectra of the filtered sections with the noise-free input data. It proves significant damage of signal energy which is undesirable especially that recovered sections demonstrates relatively noisy background. This examples explain and justify the need for using a more complex denoising method we present in this paper. 4.2 2D real data For 2D real data study we have used three different datasets characterized by different acquisition parameters, source frequency and processing steps. All sections represent pre-stack time migration (PreSTM) stacks with mild AGC scaling applied and the following spatial and frequency characteristic:

ACCEPTED MANUSCRIPT 

Section 1 – CDP spacing of 3.125 m, 4 ms sample rate, final band-pass filter of 8-12-7090 Hz Section 2 – CDP spacing of 12.5 m, 2 ms sample rate, final band-pass filter of 4-8-115-

PT



Section 3 – CDP spacing of 15 m, final band-pass filter of 8-12-65-85 Hz

SC



RI

130 Hz

NU

From the described PreSTM stacks, representative 1024 x 1024 pixels areas were extracted and served as input data in further processing. As can be seen in Fig. 17a – Fig. 17c, the chosen

MA

sections differ in terms of formation and intensity of reflectors as well as frequency of the noise. Section 1 (Fig. 17a) mostly consists of straight, parallel interfaces, their dip increasing slightly in

TE

D

deeper parts, corrupted with noise of broad frequency spectrum. Layered structure of the section is especially affected by noise in areas of poorly highlighted amplitudes. In Section 2 (Fig. 17b)

AC CE P

some aligned reflectors of different amplitudes can be observed and the image is damaged by high frequency noise. Section 3 (Fig. 17c) presents more complicated structure with varying dips and visible faults. Volume of the noise in this case seems to be relatively small, however weaker parts of reflectors seem to be damaged by somewhat low frequency energy.

We started the processing with the same FDCT settings that were used for the synthetic sections. Initially the highest 1% of all the curvelet coefficients were selected for the inverse transform. Afterwards, threshold levels were adjusted based on the restored images and the analysis of the vectors of coefficients and the frequency spectra. Taking into account the frequency characteristic of noise and the distribution of coherent energy in curvelet domain for Sections 1 and 3 we applied higher threshold levels for coarser scales and relatively lower values for finer scales,

ACCEPTED MANUSCRIPT whereas in case of Section 2 reversed approach was applied. The threshold values are summarized in Table 4. As in the case of synthetic data here we have introduced additional adjustment of the

PT

threshold levels according to dip. Fig. 18a-c presents polygons extracted from curvelet grid within

RI

which attenuation was performed with threshold level scaled by a factor c =1.5. Filtered sections

SC

are presented in Fig. 17d-f. Improvement of the recovered signal in comparison to the input data is significant. Most of the incoherent energy was suppressed without obscuring useful information

NU

(e.g. location of faults in Fig. 17f) or introducing artifacts. Improvement of data quality is most

MA

visible in Section 2 where the noise amplitudes were significant in comparison to the coherent energy. Fig. 17g-i presents the difference between the input data and the recovered sections, i.e.

D

the energy we consider to be the noise. The noise energy is incoherent and, except for the

TE

strongest event in Section 2, does not seem to follow any of the reflectors, which suggests that in

AC CE P

general the method does not damage the signal and preserves its amplitudes.

Information about frequency characteristic of the input sections as well as of the obtained results and residuals is also reflected in the corresponding frequency spectra and F-K spectra presented in Fig. 19 and Fig. 20. For Section 1 we have obtained relatively narrow frequency band of the signal (about 20-60 Hz, Fig. 19d) whereas noise energy is almost uniformly distributed in the spectrum (Fig. 19g). In case of Section 2 we notice the most significant reduction of energy after noise attenuation as compared to the input section (Fig. 19e). Analyzing residual frequency spectrum (Fig. 19h – green plot) we observe trend of growing noise energy as frequency increase up to about 65Hz. After this point nearly whole energy is treated as residual which is also reflected by the fine-scale character of the noise in the time-space domain. As we mentioned before, in case of

ACCEPTED MANUSCRIPT Section 3 volume of the noise appeared to be relatively low which is reflected by the smaller reduction of the resulting energy in comparison to the input (Fig. 19f). Residual frequency

PT

spectrum of this section (Fig. 19i) contains mainly lower frequencies and its energy consequently

RI

decreases with higher frequencies. Fig. 20 demonstrates F-K spectra of the inputs, results and

SC

residuals respectively. It illustrates the distribution of signal and noise energy for each of the sections introducing more details of the results of our denoising approach and supports the

MA

NU

interpretation.

The values of SNR (Fig. 21) confirm the effect described above. The most significant improvement

D

was observed in case of Section 2. The volume of high frequency noise of the input section in this

TE

case is reflected in a relatively low SNR before the attenuation. Section 1 shows almost twice as

AC CE P

high level of SNR after denoising in comparison to the input data. As we mentioned before, Section 3 did not contain much noise therefore the increase of SNR is not as big but still significant.

4.3 3D real data

Here we demonstrate the application of our curvelet-based denoising algorithm to a real 3D data set acquired in a complex geological setting. The data comprise a 3D partially-migrated (DMO) volume from a 17 km2 survey. It was acquired in Flin Flon mining camp (Manitoba, Canada). Due to high cultural noise associated with production and processing of ore, Flin Flon mining camp was a challenging environment for geophysical exploration. Also variations in natural terrain comprising exposed bedrock, swamps, pockets of shallow glacial deposits and lakes made it difficult to carry out the acquisition. All those inconveniences were reflected in relatively noisy

ACCEPTED MANUSCRIPT background for seismic recording and limited strength of seismic sources that could be employed as well as the irregularities in receiver and shot locations (White et al., 2012). Although the final

PT

DMO stack was post-processed with a F-XY deconvolution, it remains relatively noisy.

RI

Despite the fact that the 2D FDCT can be extended to the third dimension (Ying et al., 2005), its

SC

implementation exhibits a serious disadvantage for practical application. Curvelet transform in 3D produces increased volume of high-redundancy data that is difficult to store for most of regular

NU

computers especially when processing large data sets. Facing these difficulties, slicing the data

MA

cube into 2D sections is required, so that 2D FDCT can be used. Our approach to 3D data noise attenuation was a consequence of initial experiments with applying curvelet based denoising to

D

time slices. Visible and measurable improvement of SNR in inline and crossline sections obtained

TE

from filtered horizontal slices encouraged us to employ 2D FDCT to denoise 3D data. For the

AC CE P

presented Flin Flon dataset (171 inlines spaced every 25 m, 340 crosslines spaced every 12.5 m and 1000 time slices every 2 ms) we started by filtering the inline sections. Afterwards, data were denoised in the crossline sections. This was done due to the fact that signal may be more difficult to capture on the time slices than on the vertical sections especially in case of hard rock data lacking clear coherent events. Such approach allowed to boost coherent energy before final filtering of the time slices. Another problem in this case was the estimation of the threshold level. The number of sections and slices would make individual adjustment of the threshold for each portion of data very timeconsuming. To solve this problem we took the advantage of the fact that samples from the same data set present approximately similar characteristic in frequency and curvelet domain. After reviewing the frequency spectra and vectors of curvelet coefficients of selected samples, we opted

ACCEPTED MANUSCRIPT for partitioning into 6 scales for inlines and crosslines and 5 scales for time slices. The number of angles was set to default which is 16 at 2nd coarsest scale and doubling every second scale.

PT

Threshold level was initially set to preserve only the highest 2% of all curvelet coefficients

RI

obtained from each input section or slice. After a few tests we concluded that the inline and

SC

crossline sections are denoised efficiently when stronger attenuation of the coarsest and the finest scale coefficient is applied. For scales 1, 2, 6 (inlines filtering) and 1, 5, 6 (crosslines filtering)

NU

threshold levels were multiplied by c =1.5. In case of time slices we noticed that filtering

MA

performed better while preserving more coefficients from the coarser scales, hence threshold levels for scales 4, 5 were multiplied by the same constant c =1.5. Final increase of the mean SNR

D

after filtration with the scale-adaptive thresholds compared to the global thresholding was equal

TE

to 12% when measured for inlines and 27% when estimated for crosslines. It validates our earlier

AC CE P

conclusions that introducing more than one threshold level translates into improved results.

Fig. 22 shows a fence-plot of the 3D volume used for filtering. Comparing data before (Fig. 22b) and after (Fig. 22c) noise attenuation we observe high improvement of signal quality in all directions. The recorded data are highly contaminated with noise so that the mapping of seismic reflectors, which are discontinuous and crooked due to the geology of the area, is problematic. Curvelet noise attenuation enhances the strong amplitude reflections and makes it easier to interpret the data. Moreover, despite the data cube slicing and separate processing in each direction, signal continuity at intersections of surfaces is preserved, making interfaces easy to identify regardless of their location in the inline, crossline or time slice. Difference between input and output data (Fig. 22d) shows the amount and the generally incoherent character of the energy

ACCEPTED MANUSCRIPT separated from the signal. Fig. 23 presents the summary of SNR calculated for each inline and crossline section. Mean SNR values for inlines/crosslines were estimated at 1.36/0.63 before and

PT

2.69/1.76 after denoising. Since our SNR estimation is based on cross-correlation of traces,

RI

calculating it in the time slices is not applicable. SNR for the input sections is very regular, both for

SC

inline and crossline sections. It proves that the signal was uniformly corrupted by noise within the whole data cube. However, larger deviations of SNR in recovered data suggests that at some

NU

regions signal was either easier or more difficult to catch. It also might be caused by non-uniform

MA

distribution of strong reflections along data set. The SNR for both raw and denoised data seem to follow a similar section-to-section trend. We believe that the presented approach of running 2D

TE

D

DCT for 3D data might be sufficient substitution for 3D DCT.

AC CE P

4.4.Application for data conditioning for full-waveform inversion

Motivated by the results of 2D DCT applied to the time slices from 3D seismic data, we explore the application of curvelet denoising to enhance the signal in frequency maps (panels). By frequency map we understand the real or imaginary part of the Fourier transform for one particular frequency extracted from 2D seismic data, presented in source-receiver coordinates. It represents the data that are subjected to frequency-domain full-waveform inversion (FWI) (Pratt, 1999). The image of a frequency map is eligible for curvelet noise attenuation, since the signal has clear directional character (stripes along the diagonal of the image) and the noise is unidirectional (Fig. 24a). Here we use the data from the regional deep reflection seismic profile (POLCRUST-01 profile, see Malinowski et al. (2013) for details), acquired with the Vibroseis technique (6-72 Hz sweep

ACCEPTED MANUSCRIPT frequency was used). Mono-frequency subsets of the data were extracted in order to run frequency domain FWI. Our curvelet denoising scheme was applied separately on the real and

PT

imaginary parts of the mono-frequency data. Limited frequency band allowed for easy weighting

RI

of coefficients according to scales. Moreover, because frequency maps exhibit highly directional

SC

nature, they are typical example of data suitable for angle dependent thresholding. Fig. 24 shows the results of curvelet denoising for the real part of the data at 6 Hz. It is clear from Fig. 24b, that

NU

only the continuous events were preserved in filtration. Horizontal stripes of coherent energy

MA

visible in Fig. 24c (i.e. the difference between input and filtered panel) represent noisy receiver positions. Vertical stripes are related to the shots with amplified amplitudes, e.g. noisy shot

D

record. Both were filtered out by our method. Bottom panels of Fig. 24 shows the reconstructed

TE

phase, calculated on the filtered real and imaginary parts. Comparison of the input, filtered and

AC CE P

the residual (Fig 24 a, b, c respectively) section suggests that the phase was preserved during our denoising and it is easier to follow its changes. Initial tests with running FWI with this dataset indicate that the inversion with curvelet-denoised data outperforms the one with the normally preconditioned data (e.g. band-pass filtered). This is because the attenuation of the noise in the curvelet domain significantly improves the quality of the data and allows the FWI to fit the signal not the noise during the inversion.

5. Discussion Thresholding is not a trivial subject considering any type of filtering. It has the largest and most immediate impact on the results. We believe that the approach in which thresholding is locally adjusted at certain scales or angles leads to most efficient noise attenuation. We have already

ACCEPTED MANUSCRIPT mentioned that it is good practice to set the first approximation of the threshold level as the largest 1-2% of the curvelet coefficients. Our observations indicate that the smaller and more

PT

noisy the section, the more curvelet coefficients need to be inverted back to the x-t domain

RI

(assuming the same partitioning of frequency domain). However, this approach still fails when

SC

employed to global thresholding. Adjusting the global threshold always requires additional iterations and in the end never provides results comparable to those obtained with scale-adapted

NU

threshold levels. The weakness of global thresholding manifests itself i.e. in the form of artifacts

MA

that can be removed without damaging the signal only when certain frequency bands (via curvelet scales and/or angular wedges) are directly accessed and attenuated.

D

The direct extension of traditional thresholding approach is to set the weighting levels in the

TE

same way for specific scales. This approach can be effective even in the presence of white noise

AC CE P

and leads to the results which are much more precise then the results obtained from the inversion of globally thresholded coefficients. The access to particular frequency bands gives more leeway in separating signal from noise and avoiding artifacts. When setting the threshold levels we use the vector of coefficient and 2D frequency spectra. Thanks to better understanding and interpretation of signal and noise characteristics as well as their distribution in curvelet and F-K domain we are immediately closer to the desired, optimal solution. Hence knowing basic features of curvelet transform one is able to achieve far better results in very few iterations, while global thresholding is unable to provide comparable, artifactfree results. Transferring this scale-localized approach to angle-adjustable thresholding would be much more difficult because of the number of angular wedges. We demonstrated such method in

ACCEPTED MANUSCRIPT practice, however, its application and decision-making process should be based rather on

this issue would translate into even better noise suppression.

PT

statistical calculations than on manual selections. Nevertheless, one may expect that addressing

RI

Comparison with simple F-X filtering shows that although curvelet noise attenuation is much

SC

more costly, the results are significantly better.

The emphasis on noise reduction with a minimum number of iterations is due to the

NU

computational complexity of FDCT. Table 5 shows times for direct and inverse FDCT with wavelets

MA

and curvelets at the finest scale (spectrum partitioned into 7 scales) depending on the size of the input section. From the obtained results it can be concluded that for large datasets (i.e. 3D data)

D

complete noise attenuation might be very time-consuming, inverse FDCT being the most costly

TE

step of the algorithm. Careful choice of the threshold levels, based upon the qualitative

AC CE P

comparison of denoised images along with the analysis of frequency spectrum and vectors of curvelet coefficients allows to minimize the number of iterations and hence the computational cost and the time of calculations.

6. Conclusions

Presented paper illustrates the efficiency of curvelet noise attenuation in enhancing coherent signal in seismic data from experiments done in various environments. Applying proposed algorithm to synthetic data contaminated with noise of various characteristics, we managed to find a practical denoising scheme. The approach we developed working on synthetic data proved its effectiveness and functionality on both 2D and 3D real data. Introduced scheme is intuitive and, once the structure of curvelet representation is understood, it is simple to apply in practice. The

ACCEPTED MANUSCRIPT results of scale-adjusted thresholding demonstrate significant improvement over the images obtained with the global threshold version of the algorithm, which justifies this procedure, even

PT

though at the current stage of development it requires a lot of human interaction.

RI

Demonstrated examples show that even complex structures of crooked reflectors or varying

SC

dips can be recovered from very noisy data. Denoised sections are clearer and easier to interpret. Proper use of scale-adjusted thresholding ensures that obtained images are artifact-free.

NU

Furthermore, is was proven, that noise suppression in curvelet domain is effective not only for

MA

stacked 2D and 3D seismic data or shot gathers (e.g. Kumar 2011), but also for frequency maps (i.e. the form of data input to full-waveform inversion).

D

The thresholding algorithms presented in this paper are not the final stage of developing a

TE

curvelet-based tool for seismic data denoising. The scale-adjusted thresholding involves too much

AC CE P

human interaction to be used for large data sets consisting of numerous sections differing in noise and signal characteristics. The main goal is to automate the choice of the threshold levels based on the statistical analysis or pattern recognition methods. Extending the algorithms with angleadjusted thresholding was also only mentioned in the text, applied to a very simple example. Manually adjusting the threshold values for each angle is too laborious to be performed for more complicated data sets. However, the idea is very appealing and making the effort to automate this process is certainly worthwhile. Taking into account how promising the presented examples are, we expect that with more sophisticated methods of setting the thresholds, the results of curvelet denoising can be even more rewarding.

ACCEPTED MANUSCRIPT 7. Acknowledgments

PT

This work is funded by the Polish National Science Centre grant no 2011/03/D/ST10/05128.

SC

RI

Comments made by A. Malehmir helped us to improve our paper.

AC CE P

TE

D

MA

NU

DCT algorithm was taken from the CURVELAB project (www.curvelab.org).

ACCEPTED MANUSCRIPT 8. References

PT

Bednar, J., 1983. Applications of median filtering to deconvolution, pulse estimation, and statistical editing of seismic data. Geophysics 48, 1598–1610.

RI

Canales, L.L., 1984. Random noise reduction, in: SEG Technical Program Expanded Abstracts 1984. pp. 525–527.

SC

Candès, E.J., Demanet, L., Donoho, D.L., Ying, L., 2006. Fast Discrete Curvelet Transforms. Multiscale Model. Simul. 5, 861–899.

NU

Candès, E.J., Donoho, D.L., 1999. Ridgelets: a key to higher-dimensional intermittency? Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 357, 2495–2509.

MA

Candès, E.J., Donoho, D.L., 2000. Curvelets – A Surprisingly Effective Non adaptive Representation For Objects with Edges. Vanderbilt Univ. Press 105–120.

D

Dash, B.P., Obaidullah, K.A., 1970. Determiantion of signal and noise statistics using correlation theory. Geophysics 35, 24–32.

AC CE P

TE

Hennenfent, G., Cole, J., Kustowski, B., 2011. Interpretative noise attenuation in the curvelet domain, in: SEG Technical Program Expanded Abstracts 2011. Society of Exploration Geophysicists, pp. 3566–3570. Herrmann, Felix J., Moghaddam, P., Stolk, C.C., 2008. Sparsity- and continuity-promoting seismic image recovery with curvelet frames. Appl. Comput. Harmon. Anal. 24, 150–173. Herrmann, F. J., Wang, D., Hennenfent, G., Moghaddam, P., 2008. Curvelet-based seismic data processing: A multiscale and nonlinear approach. Geophysics 73, A1–A5. Ikelle, L., Amundsen, L., 2005. 5. Characterization of Seismic Signals by Statistical Averages. Introd. to Pet. Seismol. 181–232. Kappus, M.E., Harding, A.J., Orcutt, J.A., 1990. A comparison of tau‐ p transform methods. Geophysics 55, 1202–1215. Kumar, V., Oueity, J., Clowes, R.M., Herrmann, F.J., 2011. Enhancing crustal reflection data through curvelet denoising. Tectonophysics 508, 106–116. Malinowski, M., Guterch, A., Narkiewicz, M., Probulski, J., Maksym, A., Majdański, M., Środa, P., Czuba, W., Gaczyński, E., Grad, M., Janik, T., Jankowski, L., Adamczyk, A., 2013. Deep seismic reflection profile in Central Europe reveals complex pattern of Paleozoic and Alpine accretion at the East European Craton margin. Geophys. Res. Lett. 40, 3841–3846.

ACCEPTED MANUSCRIPT Martin, G.S., Wiley, R., Marfurt, K.J., 2006. Marmousi2: An elastic upgrade for Marmousi. Lead. Edge 25, 156–166.

PT

Neelamani, R., Baumstein, A.I., Gillard, D.G., Hadidi, M.T., Soroka, W.L., 2008. Coherent and random noise attenuation using the curvelet transform. Lead. Edge 27, 240–248.

RI

Pratt, R.G., 1999. Seismic waveform inversion in the frequency domain, Part 1: Theory and verification in a physical scale model. Geophysics 64, 888–901.

SC

Starck, J., Candès, E.J., Donoho, D.L., 2002. The curvelet transform for image denoising. IEEE Trans. Image Process. 11, 670–84.

MA

NU

Starck, J., Donoho, D.L., Candès, E.J., 2001. Very high quality image restoration by combining wavelets and curvelets, in: Laine, A.F., Unser, M.A., Aldroubi, A. (Eds.), Wavelets Applications in Signal and Image Processing IX. SPIE, pp. 9–19. Stewart, R.R., Schieck, D.G., 1989. 3‐D F‐K filtering, in: SEG Technical Program Expanded Abstracts 1989. Society of Exploration Geophysicists, pp. 1123–1124.

TE

D

Turner, G., 1990. Aliasing in the tau- p transform and the removal of spatially aliased coherent noise. Geophysics 55, 1496–1503.

AC CE P

White, D.J., Secord, D., Malinowski, M., 2012. 3D seismic imaging of volcanogenic massive sulfide deposits in the Flin Flon mining camp, Canada: Part 1 — Seismic results. Geophysics 77, WC47–WC58. Yarham, C., Herrmann, F.J., 2008. Bayesian ground‐roll separation by curvelet‐domain sparsity promotion, in: SEG Technical Program Expanded Abstracts 2008. Society of Exploration Geophysicists, pp. 3662–3666. Ying, L., Demanet, L., Candès, E.J., 2005. 3D Discrete Curvelet Transform, in: Papadakis, M., Laine, A.F., Unser, M.A. (Eds.), pp. 591413–591413–11. Zhou, B., Greenhalgh, S.A., 1994. Linear and parabolic τ- p transforms revisited. Geophysics 59, 1133–1149.

ACCEPTED MANUSCRIPT Figure captions Fig. 1. Examples of few curvelets of different scales, angles and location in spatial domain (a) and

PT

their representation in frequency domain (b) where they are positioned in the orthogonal

RI

direction; (c) Partitioning of the 2D frequency plane into concentric, rectangular annuli and their

SC

inner segmentation into angular wedges which obey the parabolic scaling. Number of wedges

NU

doubles each second scale yielding curvelets of higher frequency more anisotropic. (d) Zoom of curvelet nr 3 from (a); Each of the curvelets is described by its frequency, angle (dip) and position.

D

region indicated by dashed line ellipse.

MA

Localization (x1,x2) is possible because of rapid amplitude decay in spatial domain, outside the

TE

Fig. 2. Influence of densely partitioned frequency domain on extracting coherent energy at certain scales. The same synthetic section was transformed into curvelet domain with 4 (a) and 7 (b)

AC CE P

scales. Threshold levels for each transform are displayed by dashed lines. In consequence higher number of scales allowed thresholding at different levels (d) and better selection of coherent energy as compared to global threshold level in case of 4 scales (c). Fig. 3. Difference between choosing wavelets or curvelets at finest scale. (a) Demonstrative, noise-free synthetic section from the Marmousi2 model. (b) Section (a) with white Gaussian noise added. Recovered sections at (c) and (d) are results of global thresholding after transform using 6 scales with wavelets and curvelets at finest scale respectively. In both cases visible artifacts of different shapes arises. Introducing few threshold levels allow higher attenuation of noise and better signal passing as shown in (e).

ACCEPTED MANUSCRIPT Fig. 4. Threshold level adjusting using 2D Fourier spectrum analysis. (a) Synthetic section presented on Fig. 3a. has been highly corrupted with colored (0-50 Hz) noise. Section recovered in

PT

consequence of initial global thresholding level (b) still contains much of low frequency noise.

RI

From the 2D spectra shown as insets one can localize scales which contain noise coefficients and

SC

better adjust threshold levels for them. (c) Section recovered after increase of threshold level at scales containing low frequency curvelets still containing some residual artifacts. (d) Final section

NU

after additional attenuation of certain angles (symmetric polygons bounded by dashed lines

MA

shown in the inset (c)).

Fig. 5. Synthetic section extracted from the Marmousi2 model (a) with its frequency (b) and F-K

TE

D

spectrum (c). Note the complexity of structure and variety of reflectors amplitudes. Fig. 6. Frequency characteristic of the noise applied to the extracted part of Marmousi2 model.

AC CE P

Examples with the lowest value of norm N1 =6.386 were selected. Fig. 7. Synthetic section extracted from the Marmousi2 model (same as Fig. 5a) contaminated with noise of various characteristics. (a)-(c) correspond to Sections 1-3, (d)-(f) and (g)-(i) present Sections 4-5 and 6-9 respectively.

Fig. 8. Comparison of the frequency spectra before (green plot) and after (red plot) applying each noise to respective sections from Fig. 7. Note increasing absolute value of noise amplitude in comparison with signal energy. Each pair of the spectra were normalized to 1. Fig. 9. F-K spectra of input sections presented in Fig. 7. For each section we can observe distribution of the noise relatively to the location of signal in the 2D frequency plane. All noise

ACCEPTED MANUSCRIPT types have zero-centered Gaussian distribution but differ in respect to frequency characteristic, norm and standard deviation.

PT

Fig. 10. Partitioning of the F-K spectrum we have used during DCT. After setting threshold levels

RI

according to scales we introduce additional threshold weighting factor c for coefficients localized

SC

in certain angular wedges. Dark polygons (a) – white noise, Sections 1-3; (b) – low-frequency noise,

NU

Sections 4-6; (c) – high-frequency noise, Sections 7-9, subtracted from the curvelet grid contain coefficients thresholded with additionally weighted threshold levels. Compare those polygons with

MA

the F-K spectra presented in Fig. 9.

D

Fig. 11. Sections from Fig. 7 after denoising. Differences between obtained results reflect influence

TE

of applied noise on the results. Slight damage of weaker events amplitudes which increase with

AC CE P

standard deviation of noise is observed.

Fig. 12. Comparison of noise-free synthetic section frequency spectrum (red plot) and recovered sections frequency spectra (green plot). We can observe preserving of spectrum shape, however reduction of signal energy is visible for data contaminated with noise of higher norm. Fig. 13. F-K spectra of the filtered sections presented in Fig. 11. Comparing with the F-K spectrum of the clean synthetic section (Fig. 5c) we can observe impact of filtration on the shape of spectra. Fig. 14. Summary of the SNR before and after noise attenuation for each section. Dashed line represents SNR level of noise-free initial data. Fig. 15. Results of F-X deconvolution denoising for Section 3 (a), Section 6 (b) and Section 9 (c). In case of each section we observe significant share of unwanted energy which survived denoising.

ACCEPTED MANUSCRIPT Fig. 16. Frequency spectra of sections presented in Fig. 15 (black line – section (a), blue line – section (b), red line – section (c), green polygon – clear synthetic section). Fitting of frequency

PT

bands for obtained results proving simultaneous visible damage caused to the signal.

RI

Fig. 17. 2D real data example using PreSTM stacks. (a)-(c) Input Sections 1-3. (d)-(f) Corresponding

SC

sections recovered after noise attenuation. (g)-(i) Energy taken as noise and extracted from each

NU

section.

Fig. 18. Partitioning of the F-K spectrum we have used during DCT of Sections 1-3 (panels a-c

MA

respectively). After setting threshold levels according to scales we introduce additional threshold

D

weighting factor c for coefficients localized in certain angular wedges. Dark polygons subtracted

TE

from curvelets contain coefficients thresholded with additionally weighted threshold levels.

AC CE P

Fig. 19. Frequency spectra of the input sections (a) Section 1, (b) Section 2, (c) Section 3. (d-f)/(g-i) Comparison of the input sections frequency spectra (red plots) and recovered signal/noise frequency spectra (green plot).

Fig. 20. F-K spectra of the input sections (a) Section 1, (b) Section 2, (c) Section 3. (d-f)/(g-i) corresponding F-K spectra of the recovered signal/noise. Fig. 21. SNR values before and after noise attenuation for each section from Fig. 16. Fig. 22. (a) Cross-section of the partially-migrated (DMO) stack volume from the Flin Flon 3D survey (White et al., 2012). (b)(c) Zooms of the section obtained from (a) presenting chosen inline, crossline and time slice before and after noise attenuation respectively. (d) Energy considered as noise end extracted from input data.

ACCEPTED MANUSCRIPT Fig. 23. Comparison of the SNR for (a) inlines and (b) crosslines before and after noise attenuation. Mean values of SNR for each set of sections were marked by dashed lines.

PT

Fig. 24 (a) Plot of the real part of the 6-Hz frequency data extracted from the POLCRUST-01 deep

RI

reflection seismic profile in the source-receiver coordinate system (frequency map); (b) same

SC

panel after curvelet denoising; (c) difference between panels (a) and (b); (d) phase of the 6-Hz

NU

data; (e) phase calculated from the curvelet-denoised real and imaginary parts of 6-Hz data; (f)

AC CE P

TE

D

MA

difference between panels (d) and (e).

AC CE P

TE

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 1

AC CE P

TE

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 2

AC CE P

Fig. 3

TE

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 4

AC CE P

TE

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

AC CE P

TE

D

Fig. 5

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

AC CE P

TE

D

Fig. 6

Figure 7

AC CE P

TE

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Figure 8

AC CE P

TE

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 9

AC CE P

TE

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

AC CE P

TE

D

Figure 10

Figure 11

AC CE P

TE

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Figure 12

AC CE P

TE

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 13

AC CE P

TE

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

AC CE P

TE

D

Fig. 14

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

AC CE P

TE

D

Fig. 15

Fig. 16

AC CE P

TE

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 17

AC CE P

TE

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

AC CE P

TE

D

MA

Fig. 18

Fig. 19

AC CE P

TE

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Figure 20

AC CE P

TE

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

AC CE P

Fig. 21

TE

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

AC CE P

TE

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 22

Fig. 23

AC CE P

TE

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

AC CE P

Fig. 24

TE

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

ACCEPTED MANUSCRIPT Table 1. Threshold levels used during denoising of Sections 1-9 according to scales. Factor c denotes multiplier of threshold for angular wedges shown on Fig. 8. thresholds according to scales

threshold scaling according to angles

frequency characteristic

standard deviation

norm

Scale 1

Scale 2

Scale 3

Scale 4

Scale 5

Scale 6

1

white noise

σ1=0.1

N1 =6.386

0,14

0,12

0,12

0,12

0,16

0,16

c

= 1.5

2

white noise

σ2=0.2

N 2 =12.773

0,28

0,22

0,22

0,22

0,32

0,32

c

= 1.5

3

white noise

σ3=0.3

N 3 =19.159

0,42

0,32

0,32

0,32

0,48

0,48

c

= 1.5

σ4=0.06

N1 =6.386

0,20

0,15

0,15

0,10

0,10

0,10

c

=3

σ5=0.12

N 2 =12.773

0,40

0,30

0,30

0,22

0,15

0,15

c

=3

σ6=0.18

N 3 =19.159

0,60

0,45

0,45

0,36

0,20

0,20

c

=3

σ7=0.09

N1 =6.386

0,05

0,05

0,05

0,07

0,20

0,20

c

= 0.5

σ8=0.18

N 2 =12.773

σ9=0.27

N 3 =19.159

8 9

RI

SC

NU

MA

7

0,05

0,05

0,05

0,14

0,35

0,35

c

= 0.5

0,05

0,05

0,05

0,21

0,45

0,50

c

= 0.5

D

6

TE

5

low freq. noise (0-25Hz) low freq. noise (0-25Hz) low freq. noise (0-25Hz) high freq. noise (50+ Hz) high freq. noise (50+ Hz) high freq. noise (50+ Hz)

AC CE P

4

PT

section

noise parameters

ACCEPTED MANUSCRIPT Table 2. Applied and recovered noise parameters. white noise

low frequency noise

high frequency noise

PT

norm rms norm rms norm rms standard standard standard recovere recovere recovere recovere recovere recovere deviation applied deviation applied deviation applied applied applied applied d d d d d d 6,386

6,465

0,100

0,101

σ4=0.06

6,386

6,275

0,052

0,053

σ7=0.09

6,386

6,403

0,094

0,095

σ2=0.2

12,773

12,742

0,200

0,200

σ5=0.12

12,773

12,605

0,104

0,106

σ8=0.18

12,773

12,687

0,189

0,189

σ3=0.3

19,159

19,070

0,300

0,300

σ6=0.18

19,159

18,935

0,156

19,159

18,991

0,283

0,283

RI

σ1=0.1

σ9=0.27

SC

0,158

NU

Table 3. Recovered signal parameters. white noise

low frequency noise

high frequency noise

section 1

section 2

section 3

section 4

section 5

section 6

section 7

section 8

section 9

norm

16,484

15,710

15,014

14,448

15,781

14,675

13,731

16,039

15,743

15,425

rms

0,080

0,075

0,070

0,067

0,075

0,070

0,063

0,076

0,073

0,071

AC CE P

TE

D

MA

input section

ACCEPTED MANUSCRIPT Table 4. Threshold levels used during denoising of Sections 1-3 according to scales. Factor c denotes multiplier of threshold for angular wedges shown on Fig. 18. Scale 3 0,25 0,2 0,3

Scale 4 0,15 0,3 0,2

Scale 5 0,15 0,3 0,15

PT

Scale 2 0,25 0,2 0,3

Scale 6 0,15 0,3 0,15

c 1,5 1,5 1,5

RI

Section 1 Section 2 Section 3

Scale 1 0,25 0,2 0,3

NU

MA

D TE

512x512 1024x1024 2048x2048 4096x4096 8192x8192 16384x16384

Curvelets at finest scale Direct transform Inverse transform 0,5 s 1,3 s 1,9 s 2,6 s 6,3 s 7,8 s 21,4 s 25,7 s 90,3 s 105,8 s 403,3 s 695,3 s

AC CE P

Section size

SC

Table 5. Comparison of computation times regarding to type of DCT and data size. Wavelets at finest scale Direct transform Inverse transform 0,2 s 0,5 s 0,7 s 1,4 s 1,9 s 4,1 s 7,2 s 13,1 s 26,7 s 48,5 s 114,3 s 229,7 s

ACCEPTED MANUSCRIPT

TE

D

MA

NU

SC

RI

Strategy for complex scale- and angle-adaptive thresholding. Denoising both synthetic and real type of seismic data. Attenuation of various types of noises. Application of 2D DCT for 3D seismic data denoising.

AC CE P

   

PT

Highlights