Computational methods for seismic imaging and monitoring

(1)

Dissertations in Forestry and Natural Sciences

KENNETH MUHUMUZA

Computational methods for seismic imaging and monitoring

PUBLICATIONS OF

THE UNIVERSITY OF EASTERN FINLAND

(2)

(3)

PUBLICATIONS OF THE UNIVERSITY OF EASTERN FINLAND DISSERTATIONS IN FORESTRY AND NATURAL SCIENCES

N:o 393

Kenneth Muhumuza

COMPUTATIONAL METHODS FOR SEISMIC IMAGING AND MONITORING

ACADEMIC DISSERTATION

To be presented by the permission of the Faculty of Science and Forestry for public examination in the Auditorium MD100 in Mediteknia Building at the University of Eastern Finland, Kuopio, on October 28th, 2020, at 12 o’clock.

University of Eastern Finland Department of Applied Physics

(4)

Grano Oy Jyväskylä, 2020

Editors: Pertti Pasanen, Jukka Tuomela, Matti Tedre, and Raine Kortet

Distribution:

University of Eastern Finland Library / Sales of publications julkaisumyynti@uef.fi

http://www.uef.fi/kirjasto

ISBN: 978-952-61-3578-6 (print) ISSNL: 1798-5668

ISSN: 1798-5668 ISBN: 978-952-61-3579-3 (pdf)

ISSNL: 1798-5668 ISSN: 1798-5676

ii

(5)

(6)

Author’s address: University of Eastern Finland Department of Applied Physics P.O.Box 1627

70211 KUOPIO FINLAND

email: kenneth.muhumuza@uef.fi

Supervisors: Docent Timo Lähivaara

University of Eastern Finland Department of Applied Physics P.O.Box 1627

email: timo.lahivaara@uef.fi Dr. Teemu Luostari

University of Eastern Finland Department of Applied Physics P.O.Box 1627

email: teemu.luostari@outlook.com Professor Jari Kaipio

University of Auckland Department of Mathematics 38 Princess Street

AUCKLAND 1010 NEW ZEALAND

email: jari@math.auckland.ac.nz Reviewers: Professor Elena Kozlovskaya

University of Oulu Oulu Mining School P.O.Box 8000 90570 OULU FINLAND

email: elena.kozlovskaya@oulu.fi Professor Kristopher A. Innanen University of Calgary

Department of Geoscience 2500 University Drive NW CALGARY, AB T2N 1N4 CANADA

email: k.innanen@ucalgary.ca

Opponent: Professor Børge Arntsen

Norwegian University of Science and Technology Department of Geoscience and Petroleum

S.P. Andersen’s vei 15a N-7491 TRONDHEIM NORWAY

email: borge.arntsen@ntnu.no iv

(7)

Kenneth Muhumuza

Computational methods for seismic imaging and monitoring Kuopio: University of Eastern Finland, 2020

Publications of the University of Eastern Finland Dissertations in Forestry and Natural Sciences 2020; 393

ABSTRACT

Seismic imaging methods as tools for geophysical exploration provide excellent op- portunities for delineating spatial variations in subsurface properties and monitoring their temporal changes using active and passive-source seismic data. Traditional seismic imaging techniques such as travel-time tomography and ray-based migration do not properly correct for the conversions and bending the seismic waves can undergo while travelling through complex subsurface structures. Waveform tomography, which aims to fit the entire signal waveform during the inversion process, improves the quality of the image reconstruction of the subsurface properties. In addition, the use of Bayesian-based approaches enable the incorporation of uncertainties, which is useful for modelling errors and assessing the reliability of the reconstructions.

In this thesis, we explore geophysical subsurface imaging and monitoring using active and passive seismic methods. We seek to enhance the estimates of the subsurface properties with computationally efficient waveform-based methods and monitor temporal changes in these properties. By utilising a Bayesian-based approach, we explore uncertainty modelling in seismic imaging and propose an improvement on the Born-based approximative approaches for acoustic frequency-domain waveform inversion. Additionally, we apply a noise-based seismic approach that uses passive-source seismic data for imaging of subsurface features and monitoring of changes in subsurface physical properties. The computational implementations of the different imaging and monitoring methods are discussed and evaluated through synthetic tests of active-source seismic experiments and field measurements of seismic ambient noise.

In the first part of this thesis, we consider the problem of time-lapse seismic monitoring of CO₂injection using full-waveform inversion. The objective here is to reconstruct the difference image of the subsurface between repeated active-source surveys acquired at different times to monitor the CO₂ plume distribution. Be- cause the problem is computationally challenging, we explore efficient modelling approximations by employing the distorted Born approximation, which is based on scattering integral equation.

In the second part of the thesis, we apply a Bayesian-based approach to the problem of reconstructing the subsurface properties from active-source seismic data. The Bayesian framework provides a methodology for the formal incorporation of prior information and parameter errors into the inversion, thereby providing a measure for parameter uncertainty. By considering two choices for the prior, we study the potential for improving acoustic Born waveform inversion of seismic data for viscoelastic media using a Bayesian approximation error approach.

In the third part of the thesis, we use passive-source seismic data (ambient seismic noise) to infer the crustal structure of the Bransfield Strait. With this methodology, cross-correlating the ambient seismic noise traces results to waveform data that

(8)

are identical to actively collected reflection seismic data. The aim is to obtain information by inverting surface-wave dispersion curves from seismic ambient noise, to probe the crustal structure.

Universal Decimal Classification:53.05, 550.8, 550.344.2, 004.02

LCSH: Seismic tomography, Inversion (Geophysics), Seismology — Computer programs, Imaging systems in geophysics, Seismic waves, Seismic reflection method — Data processing, Bayesian field theory

Keywords:waveform inversion, computational seismology, inverse theory, Bayesian framework, seismic ambient noise, time-lapse seismic, waveform tomography, geophysical subsurface imaging, seismic data

vi

(9)

ACKNOWLEDGEMENTS

There is a less obvious debt to the Department of Applied Physics at the University of Eastern Finland where, under the computational Physics and Inverse Problems research group, I have broadened my understanding of inverse theory. With several academics working on different aspects of inverse problems, I have been uniquely placed to undertake my doctoral research.

I am indebted to my main supervisor Docent Timo Lähivaara for his devoted help and professional advice that enabled me to learn a lot throughout every phase of my doctoral research. I appreciate his vast knowledge and skills in the areas of geophysical inversion. His great efforts to explain things clearly and simply helped to make geophysical inverse problems fun for me. In addition, I extend my gratitude to my other supervisors Professor Jari Kaipio and Dr. Teemu Luostari for their positive and constructive guidance during my doctoral research.

I am grateful to my co-authors Professor Morten Jakobsen, Associate Professor Lassi Roininen, and Dr. Janne M. J. Huttunen. This thesis would not have been possible without their contributions to the publications and the interesting discussions regarding the thesis work. In particular, Professor Morten Jakobsen introduced me to the world of scattering theory during my master’s studies and inspired my work on non-linear inversion.

Special thanks goes to the official pre-examiners Professor Elena Kozlovskaya and Professor Kristopher A. Innanen for taking the time and effort necessary to review and assess my thesis. I sincerely appreciate all your valuable comments and input to improve this thesis.

Great thanks are due to my family members for the moral support, and the doctoral school of the University of Eastern Finland and the Academy of Finland for the financial support during the entire period of my doctoral studies.

Lastly, I would like to thank my parents for the kind of upbringing they have provided. My mum taught me to never give up on my endeavors and to always believe that my future could be as bright as my faith. This mindset instilled in me is so impactful that I don’t limit my challenges but challenge my limits.

Kuopio, September 8, 2020 Kenneth Muhumuza

(10)

(11)

LIST OF PUBLICATIONS

This thesis consists of the present review of the author’s work in the field of seismic imaging and the following selection of the author’s publications:

I K. Muhumuza, M. Jakobsen, T. Luostari, and T. Lähivaara, "Seismic monitoring of CO₂ injection using a distorted Born T-matrix approach in acoustic approximation,"Journal of Seismic Exploration27, 403–431 (2018).

II K. Muhumuza, L. Roininen, J. M. J. Huttunen, and T. Lähivaara, "A Bayesian- based approach to improving acoustic Born waveform inversion of seismic data for viscoelastic media,"Inverse Problems,36(2020) 075010 (19pp).

III K. Muhumuza, "A feasibility study on monitoring crustal structure variations by direct comparison of surface wave dispersion curves from ambient seismic noise,"International Journal of Geophysics2020(2020). Article ID 5269537 (11pp).

Throughout the overview, these papers will be referred to by Roman numerals.

AUTHOR’S CONTRIBUTION

The research presented in this thesis develops computational methods for seismic imaging and monitoring, applying them to both active and passive seismic data.

The publications contained herein are results of joint work with the thesis author, the supervisors, and the co-authors. The thesis author has been the principal (and corresponding) author of all publications, but the guidance and support of the supervisors and co-authors has been paramount. The author has been responsible for the results of the numerical experiments and implemented the numerical algorithms for the forward and inverse problems with the help of co-authors. The distorted Born T-matrix non-linear inversion method employed in PublicationIwas inspired by Professor Morten Jakobsen. The Cauchy prior for edge-preserving Bayesian inversion used in PublicationIIwas developed by co-author Lassi Roininen. The Com- putational Seismology Group at the Department of Mathematics and Geosciences, University of Trieste provided the data and frequency-time analysis (FTAN) software used in PublicationIII.

(12)

(13)

1 Introduction

Seismic imaging—which uses seismic waves to probe the earth structure—is an inverse problem in the sense that we aim to obtain a subsurface model that explains a set of observational data. Seismic waves are energy-transmitting vibrations that propagate through the earth from point to point away from a source (Figure 1.1).

The source may be either man-made, such as mechanical vibrators and explosives, or natural, such as earthquakes and ambient noise [1]. There are several types of seismic waves, and these can generally be categorised as body waves, which propagate through the earth’s interior, or surface waves that propagate along the earth’s surface. Body waves generated by the source propagate in all directions and are mainly of two types: compressional waves, known as (Primary) P-waves, and shear waves, known as (Secondary) S-waves. Surface waves travel slower than body waves and are mainly of two forms: Rayleigh waves and Love waves.

x

z

Seismic source

Seismic raypaths

Geophones

Figure 1.1:Sketch of a typical seismic exploration experiment showing a source, an array of geophones (olive triangles), and ray-traced paths: The layer cake model con- sist of three different subsurface layers having diffent material properties. Seismic waves (black arrows) from a seismic source (a vibrator) are transmitted or reflected at interfaces between the subsurface layers with different material properties and recorded by geophones.

Seismic waves carry with them information about the physical properties of the subsurface media they travel through. This information can be recorded by seismic sensors, such as geophones and seismometers on land, or hydrophones and ocean bottom seismometers onshore. The recorded seismic data can be used to delineate spatial variations in subsurface earth properties for a broad range of applications ranging from exploration of oil, groundwater, and other resources to urban applications. There are several different methods which can be employed to create subsurface structural images across a range of scales. Such methods can utilise active- source seismic data, in which the illumination is provided by man-made sources, or passive-source seismic data, in which the illumination is made by natural uncon-

(18)

trolled sources.

Monitoring subsurface properties by performing time-lapse experiments or using continuous seismic measurements is also crucial for assessing time-varying pro- cesses in the subsurface. For example, monitoring the extent of the injected CO₂ plume is essential for the characterization of the reservoir as well as the detection of potential leakage. The use of active and passive seismic methods has proven to be useful in monitoring changes in the subsurface [2–4]. Passive seismic monitoring with ambient noise has recently received considerable attention as an alternative to active-source seismic monitoring. This is because the method that does not require an active source to image or monitor the subsurface, making it low-cost and envi- ronmentally friendly. It has been shown that cross-correlating the traces of a passive experiment can generate data that are identical to actively collected seismic data [5].

The main methods commonly used in the geophysical community for seismic imaging include seismic migration, least squares migration, full-waveform inversion (FWI), phase-like inversion, and migration velocity analysis [6]. In this thesis, we consider waveform-based imaging methods to obtain quantitative estimates of subsurface properties. The subsurface properties can be represented by the spatial distribution of, e.g., P-wave velocity, S-wave velocity, porosity, and density. Full- waveform inversion has the potential to be used to reconstruct multiple parameters.

However, it is more often used to reconstruct P-wave velocity models only due to several challenges associated multi-parameter FWI, such as high computational burden and parameter cross-talk problem [7]. As with all data-fitting procedures, FWI uses a forward model to predict the seismic data (provided the model parameters are known) as part of the inversion algorithm.

At the core of any successful waveform inversion is a forward model that accurately represents the realities of seismic wave propagation in the earth. However, accurate solvers are not often applied due to the high computational cost of solving a multidimensional minimization problem. In most cases, fast and less accurate (approximate) forward models are used in the inversion procedure. For example, seismic forward modelling is typically performed using the acoustic wave equation, and approximations such as Born and distorted Born modelling [8, 9] underlie most seismic imaging techniques. These assumptions may result in acceptable recon- structed models of subsurface properties under certain circumstances, but they can lead to inaccurate models in others.

Waveform inversion can be performed in the time-space domain [10–12] or in the frequency-space domain [13, 14]. The frequency domain is considered less computationally expensive, especially when relatively few frequencies are used to obtain inversion results. This can be achieved by strategically selecting the frequencies used in the inversion [15,16]; finding a compromise between accuracy and computational cost.

In solving the waveform inversion problem of estimating subsurface properties from observed seismic data, we have applied deterministic and Bayesian approaches. Deterministic approaches to seismic inversion [17] provide a single “best”

solution to the problem, without giving explicitly any probabilistic model of uncer- tain parameters. In Bayesian approaches [18, 19], all parameters are thought of as random variables with a certain probabilistic structure. This provides, in principle, a methodology for incorporating parameter errors into the inversion, giving many equally plausible solutions with a measure of uncertainty. There are several sources which may contribute to parameter uncertainty including measurement errors, and model uncertainty arising from the lack of sufficiently accurate knowledge of struc- 2

(19)

tural model parameters.

Current research efforts within waveform inversion and monitoring strive to reduce the computational cost of the seismic inverse problem by developing efficient approximative methods. Along the same line, several authors have suggested strategies to improve acoustic FWI of seismic data for elastic or viscoelastic media for better reconstruction of P-wave velocity models [20–22]. In this thesis, the presented acoustic frequency-domain waveform methods in Publications I and II are based on computationally efficient solvers: distorted Born and Born approximations. These imaging and monitoring methods are discussed through synthetic tests of active-source seismic experiments. To address elastic and viscous effects in acoustic FWI and partially dispense with the assumptions of using an approximate forward model, we have adopted the Bayesian approximation error (BAE) approach [19] in PublicationII.

In Publication III, we applied a passive seismic monitoring method based on detecting temporal changes of velocities by computing cross-correlations of ambient seismic noise. By extracting Rayleigh surface waves from noise cross-correlations, we investigated the feasibility of the direct use of Rayleigh-wave group dispersion curves to monitor temporal variations of crustal seismic velocities. We studied the continuous ambient noise data recorded by a seismic array in West Antarctica and explored in detail the cross-correlation results for a path that crosses the Bransfield Strait. By solving an inverse problem, the retrieved dispersion curves were also used to estimate a velocity-depth profile along the studied path, gaining insight into the crustal structure of the Bransfield Basin.

The rest of the thesis is structured as follows: Chapter 2 describes wave propagation in isotropic elastic, acoustic, and viscous media, and reviews the numerical methods used for the solution of the forward problem. In Chapter 3, deterministic and Bayesian inversion frameworks are discussed, and the utilisation of the BAE approach in modelling of errors is described. Chapter 4 discusses active and passive seismic methods for imaging the subsurface and reviews the results of the publications. Finally, a summary of the various conclusions of this thesis publications is given in Chapter 5.

(20)

(21)

2 Modelling of seismic wave propagation

A thorough understanding of seismic wave propagation is essential for exploration seismology. It provides a way of understanding the behaviour of seismic waves and fundamentally facilitates an apprehension of the features observed in recorded seismic data. Modelling of seismic wavefields is one of the cornerstones of seismic data processing and imaging techniques. In addition to guiding seismic data interpretation, seismic modelling can also aid in the design and evaluation of seismic surveys in order to choose optimum acquisition methods.

It’s important to note that the exact modelling of wave propagation in the earth is often infeasible due to its extreme complexity. But if we consider various assumptions and approximations, there exist several analytical and numerical techniques for modelling seismic wave propagation. Analytical methods can only handle simple canonical models as they often fall short at explaining complex wave phenomena such as converted waves, scattering, and generation of multiples in complex heterogeneous earth models. Numerical methods are, therefore, the major tools for modelling wave propagation in different rheologies: e.g. acoustic, viscoacoustic, elastic, viscoelastic, poroelastic and poro-viscoelastic.

In this chapter, we begin by reviewing the governing equations of seismic wave propagation in linear-elastic, acoustic and viscous media, and then describe the numerical methods we used for the computation of seismic wavefields.

2.1 GOVERNING EQUATIONS

The wave equation is the starting point for describing the propagation of seismic waves in a medium. It can be formulated mathematically as a set of (partial) differential equations or an integral equation–both of which aim to describe the evo- lution of wavefields. The conservation of momentum—which follows directly from the Newton’s laws of motion—provides the basis for the formulation of the wave equation. If the medium through which seismic waves propagate is treated as a continuum, the wave equation is expressed by a partial differential equation [1, 23]

∂_jσ_ij+ f_i=ρ∂²u_i

∂t² , (2.1)

where the∂_jσ_ijterm is the divergence of the stress tensorσ_ij(i,j=1, 2, 3 respectively representingx,y,zon the Cartesian plane, and we adopt the summation convention for the repeated indices), f_i is the body force vector per unit volume,ρ is the mass density, u_i is the displacement vector and tis time. Each of the terms, σ_ij, f_i and u_i is a function of position x = (x,y,z)and timet. The stress tensor describes the internal forces that neighbouring particles in a deformable material exert on each other. The conservation of angular momentum implies that the stress at a given point in a continuum is characterized by the symmetric stress tensor, i.e,σ_ij=σ_ji.

Let us next introduce the concept of the strain tensor ε—which describes the deformation due to an applied stress σ and intrinsically involve spatio-temporal variations of the displacement field u = (ux,uy,uz). Under small deformations,

(22)

such as in most seismic experiments, the strain tensor is expressed as

ε_ij =∂_iu_j+∂_ju_i. (2.2) For a general anisotropic elastic medium—for which the material properties are direction-dependent—the relation between stress and strain tensors is given by the Hooke’s law of linear elasticity [1]

σ_ij=_C_ijklε_kl, i,j,k,l=_x,_y,_z, _(2.3) where the 3⁴ = 81 components of the elastic stiffness tensor Care the elastic constants that describe elastic properties of the medium. As a consequence of the symmetry of the stress and strain tensors, the elastic stiffness tensor satisfies the minor symmetryC_ijkl=C_jikl =C_ijlkand the major symmetryC_ijkl=C_klij. Using the minor symmetry, the number of independent components ofCreduce from 81 to 36 [1,24].

The major symmetry condition further reduces the number of independent components to 21 [1,23]. This is the maximum number of independent elastic constants for any medium. The number of elastic constants that describe a material is reduced as the medium symmetry increases. For example, transversely isotropic materials that have one axis of symmetry (5 elastic constants), and isotropic material that have an infinite number of planes of symmetry (only 2 elastic constants).

To simplify the notation, from now on, we will drop the summation convention for repeated indices. Depending on the type of anisotropy considered for the earth model, substitution of the constitutive law (2.3) into Equation (2.1) leads to the set of equations describing wave propagation in elastic media. These equations can be transformed into the frequency domain, which is obtained by performing Fourier transform. Thus, the time-domain wave equations (2.1) and (2.3) can be Fourier transformed into the following frequency-domain elastic wave equations

∇ ·σ(x,ω) +f(x,ω) =−ω²ρ(x)u(x,ω) (2.4) σ(x,ω) =_C(_x)ε(x,ω), (2.5) where ω is the angular frequency and ∇ is the divergence operator. The stress tensorσ, displacementu, body forcef, and the strain tensorεare now functions of positionxand frequencyω.

A more detailed discussion of the theory of seismic wave propagation in various idealized earth models can be found in the literature [25–27].

2.1.1 Isotropic elastic and acoustic (lossless) media

In the following, we assume the material is isotropic and has constant density. The assumption of isotropy implies that material properties at a given point are identical in all directions. For isotropic media, the stiffness tensorC(x)depends upon only two independent elastic constants, the Lamé parameters,λ(_x)_andµ(_x). In this case, the stress tensor in Equation (2.5) can be represented in the form

σ(x,ω) =λ(∇ ·u)I+2µε(x,ω), (2.6) where I is the identity matrix. Isotropic elastic media can support the propagation of both P and S body waves. The propagation velocities of these waves are determined by the density and elasticity of the material through which they travel. Hence, 6

(23)

isotropic P-wave and S-wave velocities Vp and Vs are related to the elastic Lamé parametersλandµthrough

λ+2µ=ρV_p², (2.7)

µ=ρV_s². (2.8)

Acoustic-media approximations to the elastic wave equation are often considered in exploration geophysics to make the models of wave propagation more tractable.

Acoustic media supports only P-waves, implying that the descriptive equation can be derived as a special case of Equation (2.4) by assuming that the S-wave velocity Vs is zero [28, 29]. In such a case, only the diagonal elements of the stress tensor σ are non-zero and equal to the negative pressure−p. Hence, in isotropic acoustic media, Equation (2.4) becomes

−ω²ρ(x)u(x,ω) =∇p(x,ω) +f(x,ω). (2.9) Substitutingµ=0 into (2.6) and (2.7) yield, respectively, the constitutive relations:

∇ ·u=−p/λ,

Vp=^p_λ/ρ. ^(2.10)

Thus, taking the divergence of Equation (2.9) and substituting relations (2.10) leads to the Helmholtz equation,

−∇²p(x,ω)−k²(x)p(x,ω) = f(x,ω), (2.11) where f =−∇ ·fis the dipole source term, and the wavenumberk(x)of a lossless medium is related toωby the usual formula

k(x) = ^ω

Vp(x)^. ^(2.12)

The acoustic approximation greatly simplifies the wave equation in elastic media. Nonetheless, Equation (2.11) can have the same solution as Equation (2.1) only if both of the following two conditions hold: (1) an infinite homogeneous isotropic medium; and (2) sources that generate only P-waves. But since the earth is not homogeneous, velocity contrasts between layers of rock within the earth lead to reflections at each interface. The generation of S-wave reflections—due to an incident P-wave striking the interface—breaks the assumption that the acoustic wave equation should not yield S-waves. In addition, the acoustic approximation does not account for P-S (and S-P) mode-converted reflections, which implies that the ampli- tudes for the elastic P-waves are incorrectly predicted by the acoustic model. Among the estimation frameworks presented in PublicationII, we adopt the Bayesian approximation error approach (discussed in Section 3.4), which takes into account the errors and uncertainties related to using an acoustic approximation to elastic waves.

2.1.2 Wave propagation in viscous media

The subsurface behaves like viscous media—as seismic waves propagate through the earth, they experience energy loss due to several dissipation mechanisms [1].

The anelastic attenuation of wave energy due to the viscosity of earth rocks af- fects the signal amplitude and phase of propagating waveform. This attenuation

(24)

phenomenon needs to be incorporated into the wave equation for more realistic modelling of seismic wave propagation. A common practice is to quantify seismic attenuation by the quality factor Q (or its inverse Q⁻¹, known as the dissipation factor), which is often defined as the ratio of maximum energy stored to the energy dissipated during a stress cycle. The assumption of a constant Q, that is Qinde- pendent of frequency, is usually used because it holds in the typical frequency band for exploration geophysics, and it considerably simplifies the analysis of attenuation [30, 31].

For linear viscous media, the stress depends linearly on the time history of the strain, such that the stress-strain relation is expressed as a convolution [27]:

σ=ψ∗^dε dt = ^dψ

dt ∗ε, (2.13)

whereψ is the relaxation tensor function and the symbol∗denotes time convolution. In practice, relaxation functions based on constant-Qmechanical models are commonly used, and time convolutions are avoided by introducing the so-called memory variables and fractional derivatives [32, 33]. The mechanical models, such as Kelvin-Voigt model, Maxwell model, and standard linear solid model (SLS) usually contain weightless springs, which store strain energy, and dashpots, which dis- sipates energy [27]. By taking into account these classical mechanical models, the corresponding viscoacoustic and viscoelastic wave equations can be derived both in time and frequency domain. [27, 32].

In Publication II that involves the study of viscoelastic media, we adopt the Kelvin-Voigt model, which is a linear model and well justified in modeling the viscoelastic behavior of solid rocks in the typical seismic band [34–37]. Kelvin-Voigt viscoelasticity is introduced to build a viscoelastic medium based on an isotropic elastic material. In the context of Kelvin-Voigt model, the quality factor Q that char- acterizes energy dissipation in a material is given byQ(ω) = (ωτ)⁻¹, where τ is the relaxation time —which is defined as the ratio of the viscosity of the dashpot to the elasticity constant of the spring. A more detailed explanation of wave propagation in viscoelastic media and the implementation of the Kelvin-Voigt model of viscoelasticity in the frequency domain can be found in the literature [27, 34, 38].

2.1.3 Boundary conditions

In order to accurately simulate seismic wavefields, undesired reflections from the boundaries of the computational domain should be minimised by implementing two typical boundary conditions. First, the free surface condition to represent the interface between the ground and the atmosphere, where there is the strongest dis- continuity for the material properties. On the free surface, the boundary condition is that the traction is zero [24].

Second, the absorbing boundary condition that mimic an infinite medium by absorbing any spurious reflections from the sides of the computational domain. If these artificial reflections are not taken care of adequately, they reflect back into the model domain and cause numerical errors. Several methods, such as parax- ial approximations [39], Higdon artificial boundary [40], sponge boundary condition [41, 42], and the perfectly matched layer (PML) [43, 44] have been proposed to suppress the boundary reflections. Currently, PML has become the method of choice because it effectively removes unwanted numerical reflections from the model boundaries, under any incidence angle and frequency.

8

(25)

2.2 THE SEISMIC FORWARD PROBLEM

The forward problem refers to the computation of seismic wavefields when one knows the model of the subsurface, and assuming that all input parameters are known. Thus, we solve numerically the wave equation—to determine the seismic datad—given the subsurface model parametersm, the source function, the source- receiver configuration, and the related initial-boundary conditions. The mathemati- cal formulation of the seismic forward problem can be written as [18, 45]

d=G(m), (2.14)

where G is a linear or non-linear operator that describes the explicit relationship between the seismic data and the model parameters. For example, in the acoustic case, G represents the acoustic wave equation based algorithm that simulates the pressure wavefield everywhere in the subsurface model given the P-wave velocity, density, the source and the source–receiver configuration.

Solutions of the forward problem (2.14) for realistic subsurface models is at the core of full-waveform inversion. Such solutions can be achieved by several numerical methods that can be categorised into three groups: (i) direct methods, such as finite difference and finite element [46, 47]; (ii) integral-equation methods, such as volume and boundary integral methods [48–50]; (iii) asymptotic methods, such as ray tracing methods [51]. For a given source position, the solutions to the wave equation give the predicted ground motions (displacement, velocity, and acceler- ation traces) at specified receiver locations and are often referred to as synthetic seismograms. The generation of synthetic data is an integral part of inversion algorithms and forms the foundation for benchmarking and testing new seismic processing algorithms [52–54]. Below, a brief discussion of three methods (finite element, Born approximation, and T-matrix approach)—which have been used in this thesis to solve the seismic forward problem—is presented.

2.2.1 Finite element method

The finite element method (FEM) is increasingly becoming attractive for geophysical applications that require solving wave propagation problems. The popularity of this method stems in part from its efficiency in handling complex topographies and geological interfaces by discretizing the earth model using simple geometries called elements [47, 55, 56]. The wavefield inside each element is then approximated by interpolation functions, such as polynomials of certain degrees. In the FEM con- structions, certain points (known as nodes or nodal points) are assigned to elements, and the the term “finite element mesh” is used to refer to the collection of elements and nodal points making up the computational domain of the model. The FEM offers great freedom in the selection of discretization, in terms of both the type of elements and basis functions, and it naturally satisfies the free-surface boundary condition.

Several methods that mimic FEM have been developed and successfully applied to seismic wave simulation problems. Some of the more popular ones include:

Discontinuous Galerkin methods [57–59]; Spectral element methods [60, 61]; and mixed finite element methods [62, 63, 63]. In this thesis, we use the FEM-based software, COMSOL Multiphysics^R 5.3a to numerically solve the forward problem for acoustic and viscoelastic media in the frequency domain.

(26)

2.2.2 T-matrix approach

The T-matrix (or transition operator) approach is an integral equation method based on the multiple scattering theory. It can be used to obtain approximate or exact solutions of the wave equation. We shall consider the scalar Helmholtz equation (2.11), and in the context of the scattering theory, the total wavefield p(x) in the heterogeneous earth media is treated as a reference wavefieldp0(_x)plus a perturbation which is nonlinearly related to the earth media properties. Let us decompose the velocity Vp(_x) into a constant reference velocity V₀(_x) and a perturbation so that we characterise the wavenumber in the actual medium ask²(_x) =k²₀+δL(_x)_, wherek₀ = ω/V₀ and δL(_x)is the perturbation in the squared wavenumber. The Lippmann-Schwinger solution to the Helmholtz equation (2.11) is given by [64, 65]:

p(x) =p0(x) + Z

DG⁽⁰⁾(x,x⁰)δL(x⁰)p(x⁰)dx⁰, (2.15) where D is the scattering domain in which δL(x) is non-zero and the reference medium Green’s function G⁽⁰⁾(x,x⁰) is the wavefield at position x due to a point source at positionx⁰in the reference medium. Analytical formulas for Green’s functions for homogeneous media can be found in the literature [51], and for smoothly varying heterogeneous media, the reference Greens function can be approximated using ray theory [51]. The computation of the reference Green’s function for an arbi- trary heterogeneous medium can be done using methods such the finite-difference [66] and T-matrix approach [14].

We can rewrite the Lippmann–Schwinger equation (2.15) in a more symmetric form that is convenient for discretization [67]:

p(x) =p0(x) + Z

Ddx1 Z

Ddx2G⁽⁰⁾(x,x1)V^˜(x1,x2)p(x2), (2.16) where ˜V(x1,x2) =δL(x1)δ(x1−x2)is the (two-point) scattering potential field. The main task is to compute the values of the wavefield at the different positions (represented by the vector x2) inside the scattering domain due to virtual scattering sources at different positions (represented by the vectorx1).

If we follow the discretization scheme in [67] and arrange the total (and reference) fields at the scattering domain positions into vectorsP(P0), equation (2.16) can be written in the form

P=_P₀+_G⁽⁰⁾_VP, _(2.17)

whereVis a diagonal matrix with the scattering potential field ˜V.

The T-matrix approach seeks to numerically solve for the field inside the dis- cretized scattering domain by introducing the T-matrixTthat contains all nonlinear effects of multiple scattering:

VP=TP0. (2.18)

We can then generate an integral (or matrix) equation for T-matrix from the Lipmann- Schwinger equation, which gives

T=_V+_VG⁽⁰⁾T. (2.19)

Equation (2.19) has the following exact solution [14, 67, 68]:

T=I−VG⁽⁰⁾−1

V. (2.20)

10

(27)

The T-matrix does not, therefore, depend on the source-receiver configuration and contains the whole Born Neumann series [67] when it is not used in approximation form. The independence on source-receiver geometry is useful especially in applications involving many sources and receiver positions. With this method, the full forward numerical solution only requires the knowledge about the scattering potential field and the Green’s function for the reference medium. The exact T-matrix solution can be computationally expensive, hence, approximations based on Born (Neumann) series for the T-matrix are usually employed. However, the exact solution is convenient especially in cases of time-lapse seismics [67] when the models are not too large.

2.2.3 Born approximation

The Born approximation assumes single scattering between the source and the receiver. It results from linearization of Lippmann-Schwinger integral equation (2.15) by approximating the total field under the integral with the reference field [65, 69]:

p(x) =p0(x) + Z

DG⁽⁰⁾(x,x⁰)δL(x⁰)p0(x⁰)dx⁰. (2.21) In the traditional Born approximation, the reference medium is homogeneous, implying that the computation ofG⁽⁰⁾ and p₀ is analytic. It is among the simplest solutions to the wave equation and widely used because it provides a linear approximation for inverse problems. Although its solution is quite inaccurate in general, it gives decent approximations when the difference between the actual and reference medium is not large, i.e., the values of perturbation are small. The validity of this approximation has also been widely explored in the literature [70–72]. As a weak scattering approximation, the Born approximation ignores multiple scattering effects; it is not suitable for large perturbations; and it is only valid in the low frequency regime with respect to the scattering domain. An extension of the Born approximation to nonhomogenous reference media (the so-called Distorted Born approximation), which performs better than the traditional Born approximation [24]

has been used in PublicationI.

(28)

(29)

3 Deterministic and Bayesian inversion frameworks

Numerical inversion methods are used to solve inverse problems—i.e., inferring physical parameters from observation data. There are generally two types of approaches that can be employed to obtain the parameters of interest from observation data: deterministic and Bayesian (probabilistic) inversion. In principle, deterministic frameworks provide a single “best” solution to the parameter estimation problem, but Bayesian methods can provide several plausible solutions that may explain equally well (up to data errors) the observation data.

In this chapter, we recall the formulation and notation for the seismic inverse problem and provide a brief overview of the deterministic and Bayesian approaches to inversion. We conclude the chapter by reviewing the Bayesian approximation error approach, which can be used to handle model uncertainty and modelling errors.

3.1 THE SEISMIC INVERSE PROBLEM

Formally, the inverse problem can be viewed as the inversion of the forward problem (2.14). Hence, the task is to estimate the physical parameters of interestmgiven the measured seismic data d. Since data are noisy (contain errors) due to several possible reasons, such as inaccurate measurement and uncertainty, the inverse problem entails determining the parameters from given noisy data. Here, we consider an additive noise model for measurement errors. Then, we can write the observation model that links the parameters to the noisy data, through a set of geophysical equationsG:

d=G(m) +e, (3.1)

whereerepresents the noise in the data (observation noise).

Typical seismic inverse problems are large-scale, as they involve large amounts of data and thousands to millions of parameters that describe the subsurface properties. While the forward problem is generally well-posed, the inverse problem can be ill-posed in the sense that the estimates are highly sensitive to even small perturbations in the data [45, 73]. Moreover, it is often the case that there is no explicit formula for the inverse mapping from the data space to the parameter space.

Thus, most of the techniques employed to solve inverse problems rely on solving the forward problem with several plausible input models to compare the output and choose the best model, or a distribution of models that is consistent with the data. A further difficulty is that the forward operator may depend on certain unknown parameters; for example, it may depend on a velocity field that is known only partially, or known with a large uncertainty.

Due to the aforementioned difficulties associated with seismic inverse problems, there is a need for efficient computational inversion methods that incorporate prior information about parameters and account for different sources of uncertainties. In addition, such methods should be designed to provide measures of the reliability of the inversion estimates.

(30)

3.2 DETERMINISTIC INVERSION

The deterministic framework seeks a single solution to observation model (3.1) typically using optimization algorithms without giving explicitly any probabilistic model for the errors or for the parameters. The inverse problem can be formulated within a general framework of output error minimization (least-squares optimization):

Findmsuch thatE(_m) =kG(_m)−dk²₂ is minimised, (3.2) where E(_m) is the objective function—the L2 norm of the misfit between the observed data and the forward prediction. Recall that most seismic inversion problems are ill-posed; hence, the solution by least-squares inversion must be stabilised via regularization algorithms. Tikhonov regularization is the most widely used scheme, where the regularized objective function can be written, for example, in the form

E(_m) =kG(_m)−dk²₂+γkR(_m)k²₂, (3.3) whereγis the regularization parameter andRdenotes the regularization operator, such as a smoothing operator or identity matrix.

In principle, the regularization improves the conditioning of the (deterministic) inverse problem by incorporating additional prior information in the inversion process through the regularization termγkR(m)k²₂. Besides Tikhonov, there are several regularization schemes, such as total variation regularization [74] and least absolute shrinkage and selection operator or LASSO [75]. A crucial aspect of the regularization process is deciding on a suitable value of the regularization parameter γ that provides an optimal trade-off between misfit and the conditioning of the reconstruction. A comparison of the most used techniques (generalized cross-validation, discrepancy principle, and L-curve criteria) are presented in [76] with a focus on non-linear inverse problems. In this thesis, we used the L-curve criterion whose implementation is discussed in [76] and PublicationI.

Several optimization algorithms can be used to solve the regularized or non- regularized form of the least squares problem (3.2). For linear systems, the three standard ways to solve the least squares problem include the normal equations, the QR decomposition, and the singular value decomposition (SVD) [77, 78]. For non-linear systems, these methods include steepest gradient [79], nonlinear conju- gate gradient [80], Gauss–Newton and full Newton method [81], and quasi-Newton [82, 83], to mention a few. Typically, the solution is iterative when dealing with non-linear inverse problems. In PublicationI, the deterministic non-linear seismic inversion problem was solved using a distorted Born iterative method, which has demonstrated to be similar to the Gauss-Newton method of optimization [84, 85].

Though deterministic inversion methods can provide reliable parameter estimates, especially when the number of parameters is small, they are susceptible to the numerical disadvantages of the classical optimization algorithms—their convergence is sensitive to the initial guess of the solution, and global convergence is not guaranteed. Additionally, these methods cannot be used to directly quantify the uncertainty associated with the parameter estimates. Nonetheless, deterministic methods are commonly used because they are often computationally cheaper compared to their counterparts, Bayesian methods (explained in the next section).

14

(31)

3.3 BAYESIAN INVERSION

The Bayesian framework uses a probabilistic formulation based on Bayes’ theorem [86], generating several possible model realizations that satisfy the observational data. All unknown parameters and data are modelled as random variables.

In principle, this framework provides a means to incorporating parameter errors into the inversion, giving rigorous probabilistic estimates with a measure of uncertainty. Furthermore, Bayesian inversion methodologies allow for incorporation of much geophysical information with seismic and a more flexible integration of data than deterministic ones. These methodologies, therefore, solve a long-standing challenge in the seismic industry; how to use seismic attributes in a quantitative manner for effective inversion and uncertainty quantification. However, they remain computationally expensive and difficult to be used routinely, or to be easily integrated in a multi-physical inversion procedure involving other geophysical methods. With the development of computer capabilities and sampling algorithms, Bayesian methods are increasingly being used in geophysical inversion problems.

In the Bayesian framework [18, 19], the solution to the inverse problem (3.1) is based on combining information coming from the observed datadand the model parameters of interestmof the examined medium. There are two main tasks here:

constructing the likelihood model (conditional density of the data given the parameters)π(d|m), and determining the prior density (distribution of the parameters in the absence of any data)π(m). The posterior probability density π(m|d)is given by the Bayes’ theorem,

π(m|d) = ^π(d|m)π(m)

R π(d|m)π(m)dm ∝π(d|m)π(m), (3.4) where the denominator R

π(d|m)π(m)dm is a normalizing constant, which can usually be ignored.

The construction of the likelihood (data) model assumes that the measurement process, i.e.,G(m), linking the data to the model parameters is well-known, and the observation noise properties can be measured. In this thesis, we have carried out Bayesian inversion using a linear observation model, i.e.,

d=Gm+e. (3.5)

Assuming that the observation noiseeis mutually independent with the model parametersm, the likelihood model takes the form

π(_d|_m) =πe(_d−_Gm), (3.6) whereπeis the is the probability density of the observation noisee.

Let us denote by θ ∼ N(θ∗,C_θ)the multivariate joint normal distribution with meanθ∗=E(θ)and covarianceC_θ =Cov(θ). If we assume that the noiseeis Gaus- sian distributed with meane∗ and covarianceCe, i.e.,e ∼ N(e∗,Ce), the likelihood model can be written as a Gaussian function

π(d|m)∝exp

−¹

2kLe(d−Gm)−e∗k²

, (3.7)

where Le is computed from the factorization of the inverse of the noise covariance matrixC⁻¹_e = L_e^TLe. Factorization can be done, for example, with Cholesky decomposition. In the next section, we discuss the task of constructing the prior model.

(32)

3.3.1 Prior modelling

Information about the subsurface geology or knowledge of the subsurface conditions—

which can be used to inform the process of parameter estimation—often exists in the absence of seismic data [87]. A prior modelπ(m)is used to describe this preexist- ing information regarding the parameter of interest in terms of a probability density.

Since the importance of prior information is to constrain the parameters to feasible values and how these values are distributed, it plays the role of regularization in Bayesian inversion. Thus, prior parameters influence the posterior distribution, and this influence is even greater for parameters which are not well-defined by the data.

Prior modelling remains a challenge because it can be difficult to translate the prior information (which is often in qualitative form) into a quantitative form [88, 89].

The widely used prior models in many imaging contexts are the Gaussian process priors [90]. The Gaussian prior model assumption of the geophysical parameters coupled with a Gaussian likelihood is commonly used for inversion of seismic data because it yields an analytical expression for the posterior model. However, Gaussian priors often fail to reconstruct sharp interfaces, for instance, in reconstructing the velocity structure. To potentially overcome this problem, non-Gaussian prior models have been developed in Bayesian reconstruction methods for their good edge preservation capabilities. In this thesis, we consider two choices for the prior: the Gaussian anisotropic smoothness prior [91] and the non-Gaussian Cauchy prior [92, 93].

First, we consider the Gaussian prior: we assume thatm ∼ N(_m_∗_,Cm)_{, where} m∗ is the prior mean and Cm is the prior covariance; these are chosen based on experience and prior knowledge about the parameters of interest. In this case, the prior probability density is written as

π(_m)_∝exp

−¹

2kLm(_m−_m∗)k²

, (3.8)

whereLm is computed from the factorization of the inverse of the prior covariance matrixC_m⁻¹=L^T_mLm.

Second, we consider a non-Gaussian prior, which is in the form of Cauchy distribution. In this case of the edge-preserving Cauchy prior, the probability density can be expressed as a product

π(_m)_∝

N_x

∏

j=1 N_z j

∏

⁰=1

γh

(γh)²+ (m_j,j⁰−mj−1,j⁰)²

γh⁰

(γh⁰)²+ (m_j,j⁰−m_j,j⁰−1)²^, ^(3.9) where(j,j⁰)∈I²⊂Z², whereNxandNzare numbers of grid cells inxandzdirec- tions, respectively,γis the regularization parameter, andh,h⁰ >0 are discretization steps in thex- andz-directions. For more details on the numerical implementation of Cauchy priors, see [92, 93].

3.3.2 Computation of posterior quantities

In principle, the posterior distribution of the parameters can be (mathematically) considered as the solution of the inverse problem, but the fundamental goal is to compute posterior quantities of interest. These quantities, which are computed based on the posterior distribution include point and interval estimates, such as posterior means and modes, posterior covariances, and credible intervals. There are 16

(33)

two commonly used point estimates: maximum a posteriori (MAP) estimate that maximizes the posterior distribution, and conditional mean (CM) estimate that is defined as the expected value of the unknown parameter, given the data.

In case where the likelihood and prior information are both Gaussian (and we have considered a linear model Gm), the posterior distribution of the model vec- torm conditioned by the seismic datad is Gaussian [18]. The resulting Gaussian posterior density can be written as

π(m|d)_∝exp

−¹

2kLe(d−Gm−e∗)k²−¹

2kLm(m−m∗)k²

. (3.10) The conditional mean, which also corresponds to the maximum point of the posterior distribution (in the pure Gaussian case) is then given as

mCM=arg min

m

kLe(_d−_Gm−e∗)k²−¹

2kLm(_m−_m_∗)k². (3.11) The optimization problem (3.11) has a practical closed form solution similar to the one of linear squares estimation,

mCM=C_m⁻¹+_G^TC_e⁻¹G−1

G^TC⁻¹_e (_d−e∗)C⁻¹_m m∗

. (3.12)

The covariance of the posterior reflecting the uncertainty of the estimate can also be deterministically calculated as [94]

Cm|d=G^TC_e⁻¹G+C⁻¹_m −1

. (3.13)

It is important to note that the posterior density has a practical closed form solution only in the context of linear inverse problems with additive Gaussian noise and Gaussian prior considerations. In more general situations, for example, in cases of a non-linear observation model and/or non-Gaussian distributions, it is not possible to obtain the posterior and its quantities analytically. Hence, an approximate posterior distribution and posterior quantities are obtained by utilising appropriate sampling methods. In principle, these sampling methods generate representative samples from the posterior probability density, and the generated samples can be used to approximate the posterior quantities.

One of the most popular sampling methods is the Markov chain Monte Carlo (MCMC), which allow sampling from several classes of distributions by simulat- ing random walks in the space for which the distribution is to be evaluated [95].

The widely used MCMC algorithms with generic sampling mechanisms are the Metropolis–Hastings sampler and the Gibbs sampler. Several key issues must be considered to ensure that the MCMC algorithm creates a reliable picture of the posterior distribution of interest: measuring the efficiency of the algorithm; diagnosing convergence; establishing a stopping criteria; and estimating the Monte Carlo error [96–98].

For practical geophysical inverse problems such as high-dimensional non-linear seismic waveform inversion, MCMC methods can be computationally daunting.

Several different MCMC algorithms and approaches have been developed to improve the accuracy and efficiency of generating samples from high-dimensional posterior distributions [99–101]. In PublicationII, in which a non-Gaussian Cauchy

(34)

prior was employed, we use Metropolis-within-Gibbs algorithm following the practice in [92] to explore the posterior for parameter estimation as well as uncertainty analysis.

The posterior analysis also gives Bayesian interval estimates based on quantiles of the posterior for the parameter. The interpretation of these interval estimates, termed as credible intervals can help to quantify the uncertainty, which is useful in evaluating the quality of the inversion result. In general, ifmi is the ith element of the estimated parameter vector, then a 100(1−α)% credible interval for variable mi ∈Ris an interval[a,b]that satisfies the statement:

Z a

−∞p(m_i|d)dm_i = ^α 2 =

Z _∞

b p(m_i|d)dm_i, 0≤α≤1. (3.14) Herep(_m_i|_d)is the posterior marginal density of the variablemi. An interval estimate form_i can also be constructed based on the standard deviation. For example, given a Gaussian posterior distribution, the 95% credible interval form_iis computed

as h

m_∗|d,i−1.96σ_∗|d,i,m_∗|d,i+1.96σ_∗|d,ii

, (3.15)

wherem∗|d,i represents theith component ofmi|dand σ_∗|d,i is the standard deviation of the estimate. In the case of a Gaussian prior,σ_∗|d,i can be obtained from the diagonal of the covariance estimateC_m|d, and, in the case of a non-Gaussian prior, it can be estimated as the standard deviation of the MCMC samples.

3.4 BAYESIAN APPROXIMATION ERROR APPROACH (BAE)

At the core of any successful waveform inversion is a forward model that accurately represents the realities of seismic wave propagation. Yet accurate forward models are not often employed due to high computational burden of solving a large scale multidimensional non-linear inverse problem. Hence, in most cases, fast and less accurate (approximate) forward models are used in the inversion procedure, which may adversely affect the accuracy of the solution of the inverse problem. In this thesis, we adopt the Bayesian approximation error (BAE) approach [19], which allows the incorporation of all modelling errors resulting from the approximate forward model.

In the BAE approach, two computational forward models are considered: (1) an accurate model, which simulates the underlying physical process so that modelling errors are negligible compared to other errors such as observation noise, and (2) an approximate fast model. In this thesis, we aim to apply acoustic Born inversion to viscoelastic data. Therefore, we consider the elastic-viscous modelGas the accurate model and the acoustic Born modelGas the approximate fast model.

The observation model for the BAE approach can be written as

d=G(m) +e=Gm+e+e, (3.16) wheremrepresents the vector of accurate model parameters. Heree=G(m)−Gm is the approximation error term, which describes the difference between the wavefields predicted by the accurate elastic-viscous model G(m) and the approximate acoustic Born modelGm.

Thus, the basic idea behind the BAE approach is to utilise an approximation to the accurate physical model, and take into account the related modelling errors 18

Computational methods for seismic imaging and monitoring

Dissertations in Forestry and Natural Sciences

KENNETH MUHUMUZA

Computational methods for seismic imaging and monitoring

Kenneth Muhumuza

COMPUTATIONAL METHODS FOR SEISMIC IMAGING AND MONITORING

TABLE OF CONTENTS

1 Introduction

2 Modelling of seismic wave propagation

3 Deterministic and Bayesian inversion frameworks

∏

∏