Comparing three feedback internal multiple elimination methods

Journal of Applied Geophysics 95 (2013) 66–76

Contents lists available at SciVerse ScienceDirect

Journal of Applied Geophysics journal homepage: www.elsevier.com/locate/jappgeo

Comparing three feedback internal multiple elimination methods Jiawen Song a,⁎, Eric Verschuur b, Xiaohong Chen a a b

China University of Petroleum, State Key Laboratory of Petroleum Resources and Prospecting, Beijing, China Delft University of Technology, P.O. Box 5046, 2600 GA Delft, The Netherlands

a r t i c l e

i n f o

Article history: Received 1 March 2013 Accepted 21 May 2013 Available online 30 May 2013 Keywords: Internal multiple elimination Feedback model Prediction-and-subtraction Sparse inversion

a b s t r a c t Multiple reflections have posed a great challenge for current seismic imaging and inversion methods. Compared to surface multiples, internal multiples are more difficult to remove due to poorer move-out discrimination with primaries and we are left with wave equation-based prediction and subtraction methods. In this paper, we focus on the comparison of three data-driven internal multiple elimination (IME) methods based on the feedback model, where two are well established prediction-and-subtraction methods using back-propagated data and surface data, referred to as CFP-based method and surface-based method, respectively, and the third one, an inversion-based method, has been recently extended from estimation of primaries by sparse inversion (EPSI). All these three methods are based on the separation of events from above and below a certain level, after which internal multiples are predicted by convolutions and correlations. We begin with theory review of layer-related feedback IME methods, where implementation steps for each method are discussed, and involved event separation are further analyzed. Then, recursive application of the three IME methods is demonstrated on synthetic data and field data. It shows that the two well established prediction-andsubtraction methods provide similar primary estimation results, with most of the internal multiples being removed while multiple leakage and primary distortion have been observed where primaries and internal multiples interfere. In contrast, generalized EPSI provides reduced multiple leakage and better primary restoration which is of great value for current seismic amplitude-preserved processing. As a main conclusion, with adaptive subtraction avoided, the inversion-based method is more effective than the prediction-and-subtraction methods for internal multiple elimination when primaries and internal multiples overlap. However, the inversion-based method is quite computationally intensive, and more researches on efficient implementation need to be done in future. © 2013 Elsevier B.V. All rights reserved.

1. Introduction Most current imaging algorithms are based on primary reflection energy. Therefore, multiples need to be removed from the data in a preprocessing step. For surface multiples, one important family of methods is the surface-related multiple elimination (SRME) methodology, proposed by Verschuur et al. (1992), Berkhout and Verschuur (1997), Verschuur and Berkhout (1997) and Weglein et al. (1997), which has proven to be quite successful over the years. The SRME procedure consists of two steps: multiple prediction and adaptive subtraction, in which the subtraction process is assumed to compensate for the prediction errors. In practice, the most robust and commonly used least-squares filtering method (Abma and Kabir, 2005) may lead to suboptimal removal of multiples (e.g. Guitton and Verschuur, 2004; van Borselen et al., 2003). To circumvent the limitations of adaptive subtraction, van Groenestijn and Verschuur (2009a,b) proposed estimating primaries by sparse

⁎ Corresponding author. Tel.: +86 15101175794. E-mail address: [email protected] (J. Song). 0926-9851/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jappgeo.2013.05.010

inversion (EPSI), where primaries are estimated as unknowns in a multidimensional, large-scale inversion process such that they, together with their corresponding multiples, explain the total measured data. The main advantage of EPSI over SRME is that the adaptive subtraction is avoided. Compared to surface multiples, internal multiples usually appear as less distinct arrivals in the seismic section as the downward scattering interface may also have a spatially variable character, imposing an extra imprint on the multiples. Furthermore, they often have less velocity discrimination with the primaries compared to surface multiples, as the extra bounce does not happen in the low-velocity water layer, but inside one of the other subsurface layers (Verschuur, 2006). Therefore, one has been looking for data-driven methods where exact knowledge of the downward scattering interface is not required. Currently, there are two comprehensive approaches satisfying the above criterion, one is the fully hands-off inverse scattering method, the other is the more feasible feedback method. For the inverse scattering method (see e.g. Araujo et al., 1994; Coates and Weglein, 1996; Matson et al., 1999; Weglein et al., 1997), each point in the subsurface is considered to be a possible downward scatterer. The

J. Song et al. / Journal of Applied Geophysics 95 (2013) 66–76

advantage of this method is that it can remove all internal multiples in one integrated process. The downside, however, is that it is very computationally intensive, limiting its use in practical situations. For the feedback method, the multiple problems are considered in a layer-stripping fashion, where in each step one (group of) downward reflecting interfaces is considered (Verschuur et al., 1998). As a representative feedback method, SRME has been extended for internal multiple elimination (IME) by taking the varying properties of the internal multiple generating boundaries into account in a data-driven manner (Berkhout and Verschuur, 1999, 2005; Verschuur and Berkhout, 2005). Although some interpretation is involved, the calculation time is considerably lower than the inverse scattering approach. This IME method involves back-propagation of the surface data to a certain subsurface level, at which the convolution process is carried out in order to predict the associated internal multiples. We will refer to this method as the CFP-based IME method, as will be explained later on in the paper. Jakubowicz (1998) showed that such inverse wavefield extrapolation is not required and that all convolutions can also be directly carried out at the surface level. Therefore, we will refer to this method as the surface-based IME method. Like SRME, these two IME methods share the same procedure, i.e. prediction and subtraction. However, due to small move-out discrimination between internal multiples and primaries, the limitations of adaptive subtraction may become more serious. Recently, Ypma and Verschuur (2013) redefined internal multiple elimination as a full waveform inversion process, following the principles of EPSI (van Groenestijn and Verschuur, 2009a). In this generalized EPSI method the internal multiples can be estimated and removed from the original data by plain subtraction. In this paper, we focus on the comparison of three feedback IME methods: the CFP-based method, the surface-based method and the recently developed generalized EPSI process. We begin with a general review of the theory behind these three IME methods. Next, the IME methods are compared in terms of their specific theory and implementation aspects. Then, the effectiveness of the three methods is demonstrated on synthetic and field data. 2. Review of layer-related feedback IME methods In the detail-hiding operator notation for 2D data (Berkhout, 1982), a bold quantity represents a pre-stack data volume for one frequency. Note that the matrix notation can be generalized to the full 3D case, as described by Kinneging et al. (1989). With this notation, the forward model of primaries and multiples can be described as Pðz0 ; z0 Þ ¼ P0 ðz0 ; z0 Þ þ M0 ðz0 ; z0 Þ;

ð1Þ

where P, P0 and M0 represent the total upgoing wavefield at the surface, the primaries (containing also the internal multiples) and the surfacerelated multiples, respectively. After removing multiples generated by downward scattering interfaces 0 to n − 1, Eq. (1) can be generalized to Pn−1 ðz0 ; z0 Þ ¼ Pn ðz0 ; z0 Þ þ Mn ðz0 ; z0 Þ;

ð2Þ

where Mn represents all multiples generated by interface n, and Pn represents the seismic data that do not contain multiples generated by interfaces 0 to n. For the case n > 0, it means the result after internal multiple elimination (IME). In practice, internal multiples are often related to strong reflecting interfaces, and, thus, it is of no extra benefit to apply IME interface by interface. As an alternative, Berkhout and Verschuur (1999) introduced the layer-related scheme in which all the downward reflecting effects of the overburden can be accounted for. For the layer-related version, the subsurface is manually divided into macro layers by putting virtual depth levels in between strong reflecting interfaces. All internal multiples that pass the chosen depth level at least four times can be predicted in one processing step.

67

Eq. (2) can be redefined as a layer-related process by considering subscript n as a virtual depth level. After removing all primary events above virtual depth level zn − 1, the forward model for removing internal multiples related to depth level zn is obtained: P n−1 ðz0 ; z0 Þ ¼ P n ðz0 ; z0 Þ þ δPn ðz0 ; z0 Þ þ Mn ðz0 ; z0 Þ;

ð3Þ

where P n ðz0 ; z0 Þ represents data with all primaries and multiples generated by interfaces above depth level zn being removed. It is a special case m of a more general notation P n ðz0 ; z0 Þ, in which the subscript indicates the virtual depth level up to which multiple elimination has been applied and the overbar indicates that primary events above the depth level, denoted by the superscript, are also removed. Note, for ease of n reading, we have defined P n ðz0 ; z0 Þ ¼ P n ðz0 ; z0 Þ. δPn(z0,z0) represents primary events generated by interfaces between depth level zn − 1 and zn. Mn(z0,z0) represents internal multiples related to zn. For more details on this notation, the reader is referred to Ypma and Verschuur (2013). Note that the three datasets add up to the input data P n−1 ðz0 ; z0 Þ, which is the output data from a previous internal multiple elimination step. According to how internal multiples are predicted, there are three feedback IME methods: the CFP-based method, the surface-based method and the generalized EPSI method. Construction of internal multiples is demonstrated for these three methods (see also Fig. 1). We will explain the theory and give practical implementation steps for each method. 2.1. CFP-based method As shown in Fig. 1a, internal multiples can be predicted by convolving three datasets along virtual depth level zn. The involved process can be written as: −

n

Mn ðz0 ; z0 Þ ¼ P n ðz0 ; zn ÞδPn ðzn ; zn ÞP n−1 ðzn ; z0 Þ; n

ð4Þ

where P n−1 ðzn ; z0 Þ, denoted by the black ray-path in Fig. 1a, represents the input data with receivers back-propagated to depth level zn and all events reflected above zn being removed. Subscript n − 1 means that internal multiples related to depth level 0 to n − 1 have already been removed from the data in previous IME steps. P n ðz0 ; zn Þ, denoted by the blue ray-path, represents primaries with receivers at the surface and sources at depth level zn. They can be approximated by applying n reciprocity to P n−1 ðzn ; z0 Þ. δP− n (zn,zn), denoted by the red ray-path, represents input data with both sources and receivers inverse extrapolated to depth level zn. There, only the non-casual events reflected from the overburden (zn − 1 b z ≤ zn) are kept and subsequently have been time-reversed. Note that the prediction of internal multiples according to Eq. (4) involves two convolutions along depth level zn. Two kinds of datasets are required in this process, one is the half-redatumed data with sources at the surface and receivers below or vice versa, the other is the fully redatumed data with both sources and receivers located along the chosen depth level in the subsurface. In Berkhout and Verschuur (2005), these two kinds of datasets have been called common focus point (CFP) gathers and gridpoint gathers, respectively, therefore, we have referred to this method as the CFP-based IME method. Detailed implementation steps for this method have been described by Verschuur and Berkhout (2005). Here, we introduce the practical steps briefly. First, focusing operators are determined using a chosen time level and the corresponding NMO velocity. Next, CFP gathers are created and muted using these focusing operators. Then, gridpoint gathers are constructed by another focusing step to the CFP gathers (without mute), keeping the non-casual part and time-reversing it. These so-called gridpoint gathers have both their sources and receivers along the chosen level in the subsurface. Finally, internal multiples are predicted by convolving the muted CFP gathers and the time-reversed

68

J. Song et al. / Journal of Applied Geophysics 95 (2013) 66–76

predicted using only surface data (see Fig. 1b). The convolution process can be carried out along the surface, using the following expression: H

Mn ðz0 ; z0 Þ ¼ −P n ðz0 ; z0 Þ½δPn ðz0 ; z0 Þ Pn−1 nðz0 ; z0 Þ;

ð5Þ

with n

P n−1 ðz0 ; z0 Þ ¼ P n−1 ðz0 ; z0 Þ−δPn ðz0 ; z0 Þ:

(a)

(b)

(c) Fig. 1. Schematic illustration of layer-related internal multiple prediction. In the images, dotted lines represent virtual depth levels, wavefields and their corresponding raypaths that are indicated by the same colors. (a) CFP-based prediction, where internal multiples related to current depth level zn are predicted by convolving muted CFP gathers (black and blue ray-paths) with the time-reversed gridpoint gathers (red ray-path). (b) Surface-based prediction, where the surface data above and below zn are separated and convolved along the surface with the dashed part time-reversed. (c) Generalized EPSI, where the primary impulse response above and below zn (denoted by the fine ray-path meaning that no wavelet is included) is estimated in a full waveform inversion procedure, and multiples can be constructed along the surface using convolution.

gridpoint gathers. Note that the autocorrelation of the source signature is included in this prediction output, which results in wrong amplitudes of the predicted internal multiples. Therefore, adaptive subtraction is indispensable as an additional step in this process. Note that in Berkhout and Verschuur (2005) it was demonstrated that the used velocity model only has to be approximated, as errors will cancel during the convolution steps (because the gridpoint gathers are used in a time-reversed fashion). However, the use of the (approximate) focal operators facilitates the separation of events from above and below the chosen depth level.

2.2. Surface-based method Jakubowicz (1998) proposed a new IME implementation, in which no extrapolation is involved, meaning that internal multiples can be

ð6Þ

The structure of the convolution is very similar to CFP-based prediction, but now all sources and receivers are situated at the surface. P n ðz0 ; z0 Þ, denoted by blue ray-path, and δPn, denoted by the red dashed ray-path, represent primaries reflected from the interfaces below and above depth level zn, respectively, that have been explained in Eq. (3). Pn−1 nðz0 ; z0 Þ, denoted by the black ray-path, represents data with multiples related to depth level z0 to zn − 1 removed and also primaries above zn being removed, as expressed in Eq. (6). Note that P n ðz0 ; z0 Þ is the demultiple output as well as the multiple prediction input for the current depth level zn. In practical implementations, it can be approximated by Pn−1 nðz0 ; z0 Þ or estimated iteratively. This is similar to SRME, where the multiple-free data should theoretically be used as multiple prediction operators. In practice, an iteration process is carried out, where the multiple free estimate from the previous iteration is used as a prediction operator in the current iteration, starting with the measured data as the first multiple-free estimate (Berkhout and Verschuur, 1997). In Fig. 1b, the red dashed line refers to a time-reversed application of δPn, indicated by superscript (.)H in Eq. (5), which is used to compensate the traveltime errors of the predicted internal multiples introduced by convolving P n ðz0 ; z0 Þ with Pn−1 nðz0 ; z0 Þ. It provides the reflectivity information of the downward reflecting overburden. Correct internal multiple polarity for acoustic reflections is obtained via the minus sign. In van Borselen (2002) applications of this approach to marine field data were shown and in Hanitzsch et al. (2007) it was used for removing internal multiples on land data. Compared to the CFP-based method, where inverse wavefield extrapolation is involved, the implementation steps for the surface-based method are more simple. First, data above and below the chosen time level are separated. Then, internal multiple prediction is performed via convolution along the surface. Finally, due to the effect of the source signature autocorrelation, the predicted internal multiples are adaptively subtracted from the input data. 2.3. Generalized EPSI One main drawback of data-driven IME methods based on prediction and subtraction is that primaries are often not well preserved when internal multiples overlap with primaries. In Ypma and Verschuur (2013), surface-based internal multiple elimination has been redefined as a large-scale inversion process that leads to the generalized EPSI (estimating primaries by sparse inversion) method. This is based on the work of van Groenestijn and Verschuur (2009a) for the case of surface-related multiples, where the primary impulse responses are the unknowns, that are parameterized by spikes in the space–time domain. In addition, the source wavelet is also a parameter. For this method, Eq. (3) can be rewritten as H

n

P n−1 ¼ X n S þ δXn S−X n δXn P n−1 ;

ð7Þ

where all quantities are related to surface level z0. S represents the average source wavelet for all shots, and X n and δXn represent the impulse response related to reflectors below and above depth level zn, denoted by the blue and the red dashed ray-paths, respectively (see Fig. 1c). n The internal multiple term −X n δXn H P n−1 contains no effects of source signature autocorrelation that are inherent in multiples predicted by the CFP-based and surface-based IME methods (see Eqs. (4) and (5)).

J. Song et al. / Journal of Applied Geophysics 95 (2013) 66–76

Therefore, internal multiples constructed by generalized EPSI can be subtracted from original data directly. Impulse response X n , δXn and source signature S are unknowns that can be estimated by minimizing objective function 2 H n J ¼ ∑ ∑ P n−1 −X n S−δXn S þ X n δXn P n−1 ; ω

j;k

j;k

ð8Þ

69

the events above and below virtual depth level can be separated or partly separated (near offset). Since focusing operators and phase shift operators are allowed to contain traveltime errors which can be compensated in the subsequent process, sensitivity to the velocity errors of each IME method is further reduced, therefore, NMO velocity can be used as a good alternative for internal multiple prediction. 3. Synthetic data example

It is noted in the described internal multiple elimination methods that reflections in the data from above and below the chosen virtual depth level need to be separated. In the CFP-based method, this process becomes straightforward because the muting time function has already been defined by the involved focusing operator. Berkhout and Verschuur (2005) have shown that using the correct focusing operator or using one based on a homogeneous replacement model based on the NMO velocity gave roughly the same results. The reason is that traveltime errors of constructed CFP gathers and gridpoint gathers caused by back propagation using incorrect focusing operator are fully compensated for by the time-reversion of the gridpoint gathers. Therefore, the quality of the predicted internal multiples does not heavily rely on the choice of focusing operator, which indicates that layer-related CFP method can accommodate velocity errors. Surface-based method and generalized EPSI share the same event separation process which can be achieved by two approaches in the shot domain. The first approach is direct muting where muting time for each shot gather can be calculated using NMO velocity and the vertical traveltime between source position and chosen depth level. The rationale of direct muting is that near and medium offsets that contribute most to internal multiple prediction are easy to separate. The other approach is a more elaborate extrapolation-based method, which can be used for separating events that cross at far offsets. First, the data are inverse extrapolated to twice the chosen depth level with a phase shift operator in wavenumber–frequency domain using NMO velocity. Then, causal and non-causal events are separated by a muting around t = 0 and the two separated datasets are forward extrapolated to the surface again. In this way, events above and below the chosen depth level have been completely separated, even at far offsets. Due to the fact that extrapolation in wavenumber–frequency domain is reversible, errors caused by back propagation using inaccurate velocity can be corrected in the forward propagation process with the same velocity. Direct muting and extrapolation approach for event separation can be selected according to the condition whether primary events from above and below the virtual depth level cross each other seriously. Although velocity information is required in the process of event separation, the velocity does not have to be very accurate as long as

Lateral distance (m) 0

500

1000

1500

0

500 3500

1000 2500

1500 1500

2000

Fig. 2. Subsurface model for evaluating three IME methods.

Pressure - wave velocity (m/s)

2.4. Event separation involved in the IME methods

We present a 2D synthetic dataset to illustrate and evaluate the three feedback IME methods. The dataset is based on the subsurface model shown in Fig. 2, where the sea bottom and the top salt are two major internal multiple generating interfaces. Note that a potential reservoir is located below the salt layer. 121 shots have been modeled with the shot and receiver positions ranging from 0 to 1800 m with steps of 15 m. The surface has been chosen non-reflective, so that only primary and internal multiple reflections have been simulated. In Fig. 3a, one of the shot records is shown. Internal multiples indicated by the black arrows are generated by the downward reflection of water bottom. To address these strong internal multiples, a virtual depth level is chosen between the water bottom and the top salt. Internal multiples that pass this level at least four times can be predicted. For CFP-based method, event separation takes place after back-propagation of the data, such that events are optimally separated and a simple mute can be applied. Fig. 3b and c displays the muted CFP gather and the time-reversed grid-point gather corresponding to the input shot record. According to Eq. (4), internal multiples related to current depth level are predicted (see Fig. 3d), which show good traveltime consistency with the input data. Due to the serious overlap of the first two primaries, event separation based on extrapolation is performed for the surfacebased method and generalized EPSI. As shown in Fig. 3e and f, events above and below the virtual depth level are completely separated even at the far offsets. Based on the accurate separation results, internal multiples predicted with surface-based method and generalized EPSI are displayed in Fig. 3g and h. The surface-based method provides similar prediction results with the CFP-based method, while compared with the EPSI results whose amplitude and phase are both correct, source wavelet effects can be clearly observed. To investigate sensitivity of IME methods to event separation, internal multiples are now predicted with incomplete top reflections (see Fig. 3i) and complete lower reflections, note that there is hardly any difference between the newly obtained results (Fig. 3k, l) and the previous results, which confirms

Depth (m)

where the summation signs indicate a full summation over all sources, receivers and frequencies. For internal multiple elimination (i.e. n > 0), the initial value for one of the parameter sets needs to be known, e.g. the source signature S (as output of surface multiple removal). X n , δXn and S are updated alternately in an iterative procedure. First, the update of δXn for all frequencies is obtained by taking the derivative of J with respect to δXn. Next, the update is transformed to the time domain where sparseness is imposed by selecting the largest event(s) per trace. Furthermore, a suitable time window is applied, the window starting with a small size, which increases as the iteration number grows. Then, the sparse and time-windowed update is scaled and added to the current estimate of the impulse response. Next, the update of X n is calculated in the same way as δXn. Then, S is estimated as a filter obtained by least-square matching the impulse response X n þ n δXn to P n−1 þ X n δXn H P n−1 in the time domain for all shots simultaneously. Finally, after several tens of iterations, internal multiples related to current depth level n are constructed (see Fig. 1c) and can be subtracted from the input data without adaptation.

70

J. Song et al. / Journal of Applied Geophysics 95 (2013) 66–76

Offset (m)

Offset (m) 0

-1500 -1000 -500

Offset (m) 0

-1500 -1000 -500

0.5

0.5

0.5

0.5

1.0

1.0

1.5

Time (s)

0

Time (s)

0

1.0

1.5

(a)

(b)

(c)

Offset (m)

Offset (m)

0

-1500 -1000 -500

0

1.0

-1500 -1000 -500

(d) Offset (m) 0

-1500 -1000 -500 0

0.5

0.5

0.5

0.5

1.5

1.0

1.5

Time (s)

0

Time (s)

0

Time (s)

0

1.0

1.0

1.5

(e)

(f)

(g)

Offset (m)

Offset (m)

0

-1500 -1000 -500

0

1.0

-1500 -1000 -500

(h) Offset (m) 0

-1500 -1000 -500 0

0.5

0.5

0.5

0.5

1.5

1.0

1.5

(i)

Time (s)

0

Time (s)

0

Time (s)

0

1.0

1.0

1.5

(j)

0

1.5

Offset (m) -1500 -1000 -500

0

1.5

Offset (m) -1500 -1000 -500

Time (s)

-1500 -1000 -500

0

1.5

Time (s)

Offset (m) 0

0

Time (s)

Time (s)

-1500 -1000 -500

0

1.0

1.5

(k)

(l)

Fig. 3. Comparison of event separation in internal multiple prediction. (a) Original shot gather with internal multiples indicated by black arrows. (b) Muted CFP gather. (c) Time-reversed grid-point gather. (d) Internal multiples predicted by CFP-based method. (e) Separated top part of (a). (f) Separated lower part of (a). (g) Internal multiples predicted by surface-based method. (h) Internal multiples predicted by generalized EPSI. (i) Near offsets of (e). (j) Repeat of (f). (k) Internal multiples predicted by surface-based method using incomplete top part (i) and lower part (j). (l) Internal multiples predicted by generalized EPSI using incomplete separation results.

the important contribution of the near offsets to internal multiple prediction and the robustness of IME methods to imperfect event separation.

As shown in Fig. 3, internal multiples predicted by the CFP-based and surface-based methods are not correct in dynamics (i.e. amplitude), therefore, adaptive subtraction should be applied to remove internal

J. Song et al. / Journal of Applied Geophysics 95 (2013) 66–76

multiples from original shot records while direct subtraction can be used in the generalized EPSI. For this example, adaptive subtraction of internal multiples using L2 norm (Verschuur and Berkhout, 1997) and L1 norm (Guitton and Verschuur, 2004) is performed for the first two methods. Fig. 4b and c shows the estimated primaries by CFP-based method using L2 norm and L1 norm, compared with the original shot record in Fig. 4a, internal multiples indicated by the black arrow have been attenuated, however, no obvious advantage of adaptive subtraction using L1 norm over L2 norm has been observed for this synthetic data, which may because the amplitude of the internal multiples is higher than that of the overlapping primaries. Sharing the same IME strategy with CFP-based method, i.e. prediction and subtraction,

Offset (m) -1000

-500

-1500

-1000

-500

Offset (m) 0

-1500

0

0

0.5

0.5

0.5

1.0

Time (s)

0

1.5

1.0

1.5

(a) -1500

-1000

-1500

-1000

-500

Offset (m) 0

-1500 0

0.5

0.5

0.5

1.5

Time (s)

0

Time (s)

0

1.0

1.0

1.5

(d)

0

(c)

Offset (m) 0

-500

1.0

(b) -500

-1000

1.5

Offset (m)

Time (s)

surface-based method provides the similar primary estimation results through L2-norm subtraction (Fig. 4d) and L1-norm subtraction (Fig. 4e) respectively. Fig. 4f shows the primaries estimated by generalized EPSI, note that weak primaries indicated by the white arrows have been better recovered than that obtained by the prediction-andsubtraction methods. Internal multiples at around 1.5 s (see Fig. 4) are internal multiples generated by the top salt, which are not addressed in the current step. With the estimated primaries by each IME method as input, internal multiple elimination related to the new depth level (at 700 m in Fig. 2) is performed. After removing internal multiples for two virtual depth levels (at 300 m and 700 m), most of the visible internal multiple

Offset (m) 0

Time (s)

Time (s)

-1500

71

-1000

-500

0

1.0

1.5

(e)

(f)

Fig. 4. Comparison of three IME methods for synthetic data at prestack level. White and black arrows in the images indicate internal multiples and primaries respectively. (a) Input shot record with internal multiples. (b) Results of CFP-based IME using L2-norm subtraction. (c) Results of CFP-based IME using L1-norm subtraction. (d) Results of surface-based IME using L2-norm subtraction. (e) Results of surface-based IME using L1-norm subtraction. (f) Results of generalized EPSI.

72

J. Song et al. / Journal of Applied Geophysics 95 (2013) 66–76

energy has been attenuated. For convenience of comparison, primaries estimated by the prediction-and-subtraction methods using L2-norm and generalized EPSI are sorted into CMP gathers and compared with the aid of velocity semblance panels in Fig. 5. It shows that in the original CMP gather the primaries below 1.0 s are severely distorted by internal multiples. After recursive application of the CFP-based method, strong internal multiples at around 0.5 s and 1.0 s have been effectively suppressed in the CMP gather (see Fig. 5b) and its velocity panel (see Fig. 5f). Similar results are obtained with the surface-based method (see Fig. 5c and g). However, primary events indicated by the white arrows are not well recovered because adaptive subtraction has difficulties to handle the situation when primaries and multiples strongly interfere. The CMP gather after application of the generalize EPSI method is shown in Fig. 5d, where primary events are well preserved and they exhibit a more clear and continuous behavior. The advantage of generalized EPSI over the other two methods is also confirmed by comparing the primary energy indicated by white arrows in the associated velocity panel (see Fig. 5h). Finally, the stack sections before and after internal multiple removal by the three IME methods are displayed in Fig. 6. Note that water bottom-generated internal multiples, indicated by the top white arrows in Fig. Fig. 6a, interfere with events from the potential reservoir and leave a smearing multiple leakage in Fig. 6b and c, while in the generalized EPSI

Offset (m) 0

-1000

Next, we demonstrate the recursive application of these three feedback IME methods on a field dataset from the Mississippi Canyon area in the Gulf of Mexico, for which a stacked section is shown in Fig. 7. Since this data was acquired in deep water, we can deal with just primaries and internal multiples, without the influence of surface multiples, for our considered time window up to 3.5 s. In Fig. 7, four virtual depth levels are selected in between strong reflecting interfaces, where internal multiples (see yellow arrows) related to each level from z1 to z4 (see Fig. 7) are addressed in a layer-stripping manner. One CMP gather of the original field dataset is shown in Fig. 8a. Strong primary reflections above 2.3 s are supposed to generate noticeable internal multiples. Also note that the event labeled with the white arrow is a primary reflection from the salt bottom and related to a high NMO velocity in Fig. 8e. After recursive application of the three IME methods, removed internal multiples have been displayed in Fig. 8b, c and d for the same CMP location. The removed internal multiples exhibit a smearing characteristic for the CFP-based IME and the surface-based IME, while they appear very clear and continuous for generalized EPSI. It

0

Offset (m)

1000

-1000

0

Offset (m)

1000

-1000

0

0

0.5

0.5

0.5

0.5

1.0

1.0

1.5

2.5

1.5

2.0

2.5

(d)

Velocity (km/s) 3.0

1.5

2.0

2.5

Velocity (km/s) 3.0

1.5 0

0.5

0.5

0.5

0.5

1.0

1.5

1.5

(e)

Time (s)

0

Time (s)

0

Time (s)

0

1.0

1.0

2.0

2.5

3.0

1.0

1.5

1.5

(f)

1000

1.0

(c)

Velocity (km/s) 3.0

0

1.5

(b)

Velocity (km/s) 2.0

1.0

1.5

(a) 1.5

Time (s)

0

Time (s)

0

1.5

Time (s)

4. Field data example

Offset (m)

1000

Time (s)

Time (s)

-1000

result (see Fig. 6d), all internal multiples have been removed completely, leaving the primaries well preserved.

(g)

(h)

Fig. 5. Comparison of three IME methods for synthetic data at CMP level. White arrows indicate primaries in the CMP gathers or corresponding energy in the velocity panels. (a) CMP gather with internal multiples. (b) CMP gather after CFP-based IME. (c) CMP gather after surface-based IME. (d) CMP gather after generalized EPSI. (e) Velocity semblance panel corresponding to (a). (f) Velocity semblance panel corresponding to (b). (g) Velocity semblance panel corresponding to (c). (h) Velocity semblance panel corresponding to (d).

J. Song et al. / Journal of Applied Geophysics 95 (2013) 66–76

CMP location (m)

1000 1500

0

500

CMP location (m)

1000 1500

0

0.5

Time (s)

Time (s)

0.5

0

1.0

1.5

1.0

1.5

(a)

0

500

CMP location (m)

1000 1500

0

0.5

Time (s)

0

500

1.0

1.5

(b)

0

500

1000 1500

0.5

Time (s)

CMP location (m) 0

73

1.0

1.5

(c)

(d)

Fig. 6. Comparison of three IME methods for synthetic data at stack level. White arrows indicate the multiple artifacts. (a) Original stack section with internal multiples. (b) Stack section after CFP-based IME. (c) Stack section after surface-based IME. (d) Stack section after generalized EPSI.

is observed that primary energy indicated by the white arrow in Fig. 8f and g is distorted while in generalized EPSI (see Fig. 8h) it is well preserved. We can come to the conclusion that generalized EPSI can recover primaries even if they overlap with internal multiples, whereas leastsquares subtraction, based on the minimum energy criterion, used by the other two IME methods, cannot handle this situation properly. In Fig. 9 the stack sections before and after recursive application of the three IME methods together with their difference (i.e. the removed internal multiples) are shown. Compared with the original stack section

CMP location (km) 13

15

17

1.5

in Fig. 9a, internal multiples indicated by the yellow arrows in Fig. 9b–d have been obviously removed, which can be confirmed in the stack section of the removed internal multiples in Fig. 9f–h, representing CFP-based IME, surface-based IME and generalized EPSI, respectively. Particularly, removal of internal multiples that interfere with the salt bottom at around 2.7 s can be of great value for interpretation. Detailed comparison of the results by three IME methods can be performed using the fine grids superimposed on the removed internal multiples. After close inspection, one can see that the internal multiples removed by generalized EPSI in Fig. 9h are clearer than removed by the other two methods (see Fig. 9f and g), less distortion of primaries can be observed around 2.7 s and 3.2 s. Moreover, note that in the scaled figure of internal multiples in Fig. 9i, generalized EPSI provides more continuous events which result in an improved contrast of primary energy in the stack section in Fig. 9d.

5. Discussion

Time (s)

2.0

2.5

3.0

3.5 Fig. 7. Stack section of Mississippi Canyon field data. Dotted lines indicate chosen virtual time levels in between strong reflecting interfaces and the yellow arrows indicate evident internal multiples generated by the overburden.

For prediction of surface multiples, near offsets play an important role and the same holds for internal multiple prediction. Any missing trace will act as a noise source in the predicted internal multiples and generate artifacts. Therefore, missing near offset should be interpolated in a preprocessing step (see e.g. Kabir and Verschuur, 1995). This was done for the field data. For the synthetic data, we could use all modeled offsets. The main reason why the CFP-based and the surface-based IME methods perform less effective than generalized EPSI is that leastsquares subtraction has difficulties with separating primaries and multiples when they interfere. This may yield IME methods based on prediction and subtraction to fail in recursive application, as primary distortion or multiple leakage may lead to error accumulation. It is expected that this is especially important for input data with strong surface multiple leakage. In van Groenestijn and Verschuur (2008) it was noted that noise in the input data will always end up in the estimated primaries for the case of using a prediction and subtraction based methodology, while will not leak into the primaries estimated by EPSI. This could be another potential advantage of generalized EPSI over the CFP-based and surfacebased methods for internal multiple elimination. However, from a robustness point of view, prediction and subtraction-based methods are assumed more robust than inversion-based generalized EPSI for very noisy data, because the involved sparseness constraint imposed by picking spikes in generalized EPSI can be very sensitive to noise.

74

J. Song et al. / Journal of Applied Geophysics 95 (2013) 66–76

Offset (m) 0

Offset (m) 3000

-3000 1.5

-3000 1.5

Time (s)

2.5

0

Offset (m) 3000

-3000 1.5

2.0

2.5

2.5

3.0

3.0

3.5

3.5

3.5

3.5

(b) 2.5

Velocity (km/s) 3.0

1.5

2.0

2.5

2.5

(d)

Velocity (km/s) 3.0

1.5

2.0

Time (s)

2.0

1.5

1.5

2.0

2.5

Velocity (km/s) 3.0

1.5

2.0

Time (s)

1.5

2.0

(c)

2.5

2.5

3.0

3.0

3.0

3.5

3.5

3.5

3.5

(f)

2.0

2.5

3.0

2.5

3.0

(e)

1.5

2.0

Time (s)

Velocity (km/s)

3000

2.5

3.0

1.5

0

2.0

3.0

(a)

Time (s)

Offset (m) 3000

2.0

Time (s)

Time (s)

2.0

0

Time (s)

-3000 1.5

(g)

(h)

Fig. 8. Comparison of three IME methods for field data at CMP level. All CMP gathers are displayed in the same amplitude level, white arrows indicate distortion of primaries by three IME methods. (a) Original CMP gather with internal multiples. (b) Internal multiples removed by CFP-based IME method. (c) Internal multiples removed by surface-based IME method. (d) Internal multiples removed by generalized EPSI. (e) Velocity semblance panel corresponding to (a). (f) Velocity semblance panel corresponding to (b). (g) Velocity semblance panel corresponding to (c). (h) Velocity semblance panel corresponding to (d).

Furthermore, it is very difficult to estimate the initial source signature or impulse response in the presence of strong noise. The event separation in the surface-based and generalized EPSI methods deserves some attention in order to avoid separation artifacts leaking in the prediction process. For the CFP-based method the separation is facilitated by the involved back-propagation steps and the fact that the involved focusing operators are used as mute function in the CFP domain (Verschuur and Berkhout, 2005). In the 2D implementation of the three considered IME methods, 3D effects cannot be handled properly. Due to the absence of adaptive subtraction, generalized EPSI is potentially more sensitive to 3D effects (Ypma and Verschuur, 2013). These effects could be relieved by incorporating filters as another unknown parameters that need to be estimated in the inversion process of generalized EPSI. However, the fundamental solution is to extend the three IME methods to 3D. In Wu and Dragoset (2011), surface-based method has been extended to predict internal multiples for more than one generator simultaneously, which is able to handle 3D acquisition geometry and predict internal multiples at true azimuth. First 3D data examples of this approach look very promising and the algorithm can be adopted by generalized EPSI for 3D extension. In terms of computation costs for applying IME for one virtual depth level, the CFP-based method involves two extrapolation and

two spatial convolution steps, being similar expensive as four times a SRME process. The surface-based method includes two spatial convolutions, being approximately twice as expensive as SRME. Finally, generalized EPSI is much more computationally intensive than the former two methods because it takes several tens of iterations, with each iteration at the cost of two surface-based IME processes (i.e. four SRME processes). 6. Conclusions In this paper, recursive application of three feedback internal multiple elimination methods has been compared using synthetic data and field data. It turns out that the three IME methods are very robust to velocity errors and imperfect event separation which plays an important role in internal multiple prediction. Based on the prediction and subtraction strategy, the CFP-based method and the surface-based method can provide a similar result, although they are different in the internal multiple prediction step in which multidimensional convolutions are carried out along a virtual depth level and along the surface level, respectively. In the case of internal multiples interfering with primaries, these two methods cannot recover primaries very well, as multiple leakage and primary distortion

J. Song et al. / Journal of Applied Geophysics 95 (2013) 66–76

CMP location (km) 17

1.5

2.5

3.0

1.5

2.5

2.5

1.5

13

15

1.5

3.5

(d)

13

15

CMP location (km) 17

1.5

2.5

3.0

3.0

3.5

3.5

3.5

3.5

(f) 15

CMP location (km)

17

1.5

2.5

15

CMP location (km)

17

1.5

2.0

Time (s)

Time (s)

2.0

13

2.5

3.0

3.0

3.5

3.5

13

15

17

2.0

Time (s)

13

17

(h)

(g)

CMP location (km)

15

2.5

3.0

1.5

13

2.0

3.0

(e)

17

3.0

2.0

2.5

15

2.5

CMP location (km) 17

Time (s)

Time (s)

2.5

13

(c)

2.0

2.0

1.5

3.5

CMP location (km) 17

CMP location (km) 17

2.0

(b)

CMP location (km) 15

15

3.0

3.5

13

13

2.0

(a)

Time (s)

CMP location (km) 17

3.0

3.5

1.5

15

2.0

Time (s)

Time (s)

2.0

13

Time (s)

15

Time (s)

1.5

13

Time (s)

CMP location (km)

75

2.5

3.0

3.5

(i) Fig. 9. Comparison of three IME methods for field data at stack level. Arrows in the images indicate position of internal multiples. For ease of comparison, fine grids have been superimposed on the figures of removed internal multiples. (a) Original stack section with internal multiples. (b) Stack section after CFP-based IME. (c) Stack section after surface-based IME. (d) Stack section after generalized EPSI. (e) Repeat of original stack section. (f) Removed internal multiples by CFP-based IME. (g) Removed internal multiples by surface-based IME. (h) Removed internal multiples by generalized EPSI. (i) Scaled figure of the removed internal multiples by three IME methods. All stack sections are displayed in the same amplitude level.

have been observed in the data examples. Such leakage may result in error accumulation in a recursive application. The generalized EPSI method estimates primaries via an inversion process, where the estimated primaries, together with their corresponding multiples, explain the total input data and adaptive subtraction is avoided. Internal multiples estimated by this method are much

more clear and spatially continuous, which leads to a reduced multiple leakage and a better primary contrast than with the prediction and subtraction IME methods. From a cost perspective, the CFP-based method is twice expensive as the surface-based method, and the generalized EPSI is several tens of times more expensive than the surface-based method.

76

J. Song et al. / Journal of Applied Geophysics 95 (2013) 66–76

In general, the three feedback methods are quite effective in internal multiple elimination. For the CFP-based and the surface-based IME methods, priority should be given to developing more effective and robust subtraction methods, while for generalized EPSI, reducing computation cost is the main future consideration. Furthermore, a full 3D implementation of the three IME methods will be required for application to realistic datasets. Acknowledgments This research is supported by the National Natural Science Foundation of China (Grant No. U1262207) and the National Science and Technology of Major Projects of China (Grant No. 2011ZX05023-005-005). The authors thank the China Scholarship Council for sponsoring Jiawen Song for his one-year study at Delft University of Technology. References Abma, R., Kabir, N., 2005. 3D interpolation of irregular data with a POCS algorithm. SEG, Soc. Expl. Geophys., Expanded abstracts, Houston, pp. 2150–2153. Araujo, F.V., Weglein, A.B., Carvalho, P.M., Stolt, R.H., 1994. Inverse scattering series for multiple attenuation: an example with surface and internal multiples. SEG, Soc. Expl. Geophys., Expanded abstracts, Los Angeles, pp. 1039–1041. Berkhout, A.J., 1982. Seismic Migration, Imaging of Acoustic Energy by Wave Field Extrapolation, A: Theoretical Aspects, second edition. Elsevier. Berkhout, A.J., Verschuur, D.J., 1997. Estimation of multiple scattering by iterative inversion, part I: theoretical considerations. Geophysics 62 (5), 1586–1595. Berkhout, A.J., Verschuur, D.J., 1999. Removal of internal multiples. SEG, Soc. Expl. Geophys., Expanded abstracts, Houston, pp. 1334–1337. Berkhout, A.J., Verschuur, D.J., 2005. Removal of internal multiples with the commonfocus-point (CFP) approach: Part 1 — explanation of the theory. Geophysics 70 (3), V45–V60. Coates, R.T., Weglein, A.B., 1996. Internal multiple attenuation using inverse scattering: results from prestack 1 and 2-D acoustic and elastic synthetics. SEG, Soc. Expl. Geophys., Expanded abstracts, pp. 1522–1525. Guitton, A., Verschuur, D.J., 2004. Adaptive subtraction of multiples using the L1-norm. Geophysical Prospecting 52, 27–38.

Hanitzsch, C., van Veen, L.J., Ali, J., van Borselen, R.G., 2007. Advanced multiple elimination: application to a complex seismic land data set. SEG, Soc. Expl. Geophys., Expanded abstracts, San Antonio, pp. 2525–2529. Jakubowicz, H., 1998. Wave equation prediction and removal of interbed multiples. SEG, Soc. Expl. Geophys., Expanded abstracts, New Orleans, pp. 1527–1530. Kabir, M.M.N., Verschuur, D.J., 1995. Restoration of missing offsets by parabolic Radon transform. Geophysical Prospecting 43 (3), 347–368. Kinneging, N.K., Budejicky, V., Wapenaar, C.P.A., Berkhout, A.J., 1989. Efficient 2D and 3D shot record redatuming. Geophysical Prospecting 37 (5), 493–530. Matson, K.H., Paschal, D., Weglein, A.B., 1999. A comparison of three multipleattenuation methods applied to a hard water-bottom dataset. The Leading Edge 18 (1), 120–126. van Borselen, R.G., 2002. Fast-track, data-driven interbed multiple removal: a north sea data example. EAGE, Eur. Ass. of Geosc. and Eng., Expanded abstracts, p. F040. van Borselen, R.G., Fooks, G., Brittan, J., 2003. Target-oriented adaptive subtraction in data-driven multiple removal. The Leading Edge 22, 340–343. van Groenestijn, G.J.A., Verschuur, D.J., 2008. Towards a new approach for primary estimation. SEG, Soc. Expl. Geophys., Expanded abstracts, Las Vegas, pp. 2487–2491. van Groenestijn, G.J.A., Verschuur, D.J., 2009a. Estimating primaries by sparse inversion and application to near-offset data reconstruction. Geophysics 74, A23–A28. van Groenestijn, G.J.A., Verschuur, D.J., 2009b. Estimation of primaries and near offsets by sparse inversion: marine data applications. Geophysics 74, R119–R128. Verschuur, D.J., 2006. Seismic Multiple Removal Techniques — Past, Present and Future. EAGE Publications BV. Verschuur, D.J., Berkhout, A.J., 1997. Estimation of multiple scattering by iterative inversion, part II: practical aspects and examples. Geophysics 62 (5), 1596–1611. Verschuur, D.J., Berkhout, A.J., 2005. Removal of internal multiples with the commonfocus-point (CFP) approach: Part 2 — application strategies and data examples. Geophysics 70 (3), V61–V72. Verschuur, D.J., Berkhout, A.J., Wapenaar, C.P.A., 1992. Adaptive surface-related multiple elimination. Geophysics 57 (9), 1166–1177. Verschuur, D.J., Berkhout, A.J., Matson, K.H., Weglein, A.B., Young, C.Y., 1998. Comparing the interface and point scatterer methods for attenuating internal multiples: a study with synthetic data — part I. SEG, Soc. Expl. Geophys., Expanded abstracts, New Orleans, pp. 1519–1522. Weglein, A.B., Gasparotto, F.A., Carvalho, P.M., Stolt, R.H., 1997. An inverse scattering series method for attenuating multiples in seismic reflection data. Geophysics 62, 1975–1989. Wu, Z.J., Dragoset, B., 2011. Robust internal multiple prediction algorithm. SEG, Soc. Expl. Geophys., Expanded abstracts, San Antonio, pp. 3541–3545. Ypma, F.H.C., Verschuur, D.J., 2013. Estimating primaries by sparse inversion, a generalized approach. Geophysical Prospecting 61. http://dx.doi.org/10.1111/ j.1365-2478.2012.01095.x (in press).