Time Domain Full Waveform Inversion Using ADI Modeling

84 Pages·2013·2.05 MB·English

Time Domain Full Waveform Inversion
Using ADI Modeling
Bernd Klimm
Vom Fachbereich Mathematik
der Technischen Universit¨at Kaiserslautern
zur Verleihung des akademischen Grades
Doktor der Naturwissenschaften
(Doctor rerum naturalium, Dr. rer. nat.)
genehmigte Dissertation.
Datum der Disputation: 01. July 2013
Erster Gutachter: Prof. Dr. Axel Klar
Zweiter Gutachter: Prof. Dr. Thomas Bohlen
(D 386)

Abstract
Constructing accurate earth models from seismic data is a challenging task. Tra-
ditional methods rely on ray based approximations of the wave equation and reach
their limit in geologically complex areas. Full waveform inversion (FWI) on the
other side seeks to minimize the misfit between modeled and observed data without
such approximation.
Whilesuperiorinaccuracy,FWIusesagradientbasediterativeschemethatmakes
it also very computationally expensive. To reduce the costs of our two dimensional
time domain algorithm, we apply an alternating direction implicit (ADI) scheme to
theacousticwaveequation. TheADIschemecanbeseenasanintermediatebetween
explicit and implicit finite difference modeling schemes. Being less computationally
demanding than an implicit scheme it can handle coarser discretization compared
to an explicit scheme and can be efficiently parallelized.
With the modeling in place, we test and compare several inverse schemes. To
avoid local minima and improve speed of convergence, we use different hierarchical
approaches, for example advancing over different bands of increasing frequency con-
tent. We can then demonstrate the effectiveness of the method on the Marmousi
model and the 2004 BP model.
Acknowledgment
First of all, I gratefully acknowledge Prof. Dr. Axel Klar for his interest in this
topic, his advice and continuous support during this thesis.
I would like to express my sincere appreciation and gratitude to Dr. Norman Ettrich
for his excellent advice, support, guidance and friendship at the Fraunhofer Institute
for Industrial Mathematics ITWM.
I acknowledge Fraunhofer ITWM for my PhD scholarship and in particular Dr.
Franz-Josef Pfreundt and the HPC department for providing a great work environ-
ment as well as all the resources needed for this thesis.
Many thanks go to Statoil ASA and in particular to Frank Maaø and the seismic
imaging group at the Statoil research center in Trondheim for enabling two intern-
ships in their group. I am thankful for their warm help and support as well as many
instructive discussions.
I am thankful to my co-referee Prof. Dr. Thomas Bohlen for his interest in this
work and the discussions with his group at the Geophysical Institute at Karlsruhe
Institute of Technology.
I am indebted to Dr. Leo Neseman for his feedback and suggestions to the editing
of this thesis.
I am thankful to Maxim Illyasov for providing the RTM software, that provided
help in the initial development of the code created for this thesis.
IwouldliketoextendmygratitudetoallmycolleaguesatFraunhoferITWM,whose
friendship and support I deeply enjoyed over the past three years.
Finally, I owe very special gratitude to my wife for her unconditional support and
continuous encouragement.
4
Contents
1 Introduction 7
1.1 The Seismic Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2 Imaging and Inversion . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.2.1 Role of Full Waveform Inversion . . . . . . . . . . . . . . . . . 10
1.2.2 Probabilistic vs. Deterministic . . . . . . . . . . . . . . . . . . 10
1.2.3 Overview of FWI Algorithm in this Thesis . . . . . . . . . . . 11
1.3 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2 Wave Equation and Inverse Problem 12
2.1 The Forward Problem . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.1.1 2D Problem with Constant Velocity . . . . . . . . . . . . . . . 13
2.1.2 Inhomogeneous Velocity in 1D . . . . . . . . . . . . . . . . . . 15
2.2 The Inverse Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.2.1 Existence of a Solution in 1 D . . . . . . . . . . . . . . . . . . 17
2.2.2 Derivation of the Gradient . . . . . . . . . . . . . . . . . . . . 20
3 Modeling 23
3.1 Alternating Direction Method . . . . . . . . . . . . . . . . . . . . . . 23
3.1.1 Dispersion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.2 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.2.1 Boundary Conditions in 1D . . . . . . . . . . . . . . . . . . . 25
3.2.2 Boundary Conditions in 2D . . . . . . . . . . . . . . . . . . . 26
3.3 Validation of the Code . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.4 Comparison with Explicit Finite Difference Code . . . . . . . . . . . 29
3.4.1 The Marmousi Model . . . . . . . . . . . . . . . . . . . . . . . 29
3.4.2 Comparison of the Results . . . . . . . . . . . . . . . . . . . . 30
3.5 Checkpointing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
4 Nonlinear Minimization 35
4.1 Descent Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
4.1.1 Starting Model . . . . . . . . . . . . . . . . . . . . . . . . . . 37
4.1.2 Descent Direction . . . . . . . . . . . . . . . . . . . . . . . . . 37
4.1.3 Computing the Step Size . . . . . . . . . . . . . . . . . . . . . 40
4.1.4 Stopping Criteria . . . . . . . . . . . . . . . . . . . . . . . . . 40
4.2 Improving Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . 41
4.2.1 Gradient in Water . . . . . . . . . . . . . . . . . . . . . . . . 41
4.2.2 Alternative Error Functions . . . . . . . . . . . . . . . . . . . 41
4.3 Multiscale Approaches . . . . . . . . . . . . . . . . . . . . . . . . . . 44
4.3.1 Multiscale Approach over Frequencies . . . . . . . . . . . . . . 44
4.3.2 Multiscale Approach over Time Damping . . . . . . . . . . . . 45
5
Contents
4.3.3 Multiscale Approach over Offsets . . . . . . . . . . . . . . . . 46
5 Test Cases 48
5.1 Marmousi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
5.1.1 First Gradient . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
5.1.2 General Descent Scheme . . . . . . . . . . . . . . . . . . . . . 50
5.1.3 L Norm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
1
5.1.4 Noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
5.1.5 Multiscale Schemes . . . . . . . . . . . . . . . . . . . . . . . . 62
5.2 2004 BP Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
6 Summary and Conclusion 75
Bibliography 78
6
1 Introduction
Finding natural resources is an exceedingly challenging task which is becoming more
and more difficult as resources are getting depleted. Oil & gas companies are drilling
deeper and in geologically more complex environments than before and costs of
drilling often exceed one million dollar per well. On the other side this means
that there are high incentives in trying to obtain accurate models of the area of
interest in order to better plan for the drilling and production and to reduce the
risk of unfortunate dry wells. The process of generating these models is still a
challenging task, that usually involves many disciplines from geology and geophysics
to mathematics. In this thesis, we study full waveform inversion (FWI), which
currently is the most advanced inversion algorithm to obtain earth properties from
seismic data as part of the whole exploration process.
1.1 The Seismic Problem
The earth propagates vibrations over long distances. We observe this naturally
during earthquakes, where surface waves are transmitted on continental scales. For
oil and gas exploration the earth is excited artificially by seismic experiments and
vibrations are measured locally. A seismic source generates vibrations (e.g., pressure
waves) and receivers record the reflected data at the surface. The principle is similar
to an ultrasound unit, although at much larger scales.
In marine exploration a seismic vessel creates acoustic waves inside the water with
an array of airguns. The vessel tows a long array of hydrophones some distance
behind the source. Hydrophones are basically special microphones in water that
measure pressure and record the reflected signal at the surface as displayed in Figure
1.1. One line of hydrophones is called a streamer. The vessel then travels forward
with this configuration and activates the source (i.e., shoots) in regular intervals
(e.g., every 25 m). That means that when we speak of a seismic experiment, seismic
acquisition, or survey, we actually mean multiple single shot experiments carried out
in sequence. With the different shot and receiver positions, we observe reflections
from the underground under multiple angles.
The source is designed to (ideally) produce a wavelet, a short, essentially compact
signal that often has a peak frequency around or below 50 Hz. The distance between
the source and the receivers is called offset. In modern acquisitions the streamer
lengths are usually around 8 km or more with receivers spaced as close as every 25
m. The data is usually recorded for more than 8 seconds. Furthermore, a seismic
vessel can tow several streamers simultaneously which then cover part of the 2D
surface. All in all a single seismic survey can thereby produce terabytes of data.
There are more complex acquisitions that include multiple seismic vessels in one
survey or ocean bottom cable acquisitions, where receivers are laid out on the ocean
floor.
7
1. Introduction
Figure 1.1: The marine seismic experiment
For the this thesis, we have to limit the scope of our investigation. We will
assume that we have a single two-dimensional section, that means data from a
single streamer on a straight line.
The principals of land acquisitions are quite similar to marine acquisitions. The
receivers (geophones) are laid out on the ground and a moving seismic truck gener-
ates vibrations. However, the geometry of the acquisition is often more complex as
source and receivers need to follow the natural terrain. In this work, we do not deal
with specifics of that. Furthermore, marine acquisitions are in some sense easier as
the properties in the waterlayer are often well known or can easily be determined.
This has some advantage in the computation as we see later. Nevertheless, apart
from some rather technical differences, the general scheme and many observations
we make with the marine experiment in this work could be transferred to the land
survey as well.
The next step after the data acquisition is to use the data to obtain an image
of the underground. For this we need to make assumptions towards the correct
physical equations to use. Most general, wave propation in media can be described
with the full elastodynamic wave equation
1 ∂ ∂ ∂2 1
C δ u = F (1.1)
iklj kl l k
ρ∂i ∂j − ∂t2 ρ
(cid:18) (cid:19)
where δ is the Kronecker delta, ρ the density, F are body forces, and C is the
k
stiffness tensor with up to 21 independent variables. Unfortunately, there are not
nearly enough data available to solve the full problem. Even if there were, the
required compute power would be completely out of range today. Thus we need
to approximate the earth with a less detailed model. Research is done with the
isotropic wave equation and even with certain kinds of anisotropy (e.g., transverse
isotropy) but for the scope of this work we reduce the complexity even more and
consider only the acoustic case with constant density
1 d2
∆u u = f. (1.2)
− c2dt2
8
1. Introduction
Here, f is the source term, u denotes the scalar pressure in the medium, and c is the
velocity or acoustic wave speed (i.e., speed of sound). In light of the computational
complexity of full waveform inversion and because we are testing general properties
of the inversion scheme, it seems reasonable for us to use just one parameter c.
There are several introductory books available that explain the seismic acquisition
process in depth and also discuss possible effects of elasticity and anisotropy in
greater detail, for example [Yil01], [Rob10], [GV10], [Jon10]. Note, that we also
exclude the important step of data preprocessing from this work. Steps like noise
removal, data regularization, filtering, sorting, removal of direct wave and so on are
commonly performed before inversion and imaging tools are used. Please also refer
to the cited introductory books for details.
1.2 Imaging and Inversion
Seismic data imaging is traditionally separated into two parts: velocity model build-
ing and migration. If a velocity model is available, migration means that a reflec-
tivity model is generated based on the seismic data and that velocity model. Of
course, a correct velocity model would already display the reflectors, however build-
ing a velocity model is much more difficult. So the common workflow is to start
with a good, approximate velocity model and then use this model to migrate the
seismic data to the “correct” position in the ground. The reflectivity model is then
used for interpretation. The process is called imaging or migration. There are many
industry-standard codes available, many of them us a ray approximation of the wave
equation. Most accurate among the migration algorithms is reverse time migration
(RTM), which uses the full wave equation to compute the reflectivity and is closely
related to full waveform inversion, as we will see in chapter 2.2.2. An overview over
different migration algorithms can be found in [Rob10].
Velocity model building, the focus of this work, is the more challenging part and
is done iteratively. First, velocity models might be created based on geological
reasoning and stacking analysis. Then tomography techniques are often applied,
that are commonly based on ray approximations of the wave equation. Methods
using the wave equation like full waveform inversion, which we investigate here,
are most computational expensive and might only be used for complex parts of
the whole problem, or when the simpler methods do not yield sufficiently accurate
results. Please refer to [Jon10] for an overview of common velocity model building
tools.
The whole process of seismic data processing normally encompasses multiple cy-
cles of model building, imaging, quality control and interpretation and can take
several months. What methods are employed and what amount of detail is required
is ultimately a management decision. Common questions that guide the process
might be: “Given certain amount of time and computational resources, what is the
best image/ model that we can obtain?” or “What do we have to do in order to get
an image/ model that can be used to make a decision (e.g., is it safe to drill)?”
9
1. Introduction
1.2.1 Role of Full Waveform Inversion
In this whole process, the tomography based methods are very efficient for simple
geological targets, for which the ray approximation holds, but reach their limits
in complex situations, for example around salt domes, fold and thrust belts, and
foothills [VO09]. In those cases even the separated process of building a veloc-
ity model and subsequential migration might no longer be applicable as even high
quality imaging methods with an insufficient detailed velocity model might suffer
accuracy. In those cases our full waveform inversion can lead to improvements.
Full waveform inversion goes back to Lailly in 1983 [Lai83] and Tarantola in
1984 [Tar84]. The underlying mathematical principles of least square minimization
and adjoint state method are even older. Nevertheless it has not gained much
attention until very recently because of its computational requirements and it is
not in routine productive use up to now. Inversion schemes that use raybased
approachesorcombinewaveequationmodelingwithraytheoryweredeveloped(e.g.,
[CC88], [SD92], [BCSJ01]) as well as schemes based on finite difference modeling
(e.g., [KCL86], [IDT88], [CPN+90], [PDT90]) and finite element modeling [CMS08].
A second reason, why full waveform inversion was not considered successful in the
past was because of the lack of long offsets in the seismic data. To start the inversion
for a first (low frequency) background velocity model, diving waves are needed (i.e.,
waves that are not reflected, but slowly bended back upwards, when c increases
with depth). To record diving waves with considerable depth, long receiver offsets
are needed (or crosshole data measured from within boreholes). This data is now
availableinmodernsurveysandcanbeexploitedbyfullwaveforminversion[Mor87],
[Mor88], [PW90], [PSWW96].
1.2.2 Probabilistic vs. Deterministic
Although modern explorations cover wider offsets and might use more shots and
longer receiver data, the amount of data is always limited and the inverse problem
of recovering the earth c does not have a unique solution. It is possible to use prob-
abilistic inverse theory, where a probability function p(c) is defined over the space of
velocity models. The probability measures both the misfit between computed and
measured data U(c) and U and contains a priori knowledge about possible velocity
d
models [Tar05]. However, this probability function p is not known for real world
problems and could only be approximated by sampling the model space and solving
the wave equation for each sample. Unfortunately the number of samples required
grows exponentially with the size of the model space [Fic11]. For the seismic prob-
lem the velocity model can contain a million parameters or more, which makes this
method infeasible.
Instead, we have to use deterministic inverse theory in our work. There, we define
a misfit function E which will depend on the computed and observed data. We will
then choose a starting model c , being (presumably) close to the real model c˜, and
0
start an iterative scheme with the goal of finding a c under which E(c) becomes
minimal.
10