Robust removal of free surface and internal multiples from the reflection response

(1)

Delft University of Technology

Robust removal of free surface and internal multiples from the reflection response

Slob, E. DOI 10.3997/2214-4609.201900830 Publication date 2019 Document Version Final published version Published in

81st EAGE Conference and Exhibition 2019

Citation (APA)

Slob, E. (2019). Robust removal of free surface and internal multiples from the reflection response. In 81st EAGE Conference and Exhibition 2019 [Tu_R04_07] EAGE. https://doi.org/10.3997/2214-4609.201900830 Important note

To cite this publication, please use the final published version (if applicable). Please check the document version above.

Copyright

Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy

Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.

This work is downloaded from Delft University of Technology.

(2)

Green Open Access added to TU Delft Institutional Repository

‘You share, we take care!’ – Taverne project

https://www.openaccess.nl/en/you-share-we-take-care

Otherwise as indicated in the copyright section: the publisher

is the copyright holder of this work and the author uses the

Dutch legislation to make this work public.

(3)

81st EAGE Conference & Exhibition 2019 3-6 June 2019, London, UK

Tu_R04_07

Robust Removal of Free Surface and Internal Multiples

From the Reflection Response

E. Slob1_*

1_{Delft University of Technology}

Summary

For a long time free surface related multiple reflection were removed before further processing was done. Based on the recently developed Marchenko type redatuming theory a theory was developed that is able to remove surface related and internal multiples during the redatuming step. This requires some model information that is avoided by separating the redatuming step from the multiple elimination step. The iterative solution of this method does not always converge. I propose to use a conjugate algorithm to solve the problem. I also show that the problem can be cast in a form with a self-adjoint operator that allows a faster solution. The second feature that reduces the computational cost is to use the filters that have been computed for a certain travel time as initial estimate for a following travel time. I show with a one-dimensional example that exploiting these two new aspects reduces the computational cost dramatically.

(4)

Introduction

Surface related multiple elimination schemes exist for a long time and are being used in practice. Inter-nal multiple elimination schemes lag behind, although the theory for a method that promises to remove surface and internal multiples was given more than 20 years ago (Weglein et al., 1997). The approximate prediction and subtraction of internal multiples is possible as well (Löer et al., 2016; ten Kroode, 2002). Since several years, filtering technology has been developed based on the theory of Marchenko reda-tuming to eliminate internal and free surface multiples and ghost during the redareda-tuming step (Wapenaar et al., 2014; Singh et al., 2017; Ravasi, 2017). Combining free surface and internal multiple elimination using the proposed iterative scheme is not feasible in many situations (Dukalski and de Vos, 2018). Re-cently this theory has been adapted to separate internal multiple elimination from the redatuming step (van der Neut and Wapenaar, 2016).

Here, the elimination of free surface and internal multiples are combined and separated from the reda-tuming step. This ensures that the process is automated and unsupervised. Information on acquisition geometry and the source time signature are required. The proposed solution is to modify the equation that creates a self-adjoint integral operator for which conjugate gradient iterative methods are available. These methods are unconditionally convergent. I show how this equation is obtained and how the prob-lem can be solved. I show numerical convergence and how the iterative process can be accelerated by using the result obtained at an earlier travel time as initial estimate for a following travel time. i give a 1D numerical example to illustrate the effectiveness and efficiency of the proposed method and end with conclusions.

Theory for filtering primary reflections from the impulse reflection response

To develop the theory time is denoted t and a point in space is specified with the position vector x. The acquisition surface is located at depth level ∂ D0 defined by z = z0. The surface can be reflection

free, pressure free, or rigid and these choices are defined by the parameter r as r = 0, r = −1, or r = 1, respectively. The impulse reflection response, when the acquisition surface is reflection free, is denoted R₀(x₀, x0₀,t) and R(x₀, x0₀,t) otherwise. They are related to each other as (Wapenaar et al., 2004)

R(x0₀, x₀,t) − R0(x00, x0,t) = r

Z

∂ D0

R(x0₀, x,t) ∗ R0(x, x0,t)dx,

where ∗ denotes temporal convolution. This is an equation that can be used to find R0 from R as is

done in free surface multiple elimination (Verschuur et al., 1992). This equation can be combined with filtering theory and rewritten as

U−(x₀, x0₀,t, τ) + h−(x0₀, x₀,t, τ) − Z ∂ D0 R(x0₀, x,t) ∗ [h+(x, x₀,t, τ) − rh−(x, x₀,t, τ)]dx = R(x₀, x0₀,t), (1) −U+(x₀, x0₀, −t, τ) + h+(x0₀, x₀,t, τ) − Z ∂ D0 R(x0₀, x, −t) ∗ [h−(x, x₀,t, τ) − rh+(x, x₀,t, τ)]dx = − δ (x0− x00, y0− y00,t) + rR(x0, x00, −t). (2)

A derivation of equations (1) and (2) that combines the ideas in van der Neut and Wapenaar (2016) and Singh et al. (2017) can be found in Zhang and Slob (2019). In these equations, the filters h±are non-zero only in the interval 0 < t < τ and τ is a free time parameter. The function U−(x₀, x0₀,t, τ) in equation (1) is zero for t < τ and an unmodified part of the impulse reflection response R for t > τ. The function U+(x₀, x0₀,t, τ) in equation (2) is causal and hence zero for t < 0 for any τ > 0. We reduce the time window to 0 < t < τ to evaluate equations (1) and (2) and find

h−(x0₀, x₀,t, τ) − Z ∂ D0 R(x0₀, x,t) ∗ [h+(x, x₀,t, τ) − rh−(x, x₀,t, τ)]dx = R(x₀, x0₀,t), (3) h+(x0₀, x₀,t, τ) − Z ∂ D0 R(x0₀, x, −t) ∗ [h−(x, x₀,t, τ) − rh+(x, x₀,t, τ)]dx = 0. (4)

81stEAGE Conference & Exhibition 2019 3–6 June 2019, London, UK

(5)

The significance of these equations is that once the filters are known, the primary reflections that have a two-way travel time from x0₀into the subsurface and back to x₀are found as Rr(x₀, x0₀, τ) = h−(x0₀, x₀,t ↑

τ , τ ), and Rr(x₀, x0₀, τ) denotes the new data set containing only primary reflections at the two-way travel

times and with amplitudes that are compensated for local transmission effects between the surface points x₀, x0₀ and the subsurface reflection point. To fill the new primary reflections dataset for all values of τ, equations (3) and (4) must be solved for each value of τ in the interval 0 < t < τ. It is therefore important to find a fast way to solve the equations. When r = 0 and R = R0, these equations can be solved for h±

by a Neumann-type iterative scheme (van der Neut and Wapenaar, 2016) that has been proposed to solve Marchenko redatuming equations (Wapenaar et al., 2013, 2014), in which case the series always converges. When r = ±1 the series solution is used as well (Singh et al., 2017; Zhang and Slob, 2019), but does not always converge. This has been analysed in detail for 1D Marchenko redatuming problems by Dukalski and de Vos (2018). They showed that when the subsurface reflectors are strong and/or the travel times associated with redatuming depth levels grow, the series expansion does not converge anymore. A more robust numerical method is required to solve the problem.

Numerical solutions using conjugate gradient techniques

Integral equations are often solved using numerical techniques and the family of conjugate gradient methods have been investigated over many years. Here we briefly introduce them and start with an operator notation that summarises equations (3) and (4) in a Hilbert space framework as

L u = f ,

in which u denotes the unknown filters and f denotes the impulse reflection response. The operatorL is a bounded linear operator, acting on u, with bounded inverseL−1, which maps a Hilbert space H onto itself. The space is equipped with an inner product h·, ·i and a norm || · || defined as

hu(x, x0,t, τ), v(x, x0,t, τ)i = Z ∂ D0 Z τ 0 u(x, x0,t, τ)v(x, x0,t, τ)dtdx, ||u(x, x0,t, τ|| = hu(x, x0,t, τ), u(x, x0,t, τ)i1/2,

from which it is clear that the inner product and norm are defined for each point x0and each evaluation time instant τ. The adjoint operatorL†, associated withL is defined as the operator satisfying

hL u(x,x0,t, τ), v(x, x0,t, τ)i = hu(x, x0,t, τ),L†v(x, x0,t, τ)i.

We also define a residual, r, after n iterations, given by rn= f −L u and the quantity ||r|| is minimised.

With these definitions we can give the iterative scheme based on conjugate gradients. The general scheme that can be used is given by (van den Berg, 1991)

w0= 0, wn= wn−1+ T rn−1 hLT rn−1, rn−1i , un= un−1+ wn ||L wn||2 , r0= f −L u0, rn= rn−1−_||_{L w}L wn n||2 , for n > 0, (5) in whichT is an operator of choice that we discuss now before we discuss the options for u0. We can see

that in equation (3) time convolutions occur and in equation (4) time correlations. This implies that the operatorL is not self-adjoint. In that case we need to take T = L†and hLT rn−1, rn−1i = ||L†rn||2

and the method is similar to the method given by Foll (1971). Better performance is obtained when we can work with a self-adjoint operator and we takeT = I as the unit operator in which case we have T rn−1= rn−1 and hL†rn−1,T rn−1i = hL rn−1, rn−1i. In both cases unconditional convergence of the

scheme has been proven. A self-adjoint operator is obtained by summing equations (3) and (4) in case r= −1 and subtracting one from the other when r = 1. Let’s take r = −1, which means the surface is pressure free. Summing the equations yields

h(x0₀, x₀,t, τ) −

Z

∂ D0

[R(x0₀, x, −t) + R(x0₀, x,t)] ∗ h(x, x₀,t, τ)dx = R(x₀, x0₀,t), (6) in which h = h++ h− and because convolution and correlation each other other’s adjoint operation, equation (6) has a self-adjoint operator. When h is found, we use equation (3) for t ↑ τ to compute its contribution to Rr. This comes at almost no extra cost.

(6)

Example 0 1 2 3 4 recording time (s) -1 -0.5 0 0.5 1 reflection strength

impulse reflection response

0 0.2 0.4 0.6 0.8 1 recording time (s) -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 reflection strength primary reflectivity

Figure 1 The impulse reflection response convolved with a 30 Hz Ricker wavelet (left) and the local primary reflectivity at two-way travel time (right).

200 400 600 800 1000 1200 1400 number of iterations 10-4 10-3 10-2 10-1 100 normalised error = 1 s 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Vertical traveltime (s) 0 20 40 60 80 100 120 Number of iterations

Using previous result

T=I T=L* 0 0.2 0.4 0.6 0.8 1 recording time (s) -0.5 0 0.5 h - (t, ) = 1 s 0 0.2 0.4 0.6 0.8 1 1.2 Recording time (s) 0 0.2 0.4 0.6 Vertical traveltime (s) -0.5 0 0.5

Figure 2 The number of iterations (top row) and the upgoing par of the filter, h−(t, τ) (bottom row) for truncation time τ = 1 s (left column) and τ = 4 s (right column); for three different solution strategies using conjugate gradient technique.

We solve equation (6) with the choiceT = I and T = L†and we solve the coupled equations (3)-(4) for whichT = L†is the only option. The model has a pressure free surface and three strong reflecting boundaries. The impulse reflection response dressed with a 30 Hz Ricker wavelet is shown in the left plot of Figure 1 and the corresponding primary reflectivity is shown in the right plot as a function of two-way travel time. Figure 2 shows the number of iterations for three different solutions of the same problem in the top left plot and the corresponding result for h−(t, τ) in the bottom left plot. They are shown for truncation time τ = 1 s. The colour scheme is the same of both plots in the left column. The

(7)

solution of equation (6) withT = I in shown in black and with T = L†in red, while the solution of equation (3)-(4) is shown in blue. In the bottom row the result for h−(t, τ) is shown and equation (3) has been used after the solutions of equation (6) is found. It is clear that the solution of equation (6) is preferred withT = I and that when T = L†it does not really matter which equation is solved. The most important aspect here is that it takes many iterations and when the primaries need to be filtered out this has to be done for every time point with a time step defined by the maximum frequency in the data. An interesting feature of the filter is that it changes only when an increment in the truncation time τ leads to τ + δ τ that crosses the two-way travel time of a reflector. This is shown in the bottom right plot, where the horizontal lines indicate the vertical travel time to the reflectors and the slanted line indicates the line along which the primaries are found and stored in Rr. Every time a reflector is passed

all previous events remain unchanged and new ones need to be found. It is therefore beneficial to use h±(t, τ) for all t as an initial estimate h±₀(t, τ + δ τ) as indicated in equation (5). The advantage is shown in the top right plot. The average number of iterations for all time points is four whenT = I and it is 20 when equations (3)-(4) are solved. This makes solving the free surface and internal multiples in one step a feasible approach.

Conclusions

I have shown that the combined free surface and internal multiple removal is possible and feasible. The known and new formulations have different conjugate gradient iterative solution methods that are unconditionally convergent. The fact that previous results of the filters can be used as initial estimate of the new filter makes the scheme attractive because it reduces the number of iterations dramatically. The 1D numerical examples demonstrates the advantages of the two methods.

References

van den Berg, P.M. [1991] Iterative schemes based on minimization of a uniform error criterion. Progress in Electromagnetics Research-PIER, 5, 27–65.

Dukalski, M. and de Vos, K. [2018] Marchenko inversion in a strong scattering regime including surface-related multiples. Geophysical Journal International, 212(2), 760–776.

Foll, J.L. [1971] An iterative procedure for the solution of linear and non-linear equations. Conference on Applications of Numerical Analysis. Lecture Notes in Mathematics, vol. 228. Springer, Berlin. ten Kroode, F. [2002] Prediction of internal multiples. Wave Motion, 35(4), 315–338.

Löer, K., Curtis, A. and Meles, G.A. [2016] Relating source-receiver interferometry to an inverse-scattering series to derive a new method to estimate internal multiples. Geophysics, 81(3), Q27–Q40. van der Neut, J. and Wapenaar, K. [2016] Adaptive overburden elimination with the multidimensional

Marchenko equation. Geophysics, 81(5), T265–T284.

Ravasi, M. [2017] Rayleigh-Marchenko redatuming for target-oriented, true-amplitude imaging. Geo-physics, 82(6), S439–S452.

Singh, S., Snieder, R., van der Neut, J., Thorbecke, J., Slob, E.C. and Wapenaar, K. [2017] Accounting for free-surface multiples in Marchenko imaging. Geophysics, 82(1), R19–R30.

Verschuur, D.J., Berkhout, A.J. and Wapenaar, C.P.A. [1992] Adaptive surface-related multiple elimina-tion. Geophysics, 57(9), 1166–1177.

Wapenaar, K., Broggini, F., Slob, E. and Snieder, R. [2013] Three-Dimensional Single-Sided Marchenko Inverse Scattering, Data-Driven Focusing, Green’s Function Retrieval, and their Mutual Relations. Physical Review Letters, 110(8), 084301.

Wapenaar, K., Thorbecke, J. and Draganov, D. [2004] Relations between reflections and transmission responses of 3-D inhomogeneous media. Geophysical Journal International, 156, 179–194.

Wapenaar, K., Thorbecke, J., van der Neut, J., Broggini, F., Slob, E. and Snieder, R. [2014] Marchenko Imaging. Geophysics, 79(3), WA39–WA57.

Weglein, A.B., Gasparotto, F.A., Carvalho, P.M. and Stolt, R.H. [1997] An inverse scattering series method for attenuating multiples in seismic reflection data. Geophysics, 62, 1975–1989.

Zhang, L. and Slob, E. [2019] Free-surface and internal multiple elimination in one step without adaptive subtraction. Geophysics, 84(1), A7–A11.