Combined imaging and velocity estimation by Joint Migration Inversion

(1)

Combined imaging and velocity estimation

by Joint Migration Inversion

PROEFSCHRIFT

ter verkrijging van de graad van doctor aan de Technische Universiteit Delft,

op gezag van de Rector Magnificus prof. ir. K.C.A.M. Luyben, voorzitter van het College voor Promoties,

in het openbaar te verdedigen op woensdag 11 november 2015 om 12:30 uur

door

Xander Robin STAAL

Ingenieur in de Technische Natuurkunde, Technische Universiteit Delft geboren te Lutong, Maleisie.

(2)

Dit proefschrift is goedgekeurd door de promotor: Prof. dr. ir. L. J. van Vliet

en de copromotor: Dr. ir. D. J. Verschuur

Samenstelling promotiecommissie:

Rector Magnificus, voorzitter

Prof. dr. ir. L.J. van Vliet, Technische Universiteit Delft Dr. ir. D.J. Verschuur, Technische Universiteit Delft Onafhankelijke leden:

Prof. J. Virieux, Universit´e Joseph Fourier

Prof. S.C. Singh, Institut de Physique du Globe de Paris

Prof. dr. ir. E.C. Slob, Technische Universiteit Delft Prof. dr. ir. C. Vuik, Technische Universiteit Delft

ISBN 978-94-6295-382-6

Copyright c_{2015, by X.R. Staal. All rights reserved. No part of this} publication may be reproduced, stored in a retrieval system or transmit-ted in any form or by any means, electronic, mechanical, photocopying, recording or otherwise, without the prior written permission of the au-thor.

Printed in The Netherlands by Proefschriftmaken.nl Typesetting system: LA_TEX.

The research for this thesis was financially supported by the Delphi consortium.

Papa,

ik draag dit proefschift op aan jou,

(3)

Dit proefschrift is goedgekeurd door de promotor: Prof. dr. ir. L. J. van Vliet

en de copromotor: Dr. ir. D. J. Verschuur

Samenstelling promotiecommissie:

Rector Magnificus, voorzitter

Prof. dr. ir. L.J. van Vliet, Technische Universiteit Delft Dr. ir. D.J. Verschuur, Technische Universiteit Delft Onafhankelijke leden:

Prof. J. Virieux, Universit´e Joseph Fourier

Prof. S.C. Singh, Institut de Physique du Globe de Paris

Prof. dr. ir. E.C. Slob, Technische Universiteit Delft Prof. dr. ir. C. Vuik, Technische Universiteit Delft

ISBN 978-94-6295-382-6

Copyright c_{2015, by X.R. Staal. All rights reserved. No part of this} publication may be reproduced, stored in a retrieval system or transmit-ted in any form or by any means, electronic, mechanical, photocopying, recording or otherwise, without the prior written permission of the au-thor.

Printed in The Netherlands by Proefschriftmaken.nl Typesetting system: LA_TEX.

The research for this thesis was financially supported by the Delphi consortium.

Papa,

ik draag dit proefschift op aan jou,

(4)

2 Contents 3.1.1 Reflectivity . . . 26 3.1.2 Slowness . . . 26 3.2 Objective function . . . 27 3.3 Reflectivity update . . . 29 3.4 Slowness update . . . 33 3.5 JMI algorithm . . . 36 4 Numerical examples 41 4.1 1.5-D layered model . . . 41 4.2 2-D lens-shaped model . . . 46 4.3 Marmousi model . . . 49 5 Robustness of JMI 55 5.1 Cycle-skipping . . . 55

5.2 Sensitivity to the starting model . . . 60

5.3 Sensitivity to the offset range . . . 64

6 Field data 69 6.1 The data . . . 69

6.2 Wavelet estimation . . . 70

6.3 Preconditioning . . . 74

6.4 Initial velocity model . . . 75

6.5 Inversion results . . . 77

6.5.1 Modified initial model . . . 78

7 Conclusions and recommendations 85 7.1 Conclusions . . . 85

7.1.1 Multiple scattering . . . 86

7.1.2 Smooth velocity estimates . . . 86

7.1.3 Robustness . . . 87

7.1.4 Hands-off optimisation . . . 88

7.2 Recommendations . . . 88

7.2.1 On Full Wavefield Modelling . . . 89

7.2.2 On the parameterisation . . . 91

7.2.3 On the optimisation algorithm . . . 93

7.2.4 On the input data . . . 94

A Reflectivity operators 97 Contents 3 B Gradient descent scheme 101 B.1 Update direction . . . 102 B.2 Update scaling . . . 102 B.3 Model preconditioning . . . 103 B.3.1 Model Compression . . . 104 B.3.2 Gradient optimisation . . . 105 Bibliography 107

(7)

2 Contents 3.1.1 Reflectivity . . . 26 3.1.2 Slowness . . . 26 3.2 Objective function . . . 27 3.3 Reflectivity update . . . 29 3.4 Slowness update . . . 33 3.5 JMI algorithm . . . 36 4 Numerical examples 41 4.1 1.5-D layered model . . . 41 4.2 2-D lens-shaped model . . . 46 4.3 Marmousi model . . . 49 5 Robustness of JMI 55 5.1 Cycle-skipping . . . 55

5.2 Sensitivity to the starting model . . . 60

5.3 Sensitivity to the offset range . . . 64

6 Field data 69 6.1 The data . . . 69

6.2 Wavelet estimation . . . 70

6.3 Preconditioning . . . 74

6.4 Initial velocity model . . . 75

6.5 Inversion results . . . 77

6.5.1 Modified initial model . . . 78

7 Conclusions and recommendations 85 7.1 Conclusions . . . 85

7.1.1 Multiple scattering . . . 86

7.1.2 Smooth velocity estimates . . . 86

7.1.3 Robustness . . . 87

7.1.4 Hands-off optimisation . . . 88

7.2 Recommendations . . . 88

7.2.1 On Full Wavefield Modelling . . . 89

7.2.2 On the parameterisation . . . 91

7.2.3 On the optimisation algorithm . . . 93

7.2.4 On the input data . . . 94

A Reflectivity operators 97 Contents 3 B Gradient descent scheme 101 B.1 Update direction . . . 102 B.2 Update scaling . . . 102 B.3 Model preconditioning . . . 103 B.3.1 Model Compression . . . 104 B.3.2 Gradient optimisation . . . 105 Bibliography 107

(8)

1

Introduction

Seismic methods aim at extracting information about the earth’s interior from measurements of acoustic/elastic wavefields propagating through it. As Shearer [2009] discusses, seismic methods were originally de-veloped in order to track and understand earthquakes, leading to much increased understanding about the large-scale composition of the earth’s crust, the boundaries of its tectonic plates and the mechanics of earth-quakes. Yilmaz [2007] distinguishes three modern day applications of seismic methods, each with its own zone of interest and (due to absorp-tion) bandwidth of available data:

• earthquake seismology, depths down to 100km, bandwidth roughly 10Hz.

• exploration seismology, depths down to 10km, bandwidth roughly 100Hz.

• engineering seismology, depths down to 1km, bandwidth roughly 1000Hz.

(9)

1

Introduction

Seismic methods aim at extracting information about the earth’s interior from measurements of acoustic/elastic wavefields propagating through it. As Shearer [2009] discusses, seismic methods were originally de-veloped in order to track and understand earthquakes, leading to much increased understanding about the large-scale composition of the earth’s crust, the boundaries of its tectonic plates and the mechanics of earth-quakes. Yilmaz [2007] distinguishes three modern day applications of seismic methods, each with its own zone of interest and (due to absorp-tion) bandwidth of available data:

• earthquake seismology, depths down to 100km, bandwidth roughly 10Hz.

• exploration seismology, depths down to 10km, bandwidth roughly 100Hz.

• engineering seismology, depths down to 1km, bandwidth roughly 1000Hz.

(10)

6 Introduction

Closely related methods using acoustic/elastic wavefields include med-ical diagnostics [Szabo, 2004] and weld inspection [P¨ortzgen, 2007]. These methods typically have much smaller zones of interest and larger data bandwidth.

This thesis focuses on exploration seismology with active sources and an array of detectors, all at or near the surface. We aim at extracting from the seismic data an accurate image of the reflective boundaries in the subsurface under investigation. Oil and gas companies use this type of information to improve their understanding of the geology of the subsurface, allowing them to make informed decisions on drilling new wells, estimating reserves and controlling production rates.

1.1 Migration

The process of mapping seismic data onto an image containing the reflec-tive boundaries in the subsurface is commonly called migration. There exist many flavors of migration with varying levels of assumptions on the structure and complexity of the subsurface. For an overview of various traditional migration algorithms we refer to Claerbout [1976]; Yilmaz [2001]; Biondi [2006]; Robein [2010]. Traditional migration schemes re-quire a smooth (non-reflective) input velocity model and calculate an image based on three steps:

• Forward modelling of a source response at a target location • Reverse modelling of a receiver response at the same target location • Applying an imaging condition to the source and receiver responses. Forward modelling the source response synthesizes the wavefield ar-riving at the target location during the seismic experiment based on knowledge of the source location, its waveform and the velocity model. Reverse modelling the receiver response synthesizes the wavefield leav-ing the target location based on the measured data and the velocity model. When the measured data is pre-processed such that it contains only reflected waves (removing direct waves and turning waves) and the forward modelled source response does not contain reflections, the re-ceiver response should be equal to the source response modulated by a

1.1 Migration 7

reflection coefficient. The imaging condition is aimed at extracting the reflection coefficient from these wavefields, typically by correlating the wavefields. A popular migration tool of this type is Kirchhoff migration, where the operators used for modelling the source and receiver response are based on traveltimes and corresponding amplitude weights. These operators can be obtained by ray-tracing through the smooth velocity model and the three migration steps can be interpreted as a weighted summation of the measured data over traveltime curves. More advanced methods include Wave Equation Migration (WEM) and Reverse Time Migration (RTM) where the full amplitude and phase behavior of the wavefields are taken into account.

A more advanced approach is the concept of least-squares migration [Nemeth et al., 1999], which attempts to replicate the measured seismic data based on three similar steps:

• Forward modelling source responses for all target locations • Modulate the source responses by reflection coefficients

• Forward model the resulting wavefields toward the receiver loca-tions.

The modelled data should then match the measured seismic data in a least-squares sense, migration then becomes a large-scale linear inver-sion scheme with reflection coefficients as the unknowns. Note that this approach requires forward modelling operators that honours the travel-time and amplitude behaviour of wavefields in the subsurface and may require iterative solvers, meaning that it generally requires more com-putational effort than traditional migration schemes. However, Nemeth et al. [1999] and Duquet et al. [2000] show that least-squares migra-tion can suppress migramigra-tion artifacts while boosting resolumigra-tion of the migration results, Plessix and Mulder [2004] show that this approach provides better amplitude estimates in the migration results and Tang and Biondi [2009] show that it provides better results when confronted with more complex blended source wavefields.

The migration approach in this thesis mirrors least-squares migration, where reflectivities are estimated such that the modelled data match the measured seismic data. The method is extended to incorporate surface

(11)

6 Introduction

Closely related methods using acoustic/elastic wavefields include med-ical diagnostics [Szabo, 2004] and weld inspection [P¨ortzgen, 2007]. These methods typically have much smaller zones of interest and larger data bandwidth.

This thesis focuses on exploration seismology with active sources and an array of detectors, all at or near the surface. We aim at extracting from the seismic data an accurate image of the reflective boundaries in the subsurface under investigation. Oil and gas companies use this type of information to improve their understanding of the geology of the subsurface, allowing them to make informed decisions on drilling new wells, estimating reserves and controlling production rates.

1.1 Migration

The process of mapping seismic data onto an image containing the reflec-tive boundaries in the subsurface is commonly called migration. There exist many flavors of migration with varying levels of assumptions on the structure and complexity of the subsurface. For an overview of various traditional migration algorithms we refer to Claerbout [1976]; Yilmaz [2001]; Biondi [2006]; Robein [2010]. Traditional migration schemes re-quire a smooth (non-reflective) input velocity model and calculate an image based on three steps:

• Forward modelling of a source response at a target location • Reverse modelling of a receiver response at the same target location • Applying an imaging condition to the source and receiver responses. Forward modelling the source response synthesizes the wavefield ar-riving at the target location during the seismic experiment based on knowledge of the source location, its waveform and the velocity model. Reverse modelling the receiver response synthesizes the wavefield leav-ing the target location based on the measured data and the velocity model. When the measured data is pre-processed such that it contains only reflected waves (removing direct waves and turning waves) and the forward modelled source response does not contain reflections, the re-ceiver response should be equal to the source response modulated by a

1.1 Migration 7

reflection coefficient. The imaging condition is aimed at extracting the reflection coefficient from these wavefields, typically by correlating the wavefields. A popular migration tool of this type is Kirchhoff migration, where the operators used for modelling the source and receiver response are based on traveltimes and corresponding amplitude weights. These operators can be obtained by ray-tracing through the smooth velocity model and the three migration steps can be interpreted as a weighted summation of the measured data over traveltime curves. More advanced methods include Wave Equation Migration (WEM) and Reverse Time Migration (RTM) where the full amplitude and phase behavior of the wavefields are taken into account.

A more advanced approach is the concept of least-squares migration [Nemeth et al., 1999], which attempts to replicate the measured seismic data based on three similar steps:

• Forward modelling source responses for all target locations • Modulate the source responses by reflection coefficients

• Forward model the resulting wavefields toward the receiver loca-tions.

The modelled data should then match the measured seismic data in a least-squares sense, migration then becomes a large-scale linear inver-sion scheme with reflection coefficients as the unknowns. Note that this approach requires forward modelling operators that honours the travel-time and amplitude behaviour of wavefields in the subsurface and may require iterative solvers, meaning that it generally requires more com-putational effort than traditional migration schemes. However, Nemeth et al. [1999] and Duquet et al. [2000] show that least-squares migra-tion can suppress migramigra-tion artifacts while boosting resolumigra-tion of the migration results, Plessix and Mulder [2004] show that this approach provides better amplitude estimates in the migration results and Tang and Biondi [2009] show that it provides better results when confronted with more complex blended source wavefields.

The migration approach in this thesis mirrors least-squares migration, where reflectivities are estimated such that the modelled data match the measured seismic data. The method is extended to incorporate surface

(12)

8 Introduction

multiples and internal multiples as well as transmission effects while still being based on a non-reflecting velocity model, making it a nonlinear inversion scheme.

1.1.1 Dealing with surface multiples

The dominant form of multiple scattering in most seismic surveys are surface-related multiples, being all multiple scattering events that have reflected at the surface at least once. This class of multiples has there-fore received the most attention over the years. The traditional way of dealing with surface related multiples is to predict and remove them from the dataset before migration. For an overview of commonly used multiple removal algorithms we refer to Weglein and Dragoset [2005]; Verschuur [2006]. Notably Surface-Related Multiple Elimination (SRME) [Verschuur et al., 1992] exploits the fact that seismic data itself is ac-quired at the surface and therefore contains an implicit measurement of the Green’s function needed to model the surface multiples based on the measured data. This method is remarkably effective and is completely data-driven and deterministic, although it does require a well-sampled measurement of the wavefields in order to correctly predict multiples and remove them.

Recently, the SRME concept has been recast into an inversion process, called Estimation of Primaries by Sparse Inversion (EPSI) [van Groen-estijn and Verschuur, 2009], which aims not at removing multiples but at estimating the Green’s functions that simultaneously explains both the primary wavefields and the multiples. An interesting consequence of this emphasis shift is that the method is less demanding on having well-sampled measurements as shown by van Groenestijn and Verschuur [2009], who show that EPSI can use surface multiples to reconstruct pri-maries in a near-offset gap in a seismic dataset. The idea that multiples in a seismic dataset may contain information that is missing from the primaries is extensively discussed by Berkhout et al. [2012], who show the promise of improved illumination that may be achieved by including surface multiples in the migration process rather than removing them from the dataset. This gain in illumination from multiples is convinc-ingly shown for VSP datasets by Soni [2014].

1.1 Migration 9

Methods that include surface multiples in the migration process have ex-isted for some time now, being shown by Berkhout and Verschuur [1994]; Brown and Guitton [2005]; Whitmore et al. [2010]; Zheng and Schus-ter [2014]; Tu and Herrmann [2015]. Similar to SRME, these method use the measured seismic dataset itself to extend the original modelling operator and synthesize the source response at the target location, this time including surface multiples. Typically, due to the increased com-plexity of the wavefields involved, the traditional correlation-type imag-ing will generate artifacts and least-squares migration is to be preferred.

1.1.2 Dealing with internal multiples

Internal multiples form the collection of multiple scattering events that have not reflected at the surface. These multiples are generally much weaker in amplitude than surface multiples and primaries and do not cause significant artifacts in the migration results. Nonetheless, several methods exist to predict and remove internal multiples from seismic data, including Internal Multiple Elimination (IME) [Jakubowicz, 1998; Berkhout and Verschuur, 2005; Verschuur and Berkhout, 2005], which is an extension of SRME that requires some user-specified horizons or mute functions and removes only multiples crossing the chosen horizon but is otherwise fully data-driven, and Inverse Scattering Series (ISS) [Weglein et al., 1997; Matson et al., 1999], which removes all classes of internal multiples and is fully data-driven. As with SRME, both IME and ISS are deterministic processes and require a well-sampled seismic dataset to work correctly. The IME method has also been cast into an inversion process named generalised EPSI [Ypma and Verschuur, 2013], which again relaxes its requirements for having well-sampled input data and promises the ability to reconstruct primary data that was missing in the measured seismic data.

A promising way of including both surface multiples and internal multi-ples in migration is Marchenko imaging [Behura et al., 2012; Slob et al., 2014; Wapenaar et al., 2014], which aims to include scattering into the Green’s functions such that they can be used for proper migration of primaries and multiples simultaneously. The method requires an initial estimate of the Green’s functions, which can be provided by the tra-ditional modelling operators, and updates them in a deterministic and

(13)

8 Introduction

multiples and internal multiples as well as transmission effects while still being based on a non-reflecting velocity model, making it a nonlinear inversion scheme.

1.1.1 Dealing with surface multiples

The dominant form of multiple scattering in most seismic surveys are surface-related multiples, being all multiple scattering events that have reflected at the surface at least once. This class of multiples has there-fore received the most attention over the years. The traditional way of dealing with surface related multiples is to predict and remove them from the dataset before migration. For an overview of commonly used multiple removal algorithms we refer to Weglein and Dragoset [2005]; Verschuur [2006]. Notably Surface-Related Multiple Elimination (SRME) [Verschuur et al., 1992] exploits the fact that seismic data itself is ac-quired at the surface and therefore contains an implicit measurement of the Green’s function needed to model the surface multiples based on the measured data. This method is remarkably effective and is completely data-driven and deterministic, although it does require a well-sampled measurement of the wavefields in order to correctly predict multiples and remove them.

Recently, the SRME concept has been recast into an inversion process, called Estimation of Primaries by Sparse Inversion (EPSI) [van Groen-estijn and Verschuur, 2009], which aims not at removing multiples but at estimating the Green’s functions that simultaneously explains both the primary wavefields and the multiples. An interesting consequence of this emphasis shift is that the method is less demanding on having well-sampled measurements as shown by van Groenestijn and Verschuur [2009], who show that EPSI can use surface multiples to reconstruct pri-maries in a near-offset gap in a seismic dataset. The idea that multiples in a seismic dataset may contain information that is missing from the primaries is extensively discussed by Berkhout et al. [2012], who show the promise of improved illumination that may be achieved by including surface multiples in the migration process rather than removing them from the dataset. This gain in illumination from multiples is convinc-ingly shown for VSP datasets by Soni [2014].

1.1 Migration 9

Methods that include surface multiples in the migration process have ex-isted for some time now, being shown by Berkhout and Verschuur [1994]; Brown and Guitton [2005]; Whitmore et al. [2010]; Zheng and Schus-ter [2014]; Tu and Herrmann [2015]. Similar to SRME, these method use the measured seismic dataset itself to extend the original modelling operator and synthesize the source response at the target location, this time including surface multiples. Typically, due to the increased com-plexity of the wavefields involved, the traditional correlation-type imag-ing will generate artifacts and least-squares migration is to be preferred.

1.1.2 Dealing with internal multiples

Internal multiples form the collection of multiple scattering events that have not reflected at the surface. These multiples are generally much weaker in amplitude than surface multiples and primaries and do not cause significant artifacts in the migration results. Nonetheless, several methods exist to predict and remove internal multiples from seismic data, including Internal Multiple Elimination (IME) [Jakubowicz, 1998; Berkhout and Verschuur, 2005; Verschuur and Berkhout, 2005], which is an extension of SRME that requires some user-specified horizons or mute functions and removes only multiples crossing the chosen horizon but is otherwise fully data-driven, and Inverse Scattering Series (ISS) [Weglein et al., 1997; Matson et al., 1999], which removes all classes of internal multiples and is fully data-driven. As with SRME, both IME and ISS are deterministic processes and require a well-sampled seismic dataset to work correctly. The IME method has also been cast into an inversion process named generalised EPSI [Ypma and Verschuur, 2013], which again relaxes its requirements for having well-sampled input data and promises the ability to reconstruct primary data that was missing in the measured seismic data.

A promising way of including both surface multiples and internal multi-ples in migration is Marchenko imaging [Behura et al., 2012; Slob et al., 2014; Wapenaar et al., 2014], which aims to include scattering into the Green’s functions such that they can be used for proper migration of primaries and multiples simultaneously. The method requires an initial estimate of the Green’s functions, which can be provided by the tra-ditional modelling operators, and updates them in a deterministic and

(14)

10 Introduction

data-driven scheme. The method then uses the updated Green’s func-tions to synthesize the source and receiver response in such a way that multiples map to the correct subsurface location, meaning that they will not cause migration artifacts. The image value can then be extracted by correlation or least-squares matching.

This thesis follows the migration approach of Full Wavefield Migra-tion (FWM) [Davydenko and Verschuur, 2012; Soni, 2014; Berkhout, 2014b], which also aims to include scattering into the modelling op-erators such that they can be used for proper migration of primaries, surface multiples and internal multiples simultaneously. In contrast to the Marchenko imaging approach, FWM takes a model-driven approach to estimating its modelling operators, basing them on a reflectivity and a velocity model through a Full Wavefield Modelling (FWMod) algo-rithm [Berkhout, 2014a]. The velocity model is assumed to be given, but the reflectivity model is iteratively updated in such a way that the modelled source response will match the measured data. Note that the modelling operators are a function of the estimated reflectivity model and vice versa, making this a nonlinear migration scheme. Since the proper migration of internal multiples is the truly novel aspect of this approach, we have chosen to fully focus on those type of multiples and we will use datasets where surface multiples are not an issue.

1.2 Velocity model estimation

All aforementioned migration schemes depend on an estimate of the wavefield propagation velocity in the subsurface to some degree. Specif-ically those migration schemes that output their reflectivity maps in the depth domain will suffer from poor velocity estimates, failing to map reflectivity estimates to their correct depth or even to focus the energy resulting in blurry images. Overviews of traditional velocity model esti-mation methods are given by Yilmaz [2001]; Robein [2003]; Jones [2010]. Velocity model estimation methods can be roughly divided into two cat-egories: image-domain optimisation methods including Wave Equation Migration Velocity Analysis (WEMVA) and data-domain optimisation methods including Full Waveform Inversion (FWI).

1.2 Velocity model estimation 11

1.2.1 Wave Equation Migration Velocity Analysis

Wave Equation Migration Velocity Analysis (WEMVA) [Sava and Biondi, 2004a,b] aims to optimise the modelling operators used for migration by updating the velocity model. The method extends regular migration to estimate Common Image Gathers (CIGs), for instance by estimating reflectivity as a function of reflection angle (ADCIGs) [Sava and Fomel, 2003; Biondi and Symes, 2004]. These CIGs have some property that indicates velocity model inaccuracies, events in ADCIGs appear curved when the velocity model is incorrect. WEMVA then iteratively updates the velocity model such that the value of that property is optimised, typically by maximising semblance [Symes and Carazzone, 1991] when using ADCIGs. An overview of several types of CIGs is given by Sava and Vasconselos [2011].

An important feature of WEMVA is that it uses a one-way algorithm to generate its modelling operators for migration, meaning that the veloc-ity model cannot generate scattering in the operators. As a consequence, multiples are not taken into account properly by the algorithm, but the method is very stable because the velocity updates are fully driven by the kinematic (low wavenumber) errors in the model.

1.2.2 Full Waveform Inversion

FWI [Tarantola, 1984; Virieux and Operto, 2009] aims to match the measured seismic data with modelled data based on the wave equa-tion. Typically, FWI uses as the modelling algorithm a finite-difference algorithm based on the acoustic wave equation with constant density, meaning that velocity is the only variable to be estimated. Theoret-ically, however, the method can be easily extended to estimate more elaborate medium parameters based on the elastic wave equation, as shown by Mora [1988]; Crase et al. [1990]; Sears et al. [2010]; Rizzuti and Gisolf [2014].

Contrary to WEMVA, the velocity model in FWI does generate scat-tering in its modelling operators. As a consequence, the method can properly deal with multiple scattering, but the stability of the method suffers because the velocity updates are no longer fully driven by the

(15)

10 Introduction

data-driven scheme. The method then uses the updated Green’s func-tions to synthesize the source and receiver response in such a way that multiples map to the correct subsurface location, meaning that they will not cause migration artifacts. The image value can then be extracted by correlation or least-squares matching.

This thesis follows the migration approach of Full Wavefield Migra-tion (FWM) [Davydenko and Verschuur, 2012; Soni, 2014; Berkhout, 2014b], which also aims to include scattering into the modelling op-erators such that they can be used for proper migration of primaries, surface multiples and internal multiples simultaneously. In contrast to the Marchenko imaging approach, FWM takes a model-driven approach to estimating its modelling operators, basing them on a reflectivity and a velocity model through a Full Wavefield Modelling (FWMod) algo-rithm [Berkhout, 2014a]. The velocity model is assumed to be given, but the reflectivity model is iteratively updated in such a way that the modelled source response will match the measured data. Note that the modelling operators are a function of the estimated reflectivity model and vice versa, making this a nonlinear migration scheme. Since the proper migration of internal multiples is the truly novel aspect of this approach, we have chosen to fully focus on those type of multiples and we will use datasets where surface multiples are not an issue.

1.2 Velocity model estimation

All aforementioned migration schemes depend on an estimate of the wavefield propagation velocity in the subsurface to some degree. Specif-ically those migration schemes that output their reflectivity maps in the depth domain will suffer from poor velocity estimates, failing to map reflectivity estimates to their correct depth or even to focus the energy resulting in blurry images. Overviews of traditional velocity model esti-mation methods are given by Yilmaz [2001]; Robein [2003]; Jones [2010]. Velocity model estimation methods can be roughly divided into two cat-egories: image-domain optimisation methods including Wave Equation Migration Velocity Analysis (WEMVA) and data-domain optimisation methods including Full Waveform Inversion (FWI).

1.2.1 Wave Equation Migration Velocity Analysis

Wave Equation Migration Velocity Analysis (WEMVA) [Sava and Biondi, 2004a,b] aims to optimise the modelling operators used for migration by updating the velocity model. The method extends regular migration to estimate Common Image Gathers (CIGs), for instance by estimating reflectivity as a function of reflection angle (ADCIGs) [Sava and Fomel, 2003; Biondi and Symes, 2004]. These CIGs have some property that indicates velocity model inaccuracies, events in ADCIGs appear curved when the velocity model is incorrect. WEMVA then iteratively updates the velocity model such that the value of that property is optimised, typically by maximising semblance [Symes and Carazzone, 1991] when using ADCIGs. An overview of several types of CIGs is given by Sava and Vasconselos [2011].

An important feature of WEMVA is that it uses a one-way algorithm to generate its modelling operators for migration, meaning that the veloc-ity model cannot generate scattering in the operators. As a consequence, multiples are not taken into account properly by the algorithm, but the method is very stable because the velocity updates are fully driven by the kinematic (low wavenumber) errors in the model.

1.2.2 Full Waveform Inversion

FWI [Tarantola, 1984; Virieux and Operto, 2009] aims to match the measured seismic data with modelled data based on the wave equa-tion. Typically, FWI uses as the modelling algorithm a finite-difference algorithm based on the acoustic wave equation with constant density, meaning that velocity is the only variable to be estimated. Theoret-ically, however, the method can be easily extended to estimate more elaborate medium parameters based on the elastic wave equation, as shown by Mora [1988]; Crase et al. [1990]; Sears et al. [2010]; Rizzuti and Gisolf [2014].

Contrary to WEMVA, the velocity model in FWI does generate scat-tering in its modelling operators. As a consequence, the method can properly deal with multiple scattering, but the stability of the method suffers because the velocity updates are no longer fully driven by the

(16)

12 Introduction

kinematic (low wavenumber) errors in the model but are now dominated by the reflecting (high wavenumber) errors in the model, as discussed by Symes [2008]; Plessix [2012].

A typical solution to this problem is to apply FWI only to transmitted waves such as diving waves. Symes [2008] notes that when FWI is applied to transmitted waves one runs the risk of cycle-skipping, but this can be avoided when the data contains low enough frequencies. Using FWI on transmitted waves to estimate a velocity model requires very long offsets (to measure diving waves that have propagated through the target depths) and very low frequencies. When supplied with suitable data, this approach can produce very detailed and accurate velocity models as shown by Operto et al. [2004]; Plessix et al. [2010].

A very interesting approach was presented by Tang et al. [2013] who greatly improve the stability of FWI on reflection data by decomposing the velocity update into a tomographic term and a migration term. The tomographic term is obtained by correlating only those parts of the source and receiver responses that travel in the same direction while the migration term is obtained from those parts that travel in opposite directions. The two terms can then be appropriately balanced before they are combined in the final update.

Alternative approaches are to match only adaptively selected portions of the measured data as described by Bi and Lin [2014] or to include an adaptive filter in the matching of measured and modelled data as described by Warner and Guasch [2014].

1.2.3 Hybrid methods

Recently we have seen attempts to combine the concepts of WEMVA and FWI into a single inversion scheme by Fleury and Perrone [2012] and Biondi and Almomin [2012], where the FWI data misfit criterion is complemented by a model domain criterion similar to WEMVA. To achieve this, the FWI process is adapted to estimate velocities in an extended model space similar to the CIGs discussed before. The data-domain criterion and the model-data-domain criterion can then be weighted to achieve optimal convergence properties.

Another strategy was shown by Zhou et al. [2012]; Alkhalifah [2014]; Zhou et al. [2014], who added an explicit migration/demigration step and a corresponding reflectivity model to the standard FWI formulation. Their total modelled data then consists of the modelled data based on the wave equation estimate plus the demigrated reflectivity estimate as discussed by Symes and Kern [1994].

1.2.4 Joint Migration Inversion

This thesis presents a third strategy known as Joint Migration Inversion (JMI) [Staal and Verschuur, 2013; Berkhout, 2014c]. A data-domain optimisation method that uses the FWMod algorithm to model data based on a velocity model and a reflectivity model. The method is an extension of FWM as described above, adding a velocity update to the nonlinear migration algorithm. In the FWMod algorithm, the veloc-ity model only affects the kinematics and does not generate scattering in the modelling operators, while the reflectivity model only generates scattering (including transmission effects, surface multiples and internal multiples) and does not affect kinematics. This gives JMI the stability of WEMVA, while avoiding the nonlinearities found in FWI due to the dual role of the velocity model as described before.

(17)

12 Introduction

kinematic (low wavenumber) errors in the model but are now dominated by the reflecting (high wavenumber) errors in the model, as discussed by Symes [2008]; Plessix [2012].

A typical solution to this problem is to apply FWI only to transmitted waves such as diving waves. Symes [2008] notes that when FWI is applied to transmitted waves one runs the risk of cycle-skipping, but this can be avoided when the data contains low enough frequencies. Using FWI on transmitted waves to estimate a velocity model requires very long offsets (to measure diving waves that have propagated through the target depths) and very low frequencies. When supplied with suitable data, this approach can produce very detailed and accurate velocity models as shown by Operto et al. [2004]; Plessix et al. [2010].

A very interesting approach was presented by Tang et al. [2013] who greatly improve the stability of FWI on reflection data by decomposing the velocity update into a tomographic term and a migration term. The tomographic term is obtained by correlating only those parts of the source and receiver responses that travel in the same direction while the migration term is obtained from those parts that travel in opposite directions. The two terms can then be appropriately balanced before they are combined in the final update.

Alternative approaches are to match only adaptively selected portions of the measured data as described by Bi and Lin [2014] or to include an adaptive filter in the matching of measured and modelled data as described by Warner and Guasch [2014].

1.2.3 Hybrid methods

Recently we have seen attempts to combine the concepts of WEMVA and FWI into a single inversion scheme by Fleury and Perrone [2012] and Biondi and Almomin [2012], where the FWI data misfit criterion is complemented by a model domain criterion similar to WEMVA. To achieve this, the FWI process is adapted to estimate velocities in an extended model space similar to the CIGs discussed before. The data-domain criterion and the model-data-domain criterion can then be weighted to achieve optimal convergence properties.

Another strategy was shown by Zhou et al. [2012]; Alkhalifah [2014]; Zhou et al. [2014], who added an explicit migration/demigration step and a corresponding reflectivity model to the standard FWI formulation. Their total modelled data then consists of the modelled data based on the wave equation estimate plus the demigrated reflectivity estimate as discussed by Symes and Kern [1994].

1.2.4 Joint Migration Inversion

This thesis presents a third strategy known as Joint Migration Inversion (JMI) [Staal and Verschuur, 2013; Berkhout, 2014c]. A data-domain optimisation method that uses the FWMod algorithm to model data based on a velocity model and a reflectivity model. The method is an extension of FWM as described above, adding a velocity update to the nonlinear migration algorithm. In the FWMod algorithm, the veloc-ity model only affects the kinematics and does not generate scattering in the modelling operators, while the reflectivity model only generates scattering (including transmission effects, surface multiples and internal multiples) and does not affect kinematics. This gives JMI the stability of WEMVA, while avoiding the nonlinearities found in FWI due to the dual role of the velocity model as described before.

(18)

14 Introduction 1.3 Thesis overview

Chapter 2 discusses the Full Wavefield Modelling tool used to generate seismic data based on a velocity model and a reflectivity model. As mentioned previously, the velocity model only affects the kinematics of the modelled data without generating scattering, while the reflectivity model generates all orders of scattering without affecting the kinematics of the data. The method is depth recursive and iterative, each iteration adding a new order of multiple scattering to the modelled data.

Chapter 3 discusses the Joint Migration Inversion algorithm. First, spe-cific parameterisation of the more general FWMod modelling operators are explained. Then the data-misfit criterion is explained, followed by a derivation of the updates for the reflectivity and the velocity (or rather slowness) models. Finally, an algorithmic overview of JMI is discussed. Chapter 4 applies JMI to numerical examples. The inversion scheme is illustrated using a 1.5D example, showing the the inversion strategy including stopping criteria and frequency bandwidth. 2D examples show the method’s ability to correctly incorporate internal multiples and its ability to hande complex subsurfaces.

Chapter 5 discusses the robustness of JMI. As mentioned previously, the separation of kinematic effects (being explained exclusively by the velocity model) and the scattering effects (being explained exclusively by the reflectivity model) gives JMI the stability of WEMVA, while being able to correctly deal with multiple scattering. This idea is explored further, by investigating the possibility of cycle-skipping in JMI, by investigating its sensitivity to the starting model and by investigating its dependence on offset ranges.

Chapter 6 applies JMI to a marine field dataset from the Vøring basin in Norway, provided by Statoil. Necessary pre-processing steps are dis-cussed, including wavelet estimation, data preconditioning and initial velocity model building. Results from JMI are shown and evaluated using adCIG’s, showing that the algorithm performs well on field data. Chapter 7 discusses conclusions drawn from this research and gives some recommendations for future research on the topic of Joint Migration Inversion.

2

Full Waveﬁeld Modelling

Full Wavefield Modelling does not attempt to solve the differential form of the wave equation directly, instead it uses integral operators that evaluate different aspects of wave propagation in the subsurface. This approach allows us to forward model seismic data whose kinematics are based exclusively on a background velocity model and whose reflection amplitudes depend exclusively on a reflectivity model of the subsurface. This in turn will allow us to define an effective inversion scheme that estimates a kinematically accurate velocity model as well as a true-amplitude reflectivity map of the subsurface.

2.1 Waveﬁeld operations

2.1.1 Extrapolation

Propagation of a pressure wave in a source-free half-space of an acoustic medium can be described by wavefield extrapolation, which is governed by the Rayleigh II equation [Rayleigh, 1878; Berkhout, 1982; Gisolf and

(19)

14 Introduction 1.3 Thesis overview

Chapter 2 discusses the Full Wavefield Modelling tool used to generate seismic data based on a velocity model and a reflectivity model. As mentioned previously, the velocity model only affects the kinematics of the modelled data without generating scattering, while the reflectivity model generates all orders of scattering without affecting the kinematics of the data. The method is depth recursive and iterative, each iteration adding a new order of multiple scattering to the modelled data.

Chapter 3 discusses the Joint Migration Inversion algorithm. First, spe-cific parameterisation of the more general FWMod modelling operators are explained. Then the data-misfit criterion is explained, followed by a derivation of the updates for the reflectivity and the velocity (or rather slowness) models. Finally, an algorithmic overview of JMI is discussed. Chapter 4 applies JMI to numerical examples. The inversion scheme is illustrated using a 1.5D example, showing the the inversion strategy including stopping criteria and frequency bandwidth. 2D examples show the method’s ability to correctly incorporate internal multiples and its ability to hande complex subsurfaces.

Chapter 5 discusses the robustness of JMI. As mentioned previously, the separation of kinematic effects (being explained exclusively by the velocity model) and the scattering effects (being explained exclusively by the reflectivity model) gives JMI the stability of WEMVA, while being able to correctly deal with multiple scattering. This idea is explored further, by investigating the possibility of cycle-skipping in JMI, by investigating its sensitivity to the starting model and by investigating its dependence on offset ranges.

Chapter 6 applies JMI to a marine field dataset from the Vøring basin in Norway, provided by Statoil. Necessary pre-processing steps are dis-cussed, including wavelet estimation, data preconditioning and initial velocity model building. Results from JMI are shown and evaluated using adCIG’s, showing that the algorithm performs well on field data. Chapter 7 discusses conclusions drawn from this research and gives some recommendations for future research on the topic of Joint Migration Inversion.

2

Full Waveﬁeld Modelling

Full Wavefield Modelling does not attempt to solve the differential form of the wave equation directly, instead it uses integral operators that evaluate different aspects of wave propagation in the subsurface. This approach allows us to forward model seismic data whose kinematics are based exclusively on a background velocity model and whose reflection amplitudes depend exclusively on a reflectivity model of the subsurface. This in turn will allow us to define an effective inversion scheme that estimates a kinematically accurate velocity model as well as a true-amplitude reflectivity map of the subsurface.

2.1 Waveﬁeld operations

2.1.1 Extrapolation

Propagation of a pressure wave in a source-free half-space of an acoustic medium can be described by wavefield extrapolation, which is governed by the Rayleigh II equation [Rayleigh, 1878; Berkhout, 1982; Gisolf and

(20)

16 Full Waveﬁeld Modelling

Verschuur, 2010], that reads for the 2D case: P (xn, zn) = 2sign(zm− zn) +∞ −∞ ρ(xn, zn) ρ(x, zm) ∂G ∂z(xn, zn; x, zm)P (x, zm)dx, (2.1) where P (x, z) represents a single frequency component of the pressure wavefield, ρ(x, z) is the density of the medium and ∂G_∂z(xs, zs; xd, zd)

represents a Green’s function which has a monopole source at position (xs, zs) and is recorded by a dipole receiver at (xr, zr). In discretised

computations this integral can be expressed as a vector-vector inner product or, if we extrapolate the wavefield to all lateral positions along a depth level, as a matrix-vector product:

�

P (zn) = W(zn; zm) �P (zm), (2.2)

where �P (zn) represents the wavefield measured along depth level zn.

The matrix W(zn; zm) is a propagation operator in which each index

can be related to the dipole’s Green function as follows [Berkhout, 1982]: Wij(zn; zm) = 2sign(zm− zn)

ρ(xi, zn)

ρ(xj, zm)

∂G

∂z (xi, zn; xj, zm). (2.3) In case of a homogeneous half-space below zm, we can further relate the

ith column of the propagation operator to the phase-shift operator in the kx domain [Gazdag, 1978]:

˜

w(zn; zm) = e−jkz|zn−zm|ejkxxi, (2.4a)

where xi is the source position of the Green’s function and

kz = k2_{− k}2 x, |kx| ≤ |k| −jk2 x− k2, |kx| > |k| (2.4b) with the wavenumber k = ω

c.

In our approach, we extend the one-way wavefield decomposition as discussed by Wapenaar [1998] by splitting the wavefield into four parts at any depth level zn:

• a downgoing wavefield P+_(z

n) that propagates as if the medium is

homogeneous and source-free below depth zn− ∆z

2.1 Waveﬁeld operations 17 n z 1 n z  1 n z 

 

n1 Q z 

 

n1 Q z 

 

n P z

 

n P z



z zn; n 1



  W



z zn; n 1



  W

Figure 2.1: Schematic representation of wavefield extrapolation.

• a downgoing wavefield Q+_(z

homogeneous and source-free below depth zn

• an upgoing wavefield P−_(z

homogeneous and source-free above depth zn+ ∆z

• an upgoing wavefield Q−_(z

homogeneous and source-free above depth zn

where ∆z is the depth over which we want to propagate the wavefield, in equation 2.1 this would equate to ∆z = |zm− zn|. When the medium

is laterally homogeneous in between depth levels zm and zn, we can

reconstruct the total wavefield P = P+ _{+ Q}− _{= P}−_{+ Q}+_{. Note that}

these assumptions allow us to rewrite equation 2.2 as follows: �

P+(zn) = W+(zn; zn−1) �Q+(zn−1) (2.5a)

�

P−(zn) = W−(zn; zn+1) �Q−(zn+1), (2.5b)

where each column of the propagation operators W±_(z

n, zn ± ∆z) is

calculated using the phase-shift operators from equation 2.4a using the local velocity c(xi, zn ± ∆z₂ ). For simplicity we will use the notation

zn+ ∆z₂ = zn+ in the rest of this thesis. The extrapolation procedure is

illustrated using rays in figure 2.1. Next we will relate the wavefields P+/− _{and Q}+/− _{through reflectivity operators.}

(21)

Verschuur, 2010], that reads for the 2D case: P (xn, zn) = 2sign(zm− zn) +∞ −∞ ρ(xn, zn) ρ(x, zm) ∂G ∂z(xn, zn; x, zm)P (x, zm)dx, (2.1) where P (x, z) represents a single frequency component of the pressure wavefield, ρ(x, z) is the density of the medium and ∂G_∂z(xs, zs; xd, zd)

represents a Green’s function which has a monopole source at position (xs, zs) and is recorded by a dipole receiver at (xr, zr). In discretised

computations this integral can be expressed as a vector-vector inner product or, if we extrapolate the wavefield to all lateral positions along a depth level, as a matrix-vector product:

�

P (zn) = W(zn; zm) �P (zm), (2.2)

where �P (zn) represents the wavefield measured along depth level zn.

The matrix W(zn; zm) is a propagation operator in which each index

can be related to the dipole’s Green function as follows [Berkhout, 1982]: Wij(zn; zm) = 2sign(zm− zn)

ρ(xi, zn)

ρ(xj, zm)

∂G

∂z(xi, zn; xj, zm). (2.3) In case of a homogeneous half-space below zm, we can further relate the

ith column of the propagation operator to the phase-shift operator in the kx domain [Gazdag, 1978]:

˜

w(zn; zm) = e−jkz|zn−zm|ejkxxi, (2.4a)

where xi is the source position of the Green’s function and

kz = k2 _{− k}2 x, |kx| ≤ |k| −jk2 x− k2, |kx| > |k| (2.4b) with the wavenumber k = ω

c.

In our approach, we extend the one-way wavefield decomposition as discussed by Wapenaar [1998] by splitting the wavefield into four parts at any depth level zn:

• a downgoing wavefield P+_(z

homogeneous and source-free below depth zn− ∆z

2.1 Waveﬁeld operations 17 n z 1 n z  1 n z 

 

n 1 Q z 

 

n 1 Q z 

 

n P z

 

n P z



z zn; n 1



  W



z zn; n 1



  W

Figure 2.1: Schematic representation of wavefield extrapolation.

• a downgoing wavefield Q+_(z

homogeneous and source-free below depth zn

• an upgoing wavefield P−_(z

homogeneous and source-free above depth zn+ ∆z

• an upgoing wavefield Q−_(z

homogeneous and source-free above depth zn

where ∆z is the depth over which we want to propagate the wavefield, in equation 2.1 this would equate to ∆z = |zm− zn|. When the medium

is laterally homogeneous in between depth levels zm and zn, we can

reconstruct the total wavefield P = P+_{+ Q}− _{= P}−_{+ Q}+_{. Note that}

these assumptions allow us to rewrite equation 2.2 as follows: �

P+(zn) = W+(zn; zn−1) �Q+(zn−1) (2.5a)

�

P−(zn) = W−(zn; zn+1) �Q−(zn+1), (2.5b)

where each column of the propagation operators W±_(z

n, zn ± ∆z) is

calculated using the phase-shift operators from equation 2.4a using the local velocity c(xi, zn ± ∆z₂ ). For simplicity we will use the notation

zn+ ∆z₂ = zn+ in the rest of this thesis. The extrapolation procedure is

illustrated using rays in figure 2.1. Next we will relate the wavefields P+/− _{and Q}+/− _{through reflectivity operators.}

(22)

 

n Q z

 

n Q z

 

n P z

 

n P z

 

zn  R

 

zn  R T

 

zn

 

zn  T n z 1 n z  1 n z 

Figure 2.2: Schematic representation of reflection and transmission at a single depth level.

2.1.2 Reﬂection

When a seismic wave encounters a sharp discontinuity in the subsurface, it will scatter. Part of the energy in the wave will be reflected and part of the energy will be transmitted. These effects can be neatly described by their corresponding operators. As illustrated in figure 2.2, when a wave is incident on a discontinuity we can write:

� Q+_(z n) = T+(zn) �P+(zn) + R∩(zn) �P−(zn) (2.6a) � Q−(zn) = R∪(zn) �P+(zn) + T−(zn) �P−(zn), (2.6b) where �Q+_(z

n) and �Q−(zn) are, respectively, the downgoing and upgoing

wavefield moving away from the discontinuity. Wavefields �P+_(z n) and

� P−_(z

n) are the downgoing and upgoing wavefields incident on the

dis-continuity. Operator R∪_(z

n) is the reflectivity operator associated with

the discontinuities at depth zn, acting on waves incident from above.

Re-flectivity operator R∩_(z

n) acts on waves coming from below and T+(zn)

and T−_(z

n) are the corresponding transmission operators.

In the acoustic case, the total wavefields above and below the depth level zn should be continuous, which results in the following expressions

2.1 Waveﬁeld operations 19

for the transmission coefficients in equations 2.6a,b:

T+ = I + R∪ (2.7a)

T− = I + R∩. (2.7b)

In acoustic media, we also know that R∩_(z

n) = −R∪(zn), which means

that all reflection events can be fully described by the single operator R = R∪ _{= −R}∩_.

Notice that this description of wavefield scattering is only valid for re-flecting boundaries that are horizontally oriented. We will assume gently dipping structures can be accommodated in this framework by using an effective reflectivity operator that compensates for the dip angle. When we are dealing with steeply dipping reflectors however, a downgoing wavefield may not be transformed into an upgoing wavefield after re-flection, so the equations shown above do not provide a full description of scattering in the subsurface. In this thesis, however, we will limit our-selves to the vertically-oriented scattering equations 2.6a,b. Note that Davydenko and Verschuur [2013b] describe a methodology to extend FWM to correctly handle steeply dipping reflectors.

Also notice that we do not supply a deterministic link between the medium parameters in the subsurface and the corresponding reflection operators. In this work we merely observe that reflections will occur at discontinuities in the subsurface and we will estimate the correspond-ing operators directly. Nevertheless, a derivation for acoustic reflection operators in a 1.5D medium is given in appendix A.

2.1.3 Source waveﬁelds

Wavefields in the subsurface are generated by a source term as illus-trated in figure 2.3. When a physical source is fired at depth level zn,

we simply extend equations 2.5a,b by adding a source wavefield �S±_(z

n)

to the scattered wavefields leaving that depth level: �

Q+_(z

n) = �S+(zn) + T+(zn) �P+(zn) + R∩(zn) �P−(zn) (2.8a)

�

Q−(zn) = �S−(zn) + T−(zn) �P−(zn) + R∪(zn) �P+(zn). (2.8b)

Notice that there should be no propagation and scattering effects in-cluded in the source wavefields �S, as these effects are already taken into

(23)

 

n Q z

 

n Q z

 

n P z

 

n P z

 

zn  R

 

zn  R T

 

zn

 

zn  T n z 1 n z  1 n z 

Figure 2.2: Schematic representation of reflection and transmission at a single depth level.

2.1.2 Reﬂection

When a seismic wave encounters a sharp discontinuity in the subsurface, it will scatter. Part of the energy in the wave will be reflected and part of the energy will be transmitted. These effects can be neatly described by their corresponding operators. As illustrated in figure 2.2, when a wave is incident on a discontinuity we can write:

� Q+_(z n) = T+(zn) �P+(zn) + R∩(zn) �P−(zn) (2.6a) � Q−(zn) = R∪(zn) �P+(zn) + T−(zn) �P−(zn), (2.6b) where �Q+_(z

n) and �Q−(zn) are, respectively, the downgoing and upgoing

wavefield moving away from the discontinuity. Wavefields �P+_(z n) and

� P−_(z

n) are the downgoing and upgoing wavefields incident on the

dis-continuity. Operator R∪_(z

n) is the reflectivity operator associated with

the discontinuities at depth zn, acting on waves incident from above.

Re-flectivity operator R∩_(z

n) acts on waves coming from below and T+(zn)

and T−_(z

n) are the corresponding transmission operators.

In the acoustic case, the total wavefields above and below the depth level zn should be continuous, which results in the following expressions

2.1 Waveﬁeld operations 19

for the transmission coefficients in equations 2.6a,b:

T+= I + R∪ (2.7a)

T−= I + R∩. (2.7b)

In acoustic media, we also know that R∩_(z

n) = −R∪(zn), which means

that all reflection events can be fully described by the single operator R = R∪ _{= −R}∩_.

Notice that this description of wavefield scattering is only valid for re-flecting boundaries that are horizontally oriented. We will assume gently dipping structures can be accommodated in this framework by using an effective reflectivity operator that compensates for the dip angle. When we are dealing with steeply dipping reflectors however, a downgoing wavefield may not be transformed into an upgoing wavefield after re-flection, so the equations shown above do not provide a full description of scattering in the subsurface. In this thesis, however, we will limit our-selves to the vertically-oriented scattering equations 2.6a,b. Note that Davydenko and Verschuur [2013b] describe a methodology to extend FWM to correctly handle steeply dipping reflectors.

Also notice that we do not supply a deterministic link between the medium parameters in the subsurface and the corresponding reflection operators. In this work we merely observe that reflections will occur at discontinuities in the subsurface and we will estimate the correspond-ing operators directly. Nevertheless, a derivation for acoustic reflection operators in a 1.5D medium is given in appendix A.

2.1.3 Source waveﬁelds

Wavefields in the subsurface are generated by a source term as illus-trated in figure 2.3. When a physical source is fired at depth level zn,

we simply extend equations 2.5a,b by adding a source wavefield �S±_(z

n)

to the scattered wavefields leaving that depth level: �

Q+_(z

n) = �S+(zn) + T+(zn) �P+(zn) + R∩(zn) �P−(zn) (2.8a)

�

Q−(zn) = �S−(zn) + T−(zn) �P−(zn) + R∪(zn) �P+(zn). (2.8b)

Notice that there should be no propagation and scattering effects in-cluded in the source wavefields �S, as these effects are already taken into

(24)

 

n S z

 

n S z n z 1 n z  1 n z 

Figure 2.3: Physical sources at depth level zn generate wavefields �S+/−(zn), moving away

from that depth level.

account. The source wavefields at depth level zn then are those

wave-fields that would be leaving depth level znif there was a physical source

exclusively at that depth level and there were no reflectors anywhere else in the entire medium, although a reflector at depth level zn itself is

allowed. In this thesis, we further specify that the source wavefield is given as:

S±(ω, x, z) = q(ω)√

jωδ(x − xsrc, z − zsrc) (2.9) where q(ω) is the source wavelet and δ(x − xsrc, z − zsrc) is the delta

function with (xsrc, zsrc) being the source location. This definition

en-sures that the data has the same phase and amplitude behaviour as the source wavelet. Furthermore, the source term has a radiation pattern compatible with a dipole source (oriented along the z-axis).

2.2 Waveﬁeld modelling

2.2.1 Modelling process

Full Wavefield Modelling is performed recursively in depth. For the modelling procedure, we assume that the downgoing wavefield �P+_(z

0)

2.2 Waveﬁeld modelling 21

and the upgoing wavefield �P−_(z

N) just outside our modelling grid are

known (usually we assume they are 0) and we also assume that the source wavefield as well as the reflectivity, transmission and propagation operators are known. The modelling process then consists evaluating consecutive roundtrips, each roundtrip being composed of a downward depth recursion followed by an upward depth recursion as follows:

[1] For the downgoing wave (n=0,1,2,...,N-1 )

• Apply boundary conditions and add source wavefields at depth level zn:

�

Q+_i (zn) = �S+(zn) + T+(zn) �Pi+(zn) + R∩(zn) �P_i−1− (zn).

(2.10a) • Extrapolate the wavefield to depth level zn+1:

�

P_i+(zn+1) = W+(zn+1; zn) �Q+i (zn). (2.10b)

[2] For the upgoing wavefield (m=N,N-1,N-2,...,1 ).

�

Q−_i (zn) = �S−(zn) + T−(zn) �Pi−(zn) + R∪(zn) �Pi+(zn). (2.11a)

• Extrapolate the wavefield to depth level zn−1:

�

P_i−(z_n−1) = W−(z_n−1; zn) �Q−i (zn). (2.11b)

In all above expressions, subscript i indicates the number of roundtrips that are included in the wavefields. Notice that during the ith_downward

recursion, we use the upgoing wavefield �P_i−1− at each depth level. These wavefields are not known `a priori, rather they are updated iteratively in consecutive depth recursions, starting from an initial value of �P−

0 (zn) =

(25)

 

n S z

 

n S z n z 1 n z  1 n z 

Figure 2.3: Physical sources at depth level zn generate wavefields �S+/−(zn), moving away

from that depth level.

account. The source wavefields at depth level zn then are those

wave-fields that would be leaving depth level znif there was a physical source

exclusively at that depth level and there were no reflectors anywhere else in the entire medium, although a reflector at depth level zn itself is

allowed. In this thesis, we further specify that the source wavefield is given as:

S±(ω, x, z) = √q(ω)

jωδ(x − xsrc, z − zsrc) (2.9) where q(ω) is the source wavelet and δ(x − xsrc, z − zsrc) is the delta

function with (xsrc, zsrc) being the source location. This definition

en-sures that the data has the same phase and amplitude behaviour as the source wavelet. Furthermore, the source term has a radiation pattern compatible with a dipole source (oriented along the z-axis).

2.2 Waveﬁeld modelling

2.2.1 Modelling process

Full Wavefield Modelling is performed recursively in depth. For the modelling procedure, we assume that the downgoing wavefield �P+_(z

0)

2.2 Waveﬁeld modelling 21

and the upgoing wavefield �P−_(z

N) just outside our modelling grid are

known (usually we assume they are 0) and we also assume that the source wavefield as well as the reflectivity, transmission and propagation operators are known. The modelling process then consists evaluating consecutive roundtrips, each roundtrip being composed of a downward depth recursion followed by an upward depth recursion as follows:

[1] For the downgoing wave (n=0,1,2,...,N-1 )

�

Q+_i (zn) = �S+(zn) + T+(zn) �Pi+(zn) + R∩(zn) �P_i−1− (zn).

(2.10a) • Extrapolate the wavefield to depth level zn+1:

�

P_i+(zn+1) = W+(zn+1; zn) �Q+i (zn). (2.10b)

[2] For the upgoing wavefield (m=N,N-1,N-2,...,1 ).

�

Q−_i (zn) = �S−(zn) + T−(zn) �Pi−(zn) + R∪(zn) �Pi+(zn). (2.11a)

• Extrapolate the wavefield to depth level zn−1:

�

P_i−(z_n−1) = W−(z_n−1; zn) �Q−i (zn). (2.11b)

In all above expressions, subscript i indicates the number of roundtrips that are included in the wavefields. Notice that during the ith_downward

recursion, we use the upgoing wavefield �P_i−1− at each depth level. These wavefields are not known `a priori, rather they are updated iteratively in consecutive depth recursions, starting from an initial value of �P−

0 (zn) =

(26)

2.2.2 Modelling equations

The operations described above can be summarised in modelling equa-tions as follows: � P_i+(zn) = n−1 m=0 U+(zn; zm) �S+(zm) + R∩(zm) �P_i−1− (zm) (2.12a) � P_i−(zn) = N m=n+1 U−(zn; zm) �S−(zm) + R∪(zm) �Pi+(zm) , (2.12b) where: U+(zn; zm) = _m+1 k=n−1 W+(zk+1; zk)T+(zk) W+(zm+1; zm) (2.12c) U−(zn; zm) = _m−1 k=n+1 W−(z_k−1; zk)T−(zk) W−(z_m−1; zm). (2.12d)

Note that these equations are equivalent to the modelling equations described by Berkhout [2014a], where all the transmission and reflection operators are implicitly included in secondary source terms. Also note that the equations bear strong similarities to the (generalised) Bremmer series [Bremmer, 1951; Corones, 1975; de Hoop, 1996]. We choose this formulation because it explicitly shows the complexity of the overall operators U±_(z

n; zm) used during the modelling procedure, knowing

that the same operators will be used to calculate the gradients of the inversion scheme.

2.2.3 Surface seismic data

Using the general form of equations 2.12a,b, we can now describe sur-face seismic data. In this situation, we have a source distribution that is limited to a downgoing source wavefield �S+ _{at the surface z}

0. Figure

2.5 shows an example of the modelled upgoing wavefield at the surface �

P_i−(z0), based on the model shown in figure 2.4, for different roundtrips

i. Notice that the first roundtrip will generate all primary reflection events (including transmission effects), the second roundtrip will recal-culate the primaries and will use the existing wavefields to generate

2.2 Waveﬁeld modelling 23 velocity [m/s] 2000 2200 2400 2600 depth [m] 0 50 100 150 200 250 300 350 400 450 500 velocity amplitude -0.5 0 0.5 depth [m] 50 100 150 200 250 300 350 400 450 500 reflectivity a) b)

Figure 2.4:1.5D model used to generate data in figure 2.5. a) Velocity profile. b) Reflectivity

profile.

the first-order multiples, the third roundtrip again recalculates the pri-maries and will now generate first- and second-order multiples based on the existing wavefields, and so on. In theory then, each roundtrip will increase the order of scattering included in the wavefields by one, so we have to evaluate an infinite amount of roundtrips to model all scatter-ing effects of the wavefield in the subsurface. However, in practice the amplitudes of higher-order multiples often decay so rapidly that we do not need to model them and we can limit ourselves to a small number of roundtrips.