One-step estimation of focusing operators for 3D seismic data

(1)

One-step Estimation of

(2)

(3)

One-step Estimation of

Focusing Operators for 3D Seismic Data

PROEFSCHRIFT

ter verkrijging van de graad van doctor

aan de Technische Universiteit Delft,

op gezag van de Rector Magnificus prof. dr. ir. J.T. Fokkema,

voorzitter van het College voor Promoties,

in het openbaar te verdedigen

op 28 juni 2007 om 12.30 uur

door

Marinus Joseph VAN DE RIJZEN

(4)

Toegevoegd promotor: Dr. ir. D.J. Verschuur

Samenstelling promotiecommissie:

Rector Magnificus, voorzitter,

Prof. dr. ir. A. Gisolf, Technische Universiteit Delft, promotor

Dr. ir. D.J. Verschuur, Technische Universiteit Delft, toegevoegd promotor Prof. dr. ir. A.J. Berkhout, Technische Universiteit Delft,

Prof. dr. ir. P.M.J. Van den Hof, Technische Universiteit Delft,

Dr. E. Landa, Universit´e de Pau,

Dr. ir. D.T. van Daalen, Shell International Exploration and Production, Prof. dr. ir. P.M. van den Berg, Technische Universiteit Delft,

Prof. dr. ir. J.J.M. Braat, Technische Universiteit Delft, reserve lid

ISBN 978-90-9022021-5

Copyright c 2007, by M.J. van de Rijzen, Laboratory of Acoustical Imaging and Sound

Control, Faculty of Applied Sciences, Delft University of Technology, Delft, The Nether-lands.

All rights reserved. No part of this publication may be reproduced, stored in a retrieval system or transmitted in any form or by any means, electronic, mechanical, photocopying, recording or otherwise, without the prior written permission of the author M.J. van de Rijzen, Faculty of Applied Sciences, Delft University of Technology, P.O. Box 5046, 2600 GA, Delft, The Netherlands.

SUPPORT

(5)

1

Introduction

1.1 Short introduction to seismic exploration

1.1.1 The goal of seismic exploration

In the oil and gas industry, seismic exploration is generally the first step in the production cycle of hydrocarbons. Its goal is to determine a structural map of the earth that can be used by geologists to locate potential oil and gas resources. In most seismic experiments, a source located at the surface generates a wavefield that propagates through the earth where it is reflected by the different heterogeneities present in the earth. This reflected wavefield is measured by a set of receivers, which are generally placed at the surface, making-up a common shot-gather. Repeating this experiment for different source-positions constitutes the seismic dataset. The seismic dataset can be used to obtain the desired reflectivity image of the earth by mapping back the energy in the reflected wavefield to the positions of the reflectors (Figure 1.1). This process is called imaging, also known as migration.

In order to map the energy onto the proper reflector positions, an accurate estimate of the velocity model of the earth is needed to account for wave propagation in the earth. This velocity model is a priori not known.

1.1.2 Conventional Seismic Imaging

To overcome the problem of the unknown velocity model as discussed above, most conventional imaging techniques incorporate a velocity model estimation step in the imaging procedure, resulting in a procedure called Migration Velocity Analysis (MVA) [Stork, 1992; Kosloff et al., 1996], in which three basic steps can be distin-guished (Figure 1.2a):

(10)

−2 0 2 0 0.5 1 1.5 2 2.5 3 offset (km) time (s) −2 0 2 0 0.5 1 1.5 2 2.5 3 offset (km) time (s) 0 2 4 6 8 0 0.5 1 1.5 2 2.5 3 surface location (km) depth (km) 1500 2000 2500 3000 3500 4000 4500 5000 5500 −2 0 2 0 0.5 1 1.5 2 2.5 3 offset (km) time (s) 0 2 4 6 8 0 0.5 1 1.5 2 2.5 3 surface location (km) depth (km)

Figure 1.1: An example of seismic imaging. The input are the seismic data consisting of

common-shot gathers (left upper-part) and the velocity model (left lower-part). The arrows in top of the velocity model indicate the source positions of the sources in the common-shot gathers. When the correct velocity model is used in the imaging procedure, the seismic image (right part) will reveal the structure of the earth.

• The correctness of the velocity model is verified by looking at the flatness of Common Image Gathers (CIGs), which are a by-product of the migration. • The non-flatness observed in the CIGs as function of offset is used to update

the velocity model, which can then used as new velocity model in the first step. The steps described above are repeated in an iterative procedure until the events in the CIGs are flat.

Depending on the actual implementation, other domains than offset-CIGs can be used to analyze the correctness of the velocity model. Biondi and Symes [2004] have used angle domain common image gathers (ADCIGs) in which the offsets are transformed to angles. Also the way the velocity model is updated may differ in dif-ferent implementations. Most MVA-based approaches, however, have the structure as described above.

Whereas good results have been obtained with these methods, there are some issues in these approaches:

(11)

1.1 Short introduction to seismic exploration 3 Initial Velocity Model Seismic Data Final Velocity Model N Prestack Migration Analyse CIGs CIGs Flat ? Update Velocity Model Y of Correct Operators Tomographic Inversion Seismic Data _OperatorsInitial

Estimate Correct Operators

Velocity Model

(a) (b)

Figure 1.2: Flowcharts of two different approaches to estimate a velocity model. a) The

Migration Velocity Analysis (MVA)-approach. b) The Common Focal Point (CFP)-method.

• The non-flatness in the CIG is measured in the depth-offset domain whereas the velocity model update is performed in the time-midpoint domain. Therefore, the update may be non-optimal, thereby hampering convergence.

• In each iteration of the procedure an expensive migration is performed con-tributing to the high cost of the procedure.

1.1.3 Seismic Imaging using the CFP-method

A solution to the problems mentioned above is the Common Focal Point (CFP)-method introduced by Berkhout [1997a,b]. In this (CFP)-method the velocity model es-timation is separated from the imaging procedure (Figure 1.2b). He showed that seismic imaging can also be performed without a velocity model by introducing so-called focusing operators, where each focusing operator describes the one-way wave propagation from a point in the subsurface to a grid of points at the acquisition sur-face. Often these focusing operators are parameterized as traveltime and amplitude functions. The so called focusing operators are essentially the Green’s functions of the wave equation in an inhomogeneous medium.

(12)

the image point, hence the name Common Focal Point (CFP).

Whereas in the MVA-based methods the velocity model is unknown, in the CFP-method the traveltimes defining the focusing operators are unknown. Therefore, the problem shifts from estimating a velocity model to estimating focusing operators, for different points in the subsurface.

It can be shown that applying only one focusing step to the data transforms the data to a domain that is very suitable for estimating the focusing operators. In this new domain, called the CFP-gather, the ’Principle of equal traveltime’ can be formulated [Thorbecke, 1997]. This principle states that in case the correct focusing operator has been used to create a CFP-gather, the focusing operator and the event in the CFP-gather corresponding to the image point should be kinematically the same, because they both describe one-way wave propagation from the image point to points at the acquisition surface. The residual traveltimes observed between the target event in the CFP-gather and the focusing operator can be used to update the focusing operator. When we take another look at the issues as they appear in the MVA-based methods, we see that using CFP-technology:

• No velocity model is needed, so no biases occur due to a choice of parameter-ization of the velocity model.

• The residual traveltime error is defined in a true time domain, not a time-converted depth domain, and the operators are defined in the time-domain. This makes the updating procedure much more direct.

• The procedure of creating CFP-gathers for a number of target reflectors is cheaper than a full migration, thereby lowering the cost compared to MVA.

1.2 Problem statement

Whereas the CFP-method has proven to be a valuable technique with a number of different applications, there are still some issues.

The first problem is caused by the updating of the operators. In the iterative method proposed by Berkhout [1997a] and further described by Thorbecke [1997] and Bolte [2003], the residual traveltimes observed between the operators and the CFP-gathers are divided by two and added to the operators. This leads to a new set of opera-tors which can then be used again to create new CFP-gathers. Residual traveltimes between the new CFP-gathers and the updated operators imply that the operators do not obey the ’Principle of equal traveltime’. When they do not fulfill this re-quirement, the operators are updated again and this procedure is repeated until the ’Principle of equal traveltime’ is obeyed.

(13)

1.3 Research Objectives 5

Furthermore, because of the iterative nature of the updating procedure described above, the residual traveltimes between the CFP-gathers and the operators need to be determined a number of times, which is undesirable.

Besides the operator updating procedure, there is a more fundamental problem when the CFP-method is applied on sparse 3D seismic datasets. Sparse datasets are datasets in which for each shot only a part of the acquisition surface is covered with receivers. Furthermore, only a part of the acquisition surface may be covered with shots. Applying the CFP-method directly to data acquired by these sparse acquisition geometries will not give proper results. Due to the sparse sampling of the Fresnel-zone, the event in the CFP-gather will be severely aliased [Bolte, 2003]. Because of this aliasing the residual traveltime error can not be determined and, no updating can be done.

1.3 Research Objectives

This research has two main objectives:

• Develop a new operator updating procedure in terms of a parametric inversion procedure.

– The procedure should reduce the number of updating steps, if possible, to only a one-step procedure.

– The procedure should incorporate the proper physics to guarantee con-vergence to the true operators.

• Develop a method for 3D-CFP gather construction using sparsely sampled data, based on a data reconstruction technique.

– The method that reconstructs the required seismic data needed to perform the CFP-gather synthesis should be efficient.

– The method should be able to include a-priori knowledge to improve the results of the data-reconstruction.

1.4 Thesis Outline

In this thesis a method for determining focusing operators using both 2D and 3D seismic data is presented. The current chapter contains the introduction in which the concepts of seismic exploration, seismic imaging and the CFP-method are intro-duced. Also the problem statement and research objectives are given.

The next three chapters will present the theoretical physical framework of this re-search.

Chapter 2 gives an overview of the theory of how the CFP-method can be formu-lated in terms of Rayleigh integrals. Also the iterative operator estimation procedure will be discussed.

(14)

formulated as an inverse problem. First, a general forward model of the residual traveltimes will be formulated. Next, the parameterization of the focusing opera-tors in terms of B-splines will be discussed. Finally, we will show how to solve the obtained inverse problem using the appropriate optimization methods.

Chapter 4 addresses the second objective of this research. It discusses the as-sumptions generally made in seismic imaging methods. Next, it will present how these assumptions can be used to reconstruct full 3D seismic data required to create CFP-gathers from sparsely sampled data. Finally, the implementation of the data-reconstruction will be presented.

The next part of the thesis evaluates the different methods applied to synthetic and real data.

Chapter 5 will show the results of the new updating procedure on synthetic 2D and 2D field data. Also the results of a synthetic 3D dataset will be evaluated. Chapter 6 will present the results of the data-reconstruction for synthetic 3D-data. Finally, in Chapter 7, the advantages and limitations of the methods presented in this thesis will be evaluated. Also some recommendations will be given for further research. Case Studies Conclusions Introduction Framework

Theoretical Definition and generation of focusing operators

Evaluation of the Evaluation of the

3D CFP−gather construction Introduction to

the research and the thesis 1. 6. 5. 2. 7. Conclusions and Recommendations construction theory of the 4. Explanation and 3D CFP−gather procedure 3. Exposition of new operator updating new operator updating method

(15)

2

CFP technology revisited

In this chapter we will revisit the WRW-model which forms the basis of the CFP-method and we will introduce the concept of focusing operators. An explicit formu-lation of primary reflection data in terms of integrals as implied by the spatial con-volutions embedded in the WRW-model will be given, leading to the concept of a new dataset called the CFP-gather. By performing a stationary-phase analysis of the pre-sented integral notation for the CFP-gather, ’The principle of equal traveltime’ will be derived. Finally the iterative focusing operator estimation will be presented with its limitations.

2.1 Introduction

The primary goal of this research is to estimate the kinematic part of focusing operators, which are one-way traveltimes from points in the subsurface to a grid of points located at the acquisition surface, using the Common Focal Point (CFP)-method [Berkhout, 1997a,b; Thorbecke, 1997; Bolte, 2003]. In the remainder of this thesis the term focusing operators will be used for both the kinematic part of the focusing operators and for the focusing operators themselves.

Within the CFP-method, two-way primary reflection data can be decomposed into one-way seismic data, called a gather, and focusing operators. Both the CFP-gathers and the focusing operators have applications in seismic processing. E.g., the CFP-gathers are an excellent domain to estimate AVP-information [van Wijn-gaarden, 1998] and to perform time-lapse analysis [Winthaegen et al., 2001]. The obtained focusing operators can be used for redatuming [Kelamis et al., 2000; Tegt-meier et al., 2004] and for the estimation of a velocity model via a tomographic inversion of the one-way traveltimes [Cox, 2004].

(16)

[1997a,b] and Thorbecke [1997]. It was shown by Bolte [2003] that this iterative method works well for simple geologies but that for complex geologies the method may give ambiguous results.

First we will discuss in Section 2.2 how primary seismic reflection data conveniently can be represented in terms of spatial convolutions of forward one-way wavefield operators (W ) and reflection operators (R) using the W RW -concept introduced by Berkhout [1982].

The explicit formulation of primary seismic reflection data in terms of these forward operators allows for the elimination of the effects of forward wavefield propagation on seismic data by applying focusing operators on the data (F ), as discussed by Berkhout [1997a,b]. The focusing operators are approximate inverses of the forward one-way wavefield operators (W ). In the high-frequency limit they reduce to an amplitude factor and a set of phase shifts that correspond to traveltimes when de-scribed in the time domain, as will be discussed in Section 2.3.

Canceling one of the two forward one-way propagation legs present in primary seis-mic data by means of applying the correct focusing operator, results in a new dataset where either the sources or the receivers have been backward propagated into a point in the subsurface, called the focal point, which acts as a virtual buried source or re-ceiver. A record with receivers at the surface and a buried source, or the other way round, is called a CFP-gather (Section 2.4).

Finally, we will show how the correctness of the focusing operators used in the CFP-gather synthesis can be validated by using ’The principle of equal traveltime’ [Berkhout, 1997b]. Using this principle, an iterative procedure can be developed in which the focusing operators can be estimated out of the seismic data without knowledge of the subsurface velocity model (Section 2.5)

2.2 Decomposing seismic data in terms of forward propagation

operators

In this section we will discuss how primary reflection events in seismic data can be described in terms of forward one-way wavefield operators using the W RW -method. We will consider the response of the earth as measured by a grid of receivers located at the surface due to a wave field generated by a source also located at the surface. This response can be decomposed as cascaded operators applied to the initial source wavefield distribution, S. First, the downgoing wavefield generated by the source (Figure 2.1a) propagates forward to the interface. This is described by the operator W+_{(Figure 2.1b). Next, the downgoing wavefield at the interface is reflected and an}

upgoing wavefield is generated. This is described by the operator R (Figure 2.1c). The upgoing wavefield at the interface propagates up to the acquisition surface. This is described by the operator W− _{(Figure 2.1d). Finally the upgoing wavefield}

(17)

2.2 Decomposing seismic data in terms of forward propagation operators 9

The recorded wavefield represented in matrix notation is:

P = DW−RW+_S _(2.1)

This is, in essence, the WRW-model introduced by Berkhout [1997a,b]

a) P+_(x 0, xs, ω) b)                  W+_{(η, x} 0, ω) c) P−_(x r, xs, ω)                  d) W−(xr, ξ, ω)

Figure 2.1: Overview of WRW-model. a) A downward going wavefield is generated by a

source at xs. b) The downward going wavefield propagates from the source-position to the

interface. c) After reflection an upgoing wavefield at the interface is generated. d) The upgoing wavefield propagates from the interface to the receivers.

2.2.1 The source-characteristics

(18)

the amplitude of the wavefield operators without affecting the kinematics. Because all methods described in this thesis primarily depend on the kinematics of wavefield propagation, we can, without loss of generality, use a dipole-source in our model-description. If, however, one wants to extract amplitudes using the CFP-method, the source-character must properly be taken into account.

For a point xo at the surface the downgoing wavefield of a dipole-source at surface

location xscan be described as [Wapenaar., 1990]:

P+_(x

o, 0; xs, 0; ω) = δ(xo− xs)S(ω), (2.2)

where P+_(x

o, 0; xs, 0; ω) is the downgoing wavefield at surface position xo, δ(xo−xs)

is a spatial delta pulse and S(ω) is the source-signature.

2.2.2 Forward extrapolation from source to the interface

To compute the downgoing wavefield in an arbitrary point ra in the subsurface

caused by a downgoing wavefield distribution at the acquisition surface Ω, the Rayleigh integral for forward wavefield extrapolation will be used [Berkhout, 1982]. The Rayleigh integral states that the downgoing wavefield in an arbitrary point in the subsurface, P+_(r

a; xs, 0; ω), is a weighted sum of the downgoing wavefield at

the surface P+_(x

o, 0; xs, 0; ω) where the weights, W+(ra; xo, 0; ω), are determined

by the mass-density ρ(xo, 0) at the surface and the Green’s function G(xo, 0; ra, ω):

P+(ra; xs, 0; ω) = Z xo∈Ω 2 ρ(xo, 0) ∂G(xo, 0; ra; ω) ∂z | {z } W+_(r_a_;x_o_,0;ω) P+(xo, 0; xs, 0; ω)dxo, (2.3)

where xo∈ Ω denotes the spatial coordinates of points at the acquisition surface Ω.

In this thesis, points in 3D-space are denoted by r, and points on a surface having two independent coordinates are denoted by x. Furthermore, in the remainder of this thesis it is assumed that the acquisition surface is at depth 0 and the coordi-nate indicating depth is therefore omitted. Note that the operation performed as expressed in Equation 2.3 can be considered as a spatial convolution of the downgo-ing wavefield at the surface, P+(xo; xs; ω), with a space-varying convolution kernel

given by W+ _: P+(ra; xs; ω) = Z xo∈Ω W+(ra; xo; ω)P+(xo; xs; ω)dxo, (2.4) where W+_(r

a; xo; ω) is the wavefield operator.

In case the downgoing wavefield at the surface is given by Equation 2.2 and the point in the subsurface ra is taken on an interface Σ, the expression for the downgoing

(19)

2.2 Decomposing seismic data in terms of forward propagation operators 11

2.2.3 Reflection at the interface

On the interface the downgoing wavefield is reflected and an upgoing wavefield is generated. This resulting upgoing wavefield can be expressed in terms of the incident downgoing wavefield and a reflection operator describing the reflectivity properties of the interface. For a point ξ on the interface Σ the upgoing wavefield can be expressed as follows:

P−(ξ; xs; ω) =

Z η∈Σ

R(ξ; η)P+(η; xs; ω)dηt, (2.6)

where P−(ξ; xs; ω) is the upgoing wavefield at point ξ at the interface, R(ξ; η) is

the reflection operator, and dηt(= dηx0, dηy0) denotes an element of the interface,

which is assumed to be locally plane, but not necessarily horizontal.

The space variant convolution operator R(ξ; η) actually is the inverse spatial Fourier transform of the operator ˜R(ξ, k˜ ηx0, kηy0), which gives the angle and azimuth depen-dent reflection coefficient of the interface Σ, at location ξ. Under the assumption that ˜R(ξ, k˜ ηx0, kηy0) is azimuth independent we can write:

R(ξ; ηx0; ηy0) = r(ξ; α) (2.7)

with α = arcsine(c

q k2

ηx0+k2ηy0

ω ), where c is the local acoustic propagation velocity.

In the high-frequency limit we can have a stationary phase evaluation of the integral in Equation 2.6. This is equivalent to assuming that locally P−(ξ; xs; ω) can be

described by a single plane wave with angle of incidence α that depends on ξ and xs. In this approximation and combining Equations 2.5 - 2.7 we can express the

upgoing wavefield at the interface as follows:

P−(ξ; xs; ω) = r[ξ; α(ξ; xs)]W+(ξ; xs; ω)S(ω). (2.8)

In a more general description, the assumption of a locally interacting plane wave-field does not hold and things become much more complicated. However, evaluating Equation 2.6 in the high-frequency limit, will only effect the amplitude of the result-ing wavefield expression without havresult-ing effect on the kinematic part and is therefore justified within the scope of this thesis. For a thorough discussion on the reflection operator, see de Bruin et al. [1990].

2.2.4 Forward extrapolation from the interface to the receivers

In the final step, the reflected wavefield at the interface propagates upward to the acquisition surface. Under the assumption of a plane interface, the Rayleigh integral for upgoing wavefields may be used to describe this operation. The upgoing wavefield at point xr at the surface due to the reflected wavefield at the interface Σ is given

(20)

P−(xr; xs; ω) = Z ξ∈Σ 2 ρ(ξ) ∂G(ξ; xr; ω) ∂ξz | {z } W−(xr;ξ;ω) P−(ξ; xs; ω)dξ, (2.9)

where ξ _{∈ Σ denotes the spatial coordinates of points at interface Σ. Combining} Equations 2.8 - 2.9 the upgoing wavefield as measured on receiver position xr due

to a dipole source at xs and one interface can be expressed as follows:

P−(xr; xs; ω) =

Z ξ∈Σ

W−(xr; ξ; ω) r[ξ; α(ξ; xs)]W+(ξ; xs; ω) S (ω) dξ. (2.10)

As can be seen in Equation 2.10, the event in the seismic data corresponding to primary reflection at a single interface can be thought of as a combination of one-way wavefield operators. To represent the primary reflections of all interfaces present in the subsurface the separate responses can be superimposed:

P−(xr; xs; ω) = X i Z ξ∈Σi W−(xr; ξ; ω) ri[ξ; αi(ξ; xs)]W+(ξ; xs; ω) S (ω) dξ, (2.11) where the i-th interface has a reflection coefficient ri[ξ; αi(ξ; xs)].

2.3 Focusing operators

In the previous section the concept of forward one-way wavefield operators was in-troduced to describe the effect of forward wavefield propagation in a medium on the initial wavefield as generated by the source. In a similar fashion the effect of in-verse wavefield extrapolation can be described in terms of inin-verse one-way wavefield operators, where inverse wavefield propagation can be understood as undoing the effect introduced on the seismic data by forward wavefield propagation. Within the CFP-method these inverse one-way wavefield operators are called focusing operators because when they are applied to the seismic data, they focus either the sources or the receivers into a point in the subsurface. The focusing operators are closely related to the forward wavefield operators as will be discussed in this section. First, we will introduce the formal definition of these focusing operators. Next, an ex-plicit formulation of the focusing operators will be presented by using the modified matched filter approach. Finally, we will show how this concept allows us to express the focusing operators in terms of amplitudes and a set of one way traveltimes, which are of interest for us.

2.3.1 Formal definition of focusing operators

(21)

2.3 Focusing operators 13

expressed as seeking those operators that, after they have been applied to the forward one-way wave propagation operators, produce a spatial delta-pulse centered at the focal point. The focusing operator F+ _{that corrects for the downward propagation}

leg W+ _{has to obey the following expression:}

Z

x∈Ω

W+_{ξ; x; ω}_˜ _F+_{(x; ξ; ω)dx = δ( ˜}_ξ

− ξ), (2.12)

where ˜ξ is a point in the subsurface and ξ is the focal point. Similarly, the relation between the focusing operator F− that corrects for the upward propagation leg W− is given by:

Z

x∈Ω

F−(ξ; x; ω) W−x; ˜ξ; ωdx = δ( ˜ξ− ξ). (2.13)

2.3.2 Modified matched filter approximation of focusing operators

Whereas Equations 2.12 - 2.13 give some insight in the formal definition of the different focusing procedures, practical use of the CFP-concept requires an explicit expression of the focusing operators instead of the implicit formulation as given by Equations 2.12 - 2.13.

In Berkhout [1982] different explicit formulations of the focusing operators are re-viewed of which the modified matched-filter approach is the most desirable one, leading to the following expressions:

F−(ξ, x, ω)_{≈ [W}+_{(ξ, x, ω)]}∗_, _(2.14)

and

F+(x, ξ, ω)_{≈ [W}−(x, ξ, ω)]∗, (2.15)

where the ∗ _{denotes the complex conjugate and results in time reversal when}

con-sidered in the time domain. The modified matched filter approach is correct for propagating waves only, for evanescent waves it is incorrect. But in general the effect of evanescent waves can be neglected. As can be seen, F−_{(ξ; x; ω) can be}

considered the time reverse response of a source x at the acquisition surface and a receiver at reflector location ξ in the subsurface.

Similarly, F+_{(x; ξ; ω) can be considered the time reverse response of a source at the}

reflector and a receiver at the acquisition surface.

Due to reciprocity the locations of the sources and receivers may be interchanged and, therefore, the focusing operator for focusing in detection that brings down the receivers, can be modeled by forward modeling a source at the focal point and receivers at the acquisition grid and taking the time reverse of this response.

2.3.3 High frequency approach

(22)

obvious, when both the forward wave propagation and the focusing operators are expressed in terms of their high-frequency approximations.

In the high-frequency approach the forward extrapolation operators are expressed by an amplitude factor A and a phase factor [Bleistein, 1984]:

W+_{(ξ; x}

s; ω) = A (ξ; xs) e−iωτ1(ξ;xs), (2.16)

where τ1(ξ; xs) is the one-way traveltime from source xs to reflector point ξ.

Combining the high-frequency approach with the modified matched filter expression given by Equation 2.15 we obtain:

F+(xs; ξ; ω) = A (ξ; xs) eiωτ1(ξ;xs). (2.17)

A similar expression can be obtained for the focusing operator for focusing in detec-tion:

F−(ξ; xr; ω) = A (xr; ξ) eiωτ1(xr;ξ), (2.18)

where τ1(xr; ξ) is the one-way traveltime from reflector position ξ to receiver xr.

2.4 CFP-gathers

Here we will discuss how the focusing operators, introduced in the previous section, can be used to eliminate the forward wave propagation legs present in the seismic data. First, the formal definition of the focusing operators is used to illustrate the focusing property of the focusing operators, and it will be discussed briefly how focusing operators can be used to migrate seismic data. Next, ’The principle of equal traveltime’ will be derived that can be used to validate the correctness of the kinematic part of the focusing operators. Finally, we will show a practical implementation to create CFP-gathers.

2.4.1 Focusing seismic data

As described earlier, primary reflection events can be decomposed into two one-way wave propagation operators (Figure 2.2a-b). The effect of these forward propagation operators can be undone by applying focusing operators to the seismic data thereby inverse back-propagating the sources and receivers to points in the subsurface. De-pending on which one-way wave propagation leg we want to eliminate, two different types of focusing operations can be distinguished.

Focusing the detectors

When we focus the detectors, the upward propagation leg from the reflector to the receivers is removed by applying a focusing operator to all receivers of a common shot-gather (Figure 2.2c-d), which can be formulated as follows:

P−(ξc; xs; ω) =

Z

Ω

(23)

2.4 CFP-gathers 15 a) c) e) g) 0 0.2 0.4 0.6 0.8 1.0 time 500 offset1000 1500 2000 i) 0 0.5 1.0 1.5 2.0 time -500 0 offset 500 1000 b) 0 250 500 750 1000 1250 1500 1750 2000 0 0.2 0.4 0.6 0.8 1 Surface Location (m) Time (s) d) 0 0.2 0.4 0.6 0.8 1.0 time 750 f) 0 0.2 0.4 0.6 0.8 1.0 time 500 offset1000 1500 2000 h) -0.2 -0.1 0 0.1 0.2 time -1000 -500 offset0 500 1000 j)

Figure 2.2: Overview of different steps in the CFP-procedure. a) The raypaths corresponding

(24)

where ξcis the focal point, P−(ξc; xs; ω) is the result after focusing the common-shot

gather corresponding to a source at xs, and F−(ξc; xr; ω) is the focusing operator

corresponding to ξc. As can be seen in Equation 2.19, the result after focusing is

no longer a function of the receiver coordinate because we stack all receivers of the common-shot gather.

When we consider the target reflection only we can plug Equation 2.10 into Equation 2.19, giving : P−(ξc; xs; ω)= Z Ω dxrF−(ξc; xr; ω) Z Σ dξW−(xr; ξ; ω) r[ξ; α(ξ, xs)]W+(ξ; xs; ω)S(ω) . (2.20) Interchanging the order of integration and using Equation 2.12, both integrals are eliminated and we obtain a simple expression:

P−(ξc; xs; ω) = W+(ξc; xs; ω)r[ξc; α(ξc; xs)]S (ω) . (2.21)

The result after focusing the detectors of one common shot gather, as expressed by Equation 2.21, can be interpreted as the response of a virtual receiver located at the focal point ξc in the subsurface caused by a source at xs still located at

the acquisition surface (Figure 2.2e-f). Focusing the detectors of all common shot gathers gives thus a new dataset with a virtual receiver located at the focal point and sources at the acquisition surface (Figure 2.2g-h). Such a dataset is called a CFP-gather.

Focusing the sources

Instead of focusing the receivers in a common shot gather, we can also focus the sources in a common receiver gather. Mathematically we describe focusing the sources as:

P−(xr; ξc; ω) =

Z

Ω

P−(xr; xs; ω)F+(xs; ξc; ω)dxs. (2.22)

Using Equation 2.13 this can be written as:

P−(xr; ξc; ω) = r[ξc; α(ξc; xr)]W−(xr; ξc; ω)S (ω) . (2.23)

This can be interpreted as the response a receiver at xr to a virtual source located

at the focal point ξc.

Imaging

(25)

2.4 CFP-gathers 17

primary reflection energy in the subsurface. In terms of migration this is a full stack image because the receivers and sources are all focused at the same point. Note that each CFP-gather after imaging, yields only one point in the migrated image, i.e. the full stack reflection property at the focal point. In this thesis we will primarily deal with the data in the half-migrated domain

2.4.2 The principle of equal traveltime

Applying the high-frequency approach to the response obtained after the focusing step of the detectors, as given by Equation 2.21, and combining all amplitude factors in one factor Atot we obtain:

P−(ξc; xs; ω) = Atot(ξc; xs) e−iωτ1(xs;ξc), (2.24)

where τ1(xs; ξc) are the one-way traveltimes from the focal point to the source

loca-tions. Comparing the phase terms in the result after focusing as given in Equation 2.24 and the focusing operator used, Equation 2.18, we see that they both contain the one-way traveltimes from the focal point to points at the acquisition grid, they only differ in sign.

This observation leads to ’The principle of equal traveltime’, as defined by Berkhout [1997b], that states that when the correct focusing operator has been applied to the data, the traveltimes of the target event in the CFP-gather should coincide with the time reversed focusing operator, evaluated at the source positions (Figure 2.2i). This can be verified by shifting the traces in the CFP-gather over the operator times, yielding the Differential Time Shift (DTS)-panel, in which then the target event should be a flat event, aligned at t = 0 (Figure 2.2j).

2.4.3 Discrete implementation of focusing seismic data

In the above the focusing properties of the focusing operators were illustrated by using the analytic expression given by Equation 2.19. In practice, however, seismic data is acquired on a spatial grid, and integration over the receivers is carried out as a summation of discrete numbers:

P−(ξc; xs; ω) =

X

xr

F−(ξc; xr; ω)P−(xr; xs; ω). (2.25)

Inserting the high-frequency approximation of the focusing operator, given by Equa-tion 2.18 into EquaEqua-tion 2.25 we obtain:

P−(ξc; xs; ω) =

X

xr

A (xr; ξc) eiωτ1(xr;ξc)P−(xr; xs; ω), (2.26)

in which τ1(xr; ξc) is the one-way traveltime from the focal point to the receiver xr

at the acquisition grid and A (xr; ξc) is the corresponding amplitude factor.

(26)

and traveltimes, the actual computation of one trace in a CFP-gather is a two-step procedure:

• Correct all receiver recordings of a common-shot gather with the operator times.

• Stack all corrected receiver recordings weighted by the amplitude terms, yield-ing one trace.

In case of a coarse spatial sampling, however, the discrete implementation given by Equation 2.26 can be spatially aliased and does, therefore, not resemble the true value of the integral in Equation 2.13. How this will affect the CFP-method in case of sparse 3D data will be discussed in Chapter 4. In the remainder of this chapter it will be assumed that the seismic data is properly spatially sampled and that, therefore, Equation 2.26 may be used.

0 250 500 750 1000 1250 1500 1750 2000 0 0.2 0.4 0.6 0.8 1 Surface Location (m) Time (s) a) 0 0.2 0.4 0.6 0.8 1.0 time 500 offset1000 1500 2000 b) -0.2 -0.1 0 0.1 0.2 time -1000 -500 offset0 500 1000 c) 0 0.2 0.4 0.6 0.8 1.0 time 500 offset1000 1500 2000 d)

Figure 2.3: Construction of a DTS-panel in case an erroneous operator is used. a) The

(27)

2.5 Data-driven focusing operator updating 19

2.5 Data-driven focusing operator updating

’The principle of equal traveltime’ was introduced as a tool to validate the cor-rectness of the applied focusing operators. In case the correct velocity model is known, the focusing operators can be computed by forward modeling in this veloc-ity model. In general, however, only an estimate of the velocveloc-ity model is known and incorrect estimates of the focusing operator will be used in the construction of a CFP-gather. Consequently, the corresponding event in the DTS-panel will in general not be aligned and flat at t = 0, but will exhibit residual time shifts (Figure 2.3a-c). Here we will show how these residual time shifts can be used in an iterative focus-ing operator estimation procedure. An iterative method to use these residual time shifts observed in the DTS-panels, to obtain a better estimate of the focusing op-erators, without updating the velocity model, was presented by Berkhout [1997a,b], Thorbecke [1997] and Bolte [2003]. Their method is based on the following. Consider an incorrect focusing operator used for focusing in detection:

ˆ

F−(xr; ξc; ω) = A (ξc; xr) eiω ˆτ1(ξc;xr), (2.27)

where ˆτ1(xr; ξc) denotes the incorrect traveltimes from the focal point to the

re-ceivers at the acquisition grid. Applying this focusing operator on a common shot gather, focusing the detectors is performed, aiming to backward propagate down the receivers to the focal point. The resulting CFP-gather is generally described by:

P−(ξc; xs; ω) = ˆAcomb(ξc; xs) e−iω˜τ1(ξc;xs), (2.28)

where ˜τ1(ξc, xs) is the time at which the target event will appear in the CFP-gather,

and ˆAcomb(ξc; xs) is the corresponding amplitude factor. Shifting the resulting

traces in the CFP-gather over the traveltimes of the focusing operator evaluated at the source location, ˆτ1(ξc; xs), gives the DTS-panel. Because the focusing operator

was incorrect, the event will not appear at t = 0 but at residual times given by: τDT S(ξc; xs) = ˜τ1(ξc; xs)− ˆτ1(ξc; xs). (2.29)

Next, the assumption is made that the correct operator is in between the incorrect operator and the target event in the CFP-gather that was created with the incorrect operator. This assumption can be motivated by looking at Figure 2.3d, where we see that this is actually the case for this simple model. This assumption leads to the following pragmatic approach to get a better estimate of the focusing operator:

ˆ

τ₁upd(ξc; xs) = ˆτ1(ξc; xs) +

τDT S(ξc; xs)

2 , (2.30)

where ˆτ1upd(ξc; xs) is supposed to be a better approximation of the true operator,

than the original incorrect estimate ˆτ1(ξc; xs). Using the updated operator a new

(28)

• Use current estimate of the focusing operator to create a CFP-gather. • Timeshift the CFP-gather using the focusing operator, yielding the DTS-panel. • Pick the residual times observed in the DTS-panel.

• Update the operator using the residual times if necessary and repeat this pro-cedure.

Each time the event in the DTS-panel is picked we determine whether the picked residual time shifts are below a user-specified threshold. If not, the operator is up-dated and the procedure is repeated until the residual times are small enough. Note that in this approach the focusing operators corresponding to neighboring re-flector points are updated independently, as the update of a certain focusing operator is only determined by the residual time shifts observed in the DTS-panel correspond-ing to this focuscorrespond-ing operator.

Because the updating is directly performed on the focusing operators, no velocity model is needed. In addition, ’The principle of Equal Traveltime’, used to verify the correctness of the focusing operators, is very straightforward.

2.5.1 Issues of the iterative procedure

This iterative operator estimation has shown good results on some relatively simple datasets [Bolte et al., 1999a,b]. There are, however, some issues with this approach.

Non-uniqueness of solution

As already recognized by Bolte [2003], ’The principle of Equal Traveltime’ is a nec-essary but not a sufficient condition to guarantee the correctness of the focusing operators. In other words: there may be multiple sets of focusing operators that all give flat DTS-panels but of which only one is correct. Therefore, validating ’The principle of equal traveltime’ alone does not guarantee the correctness of the focusing operator.

No convergence is guaranteed

Another important issue with the iterative approach is that the pragmatic updating procedure given by Equation 2.30 does not guarantee that the focusing operator obtained after updating is closer to the true operator than before updating. So convergence to the (non-unique) solution is not guaranteed.

Picking is laborious

(29)

2.5 Data-driven focusing operator updating 21

(Initial)

Focusing

Operator

Seismic Data

CFP−gather

Create

DTS−panel

Create

DTS−panel

Pick

PET met?

no

yes

Final

Focusing

Operator

Update

Focusing Operator

Figure 2.4: Flowchart of the iterative focusing operator updating. ’The principle of equal

traveltime’ is abbreviated to PET.

2.5.2 Solving the issues

In his PhD-thesis Bolte [2003] proposes some additional constraints to ’The prin-ciple of Equal Traveltime’ and introduces alternative updating methods that use knowledge of the subsurface and thereby reduce the null space of the solutions. Unfortunately, these alternative methods only work for some simple geologies. Fur-thermore, the methods he proposes may be even more laborious than the iterative updating method.

(30)

(31)

3

One-step focusing operator updating

In this chapter we present a new focusing operator updating procedure via a parametric inversion. The core of the parametric inversion algorithm is a forward model of the kinematics of the target event in the DTS-panels which expresses the mismatch between the focused seismic data and the focusing operators. In this forward model, the focusing operators are assumed to exhibit a certain degree of smoothness, allowing an efficient parameterization in terms of splines. Through an optimization procedure, the spline parameters that enforce the forward DTS model to match the events in the DTS-panels are retrieved in one step.

3.1 Introduction

(32)

model can be re-formulated in such a way that it does not contain any velocity model information, using a procedure called ZO-ray Mapping as will be the topic of Section 3.3. An inversion in which we directly invert the residual DTS-times for the correct operator times, using the forward model, is not advisable. First of all the number of unknowns equals the number of data-points, resulting in a potentially ill-posed inverse problem. Secondly, the forward model contains a set of constraints that depend on spatial derivatives of the focusing operators that need to be evaluated, which is not straightforward in the above mentioned forward model. How these issues can be resolved by parameterizing the focusing operators in terms of multi-variate B-splines will be discussed in Section 3.4. Finally, in Section 3.5 it will be shown how the resulting non-linear inverse problem, formulated in terms of multi-variate B-splines, can be solved, using a Non-linear Conjugated Gradient (NLCG)-scheme that allows a one-step focusing estimation procedure in which the CFP-gathers and DTS-panels only need to be created once.

3.2 Forward model of the CFP-gather synthesis

In this section the forward model of the residual times in the DTS-panel will be formulated. It consists of two key concepts. Firstly, the two-way traveltimes (twtt’s) of primary reflection events present in seismic data will be decomposed as a sum of one-way traveltimes (owtt’s) subject to constraints based on Fermat’s principle. Secondly, the residual times in the DTS-panel can be decomposed into the twtt’s of the primary event in the data and the incorrect operator times, subject to an additional constraint that is also based on Fermat’s principle. All combined, these concepts form a non-linear forward model of the observed DTS-times in terms of the known incorrect initial operators times and the correct, but unknown, operator times.

3.2.1 Decomposing twwt’s as sum of owtt’s

The point of departure in the derivation of the forward model is the forward wavefield representation of primary reflection data presented in Chapter 2 that will be repeated here for convenience:

P−(xr; xs; ω) =

Z ξ∈Σ

W−(xr; ξ; ω) r[ξ; α(ξ; xs)]W+(ξ; xs; ω) S (ω) dξt, (3.1)

where W−_(x

r; ξ; ω) and W+(ξ; xs; ω) represent forward one-way wavefield

(33)

3.2 Forward model of the CFP-gather synthesis 25

wavefield representation of primary reflection data given in Equation 3.1: P−(xr; xs; ω)=

Z ξ∈Σ

A (xr; ξ) e−iωτ1(xr;ξ)r[ξ; α(ξ; xs)]A (ξ; xs) e−iωτ1(ξ;xs)S (ω) dξt,

(3.2) in which τ1(xr; ξ) represents the owtt from reflector point ξ to the receiver location

xr and τ1(ξ; xs) represents the owtt from reflector point ξ to the source location

xs. Combining all amplitude factors, the wavelet and the reflection coefficient into

a single factor Atot, we can write:

P−(xr; xs; ω) =

Z ξ∈Σ

Atot(xr; xs; ξ; ω) e−iωτ2(xr;ξ;xs)dξt, (3.3)

where τ2(xr; ξ; xs) is the sum of the one-way traveltimes:

τ2(xr; ξ; xs) = τ1(xr; ξ) + τ1(ξ; xs). (3.4)

The integral expression in Equation 3.3 consists of a slowly varying amplitude factor, being Atot(xr; xs; ξ; ω) and a rapidly varying phase factor and can, therefore, be

approximated using the stationary phase approach, [Bleistein, 1984], which states that the main contribution of the integral is determined by those points for which the phase of the integrand is stationary:

∂τ2(xr; ξ; xs)

∂ξ

_ξ₌_ξst = 0, (3.5)

where ξstis a stationary reflector point. Note that Equation 3.5 actually represents

Fermat’s principle that states that the raypath connecting two points has a station-ary traveltime when comparing all possible raypaths. In case, for a given xr and

xs, multiple reflector points are stationary, all these points contribute, resulting in

multi-valuedness of the two-way travel times. The contribution of one stationary point is given by:

P−(xr; xs; ω) = 2π ω m 2 _A tot(xr; xs; ξst) p |detH| e −iωτ2(xs;ξst;xr)_eiπ4sgn(detH), (3.6) where m = 2 in case of 3D and m = 1 in case of 2D, H stands for the Hessian matrix of τ2 containing the second-order derivatives of τ2 with respect to the horizontal

components of ξ, and ξst is the stationary phase point. As follows from Equation

(34)

ξ

_st

xr,st

xs

Figure 3.1: Of all possible raypaths between source and receiver, the raypaths for which the

twtt is stationary (the thick line) will contribute most to the reflection event.

3.2.2 Stationary phase analysis of the Fresnel stack

After having shown that the twtt’s of primary reflection data can be decomposed into the correct owtt’s of the focusing operators, we can now analyze the impact that application of an incorrect focusing operator has on the resulting CFP-gather. For a focusing in detection step, the focusing operator is applied to the receivers of a common-shot gather, such that they are downward continuated to the focal point. Considering the case that an incorrect estimate of the focusing operator is used this procedure can be written as:

ˆ P−(ξc; xs; ω)= Z Ω ˆ F−(ξc; xr; ω) Z Σ W−(xr; ξ; ω) r[ξ; α(ξ; xs)]W+(ξ; xs; ω) S (ω) dξtdxr, (3.7) where ˆF−(ξc; xr; ω) denotes the incorrect focusing operator, and ˆP−(ξc; xs; ω) is

the result after focusing the receivers. Because an incorrect focusing operator has been used, the focusing operator will not cancel correctly the upward propagation leg from the reflector to the acquisition surface. As a consequence the receivers are not correctly focused at the focal point. Nevertheless, is is still possible to obtain a physical interpretation of focusing seismic data with an incorrect focusing operator. To see this we look at the mathematical expression for the focusing operator:

ˆ

F−(ξc; xr; ω) = ˆA(xr; ξc)eiω ˆτ1(xr;ξc), (3.8)

(35)

3.2 Forward model of the CFP-gather synthesis 27

owtt’s of the focusing operator. In the remaining we assume that the amplitude terms are correct.

Using Equation 3.6 and inserting the high-frequency approximation of the focusing operator given by Equation 3.8 into Equation 3.7 and by making use of the twtt decomposition of primary reflection data we can write the result after focusing in detection as: ˆ P−(ξc; xs; ω) = Z Ω ˆ

Atot(xr; xs; ξst; ξc; ω)e−iω(τ2(xs;ξst;xr)−ˆτ1(xr;ξc))dxr, (3.9)

where all amplitude components are combined in ˆAtot. Note that for each receiver

location xr of the common-shot gather in Equation 3.9 over which the integration

is carried out, the stationary reflector point ξst can differ and has to obey Fermat’s

principle as formulated in Equation 3.5.

The result after focusing the receivers can also be analyzed using a stationary phase approach which gives the following result:

ˆ P−(ξc; xs, ω) = 2π ω m 2 _A tot,2 p |detH|e −iω(τ2(xs;ξ_st;xr,st)−ˆτ1(xr,st;ξc))_eiπ₄sgn(detH)_, (3.10) where the H is the Hessian of τ2− ˆτ1 w.r.t. to xr, and the stationary receiver xr,st

has to obey the following constraint:

d(τ2(xs; ξst; xr)− ˆτ1(xr; ξc)) dxr _x r=xr,st = 0. (3.11)

Equations 3.10 - 3.11 tell us that the major contribution in the construction of a trace in the CFP-gather comes from those receivers in the shot-gather for which the event in the data is tangent to the operator, and appears at the twtt in the data corrected with the incorrect operator time at those stationary receivers (Figure 3.2). To obtain the spatial derivative of the twtt in the data w.r.t. to the receiver coordi-nate as expressed in Equation 3.11 we take a closer look to the total differential of the twtt which is defined as follows:

dτ2 = ∂τ2 ∂xr dxr+ ∂τ2 ∂ξ dξ + ∂τ2 ∂xs dxs, (3.12)

in which the second term on the r.h.s. vanishes because the twtt is stationary w.r.t. the reflector coordinate as expressed by Equation 3.5. Furthermore, the third term does not contribute because in the construction of a trace in the CFP-gather for focusing in detection, xs is constant and therefore dxs= 0. The only term that

contributes to the total differential is the partial derivative of the twtt w.r.t. to the receiver coordinate, which can be written using Equation 3.4 as:

(36)

0 0.5 1.0 1.5 2.0 time -500 0 offset 500 1000 xr,st

Figure 3.2: Illustration of stationary phase analysis of Fresnel stack. The main contribution

of the Fresnel stack comes from those receivers, xr,st, for which the operator (dashed white

line) and the target event in the shot-gather (black dashed line) are tangent (solid lines).

Because τ1(ξst; xs) does not directly depend on xr, the second term on the rhs in

Equation 3.13 also vanishes. Inserting Equation 3.13 with the above mentioned consideration into the Fresnel constraint as given by Equation 3.11, the constraint that determines the stationary receiver reduces to :

∂ ˆτ1(xr; ξc) ∂xr _x r=xr,st− ∂τ1(xr; ξst) ∂xr _x r=xr,st = 0. (3.14)

This expression tells us that the main contribution of a trace in a CFP-gather is determined by the receiver for which the initial incorrect operator corresponding to ξc is tangent to the correct but unknown operator corresponding to ξst. The time

at which the event appears in the CFP-gather is given by evaluating Equation 3.10 at the stationary receiver and is equal to :

ˆ

τCF P(ξc, xs) = τ1(xr,st; ξst) + τ1(ξst; xs)− ˆτ1(xr,st; ξc). (3.15)

As can been in this expression, the time at which the target event appears in the CFP-gather is a mixture of the incorrect initial focusing operator corresponding to ξc and the correct but unknown operator corresponding to ξst.

3.2.3 Forward model of the DTS-times

(37)

3.3 Zero-offset ray mapping 29

appears in the CFP-gather, given by Equation 3.15, is corrected with the operator time from the focal point to the source-location:

ˆ

τDT S(ξc; xs) = τ1(xr,st; ξst) + τ1(ξst; xs)− ˆτ1(xr,st; ξc)− ˆτ1(ξc; xs), (3.16)

where the stationary reflector point xr,st, and the stationary receiver point ξst are

implicitly given by the constraints: ∂τ1(xr,st; ξ) ∂ξ _ξ₌_ξst+∂τ1(xs; ξ) ∂ξ _ξ₌_ξst = 0, (3.17) and ∂τ1(xr; ξst) ∂xr _x r=xr,st− ∂ ˆτ1(xr; ξc) ∂xr _x r=xr,st = 0. (3.18)

Together Equations 3.16 - 3.18 constitute a forward model of the DTS-times ex-pressed in terms of the incorrect initial focusing operators and the correct but un-known focusing operators subject to two physical constraints. The expressions in this forward model depend on the reflector coordinates ξst and ξcwhich is

undesir-able for two reasons. Firstly, the CFP-location corresponding to the initial focusing operator ξc, is not known and can, therefore, not be used in an evaluation of the

forward model. Moreover, by mapping the reflector coordinates on the acquisition surface, as will be explained in the next section, it is possible to express the focusing operators as function of offset which, combined with the use of B-splines as basis functions, allows an efficient parameterization of the focusing operators as will be discussed in Section 3.4.

3.3 Zero-offset ray mapping

In the forward model presented above it was stated that the focal point locations that correspond to the incorrect operators, should be known, which can be considered as a weak form of dependency on the velocity model. In most situations, however, no knowledge on the locations of these focal points is available. Nevertheless, it is possible to formulate the forward model in such a way that the dependency on the velocity model vanishes and the forward model can be inverted without explicit velocity model information.

3.3.1 Forward model after zero-offset ray mapping

The key to the reformulation of the forward model is the assumption that for each point at the acquisition surface xo, there is a unique ZO-ray that connects the

acquisition surface at xo with the target reflector at a unique reflector location ξ

(Figure 3.3a). As a consequence, there is a monotonic relation between xo and ξ

that maps each reflector point to a point at the acquisition surface along the ZO-ray, see Figure 3.3b, hence the name ZO-ray mapping. Due to the monotonic relation the following relation is valid:

∂xo

(38)

0 x ξ

_ξ

0

x

a) b)

Figure 3.3: a) ZO-rays in a configuration that obeys the monotonic relation between reflector

locations and the ZO-ray intersection points at the acquisition surface. b) The (ξ, xo)-graph

should be monotonic so only one reflector point ξ will map into a point at the acquisition surface xo.

Because of the monotonic relation, we can uniquely express the focusing operators as function of the surface location xo instead of the reflector coordinate ξ. Doing

this for the forward model, given by Equations 3.16 - 3.18, leads to: ˆ

τDT S(xo,c; xs) = τ1(xr,st; xo,st) + τ1(xo,st; xs)− ˆτ1(xr,st; xo,c)− ˆτ1(xo,c; xs), (3.20)

subject to the constraints: ∂τ1(xr,st; xo) ∂xo ∂xo ∂ξ _x o=xo,st +∂τ1(xs; xo) ∂xo ∂xo ∂ξ _x o=xo,st = 0, (3.21) and ∂τ1(xr; xo,st) ∂xr xr=xr,st− ∂ ˆτ1(xr; xo,c) ∂xr xr=xr,st = 0. (3.22)

The first constraint as given by Equation 3.21 can be simplified by dividing by ∂xo

∂ξ,

which is justified because of Equation 3.19, which gives: ∂τ1(xr,st; xo) ∂xo _x o=xo,st +∂τ1(xs; xo) ∂xo _x o=xo,st = 0. (3.23)

In order for xo,st to correspond to a ZO-ray at xo,st, the owtt along the ZO-ray

needs to be stationary w.r.t. the reflector coordinate. This can be expressed as an additional constraint that reads:

(39)

3.3 Zero-offset ray mapping 31 0 ξ x 0

x

ξ

a) b) 0 x ξ I II III

x

0

t

I

III

II

c) d)

Figure 3.4: a) ZO-rays in a configuration that does not obey the monotonic relation between

reflector locations and the ZO-ray intersection points at the surface position. b) The (ξ; xo

)-graph is non-monotonic so multiple reflector points ξ can map into one point at the surface

xo. c) Three monotonic parts can be distinguished in the non-monotone graph. d) The

non-monotonic behavior results in multi-valuedness in the ZO-section.

The new forward model defined by Equations 3.20, 3.22, 3.23, and 3.24 no longer depends on knowledge of velocity model parameters. The only requirement is a monotonic mapping of the reflector coordinate to a point at the acquisition surface without the need of knowing the mapping function itself. The monotonic mapping assumption, however, is generally not valid.

A simple example of a case in which the assumption does not hold is given in Figure 3.4a, in which the ZO-rays for a syncline model are shown. As can be seen in Figure 3.4b, multiple reflector points can map into one acquisition point and, consequently, the corresponding mapping functions can have locations where ∂xo

∂ξ = 0, so that

dividing by ∂xo

∂ξ is not allowed. However, if we are able to distinguish between the

(40)

branches show up as multi-valuedness in the ZO-section (Figure 3.4d).

In the remainder of this thesis the term reflector point will also be used to indicate the corresponding ’zero offset point’ at the acquisition surface.

3.3.2 Direct inversion of velocity model independent forward model

In the above, a forward model of the residual times observed in the DTS-panel is presented in terms of the correct operator times evaluated on a spatial grid deter-mined by the stationary receivers and mapped stationary reflector points that are implicitly determined by the constraints.

Based on this forward model and a number of residual DTS-times, obtained by picking the target event in the DTS-panels, an inversion procedure can be designed to find the traveltimes of the correct focusing operators that explain the picked residual DTS-times. As can be seen in Equation 3.20 there are two unknown correct traveltimes per DTS-time, which suggests a system with more variables than data. However, if we realize that due to reciprocity the unknown traveltimes will also show up at other locations in the DTS-panels, the number of unknowns will be equal to the number of data points in the case that all DTS-panels corresponding to the same reflector are considered simultaneously, leading to a system that can be solved, albeit with a big risk of ill-posedness.

Furthermore, the forward model contains first-order spatial derivatives of the focus-ing operators and in case we want to use a gradient based solver, also the spatial derivatives of the constraints will be needed, giving expressions containing second-order spatial derivatives of the focusing operators. These spatial derivatives can be obtained by employing finite-difference schemes. These finite-difference schemes, however, are prone to errors and expensive to use.

In order to overcome these problems we can introduce an alternative parameteriza-tion of the traveltime funcparameteriza-tions, which will be the topic of the next secparameteriza-tion.

3.4 Spline parameterization of operators

In this section it will be shown how the traveltimes that make up the set of focusing operators that correspond to the same reflector, can efficiently be parametrized by multivariate B-splines. The use of B-splines not only reduces the number of model parameters that need to be determined in the inversion procedure, it also allows an analytic computation of the constraints and gradient information needed in the inversion procedure. All the illustrations presented in this section will be 2D and the formulation of the forward model will be 2D. The formulation of the forward model in 3D is given in Appendix C.

3.4.1 Smoothness of operator times

(41)

3.4 Spline parameterization of operators 33 0 200 400 600 800 1000 1200 1000 2000 3000 4000 5000 I II III 0 900 1800 2700 3600 4500 5400 0 900 1800 2700 3600 4500 5400 2 1.5 1 0.5 0 Reflector Location (m) Surface Location (m) Traveltime (s) III I II a) b)

Figure 3.5: Illustration of the smoothness of the focusing operator times. a) The velocity

model with the target reflector indicated by the white dashed line. Note the fault at the centre of the velocity model. b) The focusing operators corresponding to the target reflector.

I corresponds to xCF P = 600m, II to xCF P = 2700m, and III to xCF P = 4800m .

The first one is that neighboring rays in the medium in general propagate through neighboring parts in the subsurface, and will sense similar velocity structures. There-fore, the traveltimes corresponding to these rays will not vary strongly but will change smoothly along the rays.

(42)

0 0.5 1.0 1.5 2.0 -2000 -1000Offset (m)0 1000 2000 time (s)

Figure 3.6: The response of a source at the reflector at a depth z = 1050m at x = 2700m in

the model of Figure 3.5a that illustrates the wave healing effect. The source is positioned

straight below the fault in the velocity model that shows up at a depth from 8000_{− 1000m.}

The traveltimes of the first arrival correspond to the focusing operator related to the source-location and show no impact of the fault.

parameterization where the number of parameters needed to describe the surface will be significantly smaller then the number of traveltimes.

Different types of parameterizations can be thought of, that utilize the smoothness of the operators. In van de Rijzen and Verhaeghe [2005] a parameterization using orthogonal polynomials was proposed where the weights of the polynomials are the parameters. The use of these orthogonal polynomials, however, has some disadvan-tages. Firstly, the straight-forward evaluation of high-order orthogonal polynomials becomes unstable because of the large dynamical range in polynomial evaluations. Secondly, approximating a function using orthogonal polynomials is a global pa-rameterization. To get a better local approximation, the order of polynomials can become large.

(43)

3.4 Spline parameterization of operators 35

3.4.2 Spline-grids

Having decided on using a spline parameterization to approximate the traveltime surface, we can consider different types of definition grids. Here we will present the two most logical grids.

Fixed definition grid

In the first type of grid, the spatial definition grid of the acquisition coordinate is fixed for all reflector points (Figure 3.7 a) and, therefore, covers for all reflector points the same acquisition range. In this case the parameterization 2D is expressed as τ (xa, xo), where xais the acquisition location that can be either the source location

or the receiver location, and xois the mapped reflector position.

Moving window definition grid

Another parameterization grid is obtained by selecting a moving window region around the reflector coordinate xo that corresponds to the gray region in Figure

3.7a. The resulting definition grid is illustrated in Figure 3.7b, and in this case the operators are parameterized as: τ (xa− xo; xo).

x

o

x

a

x

o a

x−

x

_o a) b)

Figure 3.7: Two possible spline definition grids where the spline control points are indicated

One-step estimation of focusing operators for 3D seismic data

One-step Estimation of

One-step Estimation of

Focusing Operators for 3D Seismic Data

PROEFSCHRIFT

ter verkrijging van de graad van doctor

aan de Technische Universiteit Delft,

op gezag van de Rector Magnificus prof. dr. ir. J.T. Fokkema,

voorzitter van het College voor Promoties,

in het openbaar te verdedigen

op 28 juni 2007 om 12.30 uur

door

Marinus Joseph VAN DE RIJZEN

Contents

1

Introduction

1.1 Short introduction to seismic exploration

1.2 Problem statement

1.3 Research Objectives

1.4 Thesis Outline

2

CFP technology revisited

2.1 Introduction

2.2 Decomposing seismic data in terms of forward propagation

operators

2.3 Focusing operators

2.4 CFP-gathers

2.5 Data-driven focusing operator updating

(Initial)

Focusing

Operator

Seismic Data

CFP−gather

Create

DTS−panel

Create

DTS−panel

Pick

PET met?

no

yes

Final

Focusing

Operator

Update

Focusing Operator

3

One-step focusing operator updating

3.1 Introduction

3.2 Forward model of the CFP-gather synthesis

ξ

3.3 Zero-offset ray mapping

ξ

x

x

ξ

x

t

I

III

II

3.4 Spline parameterization of operators

x

x

x

x−

x

_ξ