Seismic Inversion Book PDF
Seismic Inversion Book PDF
Seismic Inversion Book PDF
S. P. Maurya
N. P. Singh
K. H. Singh
Seismic
Inversion
Methods:
A Practical
Approach
Springer Geophysics
The Springer Geophysics series seeks to publish a broad portfolio of scientific
books, aiming at researchers, students, and everyone interested in geophysics. The
series includes peer-reviewed monographs, edited volumes, textbooks, and
conference proceedings. It covers the entire research area including, but not
limited to, applied geophysics, computational geophysics, electrical and electro-
magnetic geophysics, geodesy, geodynamics, geomagnetism, gravity, lithosphere
research, paleomagnetism, planetology, tectonophysics, thermal geophysics, and
seismology.
123
S. P. Maurya N. P. Singh
Department of Geophysics Department of Geophysics
Institute of Sciences Institute of Sciences
Banaras Hindu University Banaras Hindu University
Varanasi, Uttar Pradesh, India Varanasi, Uttar Pradesh, India
K. H. Singh
Department of Earth Sciences
Indian Institute of Technology Bombay
Mumbai, Maharashtra, India
This Springer imprint is published by the registered company Springer Nature Switzerland AG
The registered company address is: Gewerbestrasse 11, 6330 Cham, Switzerland
Contents
v
vi Contents
1.1 Introduction
Seismic inversion is a procedure that helps extract underlying models of the physical
characteristics of rocks and fluids from seismic and well-log data. In the absence of
well data, the properties can also be inferred from the inversion of seismic data alone
(Krebs et al. 2009). In the oil and gas industry, seismic inversion technique has been
widely used as a tool to locate hydrocarbon-bearing strata in the subsurface from
the seismic reflection data (Morozov and Ma 2009; Lindseth 1979). This method
dramatically increases the resolution of seismic data and hence helps to interpret
seismic data.
The physical parameters that are of interest to a modeler performing inversion
are impedance (Z), P-wave (V P ) and S-wave (VS ) velocity and density (ρ). Lame
parameters that are sensitive towards fluid and saturation in rocks (Clochard et al.
2009) can also be derived from inverted models of impedances. The petrophysical
parameters like porosity, sand/shale ratio, gas saturation, etc. can be estimated further
with the help of inverted volumes (Goodway 2001). These petrophysical parameters
© The Editor(s) (if applicable) and The Author(s), under exclusive licence 1
to Springer Nature Switzerland AG 2020
S. P. Maurya et al., Seismic Inversion Methods: A Practical Approach,
Springer Geophysics, https://doi.org/10.1007/978-3-030-45662-7_1
2 1 Fundamental of Seismic Inversion
added strength to the seismic data interpretation which is a very crucial process for
any exploration project.
To understand seismic inversion methods, one needs to first understand forward
modeling. The seismic forward modeling uses the principle of convolution theory
which states that seismic trace can be generated by the convolution of wavelet with
earth reflectivity. The sound waves are sent to the subsurface that interacts with the
earth’s layer and returns on the surface which is recorded by the receivers. These
recorded signals are called observations (Seismic signal) and the entire process is
called forward modeling (Maurya et al. 2018). On the other hand, in seismic inversion
methods, we have observation available and we seek to find the subsurface geological
model that is a representation of observation. Both processes can be understood from
Fig. 1.1.
There are many geophysical methods used to explore oil and gas from the subsurface
but the most important technique is seismic imaging. The imaging means the visual
representation of the earth’s subsurface model. This is one of the ultimate goals
of the geophysicist. Seismic forward modeling can be implemented by using both
numerical studies in geoscience and computation technology. The forward modeling
procedure uses an elastic impedance method that generates synthetic seismograms
from velocities and densities of the subsurface layers (Connolly 1999). The elastic
impedance at each interface is calculated as a function of the offset. The resulting
impedance series is transformed into the reflectivity series and convolved with the
source wavelet to get a stacked seismic gather. The impedance (Z ) is computed from
the product of velocity (v) and density ρ.
Z = Vρ (1.1)
1.2 Seismic Forward Modeling 3
The zero offset reflectivity series can be calculated from the impedance as
following.
Z j+1 − Z j
Rj = (1.2)
Z j+1 + Z j
where Z j is the seismic impedance of jth layer, and R j is seismic reflectivity of the
interface between jth and (j + 1)th layer. The reflection coefficient for the angle-
dependent incident wave is estimated using the following formula (Bachrach et al.
2014).
1 I p Vs2 2 Is 1 Vs2 ρ
R(θ ) = 1 − tan θ
2
− 4 2 sin θ − tan θ − 2 2 sin θ
2 2
2 Ip Vp Is 2 Vp ρ
(1.3)
The synthetic seismogram is calculated from the reflection coefficient using the
following equation.
where S(t) is synthetic seismogram, W (t) is source wavelet, R(t) is the reflection
coefficient of the subsurface and N (t) is additive noise and generally assumed to be
zero for simplicity. The seismic forward modeling method gives the understanding
of seismic travel time, elastic impedance, arrival time, earth’s reflectivity, seismic
amplitude generated by the seismic wave and the other aspects.
Figure 1.2 demonstrates the use of forward modeling in faulted and anticline
geological models. Figure 1.2a shows geological model of the subsurface whereas
Fig. 1.2b depicts the corresponding seismic section generated using the forward
modeling technique. From the figure, we can see that the seismic gather (Fig. 1.2b)
shows an approximately same geological structure which is in the geological model
(Fig. 1.2a). Now our aim is to find this geological model from the observation i.e.
seismic section which can be achieved by the seismic inversion methods.
Seismic inversion methods involve mapping rock and fluid properties of the sub-
surface of the earth using seismic measurements made on the surface of the earth
as input. In fact, all inversion methods aim to estimate the geophysical properties
of the subsurface from the measurement made on the surface. In the seismic inver-
sion process, there are three main issues that need to be addressed carefully to get a
broadband spectrum with high-resolution images of the subsurface. The first prob-
lem arises due to the band-limited nature of the seismic data which means seismic
4 1 Fundamental of Seismic Inversion
Fig. 1.2 a Shows geological model and b depicts seismic section expressed geological model
data does not have a low and high-frequency component. The seismic data generally
have 10–80 Hz frequency and hence does not contains frequency less than 10 Hz and
greater than 80 Hz however these frequencies are very important for the interpreta-
tion. On the other hand, well log data has both low and high frequencies and can be
integrated with the seismic data to compensate for this. The second problem is to
use a seismic wavelet. From the convolution theory, we have seen that the seismic
wavelet is used to generate synthetic data and hence accurate wavelet estimation is
critical for successful inversion results (Russell 1988; Maurya and Singh 2018). The
third and most important problem is that the seismic inversion is non-unique. There
may be possibly more than one solution to the same problem. To reduce the number
of solution one needs other constraints to bind the inversion results. These constraints
can be prior to geological information, well log data, etc.
Figure 1.3 shows an example of seismic inversion methods. Figure 1.3a shows
seismic data from the Blackfoot field, Canada whereas Fig. 1.3b shows seismic
inversion results. From Fig. (1.3) one can notice how dramatically the image quality
has been increased. From the seismic section, one has an only amplitude and hence
cannot interpret much. However from the inverted section, one can classify the sand
formation, shale formation and hence can identify the productive zone.
To understand the seismic inversion technique, one must first understand the
physical processes which are involved in generating these data. The seismic data is
generated using a forward modeling technique that uses the convolution model for
its implementation. Initially, therefore one should look at the basic convolutional
model of the seismic trace in the time and frequency domains. This model has three
1.3 Seismic Inversion 5
components, first is earth’s reflectivity, second is seismic wavelet, and the third is the
noise. After understanding the concepts and the problems which can occur, one is in
a position to look at the methods which are currently used to invert seismic data.
The convolution model is the most common one-dimensional model for the seismic
trace. The convolution model states that the seismic trace can be generated by the
convolution of seismic wavelet with the earth’s reflectivity series along with the
addition of noises. It can be written mathematically as follows.
where * implies convolution process, S(t) is a seismic trace, R(t) is earth reflectivity,
W (t) is wavelet and N (t) is the noise component. By considering perfect case one
can consider noise component to be zero and hence Eq. 1.5 can be written simply as
Equations 1.5 and 1.6 represents the convolution model in the time domain. An
alternate form of the above equation is frequency domain, can be obtained by taking
the Fourier transform of the above equations which results in
S( f ) = W ( f ) × R( f ) (1.7)
6 1 Fundamental of Seismic Inversion
θ S ( f ) = θW ( f ) + θ R ( f ) (1.9)
Fig. 1.4 Generation of seismic data by using the convolution model. Track 1 shows the earth’s
reflectivity, track 2 shows minimum phase wavelet and track 3 shows seismic trace
1.4 The Convolution Model 7
Now, after having a good understanding of the basic convolution model which
is used by most of the seismic inversion methods. We are now in a position to look
into the seismic inversion methods and their classifications. The seismic inversion
methods are summarized in the following sections.
The seismic inversion techniques can be divided into two broad categories, Post-
stack, and Pre-stack inversion methods. The first approach is the most commonly
used where the effect of the wavelet is removed from the seismic data and a high-
resolution image of the subsurface is produced (Chen and Sidney 1997). The second
approach relies on model building from well log, seismic and geological data (Down-
ton 2005). This also generates a high-resolution image of the subsurface from which
reservoir properties are calculated. A reliable estimate of the reservoir properties is
critical in the decision-making process during the development phase (Pendrel 2006).
These inversion methods are further divided into subparts which are discussed in the
following sections.
Post-stack seismic inversion techniques are categorized as the first type. This type
of inversion results in acoustic impedance volume utilizing seismic data through the
integration of the well data and a basic stratigraphic interpretation. This impedance
volume can be used to estimate reservoir properties away from well (Russell and
Hampson 1991; Morozov and Ma 2009). Some of the advantages of post-stack
inversion are mentioned below.
1. As the acoustic impedance is a layer property; hence stratigraphic interpretation
is easier on impedance data than seismic data.
2. The reduction of wavelet effects, side lobes, and tuning enhance the resolution
of subsurface layers.
3. Acoustic impedance can be directly computed and compared to well log
measurements that serve as a link to reservoir properties.
4. Porosity can be related to the acoustic impedance. Using geostatistical methods
these impedance volume can be transformed into the porosity volume within the
reservoir.
5. Acoustic impedance can be utilized to locate individual reservoir regions.
6. It takes very less time than pre-stack inversion.
7. It does not give shear wave information to discriminate against the fluid effects.
Further, post-stack inversion can be divided into two parts namely deterministic
inversion which includes model-based inversion methods and second is Stochastic
8 1 Fundamental of Seismic Inversion
The Pre-stack inversion falls into the second category of the seismic inversion tech-
niques. The estimate of the elastic properties of the subsurface such as the S-wave
velocity of the subsurface layers which are sensitive to fluid saturation can be obtained
from Pre-stack inversion (Moncayo et al. 2012).
Pre-stack inversion transforms seismic data (angle/offset gathers) into P-
impedance, S-impedance and density volumes through the integration of well data
and horizon information from seismic data. P-impedance and V p/V s ratio are reli-
able, depending on target depth and acquisition configuration, and can be used to pre-
dict reservoir properties away from well (Carrazzone et al. 1996). Pre-stack seismic
inversion provides several benefits.
1. The P-impedance, S-impedance, and density give layer properties, whereas
seismic data is an interface property.
2. Enhanced resolution of sub-surface layers due to the reduction of wavelet effects,
tuning and side lobes.
3. Acoustic impedance can be directly compared to well log measurements which
in turn are linked to reservoir properties.
4. Compared with other inversion techniques (e.g. post-stack inversion), the data
offers additional information to distinguish between lithology and fluid effects.
The most common methods which fall in this category are simultaneous inversion,
Elastic impedance inversion, and AVO inversion methods.
Simultaneous inversion is the first type of pre-stack seismic inversion. Pre-stack
seismic gather contains additional information i.e. S-wave velocity which travels
slowly in the subsurface and contains more information about the rock properties of
the earth. This information can be estimated from the pre-stack gathers using several
seismic inversion methods. A common approach is simultaneous inversion of pre-
stack seismic data, which inverts for several rock property parameters simultaneously.
The Aki and Richards (1980) formula gives access to the approximate reflectivity at
the various offsets in the pre-stack domain and can be expressed as follows.
V P ρ VS
R(θ ) = a +b +c (1.11)
VP ρ VS
2
2
where a = 1
2 cos2 θ
, b = 0.5 − 2 VVPS sin2 θ , c = −4 VVPS sin2 θ,
10 1 Fundamental of Seismic Inversion
ρ1 + ρ2 V P1 + V P2
ρ= , ρ = ρ2 − ρ1 , V P = ,
2 2
VS1 + VS2
V P = V P2 − V P1 , VS = ,
2
θ1 + θ2
VS = VS2 − VS1 , θ =
2
These functions are used to compute the earth’s reflectivity. Thereafter, reflectivity
is convolved with the wavelet to obtain the synthetic seismic gather. Further, this
synthetic gather is compared with the original seismic gather of the area and the
misfit is computed between them. The model is subsequently perturbed and a new
comparison made to reduce the misfit and the least misfit model is selected as the
final solution.
The other inversion method used commonly is the elastic impedance inversion.
Connolly (1999) introduced the concept of Elastic Impedance and defined a function
R which is dependent on the incidence angle and is related to the relative variation
of Elastic Impedance.
1 E I 1
R(θ ) ≈ ≈ ln(E I )
2 EI 2
This function R is now called the Elastic Impedance, in analog to the acoustic
impedance concept. The angle-dependent P-wave reflectivity is also approximated
by the well-known simplified description of the Zoeppritz equations (Aki-Richards):
E I = V P(1+tan θ)
VS(8K sin θ) (1−8K sin2 θ)
2 2
ρ (1.14)
a mixture of methodologies for pre-stack and post-stack inversion. This mix allows
big data amounts to be efficiently reversed in the lack of good information.
Figure 1.5 discusses the classification of seismic inversion methods graphically
although some recent development has been done which added more classification
but above discussed inversion methods are convention techniques that are still in
use to estimate subsurface properties from seismic and well log data. These inverted
sections strengthen to seismic data interpretation and hence help to characterize the
reservoir.
The algorithm calculates a search path using a local property of the objective
function, determining an update or increase to the present model, given a starting
solution. Only a fraction of the measured update is applied, however, and a step-
length factor is used to ensure that the new objective function value is always lower
than the current value (Chunduru et al. 1997). Thus the critical factors are the update
calculations and the step length. Methods of optimization vary on the basis that how
these two parameters are calculated.
data calculated using an assumed earth model (Sen and Stoffa 2013). Physical param-
eters that characterize the characteristics of rock layers, such as compressional wave
velocity, shear wave speed, resistivity, etc., describe the earth model.
The solution of a linear equation system is the same as a quadratic function mini-
mization. Local minimization processes mentioned above are appropriate when there
is a single minimum of the cost function. In particular, however, inverse geophysical
problems are non-linear and extremely complicated. The cost function is probable
to demonstrate multiple minima in such instances. In such cases, local optimization
techniques are anticipated to fail as they tend to cover the closest local minima (Sen
et al. 1995). Local minimization algorithms are anticipated to cover a fake solution
and produce bad outcomes unless there is adequate a priori data to make an intel-
ligent guess about the original model. Gradient-based optimization systems do not
provide the means to jump from a local minimum to a worldwide solution (Maurya
and Singh 2019). However, in most situations, computationally more expensive are
the global optimization schemes compared to local optimization schemes. Figure 1.7
shows the difference between local and global solutions.
The primary motivation to develop efficient global optimization technique lies in
the fact that the technique should be practicable in obtaining better results compared
to the local optimization techniques in situations where the problem is complex
and the model dimensionality is large. There are many global optimization methods
available but Genetic algorithm and Simulated annealing are broadly used in the
seismic inversion.
The technique of the genetic algorithm (GA) is based on the analogy that the
genetic modifications that occur in the living species work towards making the species
smarter and more adaptive to the changing natural environment. One of the powerful
tools for global optimization is the genetic algorithm; it is also based on the principle
of random walking in the space parameter (Sen and Oltz 2006; Maurya et al. 2018).
GA has artificial intelligence that can handle issues that are highly nonlinear. GA
needs no derivatives or information on curvature. Therefore, once the forward prob-
lem is solved, the reverse problem can be solved automatically because the whole
exercise consists of selecting some models, calculating synthetic data, comparing
d prs with d obs , calculating the cost function or error function and selecting better and
better models through certain guidelines when deciding on the acceptance criterion
(Sen and Stoffa 2013).
14 1 Fundamental of Seismic Inversion
The simulated annealing (SA) uses techniques of random walking, but with some
artificial intelligence. It brings the concept of temperature in this technique and
which is controlling parameters even though the temperature has nothing to do with
an inverse problem. The function of energy in thermodynamics is replaced by the
function of error. This error function is also referred to in global optimization as cost
function or energy function. In thermodynamics, the probability density functions of
Gibbs can be expressed as follows.
−E
ex p K Tj
P Ej =
(1.15)
−E
ex p K Tj
where C D is the data covariance. The inverse problem begins with a very high initial
temperature of the controlling parameter. The temperature is progressively reduced
in consecutive iteration to make Gibb’s probability density function more and more
sensitive (Vestergaard and Mosegaard 1991). In this manner, the cost function is
minimized. This technique is discussed in detail in Chap. 6.
1. The spatial continuity of the well log data is quantified using variograms.
2. A statistical relationship is derived between the log and seismic data at all well
locations using cross-validation plots. The single attribute analysis, multivariate
regression uses linear relationship whereas the Neural Network uses a non-linear
relationship.
3. These linear and non-linear relationships are then used to estimate well log
property away from the boreholes.
4. The predicted well log property is evaluated for its reliability.
In single attribute analysis, first, a variety of attributes are estimated directly or
indirectly from the seismic data and are analyzed to get the best attribute. The best
attribute selected on the basis of its correlation with the target log values. Further, the
best attribute is cross plotted with target value from the well log data and a best fit
straight line is chosen which gives the desired relationship. This desired relationship
is then used to predict petrophysical parameters in the inter-well region.
Multi-attribute regression is the second type of geostatistical method used to
predict petrophysical parameters. The basic principle is similar to the single attribute
analysis. The multi-attribute analysis has only one difference from the single attribute
analysis, in the sense that this method uses more than one attribute at a time. In this
method, all the attributes are analyzed and the best combination of attributes (more
than one) is selected and cross plotted with a target value to give linear relationship
which is used for further analysis.
Till now the analysis is linear which have some limitation but now we are extending
our analysis to nonlinear. The neural network falls in this category. It is divided into
two parts namely multi-layer feed-forward neural network and probabilistic neural
network. A multilayer feed-forward neural network is an interconnection of neurons
in which data and calculations flow in a single direction, from the input data to the
outputs with intermediate one or more hidden layers. Each layer consists of nodes,
and the nodes are connected with a particular weight. These weights decide the results
of the output layer (Dubrule 2003; Hampson et al. 2001). In this book, we have used
as many input nodes as the number of attributes used for the analysis. Mostly, the
output layer consists of one node and hence one output results, since our target is to
predict one single petrophysical parameter at a time. Figure 1.8 demonstrates MLFN
architecture in which four attributes are used as four-node in the input layer, one
hidden layer with three nodes are used and finally one node of the output layer has
been used. All nodes from the input layer are connected to the hidden layer with the
weights. And all nodes of the hidden layer are connected to the output layer by the
property called summation.
An alternative type of neural network is a probabilistic neural network (Masters
1995; Specht 1990, 1991). This method is actually based on a mathematical interpo-
lation technique that uses neural network architecture for its implementation. This is
the advantage of using PNN rather than MLFN since by studying the mathematical
formulation one can understand the behavior in a better way.
Forgiven training datasets, the PNN technique assumes that the predicted log can
be written as a linear combination of the log values in the training data. Let’s we have
16 1 Fundamental of Seismic Inversion
three attributes x = A1 j , A2 j , A3 j , then the predicted log value can be estimated
as follows.
n
i=1 L i exp(−D(x, x i ))
L (x) = n (1.17)
i=1 exp(−D(x, x i ))
where
3
x j − xi j 2
D(x, xi ) = (1.18)
j=1
σj
The quantity D(x, xi ) is the distance between the input point and each of the
training points xi .
Seismic Inversion techniques have been routinely used to map the subsurface
structures both in the field of seismology and seismic exploration. In seismic explo-
ration, inversion helps extract underlying models of the physical characteristics of
rocks and fluids that are helpful for reservoir characterization. Seismic inversion has
several limitations. First, the seismic frequency band is limited to about 12–80 Hz,
and therefore, both the low- and high frequencies in the data are missing. Well-log
data provides information at these missing frequencies. Secondly, the non-uniqueness
of the solution leads to multiple possible geologic models that are consistent with the
observations. In addition, the seismic inversion methods ignore multiple reflections,
transmission loss, geometric spreading and frequency-dependent absorption. The
preferred way to reduce these uncertainties is to use additional information (mostly
1.8 Geostatistical Inversion 17
from well logs) containing low and high-frequency components which help to con-
strain the deviations of the solution from the initial-guess model. The final result,
therefore, relies on the seismic data as well as on this additional information, and
also on the details of the inversion methods themselves.
The details of all the discussed inversion methods are explained one by one in
the coming chapters. First, the mathematical background of these inversion methods
is provided and then the capacity of methods to estimate subsurface property is
presented with synthetic as well as real data examples. The effort is made to explain
every aspect of these seismic inversion methods and their application to the real data.
References
Aki K, Richards PG (1980) Quantitative seismology: theory and methods. New York, p 801
Ansari HR (2014) Use seismic colored inversion and power law committee machine based on
imperial competitive algorithm for improving porosity prediction in a heterogeneous reservoir. J
Appl Geophys 108:61–68
Bachrach R, Sayers CM, Dasgupta S, Silva J (2014) Seismic reservoir characterization for uncon-
ventional reservoirs using orthorhombic AVAZ attributes and stochastic rock physics modeling.
SEG Tech Prog Exp Abs: 325–329
Banihasan N, Riahi MA, Anbaran M (2006) Recursive and sparse spike inversion methods on
reflection seismic data. University of Tehran, Institute of Geophysics
Bortfeld R (1961) Approximations to the reflection and transmission coefficients of plane
longitudinal and transverse waves. Geophys Prospect 9(4):485–502
Carrazzone JJ, Chang D, Lewis C, Shah PM, Wang DY (1996) Method for deriving reservoir
lithology and fluid content from pre-stack inversion of seismic data. US Patent 5:583,825
Chen Q, Sidney S (1997) Seismic attribute technology for reservoir forecasting and monitoring.
Lead Edge 16(5):445–448
Chunduru RK, Sen MK, Stoffa PL (1997) Hybrid optimization methods for geophysical inversion.
Geophysics 62(4):1196–1207
Clochard V, Delépine N, Labat K, Ricarte P (2009) Post-stack versus pre-stack stratigraphic inver-
sion for CO2 monitoring purposes: a case study for the saline aquifer of the Sleipner field. In:
SEG Annual Meeting, Society of Exploration Geophysicists
Connolly P (1999) Elastic impedance. Lead Edge 18:438–452
Downton JE (2005) Seismic parameter estimation from AVO inversion. In: M.Sc. thesis. University
of California Department of Geology and Geophysics, pp 305–331
Dubrule O (2003) Geostatistics for seismic data integration in Earth models. Distinguished Instructor
Short Course, vol 6. SEG Books
Goodway W (2001) AVO and lame constants for rock parameterization and fluid detection. CSEG
Recorder 26(6):39–60
Haas A, Dubrule O (1994) Geostatistical inversion- a sequential method of stochastic reservoir
modelling constrained by seismic data. First Break 12(11):561–569
Hampson DP, Schuelke JS, Quirein JA (2001) Use of multiattribute transforms to predict log
properties from seismic data. Geophysics 66(1):220–236
Krebs JR, Anderson JE, Hinkley D, Neelamani R, Lee S, Baumstein A, Lacasse MD (2009) Fast
full-wave-field seismic inversion using encoded sources. Geophysics 74(6): WCC177–WCC188
Lancaster S, Whitcombe D (2000) Fast-track “colored” inversion. SEG Expanded Abs 19:1572–
1575
Leite EP (2010) Seismic model based inversion using matlab. Matlab-Modell Program Simul: 389
18 1 Fundamental of Seismic Inversion
Lindseth RO (1979) Synthetic sonic logs-a process for stratigraphic interpretation. Geophysics
44(1):3–26
Masters T (1995) Advanced algorithms for neural networks: a C++ sourcebook. John Wiley & Sons,
Inc.
Maurya SP, Singh KH, Singh NP (2019) Qualitative and quantitative comparison of geostatistical
techniques of porosity prediction from the seismic and logging data: a case study from the
Blackfoot Field, Alberta, Canada. Mar Geophys Res 40(1):51–71
Maurya SP, Singh KH (2019) Predicting porosity by multivariate regression and probabilistic
neural network using model-based and coloured inversion as external attributes: a quantitative
comparison. J Geol Soc India 93(2):207–212
Maurya SP, Singh NP (2018) Application of LP and ML sparse spike inversion with probabilistic
neural network to classify reservoir facies distribution—a case study from the blackfoot field,
Canada. J Appl Geophys Elsevier 159(2018):511–521
Maurya SP, Singh KH, Kumar A, Singh NP (2018) Reservoir characterization using post-stack
seismic inversion techniques based on real coded genetic algorithm. J Geophys 39(2):95–103
Maurya SP, Singh KH (2015a) LP and ML sparse spike inversion for reservoir characterization-a
case study from Blackfoot area, Alberta, Canada. In 77th EAGE Conference and Exhibition,
2015(1):1–5
Maurya SP, Singh KH (2015b) Reservoir characterization using model based inversion and prob-
abilistic neural network. In 1st International conference on recent trend in engineering and
technology, Vishakhapatnam, India
Moncayo E, Tchegliakova N, Montes L (2012) Pre-stack seismic inversion based on a genetic
algorithm: a case from the Llanos Basin (Colombia) in the absence of well information. CT&F-
Ciencia Tecnología y Futuro 4(5):5–20
Morozov IB, Ma J (2009) Accurate post-stack acoustic-impedance inversion by well-log calibration.
Geophysics 74(5):R59–R67
Pendrel J (2006) Seismic inversion- still the best tool for reservoir characterization. CSEG Recorder
26(1):5–12
Russell B, Hampson D, Schuelke J, Quirein J (1997) Multiattribute seismic analysis. Lead Edge
16:1439–1444
Russell B, Hampson D (1991) Comparison of post-stack seismic inversion methods. In: SEG
technical program expanded abstracts. Society of Exploration Geophysicists, pp 876–878
Russell B (1988) Introduction to seismic inversion methods. In: The SEG course notes, series 2
Sen R, Oltz E (2006) Genetic and epigenetic regulation of high gene assembly. Curr Opin Immunol
18(3):237–242
Sen MK, Stoffa PL (2013) Global optimization methods in geophysical inversion. Cambridge
University Press, Cambridge, p 27
Sen MK, Datta-Gupta A, Stoffa PL, Lake LW, Pope GA (1995) Stochastic reservoir modeling using
simulated annealing and genetic algorithm. SPE Formation Eval 10(01):49–56
Shuey R (1985) A simplification of the zoeppritz equations. Geophysics 50(4):609–614
Specht DF (1990) Probabilistic neural networks. Neural Netw 3(1):109–118
Specht DF (1991) A general regression neural network. IEEE Trans Neural Netw 2(6):568–576
Vestergaard PD, Mosegaard K (1991) Inversion of post-stack seismic data using simulated
annealing. Geophys Prospect 39(5):613–624
Xu S, Wang D, Chen F, Zhang Y, Lambare G (2012) Full waveform inversion for reflected seismic
data. In: 74th EAGE conference and exhibition incorporating EUROPEC 2012
Chapter 2
Seismic Data Handling
Abstract Seismic data interpretation is a very crucial step and hence needs special
care to the data. If some ambiguity will remain in the data, the interpretation will
be false and hence loss of time and money. This chapter describes some of the
processing steps of the seismic data that is necessary before proceeding towards the
interpretation of seismic data. The chapter also describes the technique of picking
seismic horizons, wavelet extraction and low-frequency model generation which are
a very important step for the seismic inversion process.
Before moving towards the seismic inversion, seismic and well log data needs to be
prepared. This preparation included seismic data conditioning, Picking of seismic
horizons, wavelet extraction from seismic as well as from well log data, time to depth
correlation and initial low-frequency model generation. These are very important as
well as very crucial steps as the inversion results depend largely on input data. These
steps are briefly described in the following sections.
Before performing the data conditioning, if the data is in the offset domain and not
in the angle domain, it needs to convert into angle gather. A full offset stack seismic
gather is not modeled by the Akki-Richards equation because it is a mix of offsets,
and hence cannot be used for pre-stack inversion (Singleton 2009). Further, this pre-
stack seismic gathers contains significant noises that can mislead the interpretation
of the inverted section. Conditioning of pre-stack seismic data is used to increase
the signal to noise ration prior to pre-stack seismic inversion. Data conditioning has
five major steps i.e. Band-pass Filtering, Muting, Super Gather, Probabilistic Radon
Transform and Trim statics. These processes are briefly explained in the following
sections along with application on real data from the Penobscot field, Nova Scotia,
Canada.
© The Editor(s) (if applicable) and The Author(s), under exclusive licence 19
to Springer Nature Switzerland AG 2020
S. P. Maurya et al., Seismic Inversion Methods: A Practical Approach,
Springer Geophysics, https://doi.org/10.1007/978-3-030-45662-7_2
20 2 Seismic Data Handling
The first step of data conditioning is bandpass filtering. In this process, a band-
pass filter is designed and applied to a subset of seismic data. This filter suppresses
unwanted low and high-frequency signals. It does not affect the wavelet shape (Mari
et al. 1999). This filter is very important because sometimes the receiver record signal
before reaching the main signal that is noise component and the same thing happened
when the main signal is stopped and still the receiver receive some component of the
signal which is again noise component and hence it needs to be separated from the
main signal. Therefore a filter is designed which removes lower frequency compo-
nent that is received before the main signal as well as removes the higher frequency
component which is recorded by receiver after end of main signal (Yilmaz 2001). An
example has been taken from the Penobscot region, Canada. Figure 2.1 shows the
effect of bandpass filtering. Figure 2.1a shows before and Fig. 2.1b show after the
bandpass filtering gathers. In this example, a 5/10–60/70 Hz frequency band-pass
filter is applied as the gather has frequency 10–60 Hz. From the figure, it can be
noticed that the resolution of the signal increases fairly and particularly, the bottom
reflector which was not clear in initial seismic gather becomes very clear in bandpass
filtered seismic gather.
Fig. 2.1 Shows cross-section of seismic data (xline 1001 and 1002), a before, b after bandpass
filtering
2.1 Conditioning of Data 21
Fig. 2.2 Shows cross-section of seismic data (xline 1001 and 1002), a before, b after muting
2.1.2 Muting
Muting is the second step of data conditioning. The muting is applied to the data to
eliminate noise components on the far offset traces. Muting removes bad traces by
setting their amplitude to be zero. It is also applied after NMO correction to remove
the effect of NMO stretch. Apart from this, the muting can also compensate for other
errors, such as refractions and geometric noise (Robinson et al. 1986). Figure 2.2
shows an example of before and after muting gather which is performed on Penobscot
seismic data after application of band-pass filtering. Figure 2.2a shows band-pass
gather before application of the muting process and Fig. 2.2b shows seismic gather
after muting. From these figures, one can notice that the far offset amplitude which
is not part of the main signal is kept to be zero.
The super gather is the third step of data conditioning used to improve the signal to
noise ration of pre-stack gather. The basic idea of super gather is to form average
CDPs by gathering adjacent CDPs and adding them together (Mari et al. 1999;
Singleton 2009). Supper gather also increases the fold of the data and eliminate
22 2 Seismic Data Handling
Fig. 2.3 Shows cross-section of seismic data (xline 1001 and 1002), a) before, b after applied super
gather
noise component. This is a very crucial step of data conditioning and hence applied
by expertise. Figure 2.3 describes the effect of super gather on input muted gathers.
Figure 2.3a shows muted gather which is used as input and Fig. 2.3b shows supper
gather results. From Fig. 2.3, one can notice that the reflector is clearer in super gather
as compared to the muted gather. Sometimes seismic trace amplitude gets distorted
due to external agent particularly near to reflector and hence super gather is used to
correct them (Chopra et al. 2006).
The Parabolic Radon Filter process is the fourth step of data conditioning. The
method eliminates multiple along with radon noise suppression. This is again a very
important step as the pre-stack gather contains a significant amount of noise and
this noisy data can mislead results and unrealistic interpretation of seismic inversion.
The basic process of the radon transform is as follows. The model parameters are set
to identify the long-period multiples or the random noise within the data. After the
model is created, the Radon Transform then subtracts the model of these multiples or
noise from the data, leaving with a data set that is greatly reduced in noise and hence
optimize the traces (Yilmaz 1990; Robinson and Treitel 2000). Figure 2.4 shows an
2.1 Conditioning of Data 23
Fig. 2.4 Shows cross-section of seismic data (xline 1001 and 1002), a before, b after application
of radon transform
example of a radon transform where Fig. 2.4a shows super gather before application
of radon transform and Fig. 2.4b depicts pre-stack gather after radon transform. The
pre-stack gather before radon transform does not reveal a consistent trend whereas
after Radon filtering, the noise, apparently due to multiples, seems to be removed.
It is also observed that the amount of seismic trace amplitude increase with offset.
Apart from this the reflector gets more flatten in the radon transform gathers.
The trim statics is the last step of data conditioning used to make the horizon more
flatten. It also fixes migration move-out problems on pre-stack seismic gathers. The
basic concept of trim statics is that it attempts to determine an optimal time shift
which has been applied to each trace in the gather. The time shift is determined
by cross-correlating each trace with a reference trace to make the input trace better
matches that reference trace (Robinson and Treitel 2000). Usually, the reference trace
is the CDP stacked trace. This step also needs the expertise to apply as in some cases
it degrades the main signal and makes data look like synthetic. An example of trim
statics is shown in Fig. 2.5. Figure 2.5a shows pre-stack gather before trim statics and
Fig. 2.5b depicts pre-stack gather after application of trim statics. From the figure,
24 2 Seismic Data Handling
Fig. 2.5 Shows cross-section of seismic data (xline 1001 and 1002), a before, b after application
of trim statics
one can see that the reflector is better aligned. The trim static was applied to the data
to reduce the residual move out, up to a maximum time shift of 10 ms.
After the conditioning of gathers, the seismic data is ready for use for seismic
inversion. The pre-stack seismic inversion methods are discussed in Chap. 4 using
trim statics gather as input. After the preparation of seismic data, the next step
involves picking seismic horizons which serve as a guide to interpolate petrophysical
parameters between the wells.
on these horizons hence they need to perform with full attention. There are many
approaches available for picking seismic horizons in 3D seismic volumes, viz. Man-
ual picking, Interpolating, Auto picking, etc. These methods are briefly explained in
the following sections. Manual picking of seismic horizons is a very old and efficient
method. Manual picking is performed in interpreting horizons on inlines, cross-lines,
time slices, and traverses by hand and by looking. The interpreter is looking for some
degree of local continuity of the large amplitude in the data along with the local sim-
ilarity of events to be picked while interpreting the horizon manually (Hildebrand
and Landmark Graphics Corp. 1992).
On the other hand, the interpolation technique is somewhat more efficient than
the traditional manual interpretation technique. However, the interpolation technique
assumes that a horizon is locally very smooth and perhaps linear between control
points. If this assumption is violated between control points then the results will
be poor (Faraklioti and Petrou 2004; Pendrel 2006). This is a comparatively fast
technique but in some areas like fault or joint, they need to perform with proper care.
The other horizon picking technique i.e. auto picking is nowadays very popular
because the method takes the least time to pick. These methods are largely used by
software developed for the interpretation of seismic data. The basic concept used in
the auto picking horizon method is quite simple. The interpreter places seed picks
on in-lines and/or cross-lines in the 3-D volume. These seed points are then used as
initial control for the auto picking operation. If it finds such a feature within specified
constraints, it picks that trace and moves on to the next trace (Dorn 1998). Simple
auto pickers allow the user to specify a feature to be tracked, an allowable amplitude
range, and a dip window in which to search. Graphically, how the auto picks works is
explained in Fig. 2.6. The auto tracker looks for similar amplitude or similar feature
in the dip window to pick and if no similar features are found then the auto tracker
stops tracking at that trace (Keskes et al. 1983; Harrigan et al. 1992).
This book uses manual picking methods to pick seismic horizons in 3D seismic
volume. The manual picking is the most efficient methods although the method is
much time taking process. Three horizons are picked in the Blackfoot seismic data
from Alberta, Canada using manual picking technique and are shown in Fig. 2.7.
Fig. 2.7 Picked seismic horizon by manual picking in Blackfoot seismic data, Canada
Two horizons are picked in the pre-stack gather of the Penobscot region, Scotian
Shelf, Canada and are shown in Fig. 2.8. It is also recommended to pick at least two
horizons for getting good results from the seismic inversion.
Fig. 2.8 Picked seismic horizon by manual picking in Penobscot seismic data, Canada
2.3 Seismic Wavelets 27
All seismic inversion methods depend on the form of forward modeling which aims
to generate synthetic gather by the convolution of earth’s reflectivity with source
wavelet. Hence, the wavelet is a very important parameter for any type of seismic
inversion. A seismic wavelet is a mathematical operator used to represent any time-
varying signal into frequency varying signals. During the acquisition of seismic
data, the signal sent by a source is some kind of wavelet. And if someone knows
this wavelet then the seismic inversion can be performed in a more sophisticated
way but unfortunately, it is unknown for most of the cases (Cheng et al. 1996).
Seismic wavelet can be divided into two broad categories i.e. Zero phase wavelet
and minimum phase wavelet.
The zero-phase wavelet has the shortest time duration and the largest peak amplitude
at time zero. The wavelet is symmetrical with a maximum at time zero (non-causal).
The phase of the zero-phase wavelet is zero for all frequency components contained
within the signal. The fact that the energy arrives before time zero is not physically
reliable but the wavelet is useful for increased resolving power and ease of picking
reflection events (peak or trough) (Russell 1988; Cheng et al. 1996). A kind of zero-
phase wavelet is a Ricker wavelet. Figure 2.9 shows three types of ricker wavelets, the
first is the 20 Hz Ricker wavelet (Fig. 2.9a), second is 30 Hz Ricker wavelet (Fig. 2.9b)
and third is 40 Hz Ricker wavelet (Fig. 2.9c) and their amplitude spectrum is shown
at the bottom of wavelet. From the figure, one can be noticed that by increasing the
frequency of wavelet, the wavelet gets sharper whereas the amplitude spectrum gets
broader.
The other wavelet which is popularly used is a minimum phase wavelet and almost
all seismic sources generate minimum phase wavelet. The minimum phase wavelet
has a short time duration and the maximum energy concentration at the start of the
wavelet near time zero. The wavelet amplitude is zero before time zero (causal)
(Russell 1988; Walden and White 1998). An example of a minimum phase wavelet
is presented in Fig. 2.10. Three minimum phase wavelet is shown with frequency
10/10–40/80, 10/20–40/60 and 10/15–60/70 and corresponding amplitude spectrums
are shown at the bottom of the figure. One can easily visualize the relative changes
in time and amplitude spectrum due to the change of frequency. From the figure, it is
noticed that all the wavelet has maximum energy concentration at near to time zero.
28 2 Seismic Data Handling
Extraction of statistical wavelet from the seismic and well log data is a very important
step as in most of the case; the source wavelet is unknown and almost all seismic
inversion method uses wavelet which is extracted either from seismic data or well
log data. Although, the wavelet extraction is a complex process, and is currently
an area of active research. There are numerous wavelet extraction methods. The
wavelet extraction in frequency domain consists of two parts, first determine the
amplitude spectrum and second determines the phase spectrum (Swisi 2009). In these
two processes, the determination of the phase spectrum is a very difficult process
and includes a major source of error in seismic inversion. There are three major
categories of wavelet extraction. The first is purely deterministic wavelet extraction
which includes measuring the wavelet directly using surface receivers and other
techniques, such as marine signatures or VSP analysis. The second approach is purely
statistical wavelet extraction which includes determination of the wavelet from the
seismic data alone. These procedures tend to have difficulty in determining the phase
spectrum reliably. And the third approach is the extraction of wavelet using well log
which uses well-log information in addition to the seismic data. In theory, this could
2.3 Seismic Wavelets 29
Fig. 2.10 Minimum phase wavelet of a 10/10–40/80 Hz, b 10/20–40/60 Hz and c 10/15–60/70 Hz
frequency
provide exact phase information at the well location (Russell 1988; Swisi 2009). The
problem is that this method depends critically on a good tie between the log and the
seismic data. The statistical wavelet extraction and extraction of wavelet using well
is most widely used and is described briefly in the following section.
The statistical wavelet extraction technique utilizes seismic data alone to extract
wavelet. The mathematical procedure used here is autocorrelation. The phase spec-
trum of the wavelet is not calculated by this method and must be supplied as a separate
parameter by the user (Edgar and Van der Baan 2011). The procedure used to extract
the statistical wavelet is as following.
• Extract the analysis window in the seismic data.
• Taper the start and end of the window with a taper length.
• Calculate the autocorrelation of the data window. The length of the autocorrelation
should be equal to half of the desired wavelet length.
30 2 Seismic Data Handling
The other type of wavelet extraction involves the use of well log data. This can
be performed in two ways. The first is to use well log data to determine both; the
full amplitude and phase spectrum of the wavelet and the second, it uses only well
log to determine the constant phase of the wavelet in a combination of the statistical
procedure described in the above section. In both approaches, the first is very popular
and most widely used one (Swisi 2009). The process of wavelet extraction from the
well log data is as follows.
• Extract the sonic, density, and seismic data analysis windows.
• Multiply sonic and density to get impedance and then transform it into reflectivity
series.
• Apply Taper on both the seismic data and the reflectivity series at the start and
end of the window.
• Next, calculate the least-squares shaping filter, W, which solves the following
equation:
Trace(S) = W ∗ Reflectivity(R)
An example of wavelet extraction is shown in Fig. 2.11 for Penobscot data, Scotian
Shelf, Canada. The figure compares three wavelets, one is extracted from well log
data, second is extracted from the raw seismic section and the third is extracted from
trim static gathers. The corresponding amplitude spectrum is shown at the bottom of
the same figure.
Further, the extracted wavelet from well log of the Blackfoot region is shown in
Fig. 2.12. The area has 13 well logs hence thirteen wavelets are extracted. Wavelet
2.3 Seismic Wavelets 31
Fig. 2.11 Extracted wavelet from a well logs, b raw seismic and c trim gathers
Fig. 2.12 Extracted wavelet from a well 29-08, b combine all wavelets together and c Blackfoot
seismic
32 2 Seismic Data Handling
from well 29-08 is shown in Fig. 2.12a. Figure 2.12b shows the average of all wavelets
extracted from all thirteen wells. Wavelet extracted from the Blackfoot post-stack
seismic section is presented in Fig. 2.12c.
Another input used by the seismic inversion method is a low frequency model. The
seismic data is band limited and contains the only frequency of ranges 10-80 Hz
and hence have a lack of low and high-frequency component. These low-frequency
components are very important in terms of interpretation as this low frequency con-
tains information of deeper reflector. During the inversion of seismic data, the low-
frequency component is supplied by the low-frequency model to get a broadband
frequency spectrum of inverted results. The low-frequency model is obtained by the
lateral interpolation and extrapolation of well log property between the wells using
a seismic horizon as a guide (Swisi 2009; Maurya and Singh 2018a, b). This inter-
polation and extrapolation are usually performed using a variety of mathematical
methods such as weighted inverse-distance, spline, kriging, and cokriging. Some of
these methods are explained briefly in the following sections.
2.4.1 Kriging
N
z0∗ = wi zi + w0 (2.1)
i=1
At each location, the estimation error is the difference between the estimated
value, z0∗ , and the true value, z0 .
2.4.2 Cokriging
We often have two separate measurements in petroleum exploration: well log curves
or key samples and a set of seismic information. The physical property of the sub-
surface, such as the acoustic and shear velocities, density, neutron porosity, etc., can
be measured straight at the well places. We can interpolate and map these read-
ings between the reservoirs using the kriging technique. However, we would like
to include the seismic information in the mapping method as they can provide the
exploration region with very excellent spatial coverage (Todorov 2000; Chambers
et al. 2000; Ahmadi and Sedghamiz 2008). Hence Cokriging is a technique when
uses both well and seismic information and predicts values away from the boreholes.
The ordinary cokriging estimator Z0∗ is written as follows.
N
M
Z0∗ = wi Zi + vj Xj (2.2)
i=1 j=1
where Z0∗ is the estimate at the grid node; Zi is the regionalized variable at a given
location, with the same units as for the regionalized variable; wi is the undetermined
weight assigned to the primary sample Zi and varies between 0 and 100%; Xj is
the secondary regionalized variable that is co-located with the primary regionalized
variable Zi , with the same units as for the secondary regionalized variable; and vj is
the undetermined weight assigned to Xj and varies between 0 and 100%.
should be remembered that this technique considers the concept of adjacent radius
and the associated strength to the distance reverse function as serious problems. A
state in which there are sufficient sample points (at least 14 points) with an appropriate
dispersion at local scale concentrations will use this technique. The primary variable
influencing inverse distance interpolator precision is the energy parameter value p
(Burrough and McDonnell 1998). Also important to the accuracy of the outcomes
are the size of the neighborhood and the number of neighbors.
N
zi .di−n
Z0 = i=1 (2.3)
N
i=1 di−n
where, Z0 is the estimation value of variable z in point I, zi is the sample value in point
I, di is the distance of the sample point to the estimated point, N is the coefficient
that determines weight based on a distance and n is the total number of predictions
for each validation case.
The low-frequency acoustic impedance (AI) model is created using the following
steps.
• Calculate the AI at well locations by multiplication of well log velocity and
density.
• Pick seismic horizons in 3D seismic volume to guide the interpolation between the
wells. These seismic horizons also provide structural information for the model.
• Use interpolation along the seismic horizons and between the well locations to
obtain the initial low AI model.
• Use the high cut filter to remove higher frequency from the model as the model
contains the only low-frequency component.
A similar approach is applied to estimate a model of velocity, density, VP /VS etc.
In this way, one can generate a low-frequency model that is used as prior information
during the seismic inversion methods to supply the low-frequency component.
The examples of low-frequency models are shown in Figs. 2.13, 2.14 and 2.15.
Figure 2.13 shows the low frequency AI model of the Blackfoot field, Canada. This
model is created using 13 well logs of the area and two picked seismic horizons.
Apart from this, a high cut filter of 10–15 Hz is used to remove higher frequency
content. Figure 2.14 shows the P-impedance model for Penobscot data, Scotian Shelf,
Canada. This low-frequency model is generated by using one well (L-30) and two
picked seismic horizons. A high cut frequency filter of 12–15 Hz is applied to remove
higher frequency from the model. Similarly, S-impedance model is generated and
shown in Fig. 2.15. The Density and VP /VS model is shown in Figs. 2.16 and 2.17
respectively.
2.4 Low-Frequency Model 35
Fig. 2.15 Low-frequency S-impedance model of Penobscot field, Nova Scotia, Canada
36 2 Seismic Data Handling
Fig. 2.16 Low-frequency density model of Penobscot field, Nova Scotia, Canada
Fig. 2.17 Low frequency VP /VS model of Penobscot field, Nova Scotia, Canada
References
Ahmadi SH, Sedghamiz A (2008) Application and evaluation of kriging and cokriging methods on
groundwater depth mapping. Environ Monit Assess 138(1–3):357–368
Burrough PA, McDonnell RA (1998) Creating continuous surfaces from point data. Principles of
Geographic Information Systems. Oxford University Press, Oxford, UK
Chambers RL, Yarus JM (2002) Quantitative use of seismic attributes for reservoir characterization.
CSEG Recorder 27(6):14–25
Chambers RL, Yarus JM, Hird KB (2000) Petroleum geostatistics for nongeostatisticians: part 1.
Lead Edge 19(5):474–479
Chang YH, Scrimshaw MD, Emmerson RHC, Lester JN (1998) Geostatistical analysis of sampling
uncertainty at the Tollesbury managed retreat site in Blackwater Estuary, Essex, UK: kriging and
cokriging approach to minimise sampling density. Sci Total Environ 221(1):43–57
Cheng Q, Chen R, Li TH (1996) Simultaneous wavelet estimation and deconvolution of reflection
seismic signals. IEEE Trans Geosci Remote Sens 34(2):377–384
References 37
Chopra S, Castagna J, Portniaguine O (2006) Seismic resolution and thin-bed reflectivity inversion.
CSEG Recorder 31(1):19–25
Dorn GA (1998) Modern 3-D seismic interpretation. Lead Edge 17(9):1262
Edgar JA, Van der Baan M (2011) How reliable is statistical wavelet estimation? Geophysics
76(4):V59–V68
Faraklioti M, Petrou M (2004) Horizon picking in 3D seismic data volumes. Mach Vis Appl
15(4):216–219
Harrigan E, Kroh JR, Sandham WA, Durrani TS (1992) Seismic horizon picking using an artificial
neural network. IEEE Int Conf Acoust Speech Signal Process 3:105–108
Hildebrand HA, Landmark Graphics Corp. (1992) Method for finding horizons in 3D seismic data.
U.S. Patent, 5:153,858
Keskes N, Zaccagnino P, Rether D, Mermey P (1983) Automatic extraction of 3-d seismic horizons.
In SEG technical program expanded abstracts. Society of Exploration Geophysicists, pp 557–559
Mari J-L, Glangeaud F, Coppens F, Painter D (1999) Signal processing for geologists and
geophysicists. Technip, Paris, p 480
Maurya SP, Singh KH (2018a) Qualitative and quantitative comparison of geostatistical techniques
of porosity prediction from the seismic and logging data: a case study from the blackfoot field,
Alberta, Canada. Mar Geophys Res 40(1):51–71
Maurya SP, Singh NP (2018b) Application of LP and ML sparse spike inversion with probabilistic
neural network to classify reservoir facies distribution—a case study from the blackfoot field,
Canada. J Appl Geophys Elsevier 159(2018):511–521
Oliver DS, Reynolds AC, Liu N (2008) Inverse theory for petroleum reservoir characterization and
history matching. Cambridge University Press
Pendrel J (2006) Seismic inversion-still the best tool for reservoir characterization. CSEG Recorder
26(1):5–12
Robinson EA, Durrani TS, Peardon LG (1986) Geophysical signal processing. Prentice-Hall
International
Robinson E, Treitel S (2000) Geophysical signal analysis. Society of exploration geophysicists
Russell B (1988) Introduction to seismic inversion methods. In: The SEG course notes series 2
Singleton S (2009) The effects of seismic data conditioning on prestack simultaneous impedance
inversion. Lead Edge 28(7):772–781
Swisi AA (2009) Post-and pre-stack attribute analysis and inversion of Blackfoot 3D seismic dataset.
Doctoral dissertation, University of Saskatchewan
Todorov TI (2000) Integration of 3C-3D seismic data and well logs for rock property estimation.
In: M.Sc. thesis, University of Calgary
Walden AT, White RE (1998) Seismic wavelet estimation: a frequency-domain solution to a
geophysical noisy input-output problem. IEEE Trans Geosci Remote Sens 36(1):287–297
Yilmaz O (1990) Seismic data processing. Society of exploration geophysicists
Yilmaz O (2001) Seismic data analysis. In: Society of exploration geophysicists, vol 1. Tulsa, OK
Chapter 3
Post-stack Seismic Inversion
Abstract Post-stack seismic inversion utilizes post-stack seismic data along with
well log data to estimate acoustic impedance. Post-stack seismic inversion is very
fast compared to other pre-stack seismic inversion methods and provides a high-
resolution subsurface image. This chapter discusses several types of post-stack seis-
mic inversion methods namely model-based inversion, colored inversion, sparse
spike inversion, and band-limited inversion. The chapter also includes the synthetic
as well as real data examples of above seismic inversion methods.
3.1 Introduction
The seismic inversion techniques outlined here utilize the post-stack seismic data
for the estimation of the physical properties of the subsurface. The efficacy of the
inversion process relies on the effective transformation of seismic amplitudes into
impedance values. Adequate precaution needs to be taken to preserve amplitude
during inversion as it ensures that the measured amplitude variations are related to
the geological bodies (Vestergaard and Mosegaard 1991). Thus, the multiples from
the seismic data need to be processed and removed. This leads to a high signal
to noise ratio and zero offset migrated section devoid of any numerical artifacts.
The seismic data is band-limited, and therefore the lack of low frequencies prevents
the transformed impedance traces from having the elementary velocity structure
important to make geological interpretation (Clochard et al. 2009). In addition, the
generated impedance volume has poor resolution incapable of identifying the thin
layers. These aspects are therefore crucial and must be taken care during seismic
inversion.
Post-stack techniques of inversion usually refer to the different workflows used
to convert stacked seismic information into parameters of quantitative rock physics.
The result is generally acoustic impedance from post-stack inversion, while pre-stack
inversion can result in acoustic impedance as well as shear impedance. Sparse spike
inversion, model-based inversion, recursive inversion, and colored inversion are the
techniques of post-stack inversion. All of these techniques are grouped together as
© The Editor(s) (if applicable) and The Author(s), under exclusive licence 39
to Springer Nature Switzerland AG 2020
S. P. Maurya et al., Seismic Inversion Methods: A Practical Approach,
Springer Geophysics, https://doi.org/10.1007/978-3-030-45662-7_3
40 3 Post-stack Seismic Inversion
There are many statistical parameters like correlation coefficient, RMS error, syn-
thetic relative error, peak signal to noise ratio, mean absolute error, sum absolute
error, etc. to compare two signals quantitatively. These parameters can be defined as
follows.
One of the most commonly used formulas in statistics is Pearson’s correlation
coefficient. Correlation between sets of data is a measure of how well they are related.
The most common measure of correlation is the Pearson Correlation.
n xy − x y
r = (3.1)
2 2 2 2
n x − x n y − y
where n is the number of data points in the signal, x is data points of the first signal
and y is data points of the second signal.
The next statistical parameter is the peak signal to noise ratio (PSNR) which can
be written mathematically as follows.
MAX f
PSNR = 20 log10 √ (3.2)
MSE
where,
m−1
MSE = |x(i) − y(i)|2 (3.3)
n i=0
and f represents the matrix data of the original signal, m and n represents number
of rows and columns in the input signal.
Further, the root means square error can be defined as follows.
n
i=1 (x i − yi )2
RMSE = (3.4)
n
One first needs to determine absolute error to calculate the relative error. The
relative error expresses how large the absolute error is in comparison to the total
size of the object being measured. The relative error is expressed as a fraction or
multiplied by 100 and expressed as a percentage.
Absolute Error
Relative Error = (3.7)
Known Value
The above discussed statistical parameters are used to calculate the quantitative
difference between the inverted and the original signals.
Band-limited inversion (BLI) is the oldest type of post-stack seismic inversion meth-
ods that utilize post-stack seismic data as input and transforms them into acous-
tic impedance. The acoustic impedance of the subsurface strengthens to seismic
data interpretation and finding prospective zone. Further, this impedance can be
transformed into many petrophysical parameters using some prediction tools like
the geostatistical approach. The mathematical equation of band-limited impedance
method can be derived from the relationship between the seismic trace and seismic
impedance (Ferguson and Margrav 1996; Maurya and Singh 2017). The impedance
is related to the velocity and density as follows.
If one knows the impedance (Z), then the earth’s reflectivity (r) can be estimated
using the following relationship.
Z j+1 − Z j
rj = (3.9)
Z j+1 + Z j
Equation 3.9 is known as normal incidence reflection formula. In Eq. 3.9, the
Z j is the seismic impedance of jth layer, and r j is seismic reflectivity of jth and
( j + 1)th interface. The above Eq. (3.9) can write as follows.
42 3 Post-stack Seismic Inversion
Or,
1 + rj
Z j+1 = Z j (3.11)
1 − rj
Similarly,
1 + r3 1 + r1 1 + r2 1 + r3
Z4 = Z3 = Z1 (3.12)
1 − r3 1 − r1 1 − r2 1 − r3
And hence one can write Eq. 3.5 for nth value as follows.
1 + r1 1 + r2 1 + rn−1
Zn = Z1 ... (3.13)
1 − r1 1 − r2 1 − rn−1
From Eq. 3.13, it is noticed that the acoustic impedance for the first layer needs to
be estimated from a continuous layer above the target area and then the impedance
of another layer can be estimated. In this method, the impedance for the jth layer can
thus be calculated as follows.
j
1 + rk
Z j+1 = Z j (3.14)
k=1
1 − rk
Now, the Eq. (3.14) is divided by Z j and then taking logarithm on both sides, we
get the equation.
j
Z j+1
1 + rk
j
ln ≤ ln ≈2 rk (3.15)
Zj k=1
1 − rk k=1
Thereafter, Eq. 3.15 can be solved for Z j+1 by considering small values of r we
get the formulae.
3.3 Band-limited Inversion 43
j
Z j+1 = Z 1 exp 2 rk (3.16)
k=1
Further, model the seismic trace as scaled reflectivity then one can write Sk =
2rk /γ , and hence Eq. 3.16 can be written as follows.
j
Z j+1 = Z 1 exp γ Sk (3.17)
k=1
Equation 3.17 integrates the seismic trace and then exponentiates the result to
provide an impedance trace. Band-limited impedance method uses this equation to
invert the seismic trace (Waters and Waters 1987; Ferguson and Margrave 1996;
Maurya and Sarkar 2016).
From Eq. 3.17, one can notice that the seismic trace gets inverted for impedance
directly although seismic data have band-limited frequency and hence inverted
impedance also have band-limited nature. To get broadband spectrum, a low-
frequency initial impedance model is added to the inverted results. The generation of
low-frequency impedance model is already discussed in Chap. 2. An important lim-
itation of the band-limited impedance inversion is that seismic data must be in zero
phases. The seismic pre-stack data could be transformed into zero phases and then
can be utilized for band-limited impedance inversion (Maurya and Singh 2015b).
Figure 3.1 depicts the flowchart of band-limited impedance inversion methods. The
figure demonstrates that the inversion method utilizes seismic and well log data as
input and acoustic impedance is estimated in the output. This method is very fast and
still in use by many oil and gas companies to understand subsurface features.
The other example of BLI to real field data is from the Blackfoot region, Alberta
Canada which is used to estimate subsurface acoustic impedance.
Post stack seismic data along with 13 well logs are used as input for BLI. The
method is applied in two steps, in the first step, composite traces near to well locations
are extracted and the algorithm is applied to these traces to get AI, and in the second
step, the entire seismic volume is inverted. This is standard procedure of performing
inversion as these methods are time-consuming and hence one needs to optimize all
the parameters of the inversion. The other benefit of using this is to cross-validate
inversion results for real data case as one knows the AI from the well log data
and hence inverted AI can be cross-checked by the original AI. Figure 3.3 shows
comparison of inverted acoustic impedance with original AI from well log along
with initial guess model for 12 wells. From Fig. 3.3, it is noticed that the inverted
acoustic impedance is following well log curve very nicely for almost all wells.
It is also noticed that some well log curves have not so nicely match of phase but
the overall trend is matching nicely.
Thereafter, quality checks of inversion results need to be performed and for that,
the cross plot of inverted and original acoustic impedance is generated which is
presented in Fig. 3.4. The figure demonstrates the cross plot of all 12 wells. The
scatter points lie near to the best fit line for all the wells which indicate the good
inversion results. A small cluster of scatter data lies away from the best fit line which
is due to the different phases of the inverted curve in some regions. These analyses
show that the method works nicely for composite trace and generates good results
for the inversion of seismic volume.
46 3 Post-stack Seismic Inversion
Fig. 3.3 Inversion analyses plot for 12 wells. The red curve shows inverted AI, the blue curve
shows the original AI and the black curve shows the initial AI model
Fig. 3.4 Crossplot of modeled and original impedance for all thirteen wells
3.3 Band-limited Inversion 47
The correlation coefficient and RMS errors are estimated between inverted and
original impedance from all the 13 wells and results are shown in Fig. 3.5. These
values indicate the quantitative capability of BLI methods. The average correla-
tion is estimated to be 0.89 and the average RMS error for all the thirteen wells is
1123 m/s*g/cc. These values justify the good performance of the algorithm.
Further, full volume inversion can be performed for Blackfoot 3D datasets to
estimate acoustic impedance volume in the inter-well region. A cross-section of
the inverted impedance (inline 28) is shown in Fig. 3.6. The variations of AI in
the subsurface vertically and horizontally are clearly visible. The figure shows that
the area have impedance variation from 5000 to 20,000 m/s*g/cc. The figure shows
Fig. 3.5 Variation of correlation coefficient (top) and RMS error for all thirteen wells of the
Blackfoot field, Canada
Fig. 3.6 Cross-section of acoustic impedance estimated using the BLI method
48 3 Post-stack Seismic Inversion
Fig. 3.7 Comparison of the amplitude spectrum of Blackfoot seismic and inverted synthetics from
BLI methods
Fig. 3.8 Process of estimation of colored operator, a using a set of wells from the area, the amplitude
spectra of the acoustic impedance for all the wells are plotted on log-log scale (I), b using a set of
seismic traces from around the wells, the average seismic spectrum is calculated (S), c operator in
time domain (I/S), and d operator in frequency domain
50 3 Post-stack Seismic Inversion
(c) The final spectrum needs to combine −90° phase shift to create the desired
operator in the time domain (Fig. 3.8c).
(d) The operator in the frequency domain is used to transform data which can be
achieved from Fourier transform (Fig. 3.8d).
Colored inversion is fast and suitable for application to 3-D datasets (Ansari
2014). Generally, traditional inversion methods are time-consuming, expensive,
require specialists and are not performed routinely by the interpretation Geophysi-
cist, whereas CI is rapid, easy to use, inexpensive, robust and does not require expert
users (Maurya and Sarkar 2016). Figure 3.9 depicts flowchart of colored inversion
methods. The figure shows that the colored inversion utilizes seismic and well log
data and generates an operator. Further, this operator convolved to the seismic trace
to get acoustic impedance as output. The method is very simple and very fast. But
some author argues (Russell 1988; Maurya and Singh 2018, 2019; Maurya et al.
2019) that the colored inversion gives an average variation of impedance and hence
one cannot use for the quantitative interpretation.
Real data example from the Blackfoot data, Alberta, Canada is presented here. 3D
post-stack seismic data along with 13 well logs are used as input for colored inversion
methods. The method is first tested on the composite trace and then applied to the
3D seismic volume. The area has thirteen wells and hence thirteen seismic trace is
extracted near to well location and colored inversion is applied one by one to all
traces. Figure 3.10 shows a comparison of inverted AI with original AI from the well
logs along with the initial AI guess model. The area has 13 wells but the results are
shown for 12 wells for simplicity. From the figure, one can notice that all the inverted
curves are matching with the well log curves very nicely. This indicates the good
performance of the algorithm.
For quality check of the inversion results, a cross plot of inverted impedance with
well log impedance is generated for all the wells and is shown in Fig. 3.11. Further, a
best fit straight line is fitted to the data and the variation of scatter points away from
the best fit line is noticed. The figure indicates that most of the scatter points lie near
to the best fit line which presents that the inverted data points are very close to the
original data points and hence indicate good performance of the CI.
For the quantitative comparison of inversion results for composite trace, some
numerical parameters are extracted and shown in Fig. 3.13. The other parameters are
given in Table 3.2. Figure 3.12 shows two histogram plots, one corresponds to the
correlation coefficient (Top) and the other corresponds to the RMS error (Bottom).
From the figure one can notices that the correlation coefficient varies from 0.86 to
0.89 with average value of 0.875 and RMS error varies from 700 to 1400 m/s*g/cc
with average value of 1200 m/s*g/cc.
Further, 3D inversion is performed and acoustic impedance in the inter-well region
is estimated. The parameters used for the colored inversion analysis are as follows.
• Colored inversion was initialized by 13 data pairs.
• Processing sample rate = 2 ms
• Threshold for inverse = 20
• Operator Length = 200
• Taper Length = 50
• Regression Line: Intercept = 5.92741, Gradient = −0.733055.
The cross-section of inverted acoustic impedance at inline 27 is shown in Fig. 3.13.
The section shows variations of AI horizontally as well as vertically in the subsur-
face and depicts very high resolution as compared with the seismic section. The
area have impedance variation from 5000 to 17,000 m/s*g/cc. The main benefit of
52 3 Post-stack Seismic Inversion
Fig. 3.10 Inversion analysis for 12 wells estimated using colored inversion methods
interpretation from the inverted AI section is that the information is given as lay-
ered whereas the seismic section gives information at the interface only. The well
log impedance is also plotted above the inverted AI section for cross verification
of results. Figure 3.13 indicates good matching between the well log and inverted
AI section. From the inverted section, one can notice that the impedance contrast
increases with increasing depth which is obvious as velocity and density increase
with increasing depth. The low impedance anomaly is also clearly visible in between
1060 and 1075 ms time intervals. One can also notice that the interpretation is easier
from the inverted section as compared with the seismic sections. This is the reason
why exploration companies use more frequently these inversion techniques. How-
ever, one seeks to find a prospective zone then the colored inversion method is the
best method as the method is less time consuming and requires less expertise.
3.4 Colored Inversion (CI) 53
Table 3.2 Quantitative comparison of seismic inversion results (AI stands for acoustic impedance,
CC stands for colored inversion, PSNR stands for peak signal to noise ratio, MAE stands for mean
absolute error and SAE stands for sum absolute error)
Colored inversion
Properties Reservoir zone CC RMS error PSNR MAE SAE
AI 6.5 × 103 –8.0 × 103 0.84 1.5 × 103 15.2 1 × 103 8 × 104
Amp. spectrum NA 0.93 0.1 70.4 0.0 10.0
Fig. 3.12 Variation of correlation (top) and RMS error (bottom) for inverted results estimated using
colored inversion methods
Further, to check frequency content in the inverted seismic, one can compare the
amplitude spectrum of input seismic i.e. Blackfoot seismic and inverted synthetic
seismic. The comparison is presented in Fig. 3.14. From the figure, one can notice
that the amplitude spectrum of inverted synthetic seismic is not matching with the
54 3 Post-stack Seismic Inversion
Fig. 3.13 Cross-section of inverted acoustic impedance estimated using colored inversion methods
Fig. 3.14 Comparison of amplitude spectrum between Blackfoot seismic and inverted synthetic
data
input seismic trace although the range of both spectra is the same. This comparison
suggests that the colored inversion method distort frequency component completely
during the process and hence dependency of the interpretation on inverted results
gets restricted to qualitative only. The colored inversion method distorts frequency
more as compared to the band-limited inversion methods. This is the drawback of
seismic colored inversion methods.
In the above section, we have presented a straightforward description and applica-
tion of colored inversion (CI). The method is called a robust process in the literature,
but it is somewhat sensitive to the chosen frequency range. Notwithstanding all this,
the process of colored inversion is simple and fast and yields informative subsurface
3.4 Colored Inversion (CI) 55
images for the interpreters. Now, we will look at the other popular post-stack seismic
inversion methods i.e. Model-based inversion.
where, S(t) is seismic trace, W (t) is wavelet, r (t) is earth’s reflectivity, n(t) is noise
component and some time it is assumed to be zero for simplicity and ∗ indicates the
convolution process. If the noise in the data is uncorrelated with the seismic signal,
the trace can be solved for the earth’s reflectivity function (Stull 1973; Russell 1988).
The main benefit of model-based inversion is that the methods involve the updating
of the geological model rather direct inversion of seismic data. This is beneficial
because the impurities in the seismic data are not incorporating in the inversion
results whereas the other inversion method transforms the main signal which has
transmission losses, spherical divergence, unwanted signals, etc. (Ferguson 1996).
The method can be classified into two parts by questioning. The first question
is; what is the fundamental mathematical relationship between the model and the
seismic data and the second question is; how to update the guessing model. To answer
this question, we have considered two different approaches namely, the Generalized
Linear Inversion (GLI) and the Seismic Lithologic (SLIM) method. These methods
are briefly discussed in the following sections.
The model-based inversion methods use the principle of generalized linear inversion
method. The basic idea behind this method is, given a set of geophysical observations,
the GLI method will find the geological subsurface model (AI in this case) which
best fits the observations in the least square sense. The generalized linear inversion
56 3 Post-stack Seismic Inversion
method is developed by Coke and Schneider (1983) and it is now very common in
use. This can be expressed mathematically as follows.
If one can write observations in the form of vector then the parameter of the k
model can be written in vector form as follows.
M = (m 1 , m 2 , . . . , m k )T (3.19)
and observation with n data points can be written in vector form as follows.
T = (t1 , t2 , . . . , tn )T (3.20)
Now, one can express the relationship between the model parameters and
observations in the functional form as follows.
If one is able to derive the functional relationship between the model parameters
and observations, any set of model parameters will be produced as output. The GLI
eliminates the need for trial and error by analyzing the error between the model
parameters and the observations and then perturbing the model parameters in such a
way as to produce an output which will produce less error (Yilmaz 2001). In this way,
one may iterate towards a solution. Mathematically, the observation can be written
by Taylor’s series expansion of the forward model.
∂ F(M0 ) ∂ 2 F(M0 )
F(M) = F(M0 ) + M + M 2 + · · · (3.22)
∂ M0 ∂ M02
F = GM (3.25)
In Eq. 3.25, the G is matrix derivatives which have n number of rows and k
number of columns. The tendency of this error is decreasing in a roughly exponential
manner with iteration, but the condition is that the initial guess model lies within the
range of convergence. The iteration is continued until the error drops below some
predetermined level. One can write solution of Eq. 3.25 as follows.
M = G −1 F (3.26)
However, one must derive the functional relationship between model parameters
and observations which is necessary to relate to each other. The simplest solution is
the convolution model which can be expressed as follows.
Equation 3.28 is not directly used by Coke and Schneider (1983) as this equa-
tion not include multiples, transmission losses, etc. in its implementation. Coke
and Schneider (1983) use a modified version of Eq. 3.28. The main advantage of
incorporating multiples and transmission loss in the solution is that they are not
included in the model parameters, although they are modeled in computing the seis-
mic response (Russell 1988). This is big advantage of GLI and hence model-based
inversion method over recursive methods discussed above since those methods incor-
porate the multiples and transmission loss into the solution if the data contains it.
Figure 3.16 displays flowchart of the GLI methods.
The workflow of model-based inversion technique is as follows (Ferguson and
Margrave 1996):
1. Calculate the acoustic impedance at well locations using the well log data and
Pick horizons in the seismic section to control the interpolation and to provide
structural information for the model between the wells in the area (Bosch et al.
2010).
2. Use the kriging interpolation technique along with the picked seismic horizons
to obtain the initial acoustic impedance model (Maurya and Singh 2015).
58 3 Post-stack Seismic Inversion
3. Extract statistical wavelet from the seismic section and convolve it with the Earth
reflectivity (generated from the initial impedance model) to obtain synthetic
seismic trace. This synthetic trace is different from the observed seismic trace.
4. The least-squares optimization is performed for minimizing the difference
between the real and modeled reflectivity section. This is achieved by analyzing
the misfit between the synthetic trace and the real trace and modifying the block
size and the amplitude to reduce the error (Brossier et al. 2015; Maurya and Singh
2017).
The Seismic Lithologic methods are described in Cooke and Schneider (1983) but
the developed methods are not commercially used as it is not promising for getting
always great results. However, Western Geophysical developed a similar technique of
Seismic Lithologic Modeling (SLIM) which is commercially used as part of model-
based inversion. Although the Western Geophysical’s developed algorithm has not
been fully published and hence the mathematics behind it is not in the public domain.
The SLIM technique involves a model being perturbed rather than a seismic segment
being directly reversed.
The flowchart of SLIM is presented in Fig. 3.15. The basic idea behind the tech-
nique is that the original geological model is produced and compared to a seismic
section, as in the GLI technique. The model is described at different control points
along the line as a sequence of layers of variable velocity, density, and thickness.
Also, either the seismic wavelet is provided or it is estimated. The synthetic model
is then compared to the seismic data and the amount of the least squared error is
calculated. The model is perturbed in such a manner that the error is reduced, and
until the convergence, the method is continued.
The user has complete power over the limitations and can integrate any source
of geological data. The main benefit of this technique over conventional recursive
techniques is that noise is not integrated into the process.
Here two examples are presented to examine the capability of MBI to estimate acous-
tic impedance, first is from synthetic data and second is real data from the Blackfoot
field, Canada. To generate a synthetic seismogram, we have utilized velocity and
density as a geological model given in Table 3.1. And then, synthetic seismograms
are generated over the 7 layer earth model by convolving 30 Hz ricker wavelet to the
reflectivity. Figure 3.16 shows inversion results for 7 layer earth model. Track 1 of
Fig. 3.16 shows the assumed geological model, track 2 shows calculated reflectivity
3.5 Model-based Inversion (MBI) 59
from given velocity and density combination, track 3 shows synthetic seismogram
generated by convolution of 30 Hz wavelet with reflectivity series and track 4 shows
a comparison of inverted impedance with true impedance. From the figure, one can
notice that the inverted impedance is following the trend of the true impedance very
nicely. The correlation coefficient between inverted and true impedance is 0.99 and
RMS error is 0.023 which indicates the good performance of the algorithm.
The real data example is applied in two steps, in the first step, composite traces near
to well locations are extracted and MBI method is applied and in the second step, it
is applied to the 3D data volume to get AI in the inter-well region. Figure 3.17 shows
Fig. 3.17 Inversion analyses for composite trace using the model-based inversion method
3.5 Model-based Inversion (MBI) 61
inversion results for composite trace. The figure shows 12 subplots each corresponds
to 12 different well logs. Each subplot has three curves, the black curves correspond
for the initial guess model, the blue curve indicates true AI from logs and the red
curve indicates inverted AI. The analysis of composite trace shows that the inverted
results are following the trend of original AI very nicely for all wells. The analysis
of composite trace is very important to cross verify inversion results as one knows
true AI from the well log data.
A cross plot between inverted and original AI has been generated for all the wells
and is shown in Fig. 3.18. Next, a best fit straight line is fitted to the data. The
deviation from this best fit line shows the quality of the inversion results. As if all
the points lying on the fitted line indicates that the inverted results are the same as
original data and if shows small deviation means small difference in between original
and inverted data points. From Fig. 3.18, one can notice that most of the data points
lie near to the best fit line and demonstrate good results.
The variation of correlation coefficient and RMS errors are shown in Fig. 3.19.
The correlation coefficient and RMS errors are very fundamental parameters to clas-
sify the difference between two similar signals. The figure shows the variation of
correlation coefficient (top) and RMS error (bottom) with well logs. From the figure,
one can notice that the correlation coefficient is very high and RMS error is very
low which is obvious. The correlation coefficients vary from 0.997 to 0.995 and the
average value is 0.99. On the other hand, the RMS error shows variation from 700
to 1100 (m/s)*(g/cc) with average value of 920 (m/s)*(g/cc). The other statistical
parameters are shown in Table 3.3. These values indicate the high performance of
the MBI algorithm.
We have seen that the inverted results are very accurate and give satisfactory
results for the composite trace and now we are in a position to apply MBI methods to
the real 3D seismic volume to estimate acoustic impedance volume. The parameters
used for model-based inversion are as follows.
62 3 Post-stack Seismic Inversion
Fig. 3.19 Variation of correlation (top) and RMS error (bottom) for results estimated MBI
Table 3.3 Quantitative comparison of seismic inversion results (AI stands for acoustic impedance,
CC stands for colored inversion, PSNR stands for peak signal to noise ratio, MAE stands for mean
absolute error and SAE stands for sum absolute error)
Model-based inversion
Properties Reservoir zone CC RMS Error PSNR MAE SAE
AI 6.3 × 103 – 8.0 0.86 2 × 103 17.9 1.6 × 103 1.1 × 105
× 103
Amp. spectrum NA 0.99 0.0 80.1 0.0 2.7
Fig. 3.20 Cross-section of inverted acoustic impedance (inline 41) estimated using the MBI
technique
Fig. 3.21 Comparison of the Amplitude spectrum of inverted synthetic and Blackfoot seismic
before application of MBI methods
amplitude spectrum of inverted synthetic is almost matching with the amplitude spec-
trum of Blackfoot seismic. This is the biggest advantage of using the MBI method
over other inversion methods. We have seen that the BLI and CI distort frequency
content but here the MBI methods have not much destroyed the frequency of the
input signal. This is obvious because the MBI method does not use seismic data
directly as input whereas the other inversion method uses seismic data directly as
input and mathematical process distort frequency content.
The model-based inversion method has a comparatively complex mathematical
background but can be easily used to the real field data. The MBI method is the
most popular post-stack seismic inversion method and the results derived from it
can be used for the quantitative interpretation. The non-uniqueness is a problem of
64 3 Post-stack Seismic Inversion
the entire inversion methods hence one uses a low-frequency model to constraint
inversion results. Comparatively, the model-based inversion methods are time taking
hence it is used when one seeks for quantitative variation of the acoustic impedance.
In the above example of the Blackfoot, one can notice that the inverted section have
very high resolution compared to other post-stack inversion methods and help to inter-
pret prospective zone in the subsurface. One cannot rely only on the interpretation
of the acoustic impedance section and hence for the confirmation of the prospective
zone, one needs to derive more petrophysical parameters. The other petro-physical
parameters can be derived using some perdition tools which are commercially avail-
able. One of these prediction tools is geostatistical methods which are largely used
in the geophysics community. This method is discussed in detail in chapter 8. In this
chapter, we will now discuss the other very important post-stack inversion method
i.e. sparse spike inversion.
The sparse spike inversion methods is another type of post-stack inversion methods
developed in the 1980s (Zhang et al. 2016) but now has become the common method
to estimate petrophysical parameters. There is an assumption in sparse spike inversion
that the earth’s reflectivity is composed of large spikes with small spikes in the
background. Lithological, these large spikes are meaningful which corresponds to
unconformity and hence the sparse spike method tries to find these large spikes
and ignore small spikes (Bosch et al. 2010; Maurya and Singh 2015; Maurya and
Sarkar 2016). The process is similar to the model-based inversion methods where one
assumes simple earth’s reflectivity model and calculate synthetic seismic trace by
the convolution of extracted statistical wavelet and thereafter, the error is estimated
between synthetic seismic and seismic data which can be minimized by adding
more number of spikes in the reflectivity series (Debeye and Riel 1990). When the
sparse spike inversion is constrained by a low-frequency model derived from acoustic
impedance well logs of geologic models, it is commonly referred to as model-based
inversion (Russell 1988).
The purpose of sparse spike impedance inversion is to obtain the acoustic
impedance volume from seismic data. The inverted impedance possesses broad-
band of spectrum, displays a blocky structure of the subsurface in the time domain,
and is directly related to the lithology. The bandwidth of seismic data is generally
10–80 Hz frequency, so lack of high and low-frequency content. To get a broadband
spectrum of the inverted impedance volume, additional information is needed which
generally comes from well log data. The Sparse Spike solution does not have any
low-frequency components due to the band-limited nature of input seismic. There-
fore, the low-frequency trend of the model is imported from the initial model. The
recursive method uses a feedback mechanism to generate a more satisfactory output.
The inversion solution may vary considerably from trace to trace, thus making the
reliability of the output weaker. A low frequency AI variation trend can be imported
3.6 Sparse Spike Inversion 65
to obtain more appropriate results and to get a better convergence for the found solu-
tions from trace to trace. The constrained option uses a low-frequency model as a
guide. The low-frequency variation is estimated from filtered well logs and this gives
much better results. The inversion replaces the seismic trace by a pseudo acoustic
impedance trace at each CDP position.
A flow chart of the sparse spike inversion method is presented in Fig. 3.22 which
includes both inversions. The sparse spike inversion methods are divided into two
broad categories on the basis of the minimization of error. The first method is called
Linear Programming inversion (LPI) that uses a l1 -norm solution for its implemen-
tation and the second method is called Maximum Likelihood inversion (MLI) which
uses l2 -norm solution for its implementation (Russell 1988; Sacchi and Ulrych 1995;
Zhang and Castagna 2011). These methods are briefly explained in the following
sections.
The maximum likelihood inversion first performs the maximum likelihood decon-
volution to estimate the earth’s reflectivity and then it is transformed into the acous-
tic impedance under the maximum likelihood inversion scheme. The fundamental
assumption of Maximum-Likelihood deconvolution is the same as the sparse spike
inversion methods which state that the earth’s reflectivity is composed of a series
of large spikes superimposed over small spikes in the background. This contrasts
with the spiking deconvolution, which assumes a perfectly random distribution of
reflection coefficients.
From the above assumptions about the model parameter, one can derive an objec-
tive function that may be minimized to yield the optimum or most likely reflectivity
(Kormylo and Mendel 1983; Mendel 2012). It is noteworthy that this method gives
an estimate of both the sparse reflectivity and wavelet.
The objective function can be written for reflectivity and noise which may be
minimized to produce the reflectivity and wavelet combination consistent with the
statistical assumption (Russell 1988; Li and Speed 2004; Maurya and Singh 2018).
The objective function E is given by
L
rk2
L
n 2k
E= 2
+ − 2mln(λ) − 2(L − m)ln(λ) (3.29)
k=1
R k=1
N2
where rk is the reflection coefficient at the kth interface, m is the total number of
reflections, L is the total number of samples presents, N is the square root of noise
variance, n k is noise at the kth sample, and λ is the likelihood that a given sample
has a reflection.
Mathematically, the expected behavior of the objective function is expressed in
terms of the parameters shown above. No assumptions are made about the wavelet.
The reflectivity sequence is postulated to be sparse, meaning that the expected number
of spikes is governed by the parameter lambda (λ), the ratio of the expected number
of nonzero spikes to the total number of trace samples. Normally, lambda is a number
much smaller than one (Hampson and Russell 1985). The other parameters needed
to describe the expected behavior are R, the RMS size of the large spikes, and N,
the RMS size of the noise. With these parameters specified, any given deconvolution
solution can be examined to see whether it is likely to be the result of a statistical
process with those parameters. For example, if the reflectivity estimate has a number
of spikes much larger than the expected number, then it is an unlikely result.
In simpler terms, we are looking for a solution with the minimum number of
spikes in its reflectivity and the lowest noise component. It is noteworthy that the
objective function for the one with the minimum spike structure is indeed the lowest
value.
Of course, there may be an infinite number of possible solutions, and it would take
too much computer time to look at each one. Therefore, a simpler method is used
to arrive at the answer. Essentially, one starts with an initial wavelet estimate of the
3.6 Sparse Spike Inversion 67
sparse reflectivity, improves the wavelet and iterate through this sequence of steps
until an acceptably low objective function is reached (Mendel 2012). Thus, there is
a two-step procedure- having the wavelet estimate, update the reflectivity, and then,
having the reflectivity estimate, update the wavelet.
After getting reflectivity from the maximum likelihood deconvolution, it is trans-
formed into the acoustic impedance which is more meaningful in comparison to the
reflectivity series (Goutsias and Mendel 1986). This conversion can be performed in
two ways; the first approach is that the acoustic impedance can be estimated by using
direct relationship between reflectivity and impedance (Velis 2006, 2007; Maurya
and Sarkar 2016).
j
1 + rk
Z j+1 = Z j (3.30)
k=1
1 − rk
But this relationship has a limitation as if the data contains significant noise then
the transformation cannot be performed properly as the relationship does not include
any noise component. A second approach is used to overcome the problem and to
generate an acoustic impedance model from the reflectivity section which included
noise component also and hence gives complete spectrum (Russell 1988; Wang et al.
2006; Zhang et al. 2016). The equation is as follows.
where Z (i) is the known impedance trend, r (i) is earth reflectivity series and n(i) is
the noise in the input trend and H (i) is a factor that would be given by
1, i < 0
H (i) = (3.32)
0, i > 0
Let us look for an example of the maximum likelihood inversion method. The
example is provided for the real data from the Blackfoot field, Alberta, Canada.
Fig. 3.23 Inversion analyses for the composite trace near to well locations for maximum likelihood
inversion methods
Thereafter, for a quality check of inverted results, a cross plot has been generated
of the inverted impedance and actual impedance for all thirteen wells and displayed
in Fig. 3.24. The distributions of scatter points show the quality check of inversion
results and hence performance of the maximum likelihood inversion. It is noticed
that the inverted data points are very close to the best fit curve which represents
the sample by sample comparison of inverted result with the original one and hence
depicts the good performance of the MLI algorithm.
Till now the results are cross-verified qualitatively but now the analysis will
be extended for quantitative comparison so that one can exactly know what is the
difference between inverted and original results.
For the composite data case, we have an estimated correlation coefficient and
RMS Error which are presented in Fig. 3.25. The figure has two curves, the top figure
shows the variation of correlation coefficient whereas the bottom figure depicts RMS
3.6 Sparse Spike Inversion 69
Fig. 3.24 Crossplot of inverted acoustic impedance for composite trace versus original impedance
from well logs
Fig. 3.25 Variation of correlation coefficient (top) and RMS error (bottom) for inverted results of
the composite trace
Error for wells. It is noticed that the correlation coefficients varies from 0.96 to 0.98
with mean value of 0.975 and RMS error varies from 1000 to 2100 m/s*g/cc with
mean value of 1300 m/s*g/cc. These values indicate that the algorithm successfully
searched the optimum solution.
Till now the analysis is limited at one point but one seeks to see the variation
of AI in the subsurface vertically as well as horizontally in the subsurface. Hence,
the maximum likelihood inversion is applied to the 3D seismic volume to estimate
AI in the inter-well region. The data used is from the Blackfoot field, Canada. The
seismic section is inverted for acoustic impedance and is shown in Fig. 3.25 at
inline 28. The special features are highlighted in the inverted impedance section.
The inverted curves show a very high-resolution image of the subsurface. The well
log impedance is also shown above it and shows nice matching between them. The
70 3 Post-stack Seismic Inversion
Fig. 3.26 Cross-section of inverted impedance (inline 41) estimated using MLI technique
3.6 Sparse Spike Inversion 71
Fig. 3.27 Comparison of the amplitude spectrum of Blackfoot seismic and inverted synthetic
seismic
expected reservoir zone horizontally which is predicted in the vertical AI section from
almost all the inversion results discussed above. It is also noticed that the reservoir
zone varies from the SW direction to the NE direction.
The Linear Programming inversion method is another type of sparse spike inversion
technique available to estimate the earth’s reflectivity series that approximates the
seismic data with minimum number of spikes (Li 2001; Russell and Hampson 1991;
Maurya and Singh 2015). The L1 norms (Linear Programming) logarithms are used
by the LPI method to estimate the acoustic impedance of the subsurface (Loris et al.
2007; Zhang et al. 2016). The LPI method tends to remove the wavelet effect from
the seismic data and hence increases the band of inversion results which maximizes
72 3 Post-stack Seismic Inversion
vertical as well as the horizontal resolution of the data and also minimizes the tuning
effects (Helgesen et al. 2000).
The mathematical background of the Linear Programming sparse spike inversion
is developed by Qing Li in (Li 2001; Sacchi and Ulrych 1995, 1996). The aim of
the seismic inversion technique is to find the subsurface model that is reflection
coefficient r (t) in this study, given that seismic data s(t) and source wavelet w(t) as
input, and then calculate acoustic impedance Z (t). There is a relationship between
data d = (x1; x2; . . .), model parameters m = (r 1; r 2; . . .) and noise n as following.
Lm + n = d (3.33)
where L is operator. Generally, the observation d is known, and then the subsurface
model can be defined by the probability p(m|d) (Li 2001; Sacchi and Ulrych 1995).
This probability can be described by the Bayes formula in the following way.
p(d|m) p(m)
p(m|d) = (3.34)
p(d)
where p(d|m) is the probability or likelihood of obtaining the data, p(m) is the prior
probability of the model, p(d) is data likelihood and enters into the problem as a
normalization factor, and p(m|d) is the posterior probability of the model (Barrodale
and Roberts 1973, 1978). One could use the MAP solution, m M A P that maximize a
posterior probability p(m|d). An objective function can be chosen as follows.
In Eq. 3.35, p(d) is a constant and hence omitted for the simplicity. The prior
knowledge of the model is generally given as a global constraint S(m). Thereafter,
one can use the principle of maximum entropy to compute the prior probability
(Debeye and Riel 1990; Sacchi and Ulrych 1995).
Consider continuous model parameters m with probability density function p(m).
The entropy h is given by the following equation.
Equation 3.36 expresses the uncertainty associated with the distribution p(m).
If the information about m is available in the form of some global constraint S(m)
then, the corresponding maximum entropy probability distribution will be given as
follows (Russell 1988; Oliveira and Lupinacci 2013):
1
p 1− p −1 |d − Lm| p
p(d|m) = ex p p (3.38)
2σ p Γ 1p p σp
Based on Eq. 3.38, one can choose the objective function with the different norms
to solve the inverse problem. In this way, one can maximize the objective function E
and can obtain a solution that gives the least error between the model parameters and
the observations (Li 2001; Oliveira and Lupinacci 2013; Maurya and Singh 2018).
The solution of the l1 − norm can be estimated as follows.
Consider the convolution theory which states that the seismograms s(t) are gener-
ated by the convolution of earth’s reflectivity series r (t) with a known wavelet w(t),
if the data does not contain type of noises (Velis 2006, 2007; Wang 2010; Maurya
and Sarkar 2016).
For the layered earth, the reflectivity function will be zero everywhere except
at those times corresponding to the two-way travel time to the jth layer. Thus, the
reflectivity function has the following mathematical form (Oldenburg et al. 1983;
Zhang and Castagna 2011).
NL
r (t) = rjδ t − τj (3.40)
j=1
where NL is the total number of layers in the earth subsurface model, r j is the
reflectivity coefficient at the jth and (j + 1)th interface. Generally, the total number
of data points in the seismogram i.e. N is much larger than the number of layers i.e.
NL in the models hence; it greatly reduces the degree of freedom of the model and
hence reduces the non-uniqueness (Brien et al. 1994; Wang 2010).
Further, recall the acoustic impedance equation and which for the jth layer is
defined as follows.
Z j = ρjvj (3.41)
where ρ j is density and v j is the velocity of the jth layer. The earth reflectivity can
be estimated from the impedance model as following.
Z j+1 − Z j
rj = (3.42)
Z j+1 + Z j
k
1 + rj 1 + rj
Z j+1 = Zj = Z1 (3.43)
1 − rj j=1
1 − rj
74 3 Post-stack Seismic Inversion
Equation 3.43 shows how the impedance is related to the reflectivity coefficients.
The reflectivity coefficients are related to acoustic impedance for continuous travel
time as follows.
1 d[ln Z (t)]
r (t) = (3.44)
2 dt
In deriving Eq. 3.44, second-order quantities have been omitted because they are
almost negligible. Equation 3.44 may be rewritten as follows.
t
Z (t) = Z (0) exp 2 ∫ r (t)dt (3.45)
0
From the above relation, we may see that the sparseness in r (t) may result in a
blocky structure in Z (t) (Oldenburg et al. 1983).
Fig. 3.29 Inversion analysis for composite trace near to well locations. The results for the four
wells are shown here
Fig. 3.30 Crossplot of inverted impedance from the composite trace versus original impedance
from well logs
Variation of correlation coefficient and RMS error between original and inverted
acoustic impedance are shown in Fig. 3.31. The top histogram represents the corre-
lation coefficient and the bottom histogram represents RMS error with wells. It can
be notice that the correlation varies from 0.987 to 0.999 with mean value of 0.993
and RMS error varies from 700 to 1400 m/s*g/cc with mean value of 950 m/s*g/cc.
These values indicate very high performance of the algorithm and represent a very
high quality of the inversion results.
76 3 Post-stack Seismic Inversion
Fig. 3.31 Variation of correlation coefficient (top) and RMS error (bottom) with wells
After getting satisfactory results for the analyses of the composite traces, now, one
can perform inversion of entire seismic volume. Figure 3.32 depicts a cross-section
of inverted impedance at inline 28. The figure shows variation of acoustic impedance
with different colors. The inverted section shows a higher resolution image of the
subsurface compare to the seismic section and highlights the thinner reflector more
clearly. It is noticed that the area has a variety of acoustic impedance from 5000 the
20,000 m/s*g/cc. The analysis of the inverted section also shows a low impedance
anomaly ranging from 6000 to 8500 m/s*g/cc in between 1060 and 1075 ms two
Fig. 3.32 Cross-section of inverted impedance (inline 28) estimated using LPI methods
3.6 Sparse Spike Inversion 77
Fig. 3.33 Comparison of amplitude spectrum between Blackfoot seismic and inverted synthetic
seismic
way travel time. This low AI zone is interpreted as sand channel zone although this is
prelim interpretation and for the confirmation of this channel needs more parameters.
The frequency content of the inverted section becomes degraded as the inversion
involves mathematical transformation which sometimes degrades frequency contents
hence needs to analyze. A comparison of the amplitude spectrum between inverted
synthetic seismic and Blackfoot seismic is shown in Fig. 3.33. The figure depicts
that both the spectrum follow each other very well and indicates that the algorithm
preserves all the frequency contents during its implementations. One more thing
that can be noticed here is that the amplitude spectrums of the inverted curves are
matching exactly to the input seismic which is not in the case of the other inversion
methods. This are great advantages of the LPI methods. Although sometimes the
degradation of the frequency content of the inverted section also depends on the
choice of appropriate parameters during the inversion process. If one does not choose
inversion parameters according to the input data, then the amplitude spectrum will
not match exactly.
The slice at 1065 ms time interval is displayed in Fig. 3.34. A horizontal display
or map view of 3D seismic data having a certain arrival time, as opposed to a horizon
slice that shows a particular reflection. A time slice is a quick, convenient way
to evaluate changes in the amplitude of seismic data. The variations of acoustic
impedance along horizontal are displayed in the figure. The interpreted sand channel
is also highlighted by the ellipse. It is noticed that the sand channel varies from NE
to SW direction.
As the seismic data shows amplitude variation along with time and cannot be
interpreted in the form of subsurface lithology whereas by using these inversion
techniques one can obtain the distribution of petrophysical parameters in the inter-
well region. This is very helpful to interpret subsurface lithology particularly in
exploration and production projects.
78 3 Post-stack Seismic Inversion
References
Ansari HR (2014) Use seismic colored inversion and power-law committee machines based on
imperial competitive algorithm for improving porosity prediction in a heterogeneous reservoir. J
Appl Geophys 108:61–68
Barrodale I, Roberts FD (1973) An improved algorithm for discrete l-1 linear approximation. SIAM
J Numer Anal 10(5):839–848
Barrodale I, Roberts F (1978) An efficient algorithm for discrete l1 linear approximation with linear
constraints. SIAM J Numer Anal 15(3):603–611
Bosch M, Mukerji T, Gonzalez EF (2010) Seismic inversion for reservoir properties combining
statistical rock physics and geostatistics: a review. Geophysics 75(5):75A165–75A176
Brien OM, Sinclair AN, Kramer SM (1994) Recovery of a sparse spike time series by l/sub 1/norm
deconvolution: IEEE Trans Signal Process 42:3353–3365
Brossier R, Operto S, Virieux J (2015) Velocity model building from seismic reflection data by
full-waveform inversion. Geophys Prospect 63(2):354–367
Brown AR (2004) Interpretation of three-dimensional seismic data. AAPG Memoir 42. SEG
Investigation in Geophysics, No. 9. AAPG, Tulsa
Clochard V, Delépine N, Labat K, Ricarte P (2009) Post-stack versus pre-stack stratigraphic inver-
sion for CO2 monitoring purposes: a case study for the saline aquifer of the Sleipner field. SEG
Annual Meeting, Society of Exploration Geophysicists
Cooke DA, Schneider WA (1983) Generalized linear inversion of reflection seismic data. Geophysics
48(6):665–676
Debeye H, Riel VP (1990) Lp-norm deconvolution: geophysical Prospecting 38:381–403
Ferguson RJ (1996) PS seismic inversion: modeling, processing and field examples. M.Sc. Thesis,
University of Calgary, Canada
Ferguson RJ, Margrave GF (1996) A simple algorithm for band-limited impedance inversion.
CREWES Res Rep 8(21):1–10
Goutsias J, Mendel JM (1986) Maximum-likelihood deconvolution: an optimization theory
perspective. Geophysics 51:1206–1220
Hampson D, Russell B (1985) Maximum-likelihood seismic inversion. Geophysics 50(8):1380–
1381
Helgesen J, Magnus I, Prosser S, Saigal G, Aamodt G, Dolberg D, Busman S (2000) Comparison
of constrained sparse spike and stochastic inversion for porosity prediction at Kristin Field. Lead
Edge 19(4):400–407
References 79
Kormylo JJ, Mendel JM (1983) Maximum-likelihood seismic deconvolution. IEEE Trans Geosci
Remote Sens 1:72–82
Lancaster S, Whitcombe D (2000) Fast-track “colored” inversion. SEG Expanded Abstracts
19:1572–1575
Leite EP (2010) Seismic model based inversion using matlab. Matlab-Modelling, Programming and
Simulations, p 389
Li Q (2001) LP sparse spike impedance inversion. Hampson-Russell Software Services Ltd, CSEG
Li LM, Speed TP (2004) Deconvolution of sparse positive spikes. J Comput Graph Statist 13(4):853–
870
Loris I, Nolet G, Daubechies I, Dahlen FA (2007) Tomographic inversion using 1-norm
regularization of wavelet coefficients. Geophys J Int 170(1):359–370
Mallick S (1995) Model-based inversion of amplitude-variations-with-offset data using a genetic
algorithm. Geophysics 60(4):939–954
Maurya SP, Sarkar P (2016) Comparison of post-stack seismic inversion methods: a case study from
Blackfoot Field, Canada. Int J Sci Eng Res 7(8):1091–1101
Maurya SP, Singh KH (2015) Reservoir characterization using model based inversion and
probabilistic neural network. Discovery 49(228):122–127
Maurya SP, Singh KH (2015b) LP and ML sparse spike inversion for reservoir characterization-a
case study from Blackfoot area, Alberta, Canada. In 77th EAGE Conference and Exhibition 2015
(1):1–5
Maurya SP, Singh NP (2017) Seismic colored inversion: a fast way to estimate rock properties from
the seismic data. Carbonate Reservoir Workshop, Nov. 30th –Dec. 1th, 2017, IIT Bombay, India
Maurya SP, Singh NP (2018) Application of LP and ML sparse spike inversion with probabilistic
neural network to classify reservoir facies distribution—A case study from the Blackfoot Field,
Canada. J Appl Geophys, Elsevier 159 (2018):511–521
Maurya SP, Singh KH (2019) Predicting porosity by multivariate regression and probabilistic
neural network using model-based and colored inversion as external attributes: a quantitative
comparison. J Geol Soc India 93(2):131–252
Maurya SP, Singh KH, Singh NP (2019) Qualitative and quantitative comparison of geostatistical
techniques of porosity prediction from the seismic and logging data: a case study from the
Blackfoot Field, Alberta, Canada. Marine Geophys Res 40(1):51–71
Mendel JM (2012) Maximum-likelihood deconvolution: a journey into model-based signal
processing. Springer Science & Business Media, Berlin
Oldenburg D, Scheuer T, Levy S (1983) Recovery of the acoustic impedance from reflection
seismograms. Geophysics 48(10):1318–1337
Oliveira SAM, Lupinacci WM (2013) L1 norm inversion method for deconvolution in attenuating
media. Geophys Prospect 61(4):771–777
Quijada MF (2009) Estimating elastic properties of sandstone reservoirs using well logs and seismic
inversion. Doctoral dissertation, University of Calgary
Russell B (1988) Introduction to seismic inversion methods. The SEG Course Notes, Series 2
Russell B, Hampson D (1991) Comparison of post-stack seismic inversion methods. In: SEG
technical program expanded abstracts, Society of Exploration Geophysicists, pp 876–878
Sacchi MD, Ulrych TJ (1995) High-resolution velocity gathers and offset space reconstruction.
Geophysics 60(4):1169–1177
Sacchi MD, Ulrych TJ (1996) Estimation of the discrete fourier transform, a linear inversion
approach. Geophysics 61(4):1128–1136
Stull RB (1973) Inversion rise model based on penetrative convection. J Atmos Sci 30(6):1092–1099
Swisi AA (2009) Post-and pre-stack attribute analysis and inversion of Blackfoot 3d seismic dataset.
M.Sc. Thesis, University of Calgary
Velis DR (2006) Parametric sparse-spike deconvolution and the recovery of the acoustic impedance.
In: SEG annual meeting. Society of Exploration Geophysicists, pp 2141–2144
Velis DR (2007) Stochastic sparse-spike deconvolution. Geophysics 73(1):R1–R9
80 3 Post-stack Seismic Inversion
Vestergaard PD, Mosegaard K (1991) Inversion of post-stack seismic data using simulated
annealing. Geophys Prospect 39(5):613–624
Wang Y (2010) Seismic impedance inversion using l1 -norm regularization and gradient descent
methods. J Inverse Ill-Posed Prob 18(7):823–838
Wang X, Shiguo Wu, Ning Xu, Zhang G (2006) Estimation of gas hydrate saturation using con-
strained sparse spike inversion: case study from the Northern South China. Sea Terr Atmos Ocean
Sci 17(4):799–813
Waters KH, Waters KH (1987) Reflection seismology: a tool for energy resource exploration. Wiley,
New York
Yilmaz O (2001) Seismic data analysis, vol 1. Society of exploration geophysicists, Tulsa, OK
Zhang R, Castagna J (2011) Seismic sparse-layer reflectivity inversion using basis pursuit
decomposition. Geophysics 76:R147–R158
Zhang Q, Yang R, Meng L, Zhang T, Li P (2016) The description of reservoiring model for gas
hydrate based on the sparse spike inversion. In: 7th international conference on environmental and
engineering geophysics & summit forum of Chinese Academy of Engineering on Engineering
Science and Technology. https://doi.org/10.2991/iceeg-16.2016.27
Chapter 4
Pre-stack Inversion
Abstract Pre-stack seismic inversion utilizes as name suggest pre-stack seismic data
along with well log data to estimate P-impedance, S-impedance, density and VP /VS
ratio away from the borehole. Pre-stack seismic inversion methods provides very
high resolution subsurface images. This method provide more detail information as
compared to the post-stack seismic data in which only P-impedance can be extracted
but in this case one can extract simultaneously P-impedance, S-impedance, density
and VP /VS ratio. The chapter includes details of Simultaneous inversion and elastic
impedance inversion. The application of these methods to the synthetic data as well
as real data are also included in this chapter.
4.1 Introduction
© The Editor(s) (if applicable) and The Author(s), under exclusive licence 81
to Springer Nature Switzerland AG 2020
S. P. Maurya et al., Seismic Inversion Methods: A Practical Approach,
Springer Geophysics, https://doi.org/10.1007/978-3-030-45662-7_4
82 4 Pre-Stack Inversion
schemes such as genetic algorithm (GA) and simulated annealing (SA) (Sen and
Stoffa 1991; Mallick 1995).
There are two distinct phases involved in the latest approach to estimating P-wave
and S-wave rock characteristics from pre-stack seismic information (Goodway et al.
1997; Ma, 2001). The first step is to obtain ordinary reflectivities of P-and S-wave
incidence through the assessment of AVO (Fatti et al. 1994). The second phase is
to introduce an inversion algorithm to the reflectivity sequence derived from AVO,
transform it into acoustic and shear impedances. Ma (2001) created a joint inversion
method that simultaneously measures the impedances of acoustics and shears from
the P-and S-wave reflectivity information derived from AVO.
Pre-stack inversion is often conducted by fitting a 3-term solution to the data,
and the reliability of the results increases with increasing incident angle. The most
accurate result of pre-stack inversion of P-wave seismic data is P-impedance, which
can be performed on short-offset data. S-impedance estimation becomes reliable as
incident angles approach 30°, whereas density evaluation (and other derived elastic
constants) becomes reliable only as incident angles approach 450 . One can invert
the CMP gathers in pre-stack inversion to achieve the impedances of compression
and shear-wave (Hampson et al. 2005). Pre-stack inversion converts seismic angle or
offsets information into P-impedance, S-impedance, and volumes of density through
seismic angle/offset gathering inclusion, well information, and fundamental strati-
graphic analysis. It is similarly feasible to combine other elastic parameters such
as P-impedance, VP /VS proportion and thickness. Typically, based on target and
acquisition, two elements (P-impedance and VP /VS ratio) is accurate and can be
used to forecast reservoir characteristics free from well control (Sherrill et al. 2008).
There are many approaches available for the pre-stack inversion methods but the
most famous approaches are simultaneous inversion and Elastic impedance inversion
methods which are described briefly in the following sections.
The simultaneous inversion method estimates simultaneously more than one param-
eter using pre-stack gather as input. The simultaneous inversion is based on three
assumptions. The first is that the linearized approximation for reflectivity holds
(Ankeny et al. 1986). The linear approximation is an approximation of a general
function using a linear function. Secondly, the reflectivity section as a function
of angle gather can be expressed by the Aki-Richards equations (Ma 2002). The
third approximation is that there is a linear relationship between the P-impedance
and S-impedance and density (Nolet 1978). Using the above three assumptions, the
P-impedance, S-impedance, and density can be modeled by perturbing an initial
guess model. Simmons and Backus (1996) used the following equations to invert the
seismic section into impedance and density sections.
4.2 Simultaneous Inversion 83
1 VP ρ
RP = + (4.1)
2 VP ρ
1 VS ρ
RS = + (4.2)
2 VS ρ
ρ
RD = (4.3)
ρ
ρ 1 VP
= (4.4)
ρ 4 VP
The relation between P-wave and S-wave velocity is given by the Castagna’s
relation (Castagna et al. 1985)
A linearized inversion approach is used to solve for the reflectivity section. Buland
and Omre (2003) have also used a similar technique which is called Bayesian lin-
earized AVO inversion. Unlike Simmons and Backus (1996), they used three terms,
VP /VP , VS /VS and ρ/ρ using the Aki-Richards approximation methods. The
small reflectivity approximation is also used by some authors to relate these parameter
changes to the original parameter. For changes in P-wave velocity, one can write
VP
≈ ln VP (4.6)
VP
where ln represents the natural logarithm. A similar equation can be written for
S-wave velocity and density as follows.
VS
≈ ln VS (4.7)
VS
ρ
≈ ln ρ (4.8)
ρ
most common incarnation is given by Aki and Richards (1980). This modification
reduced 16 equations with 16 unknown parameters to a single equation with three
unknown parameters of Knott-Zoeppritz under some assumptions. Now, first, derive
P-impedance from the post-stack seismic inversion algorithm and then this extended
to the pre-stack inversion algorithms.
In this study, we extend the work of both Simmons and Backus (2003) and Buland
et al. (1996) and build an approach that inverts directly for P-impedance (ZP = ρVP ),
S-impedance (ZS = ρVS ), and density through an approximation similar to that of
Buland and Omre (2003), using constraints similar to Simmons and Backus (1996). It
is also our aim to extend this approach to post-stack impedance inversion (Hampson
1991) so that this technique can be a generalized for pre-stack and post-stack inversion
methods. Therefore, we will first, review the principles of the post-stack inversion
method. First, combining Eqs. 4.1 and 4.6, we will get,
1 1
RPi ≈ ln Zpi = ln Zpi+1 − ln Zpi (4.9)
2 2
where i represents the ith the interface of the subsurface model, Zpi represent p-
impedance of the ith layer, Zpi+1 represent p-impedance at (i + 1)th layer and RPi is
P-reflectivity at ith interface. Now consider N sample points in the reflectivity series
then Eq. 4.9 can be written in matrix form as follows.
⎡ ⎤ ⎡ ⎤⎡ ⎤
RP1 −1 1 0 ... LP1
⎢ RP2 ⎥ ⎢ 0 −1 1 . . . ⎥⎢ LP2 ⎥
⎢ ⎥ ⎢ ⎥⎢ ⎥
⎢ .. ⎥ = ⎢ 0 0 −1 1 ⎥⎢ . ⎥ (4.10)
⎣ . ⎦ ⎣ ⎦⎣ .. ⎦
.. .. .. ..
RPN . . . . LPN
where LPi = ln(ZPi ). Further, the convolution theorem is employed which states that
the seismic trace can be generated by the convolution of the wavelet (w) with the
earth’s reflectivity (R) series (Ma 2002). Mathematically, it can be written as follows.
Si = wi ∗ Ri (4.11)
where Si represent ith sample of the seismic trace and wj represent the jth term of the
extracted seismic wavelet. Combining Eqs. 4.11 and 4.12 gives the forward model
which relates the seismic trace to the logarithm of P-impedance as follows.
1
T= W DLP (4.13)
2
where W is the wavelet matrix given in Eq. (4.12) and D is the derivative matrix
given in Eq. (4.11). If Eq. (4.13) is inverted using standard matrix inversion methods
to give an estimate of LP from knowledge of seismic trace (S) and wavelet (W ),
there are two problems. First, matrix inversion is costly and second, it is potentially
unstable (Hampson et al. 2005). Now, the post-stack inversion is extended to the
pre-stack inversion as follows.
Now, in this section, the derivation is extended for the pre-stack inversion methods.
One can redefine the Aki-Richards equation (Fatti et al. 1994) as follows.
RPP (θ ) = c1 RP + c2 RS + c3 RD (4.14)
1 1
S(θ ) = c1 W (θ )DLP + c2 W (θ )DLS + c3 W (θ )DLD (4.15)
2 2
where LS = ln(ZS ) and LD = ln(ρ). Equation 4.15 could be used for inversion,
except that it ignores the relationship between LP with LS and LD . Because we are
dealing with impedance and have taken logarithms, therefore our relationships are
different than those given by Simmons and Backus (1996) and are given as follows.
LS and LD can be estimated from the trend of the cross plot of ln(ρ) versus ln(ZP )
and ln(ZP ) versus ln(ZS ) respectively (Fig. 4.1).
Now, combining Eq. 4.15 through 4.17, one get the following equation.
86 4 Pre-Stack Inversion
Fig. 4.1 Crossplot of a ln(ZP ) versus ln(ρ) and b ln(ZP ) versus ln(ZS ) wherein both cases a best
fit straight line is added. The deviation away from this straight line LS and LD are the desired
fluid anomalies
If Eq. 4.19 is solved by matrix inversion methods, then we again have the same
problem that the low-frequency component cannot be resolved properly. So a practi-
cal approach is to initialize the solution to [LP LS LD ] = [ln ZP0 00], where ZP0 is
the initial impedance model (Larson 1999). Equation 4.19 is used for the inversion
of pre-stack seismic data. The practical of the inversion can be written as follows.
i. The following information from the seismic section is available for the input.
• One set of N angle traces;
• One set of N wavelets for each angle;
• Initial model for ZP .
ii. Next, calculate the k and m coefficient using well-log cross plot as discussed
above.
iii. Generate the initial guess model by interpolating well log property in the seismic
section.
4.2 Simultaneous Inversion 87
ZP = exp(LP ) (4.21)
The initial guess model representing the initial model of P-impedance while LS
and LD are initialized with zero values in this iteration (Maurya and Singh 2015,
2018). Further, this impedance can be transformed into lame parameters which are
discussed in detail in the following section.
The LMR method was originally proposed by Goodway et al. (1997). The lambda
rho (λρ) and mu rho (μρ) parameters are more sensitive to the fluids and saturation
compare to the velocity and density and therefore analysis in the lambda-mu-rho
domain would be more beneficial (Russell 1988). The LMR transform uses ZP and
ZS from the simultaneous inversion and transform them into lambda rho and mu rho
sections. The LMR parameters can be derived as follows.
The P-wave (VP ) and S-wave (VS ) velocity can be written in terms of lame
parameters as follows.
λ + 2μ
VP = (4.24)
ρ
And
μ
VS = (4.25)
ρ
Where λ and μ are lame parameters and ρ is density. Thereafter, Eq. 4.25 can be
rewritten as follows.
Thereafter, combining Eqs. 4.26 and 4.27, we have the following relation.
The lame parameters are derived by extracting the P-wave and S-wave reflec-
tivities from pre-stack seismic inversion and inverting them to P-wave and S-wave
impedances. Equations 4.26 and 4.28 are used in the analysis of LMR transform.
The LMR analysis is very crucial to identify gas sands (Goodway et al. 1997).
This can be done from the separation of in responses of both the λρ and μρ sections
to gas sands versus shales. In some other cases, the LMR transform is used to sep-
arate lithologies at an even finer scale so as to identify wet sands from shales. The
LMR cross plot is also used to discriminate subsurface lithology (Paul et al. 2001).
Separately neither lambda-rho (λρ) nor mu-rho (μρ) section is a good indicator of
subsurface lithology but together they can be used to reveal a great deal about lithol-
ogy. A flowchart of pre-stack simultaneous inversion and LMR transform has been
presented in Fig. 4.2.
Now we are in a position to implement the above methods to the practical data
sets. Two examples are presented, the first for the synthetic data sets where the true
output is known and hence one can verify their results and second is for the example
from the real data sets from the Penobscot field, Nova Scotia, Canada.
The first example for the pre-stack simultaneous inversion is from the synthetically
generated data sets. For the generation of synthetic seismograms, the parameters
given in Table 4.1 are used as input and forward modeling procedures are performed
to generate synthetic seismograms with angle-dependent. Figure 4.3 demonstrates
a synthetic example for pre-stack simultaneous inversion methods. The first track
of Fig. 4.3 depicts subsurface lithology model used to generate synthetics, track 2
depicts reflectivity series generated from the velocity and density combination given
in geological model, track 3 shows modeled synthetic seismic trace generated using
forward modeling for nonzero offset case (Eq. 1.10), track 4, track 5, track 6 and
track 7 demonstrate comparison of inverted ZP , ZS , ρ and VP /VS ratio with their
original values, respectively. From the figure, one can notice that the inverted curves
are following the trend of the original curve very well. The estimated correlation is
0.99 and RMS error is 0.023 for all four inverted curves which also indicate good
results and hence good performance of the simultaneous inversion methods.
The other example shown here from the real datasets of the Penobscot field,
Nova Scotia, Canada. The Penobscot survey area lies in the Scotian shelf in the
Offshore Nova Scotia, Canada. Several minor oil and gas fields are located in the
Late Jurassic to Cretaceous intervals. The seismic data used for pre-stack simul-
taneous inversion need to be conditioned to remove as many undesirable effects.
Three major undesirable effects commonly removed by data conditioning are ran-
dom noise, NMO wavelet stretch, and non-flat reflections. The conditioning of gather
is briefly explained in Chap. 3 and hence here we have used trim gather as input for
simultaneous inversion.
In a pre-stack simultaneous inversion, invert the pre-stack CMP gathers to obtain
the compressional wave impedance (ZP ), shear-wave impedances (ZS ), density (ρ)
and VP /VS ratio. These parameters are very important to interpret subsurface fea-
tures. Like other inversion techniques discussed in Chap. 3, except colored inversion,
building of an initial model is required for pre-stack inversion. The model is built by
using the P-wave, S-wave velocity and density logs at well locations. From these logs,
the P- and S-impedance logs are estimated (Z = VP ρ), and interpolated between the
wells to build the models by using horizons as guides. The models are filtered by
using a 12-Hz low-pass filter to remove the high-frequency component and well-log
heterogeneity at the seismic frequencies band (Shu-jin 2007). The details of model
building are shown in Chap. 3.
Pre-stack simultaneous inversion relies on the background relationship between
ln(ZP ); ln(ZS ), and ln(ρ) at the well locations, from which the coefficients
90
Fig. 4.4 Crossplot of ln(ZP ) versus ln(ZS ) (left) and ln(ZP ) versus ln(ρ) (right). The measured
parameters from these cross plots are shown bottom of the figures
(k; kc ; m and mc ) are calculated. Deviations of these values from the background
LS and LD are calculated from the inversion itself, and therefore they are initial-
ized equal to zero in the initial model. When the coefficients (k; kc ; m and mc )) are
calculated, they are utilized to determine the final inversion (Shu-jin 2007). These
values are given in Fig. 4.4.
After determining the coefficients of simultaneous inversion, the pre-stack inver-
sion is performed in two stages. In first stage, we apply the inversion at the well
locations for one composite trace to test the parameters and scale the seismic data.
Figure 4.5 compares the inversion result of the P-wave impedance (ZP ), S-wave
impedance (ZS ), density (ρ) and the velocity ratio VP /VS with the corresponding
parameters at well (Well-L-30). From the result, one can notice that the inverted
curves are matching with the original curves very nicely for all parameters. The
average correlation between them is 0.94 and RMS error is 1100 m/s g/cc.
Further, for the quality check of inversion results point by point, a cross plot
of inverted and original curves are generated and displayed in Fig. 4.6. The figure
has four subplots corresponding to compressional wave impedance (ZP ), shear-wave
impedances (ZS ), density (ρ) and VP /VS ratio. The distribution of scatter points from
the best fit straight line quantifies the quality of inversion results. It can be noticed
that most data points lie very close to the best fit line which indicates good quality
of inversion results.
The analysis shows that the simultaneous inversion works well for composite trace.
In the second step entire seismic volume is inverted for the petrophysical volumes.
Figure 4.7 shows a cross-section of inverted acoustic impedance. Inversion of seismic
amplitudes generated the color-coded panel, with low acoustic impedance in green
and yellow, and high acoustic impedance in pink and blue. The acoustic impedance
92 4 Pre-Stack Inversion
Fig. 4.5 Inversion analysis results for composite trace near to well location
of the area varies from 2000 m/s*g/cc to 8000 m/s g/cc. The P-impedance section
(inline 1165) shows smooth variation along with the vertical as well as horizontal
direction. From results, one can also notice that the area has no major anomaly zone.
The shear impedance section (inline 1165) is shown in Fig. 4.8, the density section
(inline 1165) is shown in Fig. 4.9 and the velocity ratio section (inline 1165) is shown
in Fig. 4.10. The shear impedance varies from 100 m/s*g/cc to 4000 m/s g/cc and
the density varies from 2 to 2.5 g/cc. These figures also show smooth variation along
4.2 Simultaneous Inversion 93
with the vertical and horizontal and hence no indication of any major anomaly zone.
If any data contains a reservoir, then low impedance and density anomaly should
appear.
Further, one also uses these sections to estimate the relationship between P-
impedance, S-impedance, and density. Cross plot among these parameters is provided
in Fig. 4.11. This cross plot also provide a derived relationship. These relationships
are better to estimate one of these parameters if one knows the other as compared with
the relationship derived from well logs which are valid only at particular locations
and not away from the well location.
In some situations, the impedance and velocity used to differentiate fluid and rock
properties and hence need other properties to discriminate subsurface lithology. The
lame parameters are found more useful in the situation discussed above. These lame
94 4 Pre-Stack Inversion
parameters are more sensitive and vary more rapidly than impedance and velocity
and hence indicate more clear fluid discrimination.
The ZP and ZS is transformed into lame parameters using LMR transform approach
discussed in the above section. The LMR feature transforms S- and P-impedance
volumes directly into Lambda-Rho and Mu-Rho volumes. This simple yet powerful
transform allows getting more physically meaningful information out of inversion
results. The workflow of the LMR transform is as follows.
• In the first step; calculate RP and RS reflectivities from pre-stack data.
• Thereafter, use the simultaneous inversion method to estimate ZP and ZS volumes
from the AVO and inversion functions.
4.2 Simultaneous Inversion 95
• Further, impedances volumes are transformed into λρ and μρ using Eqs. 4.26 and
4.28 respectively. Note that we cannot decouple density from the other terms.
• Draw cross plot of λρ versus μρ to minimize the effect of density. From this cross
plot, one can classify hydrocarbon sands from wet sands.
Figure 4.8 depicts variation of modeled λρ cross-section (inline 1165) whereas
μρ cross-section (inline 1165) is shown in Fig. 4.9. The modeled section shows very
high resolution as compared with input seismic data and describes layer properties.
One can notice that the λρ varies from 5 to 30GPa g/cc whereas the μρ varies from
0 to 10GPa g/cc. The analysis of modeled lame parameters volume demonstrates
smooth variation along vertical as well as horizontal and hence confirms that the
area does not have any major anomaly zone.
Further, a cross plot of modeled λρ and μρ is generated and is shown in Fig. 4.10.
There are no data points near to lambda-rho and mu-rho axis (no low values) which
96 4 Pre-Stack Inversion
shows that there is no major reservoir in the study area. Lambda-rho and mu-rho
parameters are very important in the gas reservoir because these parameters are
more sensitive compared to the impedances (Figs. 4.12, 4.13 and 4.14).
The Elastic impedance (EI) inversion is another approach of pre-stack inversion meth-
ods which represents a generalization of the acoustic impedance inversion of offset
dependent angle gather. Because P/S mode conversions are significant at oblique
4.3 Elastic Impedance Inversion 97
incidence angle hence EI inversion gives a more detailed picture of the subsurface
(Whitcombe 2002). The EI inversion can be derived from the small-contrast Aki and
Richards’s equation (Eq. 4.29). Connolly (1999) was the first to use range-limited
stack section to invert for elastic impedance. The reflection coefficient R(θ ) can be
calculated by using the Zeopritz equation for P-wave reflectivity for an angle θ as
follows.
VP (ti ) + VP (ti−1 )
VP = , VP = VP (ti ) + VP (ti−1 ),
2
VS2 VS2 (ti ) VS2 (ti−1 )
= − /2
VP2 VP2 (ti ) VP2 (ti−1 )
and similarly for the other variables. In Eq. 4.29, VP represents P-wave velocity, VS
represents S-wave velocity, ρ is density and ti is the time at sample i.
Further, one requires a function f(t) which has properties analogous to acoustic
impedance, such that reflectivity can be derived from the formula given for any
incidence angle θ (Tarantola 1986). Mathematically, it can be written as follows.
f (ti ) − f (ti−1 )
R(θ ) = (4.30)
f (ti ) + f (ti−1 )
Now, this function called elastic impedance (EI), and then for small to moderate
changes in impedance, the above expression can be written as
1 EI 1
R(θ ) ≈ ≈ ln(EI ) (4.31)
2 EI 2
And so, Eq. 4.29 can be rewritten as follows.
98 4 Pre-Stack Inversion
1 1 VP ρ VP VS2 VS VS2 ρ
ln(EI ) = + + −4 2 −2 2 sin2 θ
2 2 VP ρ 2VP VP VS VP ρ
1 VP
+ sin2 θ tan2 θ (4.32)
2 VP
VS2
Let K = VP2
, then Eq. 4.32 can be rewritten as follows.
We used only the first two terms of Eq. (4.29), then the above and following
expressions differ only by changing the tan2 θ to sin2 θ (Whitcombe 2002). We
substitute again ln x for x/x;
ln(EI ) = 1 + tan2 θ ln(VP ) − 8K sin2 θ ln(VS ) + 1 − 4K sin2 θ ln(ρ)
(4.36)
Now if we make K a constant we can take all terms inside the s;
1+tan2 θ 8K sin2 θ
( ) ( )
+ ln ρ (1−8K sin θ )
2
n(EI ) = ln VP − ln VS (4.37)
1+tan2 θ
( ) (8K sin2 θ ) (1−8K sin2 θ )
ln(EI ) = ln VP VS ρ (4.38)
Finally, we integrate and exponentiate (i.e., remove the differential and logarithmic
terms on both sides), setting the integration constant to zero. We have the following
equation.
4.3 Elastic Impedance Inversion 99
Equation (4.39) is the final output which is used to estimate elastic impedance
if one knows the P-wave velocity, S-wave velocity, density and angle θ (Connolly
1999; Lu and McMechan 2004; Mallick 2001; Maurya 2019).
Real data from the Penobscot field, Canada is utilized to demonstrate the capability
of EII to estimate subsurface elastic properties. Pre-stack seismic data from inline
1161–1200 and cross-line 1001–1481 along with well L-30 is used as input. Initially,
data preparation is performed which includes a partial stack, Horizon picking, depth
to time conversion, wavelet extraction and generation of the low-frequency initial
model. Figure 4.15 shows near angle stack (left) and far angle stack (right) which is
used as input for elastic impedance inversion methods. Figure 4.15 depicts that the
Fig. 4.15 Cross-section of near angle stack (top) and far angle stack (bottom) at inline 1165
100 4 Pre-Stack Inversion
signal to noise ratio is increased due to stacking which is the usual case but it is also
noticed that the resolution of the far stack section increases higher as compared to the
near stack gather (highlighted by arrow). The near stack is considered by a 0°–15°
incidence angle while the far stack is considered for 15°–30°. After preparing data,
elastic impedance inversion is applied to data in two steps, first, one composite trace
near to well a location is extracted from near as well as far angle stack gather and
inversion is performed to this trace and results are compared with well log EI.
Figure 4.16 depicts a comparison of inverted elastic impedance from the near
stack as well as a far stack with well log impedance. It is noticed from Fig. 4.16 that
the inverted EI for the near stack as well as far stack are agreeing nicely with the well
log EI. The estimated average correlation is very high (0.94) and error is 0.13 which
depict the good performance of the algorithm. For the quality check, the cross plot
is generated between original and inverted elastic impedance for near (left) as well
as far angle stack (right) data and displayed in Fig. 4.17. From Fig. 4.17, it is noticed
that the maximum scatter points lie near to the best fit line for both cases (near and
far angle stack) and depict the good performance of the algorithm.
As elastic impedance is an extension of acoustic impedance (AI) hence a compar-
ison of EI from near and far angle stack with AI is generated and plotted in Fig. 4.18.
From Fig. 4.18, it is noticed that EI is following the trend of AI very nicely. Gen-
erally, both quantities (EI and AI) are in different range as EI varies from 4000 to
8000 m/s g/cc while AI varies from 5000 to 12,000 m/s g/cc hence shows deviation
with each other particularly at larger impedance zone (Fig. 4.18). A cross plot of near
stack elastic impedance with far stack elastic impedance is generated with a color
4.3 Elastic Impedance Inversion 101
Fig. 4.17 Cross-plot of inverted and original elastic impedance for the near-angle stack and far-
angle stack
Fig. 4.19 Cross-plot of inverted near stack elastic impedance versus far angle stack impedance.
The two different zones are highlighted
bar of acoustic impedance and is shown in Fig. 4.19. From Fig. 4.19, it is noticed
that both EI (near and far angle stack) vary in a similar way. The whole cross plot
is divided into two-zone, zone 1 and zone 2. Zone 1 corresponds to maximum data
points which also follow a trend. On the other hand, zone 2 contains data points
that fall away from the trend line and indicates an anomalous zone. These small
cluster points indicate the presence of gas accumulation in the subsurface. The gas
accumulation is very minor.
Thereafter, the entire seismic section is inverted for elastic impedance in the
second step and cross-section (inline 1165) is shown in Fig. 4.9. The upper side of
the figure shows near angle stacks elastic impedance while far angle stack elastic
impedances are shown bottom of Fig. 4.9. Both sections show the variation of elastic
impedance vertically as well as horizontally in the subsurface. The improved reflector
resolution has been noticed from both the figures but it is higher in the far stack
elastic impedance section as compared with near stack elastic impedance section
comparatively.
A cross plot is always used to see any anomalous zone present in the area and
hence, near stack and far stack elastic impedance is generated for entire datasets and
shown in Fig. 4.10. The whole cross plot is divided into two-zone, zone 1 and zone 2.
Zone 1 corresponds to maximum data points which also follow a trend. On the other
hand, zone 2 contains data points that fall away from the trend line and indicates an
anomalous zone. These small cluster points indicate the presence of gas formation
in the subsurface. Although, zone 2 shows very minor accumulation of gas which
cannot be detected from the seismic section as well as inverted impedance sections.
The location of this minor accumulation of gas is shown in the seismic section and
presented in Fig. 4.15. This is standard way to present anomalous zone in the input
4.3 Elastic Impedance Inversion 103
Fig. 4.20 Cross-section of inverted elastic impedance from near angle stack (top) and far angle
stack (bottom) at inline 1165
Fig. 4.21 Crossplot of inverted elastic impedance from near angle stack versus far angle stack data.
Two different characteristics are highlighted by zone 1 and zone 2
104 4 Pre-Stack Inversion
Fig. 4.22 Anomalous zones (red) are highlighted in the seismic section
seismic data and hence one can demonstrate the capability of these inversion methods
to detect hydrocarbon zone (Figs. 4.20, 4.21 and 4.22).
References
Aki K, Richards PG (1980) Quantitative seismology. W.H. Freeman and Company, San Francisco,
CA
Ankeny L, Braile L, Olsen K (1986) Upper crustal structure beneath the Jemez Mountains volcanic
field, New Mexico, determined by three-dimensional simultaneous inversion of seismic refraction
and earthquake data. J Geophys Res Solid Earth 91(B6):6188–6198
Buland A, Omre H (2003) Bayesian linearized AVO inversion. Geophysics 68(1):185–198
Buland A, Landre M, Andersen M, Dahl T (1996) AVO inversion of troll field data. Geophysics
61(6):1589–1602
Burianyk M, Pickfort S (2000) Amplitude-vs-offset and seismic rock property analysis: a primer.
CSEG Recorder 25(9):6–16
Castagna JP, Batzle ML, Eastwood RL (1985) Relationships between compressional-wave and
shear-wave velocities in clastic silicate rocks. Geophysics 50(4):571–581
Connolly P (1999) Elastic impedance. Lead Edge 18:438–452
Demirbag E, Coruh C, Costain JK, Castagna JP, Backus MM (1993) Inversion of P-wave AVO in
offset-dependent reflectivity. Theory Practice AVO Anal Soc Expl Geophys 287–302
Fatti JL, Smith GC, Vail PJ, Strauss PJ, Levitt PR (1994) Detection of gas in sandstone reser-
voirs using avo analysis: a 3-d seismic case history using the geo-stack technique. Geophysics
59(9):1362–1376
Gardner GHF, Gardner LW, Gregory AR (1974) Formation velocity and density—the diagnostic
basics for stratigraphic traps. Geophysics 39:770–780
Goodway B, Chen T, Downton J (1997) Improved AVO fluid detection and lithology discrimination
using Lamé petrophysical parameters; “λρ” μρ and λ/μ fluid stack” From P and S inversions.
In: SEG Annual Meeting, Society of Exploration Geophysicists, pp 183–186
Gray D, Andersen E (2000) The application of AVO and inversion to the estimation of rock proper-
ties. In: SEG Technical Program Expanded Abstracts, Society of Exploration Geophysicists, pp
549–552
References 105
Hampson D (1991) AVO inversion, theory and practice. The Leading Edge 10(6):39–42
Hampson DP, Russell BH, Bankhead B (2005) Simultaneous inversion of pre-stack seismic data.
In: SEG Technical Program Expanded Abstracts, Society of Exploration Geophysicists, pp 1633–
1637
Larson RG (1999) The structure and rheology of complex fluids, vol 150. Oxford University Press,
New York
Lu S, McMechan GA (2004) Elastic impedance inversion of multichannel seismic data from
unconsolidated sediments containing gas hydrate and free gas elastic inversion of gas hydrate.
Geophysics 69(1):164–179
Ma XQ (2001) A constrained global inversion method using an overparameterized scheme:
application to post-stack seismic data. Geophysics 66(2):613–626
Ma XQ (2002) Simultaneous inversion of prestack seismic data for rock properties using simulated
annealing. Geophysics 67(6):1877–1885
Mallick S (2001) AVO and elastic impedance. Lead Edge 20(10):1094–1104
Mallick S (1995) Model-based inversion of amplitude-variations-with-offset data using a genetic
algorithm. Geophysics 60(4):939–954
Maurya SP, Singh KH (2015) Reservoir characterization using model-based inversion and
probabilistic neural network. Discovery 49(228):122–127
Maurya SP (2019) Estimating elastic impedance from seismic inversion method: a case study from
Nova Scotia field, Canada. Curr Sci 116(4):1–8
Maurya SP, Singh NP (2018) Application of LP and ML sparse spike inversion with a probabilistic
neural network to classify reservoir facies distribution—a case study from the Blackfoot Field,
Canada. J Appl Geophys 159:511–521
Mora P (1987) Nonlinear two-dimensional elastic inversion of multi offset seismic data. Geophysics
52(9):1211–1228
Nolet G (1978) Simultaneous inversion of seismic data. Geophys J Int 55(3):679–691
Pan GS, Young CY, Castagna JP (1994) An integrated target-oriented pre-stack elastic waveform
inversion: sensitivity, calibration, and application. Geophysics 59(9):1392–1404
Paul A, Cattaneo M, Thouvenot F, Spallarossa D, Béthoux N, Fréchet J (2001) A three-dimensional
crustal velocity model of the southwestern Alps from local earthquake tomography. J Geophys
Res Solid Earth 106(B9):19367–19389
Richards PG, Frasier CW (1976) Scattering of elastic waves from depth-dependent inhomogeneities.
Geophysics 41(3):441–458
Russell B (1988) Introduction to seismic inversion methods. The SEG Course Notes, Series 2
Sen MK, Stoffa PL (1991) Nonlinear one-dimensional seismic waveform inversion using simulated
annealing. Geophysics 56(10):1624–1638
Sherrill FG, Mallick S, WesternGeco L.L.C. (2008) 3D pre-stack full-waveform inversion. U.S.
Patent 7,373,252
Shuey R (1985) A simplification of the zoeppritz equations. Geophysics 50(4):609–614
Shu-jin Y (2007) Progress of pre-stack inversion and application in exploration of the lithological
reservoirs. Progr Geophys 3:032
Simmons JL, Backus MM (1996) Waveform-based avo inversion and avo prediction-error.
Geophysics 61(6):1575–1588
Simmons JL, Backus MM (2003) An introduction multicomponent. Lead Edge 22(12):1227–1262
Tarantola A (1986) A strategy for nonlinear elastic inversion of seismic reflection data. Geophysics
51(10):1893–1903
Whitcombe DN (2002) Elastic impedance normalization. Geophysics 67(1):60–62
Chapter 5
Amplitude Variation with Offset (AVO)
Inversion
Abstract In this chapter, the details of Amplitude variation with offset has been
discussed. The method utilizes CMP gathers to produce AVO attribute which is used
as good indicators of the presence of gas formation. The chapter is divided into two
parts, in the first part it discusses about fluid replacement modeling and in the second
section, it discusses about AVO inversion. Initially, the details of these methods are
provided and thereafter, synthetic as well as real data examples are discussed.
5.1 Introduction
characteristics, the rock is saturated with the fresh pore liquid, calculating the new
efficient bulk moduli and density.
The most frequently used method is to apply the equations of Gassmann. The
Gassmann equation relates the bulk modulus of the porous rock frame, mineral
matrix and pore fluids (Batzle and Wang 1992). The use of Gassmann’s equation
is based on the following assumptions. The rock is homogeneous and isotropic, all
pores are interconnected and interacting, pore pressure is balanced throughout the
rock, the pore fluid does not communicate with the solid to soften or harden the frame,
and the media is closed and no pore fluid leaves the volume of the rock (Smith et al.
2003; Wang 2001).
Changing the fluid type and its saturation value at the reservoir and creating syn-
thetic logs relating to these changes is introduced as fluid replacement modeling. In
most cases, only one well data is available in the study area which has encountered a
specified horizon (oil, gas or water). In such a situation, modeling of AVO behavior
for all conditions of the reservoir including gas, oil, and water is not possible. The
solution is the modeling of reservoir rock and existing fluids, and then the prediction
of synthetic logs with substituting fluid type using this modeling (Verm and Hilter-
man 1995). By having P-wave logs, shear wave, and density and using Zoeppritz’s
equations or one of its approximations, AVO behavior of reservoir can be estimated
in different conditions of pore fluid and comparing that with seismic data through
providing synthetic seismograph. Gassmann’s (1951) equations are used to model
reservoir rocks and fluid in almost all fluid modeling project. It is used to predict
reservoir seismic properties changes (density, compressional and shear wave veloc-
ity) caused by fluid replacement based on texture property of reservoir rock. The
Gassmann equation can be written as follows.
2
K dr y
1− Km
K sat = K dr y + φ K dr y
(5.2)
Kf
+ 1−φ
Km
+ K m2
In which K sat is the bulk modulus of the rock saturated with pore fluid, K dr y
is the bulk modulus of the dry frame (drained of any pore-filling fluid), K m is the
volumetric module resultant of forming minerals in reservoir rock and φ is reservoir
rock porosity (Kumar 2006). The supporting equation can be written as follows.
4
K sat = ρb V P2 − VS2 (5.3)
3
μ = ρb VS2 (5.4)
ρb = ρma (1 − φ) + ρ f l φ (5.7)
Gassmann’s equations relate the bulk modulus of rock to its pore, frame, and fluid
properties. Recall the Gassmann equation to estimate the bulk modulus of a saturated
rock.
2
K f rame
1− K matri x
K sat = K f rame + φ (1−φ) K f rame
(5.9)
K fl
+ K matri x
− 2
K matri x
where, K f rame , K matri x , and K f l are the bulk moduli of the saturated rock and
porous rock frame (drained of any pore-filling fluid), mineral matrix, and pore fluid,
respectively, and φ is porosity (as fa reaction). In the Gassmann formulation, the
shear modulus is independent of the pore fluid as discussed earlier and is constant
during the fluid substitutions.
In order to estimate the saturated bulk modulus (Eq. 5.9) in a given condition
of the reservoir and fluid type, it is necessary to estimate the bulk module of frame,
5.2 Fluid Replacement Modeling 111
matrix and pore fluid (Han and Batzle 2004). In the following section, we will address
formulations for calculating bulk modulus and mineral matrix density, pore fluid, and
rock frame.
Bulk module and pore fluid density (brine, oil, and gas) are estimated by an average
of fluid type values. Let’s first calculate each type of liquid (brine, oil, and gas)
characteristics.
(i) Estimation of Bulk modulus and density of brine
The bulk brine module can be estimated from known seismic velocity and brine
density. The related equation can be written as follows.
In Eq. 5.10, K brine is bulk modulus in (GPa), ρbrine is brine density in g/cm3 and
Vbrine is the P-wave velocity in brine and can be taken as m/s (Han and Batzle 2004).
The brine density can be estimated using equation suggested by Batzle and Wang
(1992). This can be written as follows.
On the other hand, the P-wave velocity of brine can also be calculated using the
formula suggested by Batzle and Wang (1992) and can be written as follows.
Vbrine = Vw + S 1170 − 9.6T + 0.055T 2 − 8.5 × 10−5 T 3 + 2.6P − 0.0029T P
− 0.0476P 2 + S 1.5 780 − 10P + 0.16P 2 − 1820S 2 (5.13)
In Eq. 5.13, Vw shows the P-wave velocity in pure water. There is a relation
between Vw , temperature (T) and pressure (P), and hence they can use to calculate
P-wave velocity as follows.
112 5 Amplitude Variation with Offset (AVO) Inversion
4
Vw = wi j T i−1 P j−1 (5.14)
i=1 j=1
In Eq. 5.14, wi j present constant and given by Batzle and Wang (1992). These
values are given in Table 5.1.
(ii) Estimation of Bulk modulus and density of oil
The oil contains some amount of gas dissolved in it and characterized by GOR (gas-
to-oil ratio). The density of oil and the bulk modulus depend on the temperature,
pressure, GOR and type of oil present. The oil density (ρoil ) values can be estimated
using the relationship provided by Batzle and Wang (1992) and Wang (2001) as
follows.
ρ S + 0.00277P − 1.71 × 10−7 P 3 (ρ S − 1.15)2 + 3.49 × 10−4 P
ρoil = (5.15)
0.972 + 3.81 × 10−4 (T + 17.78)1.175
ρ0 + 0.0012RG G
ρS = (5.16)
B0
where ρ0 is the reference density of oil measured at 15.6 °C and taken in (g/cm3 ),
RG is the gas-oil-ratio (GOR) given in litre/litre, G is the specific gravity of oil in
(API) and B0 is called the formation volume factor and can be calculated as follows.
1.175
G
B0 = 0.972 + 0.00038 2.49RG + T + 17.8 (5.17)
ρ0
Further, the P-wave velocity in oil (Voil ) can also be calculated by the relationship
given by Batzle and Wang (1992) and Wang (2001) and the equation can be written
5.2 Fluid Replacement Modeling 113
as follows.
ρ ps 18.33
Voil = 2096 − 3.7T + 4.64P + 0.0115 − 16.97 − 1 T P
2.6 − ρ ps ρ ps
(5.18)
ρ0
ρ ps = (5.19)
(1 + 0.001RG )B0
Using the relationship given above, one can estimate velocity and density of oil
and hence, the bulk modulus of oil, K oil (GPa) can be calculated as follows.
Fluid composed of brine and hydrocarbon (oil and/or gas) in the pore spaces and
hence the bulk modulus and mixed pore liquid phase density can be estimated by
reverse bulk modulus average and arithmetic density average of distinct liquid stages,
respectively (Inyang 2009). Mathematically, the bulk modulus (K f l ) can be written
as follows.
1 WS HS
= + (5.21)
K fl K brine K hyc
where W S is the water saturation (as fraction) and HS (=1 − WS) is the hydro-
carbon saturation, K hyc and ρhyc are the bulk modulus and density of hydrocarbon,
respectively. Similarly, fluid density (ρ f l ) can be written as follows.
28.8G P
ρgas ∼
= (5.23)
Z R(T + 273.15)
In Eq. 5.23, G represents the specific gravity of gas (API), R is the gas constant
taken as 8.314 and Z is the compressibility factor which can be estimated using the
relation between pseudo reduced temperature and pressure and is given as followings.
3
Z = 0.03 + 0.00527 3.5 − T pr Ppr + 0.64T pr − 0.007T pr
4
− 0.52
⎧ ⎫
⎨ Ppr1.2 /Tpr ⎬
2 1 2
+ 0.109 3.85 − T pr exp − 0.45 + 8 0.56 − (5.24)
⎩ T pr ⎭
T + 273.15
T pr = (5.25)
94.72 + 170.75G
And
P
Ppr = (5.26)
4.892 − 0.4048G
Further, the bulk modulus of gas is suggested by Batzle and Wang (1992) and can
be written as follows.
0.85 + Ppr +2 +
5.6 27.1
− 8.7 exp −0.65 Ppr + 1
P ( Ppr +3.5)2
K gas ∼
=
1 − Zpr ∂∂PZpr 1000
P
T
(5.27)
And
∂Z 3 2
= 0.03 + 0.00527 3.5 − T pr + 0.109 3.85 − T pr F (5.29)
∂ Ppr T
5.2 Fluid Replacement Modeling 115
Using the above relationship, one can calculate the bulk modulus and density
of the gas. Once these parameters are known, one can move to calculate matrix
properties.
To calculate the mineral matrix’s bulk modulus, one needs to understand the rock’s
mineral structure that can be discovered from the core sample laboratory examination.
Lithology may be presumed to be a combination of quartz and clay minerals in the
lack of laboratory information (Kumar 2006). The proportion of clay can be obtained
from the curve of volume shale (Vsh ), typically obtained from the information of the
wireline log (Gamma-ray log). It can be noted that the shale formation contains about
70% of clay along with 30% quartz. Once the mineral abundances are determined
one can calculate the bulk modulus of the matrix (K matri x ) using Voigt-Reuss-Hill
(VRH) (Fig. 5.1) averaging formula (Hill 1952).
1 Vclay Vqt z
K matri x = Vclay K clay + Vqt z K qt z + + (5.30)
2 K clay K qt z
where Vsh , K clay and K qt z is shale volume, the bulk modulus of clay and bulk modulus
of quartz, respectively. These parameters can be calculated as follows.
And
40
35
30
25
K (GPa)
20
15
10
5
0
0 0.2 0.4 0.6 0.8 1
Porosity
Fig. 5.1 Curves of Reuss (R), Voigt (V) and Hashin-Shtrickman (HS), showing bulk modulus
variation with the porosity. Using Ko = 38 GPa and Go = 44 GPa, clean sandstone parameters
(After Hashin and Shtrikman 1963)
116 5 Amplitude Variation with Offset (AVO) Inversion
The density of the mineral matrix ρmatri x can be estimated by arithmetic averaging
of densities of individual minerals as
where ρclay and ρqt z are the density of the clay and quartz minerals. The standard
value of these parameters can be found in Table 5.2 (Mavko and Mukerji 1998).
From the formula, one can notice that the K matri x and ρmatri x remain constant
during the Gassmann fluid substitution.
All the parameters in Eq. (5.34) are known from the previous formulations and
hence one can estimates K f rame which remain unchanged during fluid substitution.
Seven layers have been classified from the well 01-17 of the Blackfoot field, Canada
and the sixth layer (Layer 5) is interpreted as sandstone formation (Table 5.3). The
fluid substation is performed for the sandstone layer. Table 5.4 shows list of constantly
used as input in fluid replacement modeling.
From the analysis, it is noticed that the P-wave velocity drops steeply between 0
and 20% saturation and then increases slowly (Fig. 5.2a). The maximum decrease
in P-wave velocity is found to be 2.7%. On the other hand, the S-wave velocity
5.2 Fluid Replacement Modeling 117
Table 5.3 The layering of well 01-17 from Blackfoot field, Canada and corresponding velocity,
density and porosity values
Layer P-velocity (m/s) S-velocity (m/s) Density (kg/m3 ) Phi
Ocean (0–162 m) 1480 0 1000 NA
Layer 1 (162–510) 2978 1396 2060 0.36
Layer 2 (510–812) 3248 1628 2150 0.30
Layer 3 (812–1332) 3515 1857 2310 0.20
Layer 4 (1332–1545) 3963 2244 2520 0.08
Layer 5 (1545–1584) 3932 2217 2460 0.12
Layer 6 (1584–1600) 4690 2871 2520 0.08
increases with CO2 saturation and shows direct proportionality (Fig. 5.2b). The aver-
age increase of the S-wave velocity with respect to 100% water saturation is 2.44%.
The other parameter analysed here are density which decreases with increasing CO2
saturation (Fig. 5.2c). The average decrease of density with CO2 saturation is found
to be 4.73% as compared with 100% water saturation. Further, V P /VS the ratio is
analyses and found that the V P /VS ratio decreases with increasing CO2 saturation
(Fig. 5.2d). This can be expected as the S-wave velocity increase and P-wave veloc-
ity decrease. The maximum drop-in V P /VS ratio found to be 7.75% in comparison
with 100% water saturation. Figure 5.2d shows variation of V P /VS ratio and depicts
a sharp drop between 0 and 20% CO2 saturation and thereafter it becomes nearly
constant (Table 5.5).
Further, the effect of this fluid substitution can also be monitored from the seismic
amplitude variation. As we have already seen that the velocity and density changes
due to fluid substitution and hence these petrophysical changes are expected to be
reflected in the seismic response. In this regard, the synthetic seismograms are gen-
erated using forward modeling procedure from the velocity and density of the well
log 08-08. The forward modeling is explained briefly in the following paragraphs.
118 5 Amplitude Variation with Offset (AVO) Inversion
a 0.0
-0.5
Vp Change (%)
-1.0
-1.5
-2.0
-2.5
-3.0
0 0.2 0.4 0.6 0.8 1
CO2 saturation
b 2.0
Vs Change (%)
1.5
1.0
0.5
0.0
0 0.2 0.4 0.6 0.8 1
CO2 saturation
c 0.0
Rho change (%)
-1.0
-2.0
-3.0
-4.0
0 0.2 0.4 0.6 0.8 1
CO2 saturation
d 0.0
Vp/Vs change (%)
-1.0
-2.0
-3.0
-4.0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
CO2 saturation
Fig. 5.2 a P-wave velocity change, b S-wave velocity change, c density change, d V p/V s changes
with CO2 saturation
After fluid submission, pre-stack synthetic seismograms are generated. The angle-
dependent reflection coefficient (R pp (θ )) are calculated following Shuey’s (1985)
three-term approximation.
R pp (θ ) ≈ R(0) + G sin2 θ + F tan2 θ − sin2 θ (5.35)
5.3 AVO Modeling 119
Table 5.5 Gassmann fluid substitution results and input values (porosity, P-wave velocity, S-wave
velocity, bulk modulus, V p/V s and shear modulus)
Water CO2 VP VS ρ VP VS ρ VP /VS
Sat Sat change change change change
1 0.0 3969.9 2179.2 2.44 0.0 0.00 0.00 0.00
0.9 0.1 3875.5 2183.3 2.43 −2.4 0.19 −0.37 −2.56
0.8 0.2 3872.6 2187.5 2.42 −2.5 0.37 −0.74 −2.85
0.7 0.3 3873.8 2191.4 2.41 −2.4 0.56 −1.11 −2.96
0.6 0.4 3878.7 2195.5 2.40 −2.3 0.75 −1.48 −3.02
0.5 0.5 3884.5 2199.7 2.39 −2.1 0.94 −1.85 −3.06
0.4 0.6 3890.8 2203.8 2.38 −2.0 1.13 −2.22 −3.08
0.3 0.7 3897.5 2208.0 2.37 −1.8 1.32 −2.59 −3.10
0.2 0.8 3904.4 2212.2 2.36 −1.6 1.51 −2.96 −3.12
0.1 0.9 3911.4 2216.5 2.35 −1.5 1.71 −3.33 −3.13
0.0 1.0 3918.6 2220.7 2.35 −1.3 1.90 −3.70 −3.14
Depth of 1545 m in the well site. Reservoir thickness 39 m
And the variation of the reflectivity versus the angle, the AVO gradient, G can be
calculated as follows.
Vρ 1 2Vs2 4Vs2 Vs
G = R(0) − + 2
− (5.37)
Vρ 2 Vs V p2 Vs
And F, the reflectivity at the far angles (reflection angles higher than 30°) can be
expressed as follows.
1 V p
F= (5.38)
2 Vp
Each elastic property is defined on each side of the interface where the reflection
is happening as follows.
V p = V p2 − V p1
V p2 + V p1
Vp =
2
Vs = Vs2 − Vs1
120 5 Amplitude Variation with Offset (AVO) Inversion
Vs2 + Vs1
Vs =
2
Vρ = Vρ2 − Vρ1
Vρ2 + Vρ1
Vρ =
2
Index 1 refers to the vertical location at which the calculation of the reflection
coefficient is performed, the mean above the reflection interface; while index 2 refers
to the sample immediately below, the mean below the reflection interface (Shuey
1985). Figure 5.3 shows the process of the generation of the synthetic trace.
By using the linear approximation of Shuey within a geostatistical inversion
approach, we can obtain the AVO normal-incidence, R(0) and gradient, G, cubes
immediately and as part of the inverse solution (Rutherford and Williams 1989;
Avseth and Bachrach 2005; Castagna and Backus 1993; Maurya and Singh 2019).
The propagation of uncertainty directly from the inverse problem to the inverted vol-
umes of AVO allows for better risk assessment of interest anomalies in amplitude,
which may be related to actual accumulations of hydrocarbons. It is crucial to note,
however, that any alternative approximation can be used to calculate the coefficients
of angle-dependent reflection.
Fig. 5.3 Forward modeling procedure where the first track shows impedance, the second track
shows reflectivity, third track shows wavelet and fourth track shows synthetic seismograms
5.3 AVO Modeling 121
The fluid replacement modeling has been performed for well log data (Penobscot
L-30) from the Penobscot field, Scotian Shelf, Canada. The properties of fluid
replacement modeling are as follows (Table 5.6).
With all the above parameters, AVO modeling is performed for Penobscot data,
Canada and first, variation in rock properties like velocity, density, poison’s ratio is
monitored and then seismic amplitude variation. Three saturations have been mod-
eled, the first saturation is taken as pure brine, the second saturation is pure gas and
the third saturation has been taken as pure oil. Figure 5.4 shows variation of rock
properties due to fluid substitution. It is noticed that the average change in P-wave
velocity in the target zone is as follows. The change in velocity has been monitored
by replacing reservoir fluid by pure brine, pure gas, and pure oil case one by one
respectively. The 20% decrease has been noticed for pure brine case, a 23% decrease
has been noticed for the pure gas case and a 27% decrease has been noticed for the
pure oil case respectively.
comparatively smooth sand that the existence of gas dramatically improves the rock’s
compressibility, the velocity falls and the amplitude reduces producing a negative
polarity, which explains the word bright spot. A seismic amplitude anomaly of high
amplitude that could indicate the presence of hydrocarbons is called a bright spot.
Bright spots are the result of major changes in acoustic impedance and effect of
tuning, such as when gas sand underlies shale, but can also be caused by occurrences
other than the presence of hydrocarbons, such as a lithological change (White 1977;
Mocnik 2011). The term is frequently used as a hydrocarbon indicator. These ampli-
tude peaks can be caused by a gas resulting in increase in the reflection coefficient
in the pore space. It is used to define the rise in amplitude rather than the presence of
hydrocarbons, and in unconsolidated clastic rocks, it is generally higher. As shown
in Fig. 2.1, the reduced sand impedance generates an amplitude increase over the
crest of the hydrocarbon structure.
ii. Flat spot
The flat spot depicts a seismic reaction of hydrocarbon contact where it appears to be
flat. Such contact can be between gas and oil, between oil and water, or between gas
and water. To depict a flat spot, the hydrocarbon reservoir is much thicker than the
vertical resolution. Flat spots are often hard to locate; channel edge or base, low angle
flaws, or artifacts processing can often be misunderstood as flat spots (Luping 2008;
Mocnik 2011). Low saturated gas in a reservoir can also cause flat spots. Finding the
so-called flat spot at the hydrocarbon-water touch is very crucial. This is a difficult
reflector (increase in impedance) and should be on the same TWT as the shift in
amplitude. If both oil and gas are present, two separate flat spots should be apparent,
one at the contact between gas and oil and one closer to the contact between oil and
water (Backus and Chen 1975; Swan et al. 1993).
iii. Dim spot
Highly consolidated sands with much higher acoustic impedance than the overlying
shale are causing dim spots. As seen in Fig. 5.5, there is a strong peak at the top
of the sandstone. Adding hydrocarbons to the pore space can cause the sandstone’s
velocity and density to reduce; however, the polarity of the reflection coefficient
will not reduce sufficiently (Chen and Sidney 1997). The hydrocarbon decreases the
impedance of acoustics and the coefficient of reflection, thus producing a dim spot.
This situation can often be hard to interpret with confidence, especially if minor fault
affects the structure. We expect to see a “dim place” or a decline in amplitude at the
top reservoir for difficult brine sand relative to the shale and difficult hydrocarbon
sand relative to the rock (Fig. 5.5). A dim spot is not simply to acknowledge and
should be in line with the flat spot’s structure and TWT (Mocnik 2011; Brown 2012).
Therefore, the elements to consider during the analysis of seismic information are
several and mostly rely on the acoustic effects of a gas accumulation, the overlying
material, the porosity, the depth, the water saturation and the configuration of the
reservoir. Thus, the only observation of an amplitude anomaly may not be adequate
to attribute a hydrocarbon origin; it is very essential to consider all the impacts on
124 5 Amplitude Variation with Offset (AVO) Inversion
Fig. 5.5 Schematic model of a bright spot, polarity reversal and dim spot for different oil/gas brine
sand responses (From Bacon et al. 2007)
the seismic signal generated by gas/oil presence, including, in specific, its elements,
which is not only amplitude but also frequency and phase (Li and Mueller 1997).
iv. Polarity reversal
If the reservoir material in its water-saturated state has an acoustic impedance greater
than the surrounding material, and if the local water-to-gas replacement depresses
the acoustic impedance of the gas saturated zone to less than that of the surrounding
material, then the reflection from the top of the reservoir shows a reversal of polarity
over gas. This is not needed proof for gas as these circumstances may or may not
occur; it may also be hard to identify polarity reversal owing to interferences and
tiny fault impacts (Keys 1989; Chen and Sidney 1997; Luping 2008).
Phase change, also recognized as polarity reversal, happens when the surrounding
reservoir has a reduced reservoir rock velocity. This can happen when a partly consol-
idated sand becomes moist. This leads to mildly greater acoustic impedance than the
overlying shale. The bottom of the sand correlates with a small favorable coefficient
of reflection. The velocity and density of the sandstone decrease as hydrocarbons
are introduced to the pore space; which also results in a reduction in the acoustic
impedance to the stage where it is mildly less than the overlying shale (Tai et al. 2009;
Mocnik 2011). The reflection at the base then switches stage, from a maximum to a
soft trough, which can be seen in Fig. 5.5.
Chiburis et al. (1993) indicated that the use of AVO for identifying fluid is achieved by
using AVO comparing the actual information with a synthetic seismogram. Hampson-
Russell’s AVO is used in this study to generate AVO synthesized seismograms for
5.3 AVO Modeling 125
fluid-saturated soils (petroleum, brine, and gas) with input density and recorded
velocity (P-wave and S-wave) from well L-30 of Penobscot field, Canada. With the
use of Zoeppritz and elastic flow equations, synthetics are produced for interval 2100–
2300 ms and analytical outcomes are contrasted. While both equations calculate
the amplitudes of seismic waves, the equations of Zoeppritz recognize only plane-
wave amplitudes of reflected P-waves and disregard intermittent multiples and mode-
converted waves (Russell 1999). On the other side, the elastic wave algorithm designs
multiple waves as well as converted mode waves. Synthetics are developed using
both algorithms taking into consideration failures in transmission and distributing
geometry. These would lead to a fake Class III AVO. The seismic event of concern
is the elevated amplitude corresponding to the boundary between reservoir sand and
enclosing shales. Off the three synthetic seismograms (oil, brine, and gas) all show
large amplitudes of normal incidents, the gas synthetics show the largest amplitude,
brine shows the smallest, and oil falls between brine and gas.
These hydrocarbon indicators can be observed by the modeling of gas, oil, and
brine in the reservoir. We have taken three cases, first, the reservoir is saturated with
pure brine and corresponding changes in seismic sections are noticed. Thereafter, in
the second and third cases, pure oil and pure gas substitution have been performed
and corresponding seismic amplitude changes are noticed. Further, a comparison of
seismic amplitude with offset in all three cases is also performed. These changes are
provided in Figs. 5.6, 5.7 and 5.8.
From fluid replacement, it can be seen that the impact of gas saturation in the
reservoir is greater than the same quantity of oil. The gas-filled reservoir rocks have
the lowest impedance and the largest amplitudes of reflection, while brine rocks have
the highest impedance and the lowest amplitudes of reflection. On the other side, oil
sands have characteristics that lie between the characteristics of gas and brine.
Fig. 5.6 Variation of seismic amplitude with offset for pure brine (left) and comparison of seismic
amplitude (right)
The subsurface of the earth is quite complicated, so that distinct rocks have distinct
reactions to AVO, yet these rocks are filled with the same fluid or have the same
porosity.
Rutherford and Williams (1989) developed the classification scheme for AVO
anomalies and altered it by Ross and Kinman (1995) and Castagna and Swan (1997).
Class 1: High-impedance contrast with decreasing AVO. The layer has a higher
impedance than the surrounding shales.
Class 2: Near-zero impedance contrast between the sand and surrounding shales.
Class 2p: Same as class 2, but with polarity change;
Class 3: Low impedance with increasing AQVO compared to surrounding shales.
Class 4: Low impedance sand with decreasing AVO.
Graphically, it can be represented as in Fig. 5.9. Figure 5.9 shows plot of the
5.4 Amplitude Variation with Offset (AVO) Analysis 127
Fig. 5.7 Variation of seismic amplitude with offset for pure oil case (left) and comparison of seismic
amplitude (right)
reflection coefficient with the incident angle. The four different kinds of AVO class
is defined. These classifications are used to detect gas-bearing formation from no gas
formations.
Basic AVO theory is well understood because it is widely used as a tool in hydrocar-
bon detection (Smith and Gidlow 1987). A few of the most important ideas will be
highlighted to keep in mind when doing AVA analysis. The AVO can be understood
by study the partition of energy from an interface. Figure 5.10 shows the theoretical
energy partition at a boundary. This figure illustrates an important point that accounts
for AVA phenomena: the conversion of P-wave energy to S-wave energy. Though the
majority of seismic data is recorded simply as a single component pressure wave,
the fact that the Earth is elastic causes amplitudes of P-wave arrivals to be a function
of S-wave reflection coefficient (Rs ) (Castagna and Smith 1994). In practice, Rs is
tricky to obtain and the P-wave reflection coefficient (Rp ) is what we have in the vast
majority of cases (Smith and Sutherland 1996).
128 5 Amplitude Variation with Offset (AVO) Inversion
Fig. 5.8 Variation of seismic amplitude with offset for pure gas case (left) and comparison of
seismic amplitude (right)
Using Snell’s law, Knott in 1899, and Zoeppritz in 1919, developed general expres-
sions for the reflection of compressional and shear waves at a boundary as a function
of the densities and velocities of the layers in contact (Knott 1899; Zoeppritz 1919).
Although Zoeppritz was not the first to publish a solution, his name is associated with
the cumbersome set of formulas describing the reflection and refraction of seismic
waves at an interface (Aki and Richards 1980).
i. Zoeppritz equation
The relationship between the incident and reflection/transmission amplitudes for
plane waves at an elastic interface is described by the Zoeppritz equations. These
equations give us the exact amplitudes as functions of the incidence angle. Figure 5.10
explains partitioning of energy of a boundary.
5.4 Amplitude Variation with Offset (AVO) Analysis 129
⎡ ⎤ ⎡ sin θ2
⎤
R p (θ1 ) − sin θ1 − cos φ2 cos φ2
⎢ R S (θ1 ) ⎥ ⎢⎢ cos θ1 − sin φ1 cos θ2 − sin φ2 ⎥
⎥
⎢ ⎥
⎣ TP (θ1 ) ⎦ = ⎢ ρ2 VS2 ⎥
2
V P1 V P2 ρ2 VS2 V P1
⎣ sin 2θ1 VS1
cos 2φ1 ρ V2 V
1 P1
cos 2φ2 ρ1 VS1
cos 2φ2 ⎦
S1
ρ2 V P2 ρ2 VS2
TS (θ1 ) − cos 2φ1 VS1
sin 2φ1 cos 2φ2 − ρ1 VP1 sin 2φ2
V P1 ρ1 V P1
⎡ ⎤
sin θ1
⎢ cos θ1 ⎥
×⎢
⎣ sin 2θ1 ⎦
⎥ (5.39)
cos 2φ1
P VS ρ
R(θ ) = a +b +c , (5.40)
VP VS ρ
2 2
where a = 2cos1 2 θ , b = 0.5 − VVPS sin2 θ , and c = 4 VVPS sin2 θ . Further, V P
represents p-wave velocity and VS presents S-wave velocity, ρ is density and R is a
reflection coefficient.
iii. Wiggins approximation
Wiggins et al. (1983) separated the three terms of the Zoeppritz equations related to
perturbations in the elastic parameters of interest (Russell 1988). This can be written
as follows.
The method of gradient analysis enables parameters of the AVO function to be tested
on a chosen seismic data collection. Therefore, this feature can be used to see what
attributes to use in the volume and map functions of the AVO attribute. It also enables
one to determine an AVO anomaly category (Ramos and Davis 1997).
The chart produces the Intercept when one models the amplitude of the wave for
a reflector against the trace offset. This is where the zero-offset line meets the trend
of amplitude measurements. The gradient, which is the curve slope produced by the
plot points, is also produced (Li et al. 2007). It is then possible to use the sums or
132 5 Amplitude Variation with Offset (AVO) Inversion
differences of these gradients and intercept values to map AVO anomalies. The offset
is traditionally described as the sine squared of the angle of the event (Figs. 5.11 and
5.12).
Fig. 5.12 Analysis of AVO anomalies at around 1626 ms. The red line on the seismic display
shows the time location at which the amplitudes have been extracted. Those amplitudes are plotted
as squares on the right-hand graph. The curve which has fit through the picks is a plot of Aki-Richards
two-term equation
5.4 Amplitude Variation with Offset (AVO) Analysis 133
Fig. 5.13 Analysis of AVO anomalies at around 2353 ms. The red line on the seismic display shows
the time location at which the amplitudes have been extracted. Those amplitudes are plotted as red
squares on the right-hand graph. The curve which has fit through the picks is a plot of Aki-Richards
two-term equation
Figures 5.13, 5.14 and 5.15 show the curves of the amplitude variations with
incidence angles for several CMP gather selected from the dataset. These amplitude
curves are utilized to estimate the intercept (A) and gradient (B) values. We would
like to analyze the AVO anomaly at around 1626, 2353 and 2365 ms. The anomaly
plot is shown in Fig. 5.13. The red line on the seismic display shows the time location
at which the amplitudes have been extracted. Those amplitudes are plotted as red
squares on the right-hand graph. The curve which has fit through the picks is a plot
of Aki-Richards two-term equation. We can confirm this by the information at the
top of the graph.
Notice that we are seeing a classic class 2p AVO anomaly with amplitudes increas-
ing for both the trough at the top of the sand (red) and the peak at the base of the
sand (green). Notice also that the fit of the AVO curves is good. Mathematically, this
is expressed by the normalized correlation between the picked amplitudes and the
curves, printed at the top of the graph. The analysis for time 2353 ms and 2365 ms
are shown in Figs. 5.14, and 5.15 respectively. The graph shows class 2p anomaly
for 2353 ms time and class 4 anomaly for 2365 ms time respectively.
b. AVO attributes volume
AVO attributes volumes such as intercept (A), gradient (B) and their proportions
such as AVO product (AB), Poisson coefficient (A + B), shear reflectivity (A-B)
134 5 Amplitude Variation with Offset (AVO) Inversion
Fig. 5.14 Analysis of AVO anomalies at around 2363 ms. The red line on the seismic display shows
the time location at which the amplitudes have been extracted. Those amplitudes are plotted as red
squares on the right-hand graph. The curve which has fit through the picks is a plot of Aki-Richards
two-term equation
Fig. 5.15 Cross-section of intercept (inline 1163) estimated from Aki-Richards two-term equation
5.4 Amplitude Variation with Offset (AVO) Analysis 135
and fluid coefficient (FF) can be generated straight from pre-stack seismic NMO-
corrected CMP-sorted data. Now we use AVO Gradient Analysis to examine the AVO
anomaly, we will extend the calculation to the whole sample to see the allocation of
AVO anomalies. The method of the AVO Attribute utilizes the Aki-Richards equation
of two or three terms to obtain AVO attributes from the seismic data. The attributes are
based on intercept, gradient and curvature combinations as described by the equation
of Aki-Richards. The angle gathers volume is used as input. Different combinations
of AVO attributes can be generated as per requirement. Here, two basic types of
attributes namely AVO A and AVO B are generated.
The Aki-Richards approximation for two terms only is used because we only
have incident angles less than 30°. We need high angle data, usually exceeding 45°,
to reliably obtain three terms. The calculated attributes from the Penobscot field,
Canada are shown in Figs. 5.15 and 5.16 for AVO A and AVO B respectively. One
can notice the variation of intercept and gradient attributes and hence can utilize to
detect the hydrocarbon zone. Although, the present datasets did not have the presence
of any major hydrocarbon reservoir.
Cross-plotting of the intercept (A) against gradient (B) is an efficient interpre-
tation technique for identification of AVO anomalies (Fig. 5.17). This method was
developed by Castagna et al. (1998). It is based on two ideas: (1) the Rutherford and
Williams (1989) AVO classification scheme described below and (2) the so-called
Mudrock Line. The cross plot shows the anticipated context pattern through the ori-
gin, with anomalous occurrences in quadrants 1 and 3, corresponding with anomalies
Fig. 5.16 Cross-section of the gradient (inline 1163) estimated from the Aki-Richards two-term
equation
136 5 Amplitude Variation with Offset (AVO) Inversion
Fig. 5.17 Crossplot of intercept and gradient near to horizon 2 from Penobscot data, Canada. The
Class 3 anomaly is highlighted
in class 3 AVO (Fig. 5.17). The slices of a different combination of attributes are
plotted in Fig. 5.18. These slices show the nice variation of extracted attributes in
the subsurface.
iv. Shuey approximation
Another useful simplification of the Zoeppritz equation was proposed by Shuey
(1985), who decomposed the reflectivity into the normal-incidence term and correc-
tion terms principally depending on the Poisson’s ratio and density variations across
the boundary. The Shuey approximation can be written as follows.
σ 1 V P
R(θ ) ≈ R0 + A0 R0 + sin2 θ + (tan2 θ − sin2 θ ) (5.42)
(1 − σ )2 2 VP
V P
V P ρ
In Eq. 5.42, R0 = 1
+ ρ
, A0 = B − 2(1 + B) 1−2σ and B = V P
VP
ρ .
VP + ρ
2 VP 1−σ
The first term in Eq. 5.42 describes the amplitude at θ = 0, the second term rep-
resents an amplitude correction at intermediate angles, and the third term describes
the amplitude at wide angles. For a rock sample under unidirectional pressure, the
Poisson’s ratio σ is the ratio of the transverse expansion to the longitudinal compres-
sion, or the ratio of shear strain to principal strain (Yilmaz 2001). For isotropic rock,
σ can be expressed through the ratio of the P- and S-wave velocity as follows.
5.4 Amplitude Variation with Offset (AVO) Analysis 137
Fig. 5.18 Depicts slices from gradient (B), Intercept (A), A*B and Poisson’s ratio (From top to
bottom)
2
−2
VP
VS
σ = 2 (5.43)
2 VVPS − 2
Thus, the Poisson’s ratio increases when V P /VS increases, and vice versa, and
therefore this ratio is typically low for gas reservoirs (Ma and Morozov 2010). It
typically equals to 0.1, for gas sands. Changes in gas or fluid saturation change
the Poisson’s ratio significantly because of the changes in rock bulk modulus and
consequently the P-wave velocity. At the same time, the shear modulus changes
only slightly and therefore fluid saturation has little effect on the S-wave veloc-
ity (Gassmann 1951). Further, an increase in fluid saturation decreases the P-wave
reflectivity and hence the Poisson’s ratio.
v. Fatti approximation
An alternate form of the Aki and Richards’s equation was given by Fatti et al. (1994).
This expression can be written as follows.
R(θ ) = c1 R P + c2 R S + c3 R D (5.44)
138 5 Amplitude Variation with Offset (AVO) Inversion
2
In Eq. 5.44, c1 = 1 + tan2 θ, c2 = −8γ sin2 θ , γ = VP
VS
, and c3 =
− 21tan θ + 2γ sin θ where R P and R S are the P- and S-wave reflectivities which
2 2 2
Equation (5.44) allows us to calculate R P and R P from seismic data. The dif-
ference between the P-wave and S-wave reflectivities, (R P −R S ), can be used as an
indicator differentiating the shale over brine-sand and shale over gas-sand cases.
R P −R S values are negative for shale over gas-sand and always more negative in the
case of shale over brine-sand (Castagna and Smith 1994). The R P −R S tend to be
constant and near zero for non-pay reservoirs. The reflectivities R P and R S can also
be transformed into two new attributes i.e. the Fluid Factor (FF) and Lambda-Mu-
Rho (LMR). Figures 5.19, 5.20 and 5.21 shows the variation of R P , R S and cross
the plot between them. Figure 5.22 shows cross-section of R P −R S respectively. A
Fig. 5.19 Cross-section (inline 1165) of inverted R P attributes estimated from Fatti et al. (1994)
two-term equation
5.4 Amplitude Variation with Offset (AVO) Analysis 139
Fig. 5.20 Cross-section (inline 1165) of inverted R S attributes estimated from Fatti et al. (1994)
two-term equation
small anomaly zone can be seen from the cross plot (Fig. 5.21).
In a clastic sedimentary sequence, the Fluid Factor is defined so that it is high-
amplitude for reflectors that lie far from the mudrock line and low-amplitude for all
reflectors on the mudrock line. The equation defining the FF, according to Fatti et al.
(1994) can be written as follows.
140 5 Amplitude Variation with Offset (AVO) Inversion
Fig. 5.22 Cross-section of R P − R S at inline 1165. The anomalous zone is clearly visible
V P VS VS
F = − 1.16 (5.48)
VP V P VS
The FF equals zero when layers both above and below the reflecting boundary lie
on the mudrock line, such as shale over brine sand. By contrast, the FF is nonzero
when one of the layers lies on the mudrock line and the other one lies away from it
(Fatti et al. 1994). In cases of gas sands, the FF will be non-zero at both the top and
bottom of gas.
The Lambda-Mu-Rho attributes (LMR) are defined so that the Lame’s elastic
parameters λ and μ are combined with density ρ in the form of λρ and μρ, as was
first proposed by Goodway et al. (1997). Pre-stack seismic CMP gathers are inverted
to extract the P-impedance and S-impedance, and from these impedances, the λρ
and μρ products are extracted. Starting from the equations for wave velocities, we
have the followings relation.
λ + 2μ
VP = (5.50)
ρ
5.4 Amplitude Variation with Offset (AVO) Analysis 141
2μ
VS = (5.51)
ρ
We have
And therefore
λρ = Z 2P − 2Z S2 (5.54)
References
Aki K, Richards PG (1980) Quantitative seismology. W.H. Freeman and Company, San Francisco,
CA
Aki K, Richards PG (2002) Quantitative seismology, 2nd edn. University Science Books, Sausalito,
CA
Avseth P, Bachrach R (2005) Seismic properties of unconsolidated sands: tangential stiffness, Vp/Vs
ratios and diagenesis. In: SEG Technical Program Expanded Abstracts, Society of Exploration
Geophysicists, pp 1473–1476
Backus MM, Chen RL (1975) Flat spot exploration. Geophys Prospect 23(3):533–577
Bacon M, Simm R, Redshaw T (2007) 3-D seismic interpretation. Cambridge University Press,
Cambridge
Batzle M, Wang Z (1992) Seismic properties of pore fluids. Geophysics 57(11):1396–1408
Brown AR (2012) Dim spots: opportunity for future hydrocarbon discoveries? Lead Edge
31(6):682–683
Castagna JP, Backus MM (eds) (1993) Offset-dependent reflectivity—theory and practice of AVO
analysis. Soc Explor Geophys
Castagna JP, Smith SW (1994) Comparison of AVO indicators: a modeling study. Geophysics
59(12):1849–1855
Castagna JP, Swan HW (1997) Principles of AVO cross plotting. Lead Edge 16(4):337–344
Castagna JP, Swan HW, Foster DJ (1998) Framework for AVO gradient and intercept interpretation.
Geophysics 63(3):948–956
Chen Q, Sidney S (1997) Seismic attribute technology for reservoir forecasting and monitoring.
Lead Edge 16(5):445–448
Chiburis E, Leaney S, Skidmore C, Franck C, McHugo S (1993) Hydrocarbon detection with AVO.
Oilfield Rev 5:1–3
Coulombe AC (1993) Amplitude versus offset analysis using vertical seismic profiling and well log
data. The crewes project, consortium for research in elastic wave exploration seismology, 167p
Dey-Sarkar SK, Foster DJ, Smith SW, Swan HW, Atlantic Richfield Co. (1995) Seismic data
hydrocarbon indicator. U.S. Patent 5: 440,525
142 5 Amplitude Variation with Offset (AVO) Inversion
Dufour J, Goodway B, Shook I, Edmunds A et al (1998) AVO analysis to extract rock parameters
on the Blackfoot 3c-3d seismic data. In: 68th Annual International SEG Mtg, pp 174–177
Dufour J, Squires J, Goodway WN, Edmunds A, Shook I (2002) Integrated geological and geo-
physical interpretation case study, and lame rock parameter extractions using AVO analysis on
the Blackfoot 3C-3D seismic data, southern Alberta, Canada. Geophysics 67(1):27–37
Fatti JL, Smith GC, Vail PJ, Strauss PJ, Levitt PR (1994) Detection of gas in sandstone reser-
voirs using AVO analysis: a 3-d seismic case history using the geo-stack technique. Geophysics
59(9):1362–1376
Gardner GHF, Gardner LW, Gregory AR (1974) Formation velocity and density —the diagnostic
basics for stratigraphic traps. Geophysics 39:770–780
Gassmann F (1951) Elastic waves through a packing of spheres. Geophysics 16(4):673–685
Goodway B, Chen T, Downton J (1997) Improved AVO fluid detection and lithology discrimination
using Lamé petrophysical parameters; “λρ”, “μρ”, & “λ/μ fluid stack”, from P and S inversions.
In: SEG Technical Program Expanded Abstracts, Society of Exploration Geophysicists, pp 183–
186
Hammond AL (1974) Bright spot: better seismological indicators of gas and oil. Science
185(4150):515–517
Han DH, Batzle ML (2004) Gassmann’s equation and fluid-saturation effects on seismic velocities.
Geophysics 69(2):398–405
Hashin Z, Shtrikman S (1963) A variational approach to the theory of the elastic behavior of
multiphase materials. J Mech Phys Solids 11(2):127–140
Hill R (1952) The elastic behavior of a crystalline aggregate. Proc Phys Soc. Sect A 65(5):349
Inyang CB (2009) AVO analysis and impedance inversion for fluid prediction in Hoover Field, Gulf
of Mexico, Doctoral dissertation, University of Houston
Keys RG (1989) Polarity reversals in reflections from layered media. Geophysics 54(7):900–905
Knott CG (1899) Reflexion and refraction of elastic waves, with seismological applications. Lond
Edinburgh Dublin Philosop Maga J Sci 48(290):64–97
Koefoed O (1955) On the effect of Poisson’s ratios of rock strata on the reflection coefficients of
plane waves. Geophys Prospect 3(4):381–387
Koefoed O (1962) Reflection and transmission coefficients for plane longitudinal incident waves.
Geophys Prospect 10(3):304–351
Kumar D (2006) A tutorial on Gassmann fluid substitution: formulation, algorithm and Matlab
code. Matrix 2:1
Larsen JA, Margrave GF, Lu HX (1999) AVO analysis by simultaneous PP and PS weighted stack-
ing applied to 3C-3D seismic data. In: SEG Technical Progress Expanded Abstract, Society of
Exploration Geophysicists, pp 721–724
Li XY, Mueller MC (1997) Case studies of multicomponent seismic data for fracture characteriza-
tion: Austin Chalk examples. In: Carbonate seismology, Society of Exploration Geophysicists,
pp 337–372
Li Y, Downton J, Xu Y (2007) Practical aspects of AVO modeling. Lead Edge 26(3):295–311
Luping G (2008) Pre-stack inversion and direct hydrocarbon indicator. Geophys Prospect Petroleum,
3
Ma J, Morozov I (2010) AVO modeling of pressure-saturation effects in Weyburn CO2 sequestration.
Lead Edge 29(2):178–183
Maurya SP, Singh NP (2019) Seismic modelling of CO2 fluid substitution in a sandstone reservoir:
a case study from Alberta, Canada. J Earth Syst Sci 128(8):236
Mavko G, Mukerji T (1998) Bounds on low-frequency seismic velocities in partially saturated
rocks. Geophysics 63(3):918–924
Mocnik A (2011) Processing and analysis of seismic reflection data for hydrocarbon exploration in
the plio-quaternary marine sediments, Doctoral thesis, University of Trieste
Muskat M, Meres MW (1940) Reflection and transmission coefficients for plane waves in elastic
media. Geophysics 5(2):115–148
References 143
Ostrander W (1984) Plane-wave reflection coefficients for gas sands at nonnormal angles of
incidence. Geophysics 49(10):1637–1648
Pendrel J, Dickson T (2003) Simultaneous AVO Inversion to P Impedance and Vp/Vs. In: CSEG
Annual Meeting, Expanded Abstract
Ramos AC, Davis TL (1997) 3-D AVO analysis and modeling applied to fracture detection in
coalbed methane reservoirs. Geophysics 62(6):1683–1695
Ross CP, Kinman DL (1995) Nonbright-spot AVO: two examples. Geophysics 60(5):1398–1408
Russell B (1988) Introduction to seismic inversion methods. The SEG Course Notes, Series 2
Russell H (1999) Theory of the STRATA Program. Hampson-Russell, CGG Veritas
Rutherford SR, Williams RH (1989) Amplitude-versus-offset variations in gas sands. Geophysics
54(6):680–688
Sheriff, R.E., 1973. Encyclopedic dictionary of exploration geophysics: Tulsa. Society of Explo-
ration Geophysicists
Shuey R (1985) A simplification of the zoeppritz equations. Geophysics 50(4):609–614
Simmons JL, Backus MM (1996) Waveform-based AVO inversion and AVO prediction-error.
Geophysics 61(6):1575–1588
Smith GC, Gidlow PM (1987) Weighted stacking for rock property estimation and detection of gas.
Geophys Prospect 35(9):993–1014
Smith GC, Sutherland RA (1996) The fluid factor as an AVO indicator. Geophysics 61(5):1425–1428
Smith TM, Sondergeld CH, Rai CS (2003) Gassmann fluid substitutions: A tutorial. Geophysics
68(2):430–440
Swan HW, Castagna JP, Backus MM (1993) Properties of direct AVO hydrocarbon indicators. Offset
dependent reflectivity-theory and practice of AVO anomalies. Invest Geophys 8:78–92
Tai S, Puryear C, Castagna JP (2009) Local frequency as a direct hydrocarbon indicator. In: SEG
Technical Program Expanded Abstracts, Society of Exploration Geophysicists, pp 2160–2164
Verm R, Hilterman F (1995) Lithology color-coded seismic sections: the calibration of AVO cross
plotting to rock properties. Lead Edge 14(8):847–853
Wang Z (2001) Fundamentals of seismic rock physics. Geophysics 66(2):398–412
White RS (1977) Seismic bright spots in the Gulf of Oman. Earth Planet Sci Lett 37(1):29–37
Wiggins R, Kenny GS, McClure CD (1983) A method for determining and displaying the shear-
velocity reflectivities of a geologic formation. European patent application, 113944
Xu Y, Bancroft JC (1997) Joint AVO analysis of PP and PS seismic data. The CREWES Project
Research Report, 9
Yilmaz O (2001) Seismic data analysis, vol 1. Society of Exploration Geophysicists, Tulsa, OK
Zhu X, McMechan GA (1990) Direct estimation of the bulk modulus of the frame in a fluid-
saturated elastic medium by Biot theory. In: SEG Technical Program Expanded Abstracts, Society
of Exploration Geophysicists, pp 787–790
Zoeppritz K (1919) On the reflection and propagation of seismic waves. Gottinger Nachrichten
1(5):66–84
Chapter 6
Optimization Methods for Nonlinear
Problems
6.1 Introduction
Fig. 6.1 Schematic diagram showing the process of solving nonlinear problems using an
optimization technique
minimize cost and time and hence to decide better transportation service one can use
to move from one place to the other place. These types of problems can be solved as
follows.
First, the physical problems are transformed into mathematical modeling using
the information provided in the problem. Thereafter, this mathematical modeling is
transformed into mathematical formulation according to the optimization technique
using objective function and constraints. Further, the solution (output) is generated
by solving some kind of optimization technique. Further, the decision is made to
decide either the solution is optimal or not. Graphically, it can be represented as
follows (Fig. 6.1).
One of the geophysical inversion’s main objectives is to discover earth mod-
els that explain geophysical findings. Thus, in many geophysical applications, the
mathematics branch known as optimization has discovered important use. In this
sense, geophysical inversion includes finding the ideal value of a multiple variables
function.
The feature we want to minimize (or maximize) is a misfit (or fitness) func-
tion that characterizes the distinctions (or similarities) between observed and syn-
thetic information calculated using an assumed earth model. The earth model is
depicted by physical parameters that characterize rock layers’ characteristics, such
as compressional wave velocity, shear-wave velocity, resistivity, etc.
In the estimation of material characteristics from geophysical information, both
local and global optimization techniques are used. The objective of this chapter is
to explain the implementation of geophysical problems of several newly developed
techniques of local and global optimizations. Although we emphasize the implemen-
tation elements of these algorithms, we explain in sufficient detail several compo-
nents of the theory for readers to know the underlying basic values on which these
algorithms are based.
There are many optimization techniques, and they are selected on the basis of
objective function and constraints provided in the problem. For many geophysical
applications, the misfit surface can be extremely complex and distinguished by var-
ious mountains and valleys as a function of the model parameters outlined by the
mismatch between the expected and observed geophysical information (Jervis et al.
1996). Such a feature will, therefore, have several minima and maxima; the mini-
mum of all minima will be called the global minimum, and all other minima will be
called local minima. Note that the global minimum is one of the local minima, but
the opposite is not true, and multiple minima of nearly the same depth may also be
6.1 Introduction 147
available (Sen and Stoffa 2013). They use the misfit function’s local characteristics
to calculate an update to the present response and search downhill direction. Thus,
if the starting point is closer to one of the local minima than the global minimum,
these algorithms will miss the global minimum solution. Geophysicists have been
plagued by the local minimum syndrome for over a century now.
There are two broad categories available for the optimizations, the first is local
optimization methods and the second is called global optimization methods. The local
optimization method generally looks for a local minimum solution in the vicinity of
a startup model. On the other hand, the global optimization methods have a tendency
to reach the global minimum by jumping from the local minima (Sen and Stoffa
2013). These local and global optimization methods are further divided into several
subparts. In the following sections, they are explained briefly. Before discussing these
optimization methods, we are in a position to discuss different kinds of optimization
methods that are used by both methods.
One can characterize the fitness function as continuous or discrete. When variables
can take only finite numbers in the fitness function, they are called discrete problems
of optimization. On the other hand, if variables in the error function can take a
continuous set of real values, the continuous optimization problems arise (Dorrington
and Link 2004; Sen and Stoffa 2013). We have provided five types of fitness function
produced by alteration in the l1 and l2 norm. These fitness functions are mostly used
in all geophysical problems. They are as follows.
The first fitness function we are presenting here is based on L 2 norm. This L 2
norm is applied between observed and modeled data (Moncayo et al. 2012). It can
be represented as follows.
n
1 2
E= Sobji − Smodi (6.1)
n i=1
In Eq. 6.1, E is the error between modeled and observed data, Sobji is observed
data at ith sample, Smodi is modeled data at ith sample and n is the total number of
sample points in the data. This method is largely used in geophysics because of fast
convergence.
The second fitness equation is taken from the mean square error between modeled
and observed data (Pullammanappallil and Louie 1994). The mathematical formulae
can be expressed as follows.
1 2
n
E= Sobji − Smodi (6.2)
n i=1
148 6 Optimization Methods for Nonlinear Problems
The fitness function presented in Eq. 6.2 takes a relatively larger time to converge
as compared to the first fitness function.
The third fitness equation is borrowed from the Misra and Sacchi (2008) which
can be expressed as follows.
n
E= Sobji − Smodi 2 (6.3)
i=1
This fitness equation is very time consuming and hence before using them, proper
care should be given.
The fourth fitness equation presented here is based on l1 norm. This is also called
the least absolute deviation. This is calculated between observed and modeled data
sets sample by sample. In order to constrain the solution, one additional term is
added to the equation and this additional term is called a priori impedance model.
Mathematically this fitness function can be expressed as follows.
n
n
i=1 Sobji − Smodi i=1 Z obji − Z modi
E = W1
n + W2
(6.4)
Sobj
i=1
n Z obj
i=1
In Eq. 6.4, W1 and W2 are weights applied to the two terms, respectively. The
weighting factor can be chosen as W1 = W2 = 1 in most cases (Ma 2002). The Z prii
is the prior low-frequency impedance at ith time sample and Z modi is the model
impedance at ith time sample. The last (fifth) fitness equation is again based on root
mean square error (l2 norm). In this equation, one additional term is added which
comes from the initial impedance model. This additional term act as constraints for
the solution (Ma 2002). Using this equation the non-uniqueness of the solution can
be reduced to some level. The expression can be written as follows.
n n
1 2 1 2
E= Sobji − Smodi + Z prii − Z modi (6.5)
n i=1 n i=1
In Eq. 6.5, the Z prii is the prior low-frequency impedance at ith time sample and
Z modi is the modeled impedance at ith time sample. The convergence of this fitness
function (Eq. 6.5) is very fast.
The local optimization algorithms move from solution to solution in the search space
of candidate solutions by applying local changes, until a solution deemed optimal is
found. Most of the local optimization methods are solved in an iterative manner, and
the objective of these algorithms is to move towards the solution in each iteration
6.3 Local Optimization Methods 149
and hence guarantee a reduction in the value of the objective function (Bremermann
1958, 1962). The local optimization methods get stuck at the local minimum in the
vicinity of the initial model and therefore, the choice of the initial model is of utmost
importance. If a choice of initial guess model will be made away from the global
solution then, the algorithm will be trapped in a local minimum. The basic working
principle of local optimization methods is presented in Fig. 6.2.
The local optimization methods start with a given startup/initial model. Thereafter,
objective function is employed to estimate error and if this error is small enough,
then the initial guess model will be treated as a final solution. If the error is not
small enough, the algorithm calculates a search address using a local property of
the objective function, which determines an update or increment of the current one
model (Hejazi et al. 2013). There are many local optimization methods available but
some of them are as follows.
The steepest descent method is an algorithm used to find the nearest local minimum
of an inverse problem which presupposes that the gradient of the inverse problem
can be computed. This method is also called the gradient descent method (Sen and
Stoffa 2013; Deift and Zhou 1993). Mathematically it can be expressed as follows.
Let’s a function f : Rn → R is such that the function f is differentiable at
some point x0 , and then the direction of steepest descent is the vector which can be
estimated as −∇ f (x0 ). Now, consider the function that can be written as follows.
In Eq. 6.6, u is a unit vector; that is, ||u|| = 1. Now, partial differentiation is applied
to Eq. 6.6, and we have the following equation.
∂ f ∂ x1 ∂ f ∂ xn
ϕ (t) = + ... + (6.7)
∂ x1 ∂t ∂ xn ∂t
∂f ∂f
ϕ (t) = u1 + . . . + un
∂ x1 ∂ xn
And therefore
where θ is the angle between ∇ f (x0 ) and u. It follows that ϕ (0) is minimized when
θ = π , which yields the following equation.
∇ f (x0 )
u= , ϕ (0) = −∇ f (x0 ) (6.9)
∇ f (x0 )
From Eq. 6.9, we are noticed that the problem is reduced to the single variable
from the several variables. Now, we want to find the minimum of ϕ(t) for t > 0,
x1 = x0 − t0 ∇ f (x0 ) (6.11)
and continue the process, by searching from x1 in the direction of −∇ f (x0 ) to obtain
x2 by minimizing ϕ1 (t) = f (x1 − t∇ f (x1 )), and so on.
This is the Method of Steepest descent, given an initial guess x0 , the method
computes a sequence of iterates {xk }, where
The conjugate gradient method is again local optimization methods used to find the
optimum solution to the inverse problems. It is a mathematical technique used to
solve the linear and non-linear inverse problem in an iterative manner. Although, the
method can also be used as a direct search technique and hence produce a numerical
solution to the problem (Polyak 1969; Sen and Stoffa 2013). It is a very useful
algorithm to solve large dimensional linear and non-linear inverse problems. The
method is very popular because of the very fast convergence rate and also need less
space to store computed data (Shewchuk 1994). The inverse problems are solved by
using the conjugate gradient method in the following way.
First, compute synthetic data dn = g(m n ) and then compute the derivative of
data with respect to the current model given by G n (m n ). Thereafter, compute data
residual d = dn − dobs , and model residual m n = m n − m pr . Further, compute
the regularized objective function as follows.
1 T
E(m n ) = dn dn + m nT m n (6.14)
2
Thereafter, compute the direction of steepest ascent. This can be achieved as
follows.
γn = G nT dn + m n (6.15)
(γn − γn−1 )T γn
σn = (6.16)
γn−1
T
γn−1
152 6 Optimization Methods for Nonlinear Problems
In last, compute optimum step length μn using a linear search and update model
in the following way.
m n+1 = m n − μn ϕn (6.17)
These steps repeated until one finds an optimum solution to the inverse problem.
Newton method is again an iterative local optimization method used to find roots of
an inverse problem which is differentiable. The method is applied to the derivative f
of a twice twice-differentiable function f to find the roots of the derivative (Dennis
and Moré 1977; Sen and Stoffa 2013). The basic process of the Newton method is
to construct a sequence xn from an initial guess model x0 which have tendency to
converge towards some point x ∗ such that f (x ∗ ) = 0.
The second-order Taylor expansion f T (x) of f around xn gives the following
expression.
1
f T (x) = f T (xn + x) ≈ f (xn ) + f (xn )x + f (xn )x 2 (6.18)
2
Now, we want to find x such that xn +x is a stationary point. By differentiation
of Eq. 6.18, we will have the following equation.
d 1
0= f (xn ) + f (xn )x + f (xn )x 2 = f (xn ) + f (xn )x (6.19)
dx 2
f (xn )
x = − (6.20)
f (xn )
Equation 6.20 is the solution to the equation. In this solution, it is hoped that
xn+1 = xn + x = xn − f (xn )/ f (xn ) will move closer to a stationary point x ∗
(Kelley 1999).
Many inverse geophysical problems are nonlinear optimization problems. Local
optimization techniques, e.g. generalized linear inversion, steepest descent, etc.,
therefore do not provide a satisfactory solution as they usually converge to a local
minimum, based on the starting model’s selection (Broyden 1967). Therefore, global
methods of optimization such as genetic algorithms or simulated annealing are
appropriate to solve these issues.
6.4 Global Optimization Methods 153
The objective of the global optimization method is to find the optimum solution glob-
ally of the inverse problem in the presence of many local solutions. Recently, global
optimization methods are routinely used to solve the geophysical inverse problem
because of the availability of advanced computing systems and large memory. The
global optimization methods attempt to find the global optimum solution of the misfit
function rather than finding a local optimum solutions like local optimization meth-
ods (Sen et al. 1995; Kelley 1999). A solution is called global optimal if there are no
other feasible solutions available with better objective function values. On the other
hand, a solution is termed as a locally optimal solution if there is no other solution
with better objective function value is available in the vicinity (Kelley 1999). These
objective functions play a very important role to optimize the problem.
The main advantages of using global optimization methods over local optimization
methods are not dependent on the initial guess model. There are many global opti-
mization techniques are available but some important algorithm which is regularly
used in geophysical problems are as follows.
1. Genetic algorithm (GA).
2. Simulated annealing (SA).
These methods are discussed briefly in the following sections. These methods are
largely used in geophysical problems hence they are applied, first, synthetic data sets
and then to the real data sets to estimate subsurface rock property. Figure 6.3 depicts
flowchart of global optimization methods.
create randomness to the individual solution (Sivanandam and Deepa 2007; Morgan
et al. 2012). These methods are explained briefly in the following subsections.
i. Selection
The selection is the first genetic operation carried out in the population after determin-
ing the fitness of each individual model. Based on their fitness values, the selection
method is used to pair individual models. The model’s fitness value is used to deter-
mine whether or not they are being chosen. Compared to the low fitness value, the
models with very high value have an enormous opportunity to be chosen. This is
because the fitter models are more frequently selected for reproduction. The proba-
bility of selection is directly associated with the model’s fitness value (Goldberg and
Holland 1988). There are six major selection types available and they are as follows.
• Stochastic uniform
• Remainder
6.4 Global Optimization Methods 155
• Roulette
• Uniform
• Tournament
• Rank.
Stochastic uniform
Stochastic uniform lays out a line in which each parent corresponds to a section of
the line of length proportional to its expectation. The algorithm moves along the line
in steps of equal size, one step for each parent. At each step, the algorithm allocates
a parent from the section it lands on. The first step is a uniform random number less
than the step size (Tran and Hiltunen 2012).
Remainder
Remainder assigns parents deterministically from the integer part of each individ-
ual’s scaled value and then uses roulette selection on the remaining fractional part
(Velez 2005).
Roulette
Roulette simulates a roulette wheel with the area of each segment proportional
to its expectation. The algorithm then uses a random number to select one of the
sections with a probability equal to its area (Yang and Honavar 1998).
Uniform
Uniform functions select parents at random from a uniform distribution using the
expectations and number of parents. This results in an undirected search. Uniform
selection is not a useful search strategy, but you can use it to test the genetic algorithm
(Yang and Honavar 1998).
Tournament
Tournament selects each parent by choosing individuals at random, the number
of which you can specify by Tournament size, and then choosing the best individual
out of that set to be a parent. When Constraint parameters > Nonlinear constraint
algorithm is Penalty, GA uses Tournament with size 2 (Yan-Jun et al. 2016).
Rank selection
The fitness values of the models are evaluated in the rank selection method (Baker
1987; Whitley 1989), and then they are sorted to assign a rank to each model. A
model’s rank can range from 0 (for the best model) to n − 1 (for the worst model).
The selection is done in such a way that the best model in the population contributes a
number of copies that are an integral number of copies that the worst model receives
(the number is predetermined). Thus the choice of ranks fundamentally exaggerates
the distinction between models with fitness values that are almost identical (Yuan
and Zhuang 1996).
ii. Crossover
The crossover is the next genetic operator to be used after selecting and pairing the
models. The fundamental crossover principle is the sharing of genetic information
among paired models. In this manner, it can also be understood that the crossover
exchanges some data between paired models and thus generates fresh models (Sen
and Stoffa 2013). In geophysical problems, there are two kinds of crossover widely
156 6 Optimization Methods for Nonlinear Problems
used namely single-point crossover and multipoint crossover. The fundamental idea
of the single point crossover is that with the uniform probability distribution, the
model’s one-bit location is chosen at random. Then all the bits falling on the right
side of the selected bit will be exchanged between two models, producing a new
model. On the other hand, a bit position is randomly selected in the multipoint
crossover process and all the bits right to this selected bit are exchanged among the
paired models. In addition, another bit location is selected at random from the second
model parameter and all other bits falling straight to this bit are swapped again. For
each model parameter, this method is repeated (Goldberg and Holland 1988).
Figures 6.4 show how single- and multipoint crossovers are performed. In single-
point crossover (Fig. 6.4), one bit is selected at random and all the bits between
model m i and model m j after the crossover points are exchanged with each other.
This results in the generation of two new models. On the other hand the multipoint
crossover exchanges binary information between each model parameter by choosing
a crossover point independently and performing crossover on a model-parameter-
by-model-parameter basis (Yuan and Zhuang 1996).
Once a crossover site has been randomly selected, it is decided on the basis of
the probability of crossover whether or not crossover will be performed; i.e. the
crossover rate is controlled by a probability px specified by the GA designer. A
high probability of crossover implies that crossover between the paired models or
between the present model parameter of the paired models is likely to happen in the
event of multipoint crossover (Boschetti et al. 1996; Sen and Stoffa 2013; Maurya
et al. 2017).
iii. Mutation
The last genetic operator is mutation. The mutation is a process used to create ran-
domness to the crossover. Generally, the process of mutation is performed during the
crossover process. The mutation rate which is used to decide the number of walks
in the model space is specified by a probability that is determined by the algorithm
designer. The low mutation probability means fewer walks in the model space and
the conversion of the problem will be very fast (Holland 1975; Deb et al. 2002). On
the other hand, the high mutation probability will result in a large number of random
walks in the space but this can take more time to converge the algorithm (Mallick
1995; Aleardi et al. 2016). Figure 6.5 shows the basic concept of the mutation. One
new solution (children) is selected from the crossover and treated as parents to per-
form mutation. For this purpose, two mutation points are selected and interchanges
with each other to produce randomness to the solution (Mallick 1995).
The inversion based on genetic algorithm does not give an absolute subsurface
impedance model, it gives a relative variation of acoustic impedance. The practical
process of the genetic algorithm is as follows.
Mathematical example
The process can understand in a specific way with the examples. Let us consider a
function
cos(3π x)
f (x) = , wherex ∈ [0.11.1] (6.21)
x
This function graphically can be presented as follows (Fig. 6.6).
The function has a local and global maximum and local and global minimums
and they are highlighted. The function has a global minimum fall at x = 0.30 and the
corresponding function value is −3.1705. To demonstrate the capability of a genetic
algorithm to find the solution to the problem, the method is applied to Eq. 6.21. The
estimated results are as follows.
x = 0.296878161902085
And
f (x) = −3.17151688802542
This is very close to the real solution and describes the high performance of the
algorithm.
158 6 Optimization Methods for Nonlinear Problems
Synthetic example
The above-discussed example is presented in a mathematical case although the
GA largely used in seismic inversion purpose. Hence, an example is present
to demonstrate the use of GA for seismic inversion. In this regard, the param-
eters are optimized by applying the first synthetic data and then real com-
posite trace. A seven-layer earth model is assumed with acoustic impedance
8050, 10,000, 7260, 9120 and 11,440, 7820, 10,250 m/s ∗ g/cc. The result for a
seven-layer Earth’s model is illustrated in Fig. 6.7. The first track of Fig. 6.7 shows
assumed subsurface seven-layer earth model, second track depicts corresponding
synthetic seismogram generated by the convolution is 40 Hz Ricker wavelet with
earth reflectivity, the same trace shown 8 times to looks like data and third track
depicts comparison of modeled (dot line) impedance with observed impedance (solid
line). From the figure, it can be examined that the inverted acoustic impedance is
very close to the modeled acoustic impedance and depicts the good performance of
the algorithm. To quality check of inverted results, crossplot between inverted and
original data is generated and shown in Fig. 6.8. The red solid line in the figure rep-
resents best-fit line. It can be noticed that the scatter points are very close to the best
fit line and hence indicates that the inverted results are in the vicinity of real value
and indicates the good performance of the algorithm. The correlation is estimated to
be 0.96 between modeled and observed impedance which is very high and supports
the high performance of the algorithm.
Figure 6.9 shows visualization of error minimization. One can notice that the error
decreases exponentially with the iteration and reached to minimum value of 0.0069
after 100 iterations. The termination criterion is set to be 100 iterations by considering
6.4 Global Optimization Methods 159
Fig. 6.7 Demonstrate seismic inversion results for synthetic data case whereas a represents a geo-
logical model, b depicts synthetic trace generated over geological model and c compares inversion
results with the original value
synthetic data case. The synthetic example contains 7 layers only and hence one can
compare layer by layer the inversion results. This comparison is not easy for real
data cases as the subsurface contains a large number of layers. Table 6.1 shows a
quantitative comparison of inverted results with original value for each layer. From
the table, one can notice that the modeled impedance is quite close to the observed
impedance. The change percentage between observed and modeled data points is
also estimated and found very low (<5.967% average).
Real data example
Now, we are in a position to move towards real data example as the mathematical
and synthetic example shows satisfactory results. These analysis always important as
the global optimization methods implemented by many parameters and hence they
become very crucial to take input. One other reason is that this global optimization
160 6 Optimization Methods for Nonlinear Problems
Fig. 6.9 Visualization of errors for the inversion of synthetic data with genetic algorithm
Table 6.1 Comparison of observed and modeled AI for synthetic data case
S. No. Layers Observed AI Modeled AI Change (%)
1. Layer 1 8050 8300.93 3.12
2. Layer 2 10,000 9162.50 8.38
3. Layer 3 7260 8310.48 14.47
4. Layer 4 9120 9393.93 3.00
5. Layer 5 11,440 10,817.03 5.45
6. Layer 6 7820 8261.50 5.65
7. Layer 7 10,250 10,423.94 1.70
is much time taking and hence one needs to optimize parameters for small datasets
before moving for large data sets. The real data is taken from the Blackfoot field,
Alberta, Canada. The GA is applied to estimate acoustic impedance in the inter-well
region. The real data example is divided into two parts, first, one composite trace
near to well 08-08 is extracted and GA is performed and in the second part, the entire
seismic section is transformed into the AI using GA technique.
For real data case, a real coded GA is implemented where an individual is repre-
sented by an array of real numbers and these are velocity, density, and thickness for
each layer. In the first step, define a search width for each of the model parameters,
which is obtained from a priori information. The lower bound is decided to be used
2250 m/s * g/cc and the upper bound is 21,000 m/s * g/cc for impedance which
is a reasonable range for seismic reflection data. GA starts with the generation of
a random population of 200 chromosomes. Using forward modeling procedure the
6.4 Global Optimization Methods 161
synthetic seismograms are generated for each individual chromosomes and corre-
sponding fitness value is evaluated. On the basis of these fitness values, two best
solutions/chromosomes have been selected to perform the genetic operation.
Figure 6.10 shows comparison of inverted impedance and real impedance from
well log. Track 1 shows extracted seismic trace, the same trace plotted 8 times to
look like data, track 2 shows extracted wavelet from well log data which is used for
forwarding modeling and track 3 shows comparison of model, original and inverted
impedances. From track 3 (Fig. 6.10), one can notice that the inverted curves are
very close to the original curves. The correlation between them is noticed to be 0.88
which is quite good by considering noise in the data.
Again, for the quality check of the inversion results, a cross plot between inverted
and original impedances is generated and displays in Fig. 6.11. A best-fit line is fitted
to the data. Now, one can observe that the distribution of scatter points lies close to
the best fit line which indicates that the inverted results are close to the original
one. The statistical analyses of inverted results are presented in Table 6.2. The table
shows the minimum and maximum AI estimated for composite trace. The variation
of error with the iteration is shown in Fig. 6.12. The stopping criterion is set to be
900 iterations. One can increase the iteration number to get more close to the real
solution but here we have decided to use 900 iterations as a stopping criterion to
Fig. 6.10 Inversion results for the composite trace from the Blackfoot field, where a represents seis-
mic trace (one trace plotted 8 times), b shows extracted wavelet and c compares inverted impedance,
model impedance, and original impedance
162 6 Optimization Methods for Nonlinear Problems
optimize time and memory. From the error analysis, one can notice that it decreases
exponentially with iteration.
The analysis of composite trace shows satisfactory results and therefore, the GA
algorithm is applied to the entire seismic section to estimate AI in the inter-well
region. Figure 6.13 (top) shows a cross-section of seismic section which is used
6.4 Global Optimization Methods 163
Fig. 6.13 Cross-section of seismic data (top) from the Blackfoot field and inverted impedance
(bottom) derived using GA
as input for GA and Fig. 6.13 (bottom) shows cross-section of inverted acoustic
impedance (AI) at inline 27. The inverted section shows higher resolution as com-
pared with the seismic data and shows a clear boundary of subsurface layers. From
the inverted AI section, the thinner layer cannot be interpreted because of seismic
bandwidth and also because the algorithm does not include constraints from well log
164 6 Optimization Methods for Nonlinear Problems
information as other inversion methods used. The section shows that the area have
acoustic impedance variation from 8000 to 14,000 m/s * g/cc. The sand channel is
clearly identifiable with low AI between 1055 and 1065 ms TWT. The minimum
RMS error is achieved to be 0.12 and the mean error as 0.20 during the optimization.
One can include a low-frequency model to the inverted section to get information
on deeper reflector. The method is time-consuming and hence can be performed
on high performing computing. It is suggested that before the application of these
methods to the large datasets, the number of parameter values needs to be tested.
The Simulated Annealing (SA) is another type of global optimization method used
to estimate the global optimum solution of a function using temperature cooling
principle. The methods are very useful to the geophysical problem and some nice
examples are presented in Kirkpatrick et al. (1983), Geman and Geman (1987),
Davis (1987), Van Laarhoven and Aarts (1987), Aarts and Korst (1988), Sen and
Stoffa (1991) and Afanasiev et al. (2014).
The principle of simulated annealing is borrowed from statistical mechanics. The
method involves an analysis of a large number of atoms in samples of liquids or solids.
A physical annealing process occurs when the temperature of a solid is increased
by heating and hence all the particles are distributed randomly in a liquid phase
(Afanasiev et al. 2014). Thereafter, it is followed by slow cooling of the particle
temperature and they arrange themselves in the lower energy state and resultant
crystallization occurs (Sen and Stoffa 2013).
The SA algorithm attempts to sample the model space according to the Gibbs
distribution. The Gibbs function is largely controlled by the temperature of the model.
That’s why the simulate annealing always look towards the global minimum of an
error function without depending on the starting point (Sen and Stoffa 2013). The
Gibbs distribution can be written as follows.
exp − KETi 1 Ei
P(E i ) =
= exp − (6.22)
exp − j
E Z (T ) KT
j∈S KT
where the model state is i, the energy of the medium is represented as E i , the set S
consists of all possible configurations, K is Boltzmann’s constant, T is the tempera-
ture, and Z(T ) is the partition function/energy. This partition energy can be written
as follows.
Ej
Z (T ) = exp (6.23)
j∈S
KT
6.4 Global Optimization Methods 165
The simulated annealing algorithm starts with an initial temperature that is termed
as initial model m 0 , and the energy associated with this model can be expressed as
f (m 0 ). The SA generate new model, let say it is m n and the energy associated with
this new model can be expressed as f (m n ). Thereafter, compare the new energy
f (m n ) with the earlier energy state i.e. f (m 0 ). If the new energy state is lower than
that of the starting model, then the new model is accepted without any condition. On
the other hand, if the new energy state is larger than that of the initial model, then
the new model is accepted with the following probability (Ma 2002).
( f (m n ) − f (m 0 ))
p = exp − (6.24)
T
m j = m i + m i (6.25)
Now, lets the difference in energy between the initial model and newly generated
model is presented by E i j . This energy difference can be given by the following
equation.
E i j = E m j − (m i ) (6.26)
The new model is better or not from the previous one is decided based on the
value of E i j . If the energy difference is zero (E i j = 0), then the new model is
always get accepted. On the other hand, if the energy difference is greater than zero
(E i j > 0), then the new model gets accepted with the following probability.
166 6 Optimization Methods for Nonlinear Problems
E i j
P = exp − (6.27)
T
m imax − m imin
M= (6.28)
m i
From Eq. 6.28, it can be noted that the parameter M can take different values
for different model parameters. However, we will suppose that M is the same for
all model parameters without loss of generality. This results in a model space made
up of MN models. From this point in time, it is convenient to present the model
parameter by m i j where i = 1, 2, . . . , N and j = 1, 2, . . . , M.
6.4 Global Optimization Methods 167
The algorithm starts with an initial guess model m 0 and then each parameter of the
model is sequentially visited. The following marginal probability density function
(pdf) is assessed for each parameter of the model.
E (m̄|m i =m i j )
exp − T
P m̄|m i = m i j =
(6.29)
M E(m̄|m i =m ik )
k=1 − T
The sum is over M values permitted by the model parameter and E m̄|m i = m i j is
the energy of the model vector m̄ whose model parameter has a value m i j . Then a value
from the previous allocation is taken. Essentially, a random number is taken from
the uniform distribution U [0, 1] is mapped to the pdf given in Eq. 6.29. To perform
this, first, we need to compute the cumulative probability Ci j from P m̄|m i = m i j .
This can write as follows.
j
Ci j = P(m̄|m i = m ik ), J = 1, 2, . . . , M. (6.30)
k=1
Thereafter, one can draw a random number r from the uniform distribution. At
j = k, where Ci j = r, we select m i j = m ik = m i . The distribution is almost uniform
at high temperatures and therefore every parameter of the model is almost equally
probable to be selected. Only the peak corresponding to the model parameter that
most dominates the error function at very low temperatures.
The new value of the model parameter thus selected and replaces the old value of
the corresponding model parameter in the model vector, and for each model param-
eter, the process is repeated sequentially. After having examined each parameter of
the model once, we may have a model different from the starting model. This is one
iteration of the algorithm consisting of assessing the error function N × M times.
Unlike Metropolis SA, the heat bath does not involve testing for accepting a model
after it has been produced, but this algorithm involves significant computational cost.
When the problems have a large number of variables, it is faster than the Metropolis
SA (Vestergaard and Mosegaard 1991).
Fast simulated annealing (FSA)
Simulated annealing is a stochastic system for looking through the ground state. Fast
simulated annealing (FSA) is a semi-neighborhood search and comprises of intermit-
tent long bounces. SA’s Markov chain model can be used to prove that the asymptotic
convergence given by the Boltzmann distribution to a stationary distribution can be
achieved, and the global minimum energy state has a probability of 1 when the
temperature goes to zero. However, the speed at which the temperature is reduced
depends on whether or not the global minimum energy state is reached. It means
it relies on the cooling schedule. Geman and Geman (1984) demonstrated that the
following cooling plan provides a required and adequate condition for convergence
to the worldwide minimum for SA.
168 6 Optimization Methods for Nonlinear Problems
T0
T (k) = (6.31)
ln(k)
The search width is restricted as the lower limits (LB) are (7500, 9000, 6500,
8500, 11,000, 7000, 9500 m/s * g/cc) and the upper limits (UB) are (9500, 10,500,
7500, 1000, 12,000, 8000, 11,000 m/s * g/cc). The advantage of using LB and UB
is creating a window within that the search for optimum solution performed within
the interval to reduce time and space. The Simulated Annealing algorithm is used to
restore acoustic impedance in the synthetic datasets. The seismic inversion method
starts by using the original model parameters to calculate the objective function. From
either a random model or a macro model a starting model is traced. According to
Metropolis criteria, a fresh model is adopted or dismissed. Large numbers of iterations
are carried out at a specified temperature. Then the temperature will be reduced until
eventually the convergence criterion is met. The outputs at each sample interval are
the optimized acoustic impedance values.
Figure 6.14 shows inversion results for synthetic data. The first track of Fig. 6.14
displays the geological subsurface model used to generate synthetic data, the second
track demonstrates the respective synthetic traces, the same trace plotted 8 times to
look like data and the third track shows Simulated annealing reversal outcomes. The
inverted trace is very close to the original curves. The correlation between observed
and modeled impedance is noted to be 0.99, which is very high and depicts the good
performance of the algorithm. The estimated P-impedances can be seen to match
their real models quite nicely. From Fig. 6.14, track 3, it is evident that the original
impedance can get back quite nicely using a Simulated annealing technique.
Layer by layer comparison of inverted impedance is presented in Table 6.3. The
table demonstrates comparison of observed and modeled acoustic impedance esti-
mated by the simulated algorithm. From this comparison, it can be seen that the
modeled AI is very close to the observed AI for all 7 layers and hence the change
Fig. 6.14 Demonstrate seismic inversion results for synthetic data case whereas track 1 represents
the geological model, track 2 depicts synthetic trace generated over geological model and track 3
compares inversion results with the original value
170 6 Optimization Methods for Nonlinear Problems
Table 6.3 Comparison of observed and modeled AI for synthetic data case
S. No. Layers Observed AI Modeled AI Change (%)
1. Layer 1 8050 8436.99 4.81
2. Layer 2 10,000 10,767.11 7.67
3. Layer 3 7260 7569.52 4.26
4. Layer 4 9120 9290.67 1.87
5. Layer 5 11,440 11,819.47 3.32
6. Layer 6 7820 7382.39 5.6
7. Layer 7 10,250 9636.29 5.99
percentage is very low (<4.788, average). Figure 6.15 shows a cross-plot of modeled
and observed acoustic impedance. This cross plot shows visually how close both
results are. A best-fit line is fitted to the scatter data. From Fig. 6.15, it is visible that
the scatter points are very close to the best fit line and show excellent SA algorithm
efficiency. Figure 6.16 demonstrates annealing error assessment. The plots show that
the error decreases with iteration but exponentially as in the case of a genetic algo-
rithm. The stopping criterion is set to be 500 iterations. The simulated annealing is
comparatively faster than the genetic algorithm.
A comparison of results found in the genetic algorithm and simulated annealing
for synthetic data case has been performed. Figure 6.17 shows a cross plot of observed
and modeled acoustic impedance for both the case. It’s noticeable from Fig. 6.17 that
the scatter points are very near to the best fit line for both instances (GA and SA)
and show excellent algorithm efficiency. The close examination of the cross plot
indicates that the data points from simulated annealing are very close to the best
fit line whereas, for the genetic algorithm, the data points are scattered relatively.
The correlation coefficient for SA and GA is anticipated to be 0.987 and 0.937
respectively. This analysis indicates that the inversion results for synthetic data cases
are slightly better for simulated annealing as compared with the genetic algorithm.
This is not always the case, it largely depends on data parameters and input data.
6.4 Global Optimization Methods 171
Fig. 6.18 Inversion results for the composite trace from the Blackfoot field, where a represents seis-
mic trace (one trace plotted 8 times), b shows extracted wavelet and c compares inverted impedance,
model impedance, and original impedance
extracted wavelet used to generated modeled trace and track 3 demonstrates an esti-
mated inversion outcome using simulated annealing algorithm. By examining the
inverted outcomes, one can see that the inverted curve is quite nicely following the
original impedance from well log. The correlation between estimated and observed
AI is estimated to be 0.95 which is quite high and indicates good inversion results.
A cross plot between inverted impedance and the original impedance is produced
and displayed in Fig. 6.19 for quality check of inversion results. In addition, the best
fit line is fitted to the data. From the cross plot, one may notice that the scatter points
fall close to the best fit line and show excellent algorithm efficiency. The statistical
analysis of composite trace is given in Table 6.4.
Once acceptable outcomes have been obtained from the composite trace, the SA
technique is applied to the entire seismic segment to assess the acoustic impedance of
the subsurface in the inter-well region. Figure 6.20 (top) demonstrates seismic input
information (Inline 65) from the field of Blackfoot, Alberta, Canada and Fig. 6.20
(bottom) shows inverted results at inline 65. The inverted section demonstrates com-
paratively greater resolution and offers the layer property while seismic informa-
tion only offers the property of the interface. The impedance varies from 7500 to
13,500 m/s * g/cc. The inverted acoustic impedance section highlights sand channel
between 1055 and 1065 ms TWT.
6.4 Global Optimization Methods 173
Fig. 6.20 Cross-section of seismic data (top) from the Blackfoot field and inverted impedance
(bottom) derived using SA
174 6 Optimization Methods for Nonlinear Problems
References
Du Z, MacGregor LM (2010) Reservoir characterization from joint inversion of marine CSEM and
seismic AVA data using Genetic algorithms: a case study based on the Luva gas field. In: SEG
technical program expanded abstracts, society of exploration geophysicists, pp 737–741
Geman S, Geman D (1984) Stochastic relaxation, Gibbs distributions, and the Bayesian restoration
of images. IEEE transactions on pattern analysis and machine intelligence PAMI 6(6):721–741
Geman S, Geman D (1987) Stochastic relaxation, Gibbs distributions, and the Bayesian restoration
of images. In: Readings in computer vision. Morgan Kaufmann, pp 564–584
Gill PE, Murray W, Wright MH (1981) Practical optimization. Academic Press, London
Goffe WL, Ferrier GD, Rogers J (1994) Global optimization of statistical functions with simulated
annealing. J Econometrics 60(1–2):65–99
Goldberg DE, Holland JH (1988) Genetic algorithms and machine learning. Mach Learn 3(2):95–99
Hejazi F, Toloue I, Jaafar MS, Noorzaei J (2013) Optimization of earthquake energy dissipation
system by genetic algorithm. Comput Aided Civ Inf 28(10):796–810
Holland J (1975) Adaptation in natural and artificial systems: an introductory analysis with
application to biology. In: Control and artificial intelligence
Jervis M, Sen MK, Stoffa PL (1996) Prestack migration velocity estimation using nonlinear methods.
Geophysics 61(1):138–150
Kelley CT (1999) Iterative methods for optimization. Soc Ind Appl Math
Kirkpatrick S, Gelatt CD, Vecchi MP (1983) Optimization by simulated annealing. Science
220(4598):671–680
Ma XQ (2002) Simultaneous inversion of pre-stack seismic data for rock properties using simulated
annealing. Geophysics 67:1877–1885
Mallick S (1999) Some practical aspects of prestack waveform inversion using a genetic algorithm:
an example from the east Texas Woodbine gas sand. Geophysics 64(2):326–336
Mallick S (1995) Model-based inversion of amplitude-variations-with-offset data using a genetic
algorithm. Geophysics 60(4):939–954
Maurya SP, Singh KH, Kumar A, Singh NP (2017) Reservoir characterization using post-stack
seismic inversion techniques based on a real coded genetic algorithm. J Geophys 39(2):95–103
Metropolis N, Rosenbluth AW, Rosenbluth MN, Teller AH, Teller E (1953) Equation of state
calculations by fast computing machines. J Chem Phys 21(6):1087–1092
Misra S, Sacchi MD (2008) Global optimization with model-space preconditioning: application to
AVO inversion. Geophysics 73(5):R71–R82
Moncayo E, Tchegliakova N, Montes L (2012) Pre-stack seismic inversion based on a genetic
algorithm: a case from the Llanos Basin (Colombia) in the absence of well information. CT&F-
Ciencia Tecnología y Futuro 4(5):5–20
Morgan EC, Vanneste M, Lecomte I, Baise LG, Longva O, McAdoo B (2012) Estimation of free
gas saturation from seismic reflection surveys by the genetic algorithm inversion of a P-wave
attenuation model. Geophysics 77(4):R175–R187
Mosegaard K, Vestergaard PD (1991) A simulated annealing approach to seismic model optimiza-
tion with sparse prior information 1. Geophys Prosp 39(5):599–611
Polyak BT (1969) The conjugate gradient method in extremal problems. USSR Comput Math Math
Phys 9(4):94–112
Pullammanappallil SK, Louie JN (1994) A generalized simulated-annealing optimization for
inversion of first-arrival times. Bull Seismol Soc Am 84(5):1397–1409
Rebbi C (1987) Monte Carlo calculations in lattice gauge theories. Applications of the Monte Carlo
method in statistical physics. Springer, Berlin, pp 277–297
Sen MK, Stoffa PL (1991) Nonlinear one-dimensional seismic waveform inversion using simulated
annealing. Geophysics 56(10):1624–1638
Sen MK, Stoffa PL (2013) Global optimization methods in geophysical inversion. Cambridge
University Press, Cambridge
Sen MK, Datta-Gupta A, Stoffa PL, Lake LW, Pope GA (1995) Stochastic reservoir modeling using
simulated annealing and genetic algorithm. SPE Formation Eval 10(01):49–56
176 6 Optimization Methods for Nonlinear Problems
Shewchuk JR (1994) An introduction to the conjugate gradient method without the agonizing pain,
Technical report, Carnegie Mellon University
Sivanandam SN, Deepa SN (2007) Introduction to genetic algorithms. Springer Science and
Business Media
Szu H, Hartley R (1987) Fast simulated annealing. Phys Lett A 122(3–4):157–162
Tran KT, Hiltunen DR (2012) One-dimensional inversion of full waveforms using a genetic
algorithm. J Environ Eng Geophys 17(4):197–213
Van Laarhoven PJ, Aarts EH (1987) Simulated annealing. Simulated annealing: theory and
applications. Springer, Dordrecht, pp 7–15
Velez-Langs O (2005) Genetic algorithms in the oil industry: an overview. J Petrol Sci Eng 47(1):15–
22
Whitley LD (1989) The GENITOR algorithm and selection pressure: why rank-based allocation of
reproductive trials is best. Icga 89:116–123
Yang J, Honavar V (1998) Feature subset selection using a genetic algorithm. In: Feature extraction,
construction and selection. Springer, Berlin, pp 117–136
Yan-Jun H, Ding-Hui Y, Yuan-Feng C (2016) Reservoir parameter inversion of co2 geological
sequestration based on the self-adaptive hybrid genetic algorithm. Chin J Geophys Chin Ed
59(11):4234–4245
Yuan Y, Zhuang H (1996) A genetic algorithm for generating fuzzy classification rules. Fuzzy Sets
Syst 84(1):1–19
Chapter 7
Geostatistical Inversion
7.1 Introduction
© The Editor(s) (if applicable) and The Author(s), under exclusive licence 177
to Springer Nature Switzerland AG 2020
S. P. Maurya et al., Seismic Inversion Methods: A Practical Approach,
Springer Geophysics, https://doi.org/10.1007/978-3-030-45662-7_7
178 7 Geostatistical Inversion
Seismic attributes can be conveniently defined as the quantities that are measured,
computed or implied from the seismic data. From the time of their introduction in
the early 1970s, seismic attributes have gone along the way and they became an aid
for geoscientists for reservoir characterization and also as a tool for quality control
(Chambers and Yarus 2002; Chopra and Marfurt 2007). Seismic attributes are specific
measurements of geometric, kinematic, dynamic, or statistical features derived from
seismic data. Some are more sensitive than others to specific reservoir environments;
some are better at revealing subsurface anomalies that are not easily detectable, and
some have been used as direct hydrocarbon indicators (Taner et al. 1994; Chopra
and Marfurt 2005, 2007). Different authors introduced different kinds of attributes
and their uses. With the introduction of 3D seismic techniques and associated tech-
nologies and introduction of seismic sequence attributes, coherence technology in
the mid-1990s, and spectral decomposition in the late 1990s has changed the seis-
mic interpretation techniques and provided essential tools that were not available
for geoscientists earlier. With the introduction of 3D visualization techniques, the
use of seismic attributes has attained a new dimension. The development of a wide
variety of seismic attributes warrants a systematic classification. Also, a systematic
approach is needed to understand the use of each of these attributes and also their
limitations under different circumstances.
Though the purpose of this book is to understand different attributes that can be
used as tools in interpretation, it is useful to understand the classification of different
attributes at this stage. The more seismic attributes available, the more confusion geo-
scientists may have in selecting the appropriate ones. Improved documentation and
a greater understanding of the underlying rock physics that give rise to an anomaly,
the role that seismic processing plays in enhancing or destroying attribute anomalies,
and the value of seismic attributes in the absence of good well control are areas where
assistance is needed (Russell et al. 1997). Guidance should be available to the inter-
preter to help with the choice of the seismic attributes that will provide the greatest
insight into the current problem under investigation. The one way to classify seismic
attributes is based on their characteristics and hence geometrical and physical.
7.2 Seismic Attributes 181
Physical attributes are defined as those attributes which are directly related to the
wave propagation, lithology, and other parameters. These physical attributes can be
further classified as pre-stack and post-stack attributes. Each of these has sub-classes
as instantaneous and wavelet attributes. Instantaneous attributes are computed sam-
ple by sample and indicate a continuous change of attributes along the time and space
axis (Brown 1996; Taner 2001; Maurya et al. 2019). The Wavelet attributes, on the
other hand, represent characteristics of wavelet and their amplitude spectrum. The
magnitude of the trace envelope is proportional to the contrast between the acoustic
impedance; frequencies are related to the density of the bed, wave dispersion, and
absorption. Consequently, these attributes are mostly used for lithological classifi-
cation and reservoir characterization. The classification of physical attributes is as
follows.
(i) Post stack attributes
(ii) Pre-stack attributes
i. Post-stack attributes
Post stack attributes are derived from the stacked seismic data. Stacking is an average
method that removes data linked to offset and azimuth. It may be possible to stack or
migrate CDP input information. It should be noted that time migrated information
will retain their time relationships, thus temporal factors such as frequency will retain
their physical characteristics as well. Frequency is substituted by wave number for
depth migrated data, which is a function of velocity and frequency of propagation.
Post-stack characteristics are a more manageable strategy in original recognition
inquiries for looking big quantities of information.
182 7 Geostatistical Inversion
The Attributes are the result of the properties derived from the complex seismic
signal. The concept of complex traces was first described by Taner et al. (1979). The
complex trace is defined as follows.
where CT (t) represent complex trace, T (t) represent seismic trace, H(t) is Hilbert’s
transform of a seismic trace (T(t)). Further, post-stack attributes are classified into
many subparts based on complex trace and some of them are discussed in the
following sections.
(a) Signal envelope (E)
The Signal envelope (E) is also called Reflection Strength is calculated from the
complex trace as follows.
E(t) = SQRT T 2 (t) + H 2 (t) (7.2)
where SQRT represents square root. The signal envelope of the seismic signal has a
low-frequency appearance with only positive amplitudes. It often highlights the main
seismic features. The envelope represents the instantaneous energy of the signal and
is proportional in its magnitude to the reflection coefficient (Brown 1996).
The envelope is useful in highlighting discontinuities, changes in lithology, faults,
changes in deposition, tuning effect, and sequence boundaries. It also is proportional
to reflectivity and therefore useful for analyzing AVO anomalies. If there are two
volumes that differ by constant phase shift only, their envelopes will be the same
(Russell et al. 2001). This attribute represents mainly the acoustic impedance contrast,
hence reflectivity. This attribute is mainly useful in identifying bright spots, gas
accumulation, Sequence boundaries, major changes or depositional environments,
Thin-bed tuning effects, Unconformities, Major changes of lithology.
(b) Envelope derivative (RE)
Envelope derivative or time rate of envelope shift indicates the energy variability
of the reflected events. Events with a sharp relative increase also indicate a broader
bandwidth and therefore fewer impacts on absorption (Chopra and Marfurt 2007).
This attribute is also a physical characteristic and can be used to identify poten-
tial impacts of fracturing and absorption. The envelope derivative can be written
mathematically as follows.
dE(t)
RE(t) = (7.3)
dt
The derivative of the envelope highlights the
• Change in reflectivity and is also related to the absorption of energy.
• The sharpness of the rise time relates to absorption sharp interfaces
• Shows discontinuities,
7.2 Seismic Attributes 183
d 2 E(t)
DDE(t) = (7.4)
dt 2
The second derivative of the envelope highlights the interfaces very well. This
attribute is not too sensitive to the amplitude and can highlight even weak events. It is
used to show all reflecting interfaces visible within seismic bandwidth, the sharpness
of events, sharp changes of lithology (Chopra and Marfurt 2007). It is a very good
presentation of the image of the subsurface within the seismic bandwidth.
(d) Instantaneous phase
Due to the definition of wavefronts as lines of the continuous stage, the phase attribute
is also a physical characteristic and can be used efficiently as a discriminator for the
classification of geometric shapes. The instantaneous phase attribute is given as
follows.
H (t)
φ(t) = arctan (7.5)
T (t)
The seismic trace T (t) and its Hilbert transform H (t) are related to the envelope
E(t) and the phase φ(t) by the following relation.
The cosine of the instantaneous phase is also independent of amplitude and shows
bedding very well. This attribute is smoother than phase (which has discontinuities)
and thus is useful for automatic classification procedures.
(f) Instantaneous acceleration AC(t)
Instantaneous acceleration is the time derivative of the instantaneous frequency. It
highlights the change in instantaneous frequency which may be indicative of thin
bedding or absorption. Mathematically, the AC can be represented as follows.
d (F(t))
AC(t) = (7.7)
dt
The instantaneous acceleration used to accentuates bedding differences, higher
resolution, may have a somewhat higher noise level due to differentiation, may have
some relation to elastic properties of beds.
(g) Instantaneous frequency (IF)
Instantaneous frequency is the time derivative of the phase, i.e., the rate of change
of the phase. The instantaneous frequency attribute has been shown to be related to
the centroid of the seismic wavelet power spectrum (Cohen and Lee 1990; Barnes
1991, 1993). The instantaneous frequency attribute reacts to both impacts of wave
propagation and deposition features, so it is a physical attribute and can be used as
an efficient discriminator. Mathematically, it can be written as follows.
d (φ(t))
IF(t) = (7.8)
dt
Instantaneous frequency represents the mean amplitude of the wavelet. Instanta-
neous frequency can indicate bed thickness and also lithology parameters. It indi-
cates the edges of low impedance thin beds and also used for hydrocarbon indicator
and fracture zone indicator by low-frequency anomaly. This effect is sometimes
accentuated by the unconsolidated sands due to the oil content of the pores.
(h) Thin bed indicator
This attribute is designed to highlight the location where the instantaneous frequency
jumps or goes in the reverse direction (Tanner 1979). These jumps could be the result
of very close reflectors, i.e., thin beds. It is calculated as follows.
For this attribute, calculation window length is used to calculate the F (t). It
computed from large spikes of instantaneous frequency, indicate overlapped events,
indicates thin beds when laterally continuous indicates non-reflecting zone when
7.2 Seismic Attributes 185
it appears laterally at random, like salt and pepper and also used to fined detail of
bedding patterns.
(i) Instantaneous bandwidth B(t)
This attribute is described by Barnes 1991. It can be written as follows.
d (E(t))
B(t) = dt
(7.10)
2.E(t)
E(t)
Q(t) = ∗IF(t) ∗ ∂E(t)
(7.11)
∂t
quantities of data, however, which can be directly connected to fluid content and
fracture orientation. This class includes AVO, velocities and azimuthal variety of all
characteristics. The pre-stack attributes are further classified.
(a) RMS velocities of reflectors
This may be calculated from time migration velocity analysis, independent of the
major influence of dips. This is used for sand/shale ratios estimation, high-pressure
shale zone detection, major lithologic change detection, and etc.
(b) Discontinuity
It is again a geometrical attribute and measures the lateral relations in the data.
It is designed to emphasize the discontinuous events such as faults. High amplitude
values on this attribute correspond to discontinuities in the data, while low amplitude
values correspond to continuous features. Discontinuity varies between zero and one,
where zero is continuous and one is discontinuous. The derivation of the attribute is
prepared by considering several traces together to reveal the discontinuity (geometry)
of the beds (Soubotcheva and Stewart 2004). This attribute can be used to understand,
coherency at maximum coherency direction, minimum coherency direction, event
terminations, picked horizons, fault detection, zones of parallel bedding, zones of
chaotic bedding, non-reflecting zones, converging and diverging bedding patterns
and unconformities.
Some important attributes are estimated from the Blackfoot seismic data and pre-
sented in Figs. 7.2 and 7.3. The high amplitude anomaly is highlighted by the arrow.
One can also notice that some of the above-discussed attributes are dramatically
increased visibility of subsurface features and hence seismic amplitudes are very
important steps used in support of seismic data interpretation. Figures 7.4 and 7.5
shows slice at a 1060 ms time interval of all the attributes presented in Figs. 7.2 and
7.3. These slices represent the horizontal distribution of the attributes and also one
can notice the special features of the subsurface are increased.
This is the first type of geostatistical inversion method used by geoscientists to predict
subsurface rock properties away from the boreholes. In this methodology, the goal
is to discover an operator potentially linear that can predict well logs from nearby
information. Indeed, we choose not to evaluate the seismic information itself, but
the seismic attributes. The reason behind is that many of these attributes will be
nonlinear, thereby improving the technique’s predictive power as compared with the
raw seismic data. One other reason is that the breakdown of input information into a
component is often beneficial. This method is called preprocessing or extraction of
features and can often significantly enhance the efficiency of a pattern recognition
scheme by decreasing information dimensionality (Russell et al. 2001).
7.3 Single Attribute Analysis (SAA) 187
Fig. 7.2 Shows seismic attributes estimated from the Blackfoot seismic data sets, where a rep-
resents amplitude envelope, b amplitude weighted cosine phase, c amplitude weighted frequency,
d amplitude weighted phase, e apparent polarity and f represents cosine instantaneous phase
The single attribute analyzes all or some of the inner and external attributes to
determine the most helpful attribute for anticipating target logs. The single attribute
expression implies the attributes are used separately as opposed to the Multi-attribute
where the characteristics are used in groups.
Single attribute analysis processes in which one attribute at a time used to find
out the relation between attribute and desired well log property. In the first process
toward prediction of petrophysical parameters, seismic attributes are extracted from
the seismic data. In the second step, the best attribute is selected. Thereafter, that
particular attribute is cross plotted with the target log for deriving a linear relationship
that is utilized to predict log property.
Table 7.1 contains a list of attributes and their correlation with desired well log
property. Each line of the table contains information for a particular attribute. Usu-
ally, it is best to use all the attributes and rank them according to their accuracy
to predict the target logs. The attributes are ordered by increasing prediction error,
i.e., by decreasing suitability for predicting the target. Nonlinear transformations
will check to see if the inverse, log, square root or other easy transformations will
generate a linear connection between good log and seismic attributes. For example,
188 7 Geostatistical Inversion
Fig. 7.3 Shows seismic attributes estimated from the Blackfoot seismic data sets, where a represents
derivative, b derivative instantaneous amplitude, c instantaneous frequency, d instantaneous phase,
e quadrature trace and f second derivative instantaneous amplitude
a relationship between a log attribute and a seismic attribute may not be linear (it is
cross plot may resemble a curve), but if one or both of the attributes are transformed
into log values, a linear relationship may appear, with the points showing a straight
line on the cross plot.
From Table 7.1, one can notice that attribute 1/impedance is the best one and
hence can be used for the prediction process. Thereafter, this selected attribute is
cross-plotted with well log porosity (density porosity) to predict porosity in the
inter-well region and represented in Fig. 7.6. In this analysis, the assumption is that
the target log has been integrated to travel time at the same sample rate as the seismic
attribute. Effectively, this integration reduces the target log to the same resolution as
the attribute, which is usually significantly coarser than the log property. Each point
in the cross plot consists of a pair of numbers corresponding to a particular time
sample.
Assuming a linear relationship between the target log (density porosity) and the
attribute, a straight line may be fit by regression in the following way.
7.3 Single Attribute Analysis (SAA) 189
Fig. 7.4 Shows slices at 1065 ms time from seismic attributes estimated from the Blackfoot seismic
data sets, where a represents amplitude envelope, b amplitude weighted cosine phase, c ampli-
tude weighted frequency, d amplitude weighted phase, e apparent polarity and f represents cosine
instantaneous phase
y = qx + p (7.12)
The coefficients p and q in the above equation may be derived by minimizing the
mean-squared prediction error as following.
1
N
E2 = (yi − p − qxi )2 (7.13)
N i=1
where the sum is overall points in the cross plot. The calculated prediction error E
is a measure of the goodness-of-fit for the regression line defined by Eq. (7.12). An
190 7 Geostatistical Inversion
Fig. 7.5 Shows slices at 1065 ms time of seismic attributes estimated from the Blackfoot seismic
data sets, where a represents derivative, b derivative instantaneous amplitude, c instantaneous
frequency, d instantaneous phase, e quadrature trace and f second derivative instantaneous amplitude
where,
1
N
σxy = (xi − mx ) yi − my (7.15)
N i=1
7.3 Single Attribute Analysis (SAA) 191
Table 7.1 Single attribute list and their correlation with target porosity
S. No. Target Attribute Error Correlation
1. Porosity 1/(Impedance) 5.01322 0.381154
2. Porosity Impedance 5.041142 −0.36842
3. Porosity (Impedance) * 2 5.066357 −0.35646
4. Porosity Derivative 5.261693 0.24177
5. Porosity Quadrature trace 5.275328 −0.23145
6. Porosity Amplitude weighted cosine phase 5.275942 0.230968
7. Porosity Filter 5.306574 0.205723
8. Sqrt (porosity) Derivative 5.308363 0.249109
9. Porosity Amplitude weighted phase 5.320532 −0.19307
10. Porosity Cosine instantaneous phase 5.343901 0.169711
11. Porosity Instantaneous phase 5.36141 −0.14976
12. Porosity Derivative instantaneous amplitude 5.400413 0.090291
13. Porosity Second derivative 5.403192 −0.08445
14. Porosity Integrated absolute amplitude 5.403278 −0.08426
15. Porosity Instantaneous frequency 5.411403 0.064118
Fig. 7.6 Crossplot of density porosity with attributes. A best fit straight line is fitted to the data
(redline)
1
N
σx = (xi − mx ) (7.16)
N i=1
1
N
σy = xi − my (7.17)
N i=1
192 7 Geostatistical Inversion
1
N
mx = (xi ) (7.18)
N i=1
1
N
my = (yi ) (7.19)
N i=1
Fig. 7.7 The plot shows the target logs in black with the predicted logs in red. The average error
and well name reported at the top of the plot
7.3 Single Attribute Analysis (SAA) 193
Fig. 7.9 Cross-section (inline 244) of predicted porosity estimated using a single attribute analysis
algorithm
one can not differentiate the thin subsurface layer. This is because the prediction
depends only on one attributes hence the resolution is poor.
The use of single attribute analysis over 3D datasets is very limited. In most
cases, the single attribute analysis shows poor correlation of predicted curve with
the actual curve and hence cannot use for quantitative interpretation. Although, the
methods can be utilized for qualitative interpretation as the method is very fast and
hence within a limited time, one can generate a volume of petrophysical properties.
To overcome the problem encounter in single attribute analysis, the multi-attribute
regression needs to be analyzed.
194 7 Geostatistical Inversion
Fig. 7.10 3D view of predicted porosity estimated using single attribute analysis technique
The multi-attribute regression utilized more than one attribute at a time and predict
petrophysical parameters away from the boreholes. This can be achieved by the
extension of the conventional linear analysis (single attribute analysis) to multiple
attributes (multi-attribute regression). Assume, for simplicity, that we have three
attributes (a1 , a2 and a3 ) good enough to predict porosity. Figure 7.11 shows the
target log along with three attributes and represents how these attributes used to
predict log value.
At each time sample, the target log is modeled by the linear equation as follows.
where w0 , w1 , w2 , and w4 are weights. The weights in Eq. 7.20 may be derived
by minimizing the mean-squared prediction (Hampson et al. 2000, 2001) error, as
follows.
1
N
E2 = (Li − w0 − w1 a1i − w2 a2i − w3 a3i )2 (7.21)
N i=1
7.4 Multi-attribute Regressions (MAR) 195
Fig. 7.11 Assuming the case of three seismic attributes, each target log sample is modeled as a
linear combination of attribute samples at the same time
The solution for the four weights produces the standard normal equations which
can be estimated in the following way.
Multi-attribute linear regression is an extension of simple attribute regression
methods with M number of variables. Now we will use a M number of seismic
attributes, a1 , a2 , . . . aM , to predict the log properties (L). To do this, one must deter-
mine the M + 1 weights w1 , w2 , . . . wM , which, when multiplied by the particular
set of attribute values, give the closest result to the log in a least-squared sense. For
simplicity, assume that M = 3 (Hampson et al. 2001). If one has N samples in log
data, then, the following equations can be written.
where aij is the jth sample of the ith attribute. Notice that the above equations can be
written in matrix form as follows.
⎡ ⎤ ⎡ 1 a a21 a31
⎤⎡ ⎤
L1 11 w0
⎢ L2 ⎥ ⎢ 1 a a22 a32 ⎥⎢
⎥⎢ w1 ⎥
⎢ ⎥ ⎢ . 12
⎥
⎣ L3 ⎦ = ⎢
⎣ ..
.. .. .. ⎥⎣
⎦ w2 ⎦
(7.22)
. . .
L4 1 a1N a2N a33 w3
L = AW (7.23)
196 7 Geostatistical Inversion
Just as in the single-attribute case, the mean-squared error (7.23) calculated using
the derived weights constitutes a goodness-of-fit measure for the transform, as does
the normalized correlation defined in Eq. (7.14), where the x-coordinate is now the
predicted log value and the y-coordinate is the real log value (Russell et al. 1997;
Swisi 2009).
We derived equations in the previous sections that enable us to identify optimal oper-
ators for any set of attributes. These operators are optimal in the sense of minimizing
the mean-squared error of prediction between the actual target logs and the predicted
target logs. The next issue to be dealt with is how these attributes are selected.
There are two approaches available to select attributes. The first approach is to
exhaustive search and the second approach is stepwise regression methods. The first
approach is very much time-consuming and hence less used. The second approach
i.e. stepwise regression although less optimal, a procedure developed by Draper and
Smith in (1966). In this procedure, the assumption is that if you already know the
best combination of M attributes, then the best combination of M + 1 attributes
includes the previous M attributes as members. Of course, it is necessary to rederive
the previously calculated coefficients. The process is shown in the section below.
First, by exhaustive search, find the single best attribute and solve the optimum
coefficients and calculate the error of prediction. The best attribute is the one with
the lowest prediction error. Let say this attribute to be a1 .
Second, find the best attribute pair, assuming that attribute 1 (a1 ) is the first mem-
ber. Form all pairs for each other attribute in the list, solve the optimum coefficients
and calculate the error of prediction. The best pair is the one with the lowest error of
prediction. Let say the second attribute of the best pair is a2 .
Similarly, in the third step, find the best triplet of attributes, assuming that the
first two members are a1 and a2 . For each other attribute in the list, all triplets are
7.4 Multi-attribute Regressions (MAR) 197
formed and solve the optimum coefficients for each triplet and calculate the error of
prediction. The best triplet is a3 .
Continue as long as you desired this process. The first thing to note is that this
process’s calculation time is much shorter than exhaustive searching (Russell et al.
1997, 2001). The problem with stepwise regression is that one cannot be sure that the
optimal solution will be derived. In other words, the combination of five attributes
found may not be the best five by exhaustive search. Each additional attribute found,
however, has an error of prediction less than or equal to the combination found
previously. If the new prediction error is greater than the combination of the previous
set, all the weights for this new attribute to be zero and the prediction error is equal
to the previous set.
Table 7.2 Multi-attribute list and their error with target porosity
S. No. Target Final attribute Training error Validation error
1. Porosity 1/(Impedance) 4.963269 5.046682
2. Porosity Amplitude weighted cosine phase 4.813295 4.959127
3. Porosity X-Coordinate 4.77781 4.925902
4. Porosity Apparent polarity 4.745512 4.95469
5. Porosity Instantaneous frequency 4.710768 4.973616
6. Porosity Dominant frequency 4.684892 4.982152
7. Porosity Filter 25/30–35/40 4.659465 4.991252
198 7 Geostatistical Inversion
The desired attributes are then cross plotted against the well log porosity and a best
fit straight line provides the relation between attributes and porosity at a particular
sample point which is further utilized to predict porosity in the inter-well region. The
process is similar to the single attribute analysis but the only difference is that in multi-
attribute regression uses more than one attribute at a time. The MAR is first applied
to the composite trace near to well locations and porosity is predicted. Although the
porosity is already known at a well location to cross-verify predicted porosity, the
composite trace analysis needs to be performed. Figure 7.13 shows a comparison of
predicted porosity (red lines) with actual porosity (black lines) for well 01-17, 08-08
and 16-08. From the figure, one can notice that the predicted porosity is followed the
trend to the original porosity for all wells. The average correlation between the two is
estimated to be 0.58 and error 4.61. Further, a cross plot is generated between actual
porosity from well log data and predicted porosity and displayed in Fig. 7.14. The
distribution of scatter points shows that the predicted porosity closes to the original
one. The composite trace analysis shows quite good prediction results and hence one
can go for 3D porosity prediction.
Now that we have derived the relationship between attributes and target logs, one
can apply the analysis to the whole 3D volume. A cross-section of predicted 3D
porosity volume at inline 27 is presented in Fig. 7.15. The predicted section shows
very high-resolution subsurface images as compared to the input seismic section. The
well log porosity from well 08-08 is also plotted on the inverted section and shows
a good match between them. The close inspection of the inverted section shows
an anomaly zone with very high porosity ranging from 0.12 to 0.15 particularly in
the time interval between 1055 ms and 1070 ms two way travel time. This high
porosity is correlated with low impedance anomaly interpreted in chapter 4 and
hence interpreted as Glauconitic sand channel. A 3D volume of predicted porosity
is displayed in Fig. 7.16. The figure highlighted the anomalous zone in all directions
in the subsurface.
The benefit of the 3D section is that one can analyze the anomalous zone vertically
as well as horizontally. Till now, the discussion is performed for vertical variation of
the sand channel but now, one can generate a horizontal variation of the channel to
7.4 Multi-attribute Regressions (MAR) 199
Fig. 7.13 Comparison of predicted porosity with the actual porosity from the a well 01-17, b well
08-08, and c well 16-08
Fig. 7.14 Crossplot of predicted porosity versus actual porosity for all the 13 wells. The red solid
line indicates the best fit line
200 7 Geostatistical Inversion
Fig. 7.15 Cross-section of predicted porosity at inline 61. The anomalous zone (high porosity) is
highlighted by the ellipse
Fig. 7.17 Slice at 1065 ms two way travel time of predicted porosity
monitor subsurface features. A slice at 1065 ms time intervals is shown in Fig. 7.17
which shows the variation of porosity in the subsurface. One can notice that the
horizontal variation of porosity is nicely visible and can be used to analyze the details
of the sand channel. This slice indicates the variation of the sand channel by blue
to pink hot color. From this analysis, the width of the anomalous zone horizontally
can also be estimated. It is also noticed that the sand channel varies from north-east
(NE) to south-west (SW) direction. As the seismic data shows amplitude variation
along with time and cannot be interpreted subsurface lithology whereas by using
the above techniques one can obtain the distribution of petrophysical parameters in
the inter-well region. In this example, porosity prediction is performed but the other
petrophysical properties like permeability, saturation, shale sand ratio, gamma-ray,
etc. can also be generated which helps to decide movable hydrocarbon saturation, net
pay, etc. This technique is very helpful for the exploration and production projects.
The assessment was linear so far. A linear straight line was fitted to the cross plot
between the target logs and attribute(s) and weights are calculated for the regression
line by minimizing the mean-square error. But, you might think a greater order curve
would fit better with the points rather than fitting a straight line. There are several
choices to calculate this curve. One choice is to apply a nonlinear transformation to
either or both of the variables and to fit a straight line to the transformed data (Schultz
et al. 1994). A second choice is to match a polynomial of greater order. A third choice
202 7 Geostatistical Inversion
is also available which is to use a neural network to obtain the relationship between
target logs and attributes. In geophysics, neural networks are becoming increasingly
common (McCormack 1991; Schuelke et al. 1997) nowadays.
Neural networks have acquired popularity over the past century in geophysics.
They were effectively implemented to a verity of issues. Neural networks have been
used in the geophysical domain for waveform recognition and first-break pick-up
(Murat and Rudman 1992; McCormack et al. 1993); for electromagnetic (Poulton
et al. 1992), magnetotelluric (Zhang and Paulson 1997) and seismic inversion (Röth
and Tarantola 1994; Langer et al. 1996; Macías et al. 1998); shear-wave splitting (Dai
and MacBeth 1994), well-log analysis (Huang et al. 1996), trace editing (McCormack
et al. 1993), seismic deconvolution (Wang and Mendel 1992; Calderón-Macías et al.
1997), and event classification (Dowla et al. 1990; Romeo 1994); and for many other
problems.
Many types of neural networks exist but two popular types’ namely multi-layer
feed-forward neural network (MLFN) and Probabilistic neural network (PNN). Both
methods have their advantages and disadvantages to the application of data. Both
methods are described briefly in the following sections along with real data examples.
In this section, we limit our discussion to the static feed-forward neural networks
(MLFN). It is called a static algorithm because the weights estimated during the
process are fixed and do not change with time. The feed-forward indicates that the
output is not feedback, i.e., refed, to the network. Thus, this type of network does not
iterate to a final solution but directly translates the input signals to output independent
of previous input.
McCulloch and Pitts (1943) conceived the mathematical perceptron some 55 years
ago to mimic the behavior of the biological neuron (Fig. 7.18). The biological neuron
consists primarily of three components/layers, one input layer, one or more hidden
layers, and one output layer. Each layer of the network consists of nodes and connects
these nodes with each other. These connections represent weights. These weights are
the final element of the production decision as they play a very significant part in
the implementation of MLFN. Figure 7.18 shows the MLFN’s network architecture.
The amount of input nodes is determined by the amount of characteristics during the
MLFN phase, and one chooses to use them to forecast rock characteristics (Maurya
and Singh 2019). The number of attributes is determined by the convolution operator
process in which a specific combination of attributes has been tested and the combi-
nation with maximum similarity is selected against the desired log property (Russell
et al. 2001; Svozil et al. 1997).
MLF neural networks are the most common neural networks, trained with a learn-
ing algorithm for back-propagation. For the formal description of the neurons, we
can use the so-called mapping function Γ , that assigns for each neuron i a subset
Γ (i) ⊆ V which consists of all ancestors of the given neuron. A subset Γ −1 (i) ⊆ V
7.5 Neural Network Techniques 203
Fig. 7.18 Representation of MLFN architecture. Ai indicates desired attributes, wij indicated
weights
than consists of all predecessors of the given neuron i. Each neuron in a particular
layer is connected with all neurons in the next layer. The connection between the ith
and jth neuron is characterized by the weight coefficient ωij and the ith neuron by
the threshold coefficient vi . The weight coefficient reflects the degree of importance
of the given connection in the neural network (Svozil et al. 1997). The output value
(activity) of the ith neuron xi is determined by Equations given below.
xi = f (ψi ) (7.26)
ψi = vi + ωij xj (7.27)
j∈ i−1
where ψi is the potential of the ith neuron and function f (ψi ) is the so-called transfer
function. The threshold coefficient can be understood as a weight coefficient of the
connection with formally added neuron j, where xj = 1.
The transfer function can be written as follows.
1
f (ψi ) = (7.28)
1 + exp(−ψ)
The supervised adaptation process varies the threshold coefficients vi and weight
coefficients ωij to minimize the sum of the squared differences between the computed
and required output values (Luo et al. 2000). This is accomplished by minimization
of the objective function:
204 7 Geostatistical Inversion
1 2
E= xp − x̂p (7.29)
p
2
where xp and x̂p are vectors composed of the computed and required activities of the
output neurons and summation runs over all output neurons p.
The neural network works in two modes, first is training and the second is the
prediction. You need two information sets, the practice set and the set you want to
predict (test set) to train the neural network and to predict.
The training mode starts with arbitrary weight values that could be random num-
bers and iteratively proceeds. An epoch is called every iteration of the full training
set. The network adjusts the weights in the direction that decreases the error in each
epoch. The weights gradually converge to the locally optimal set of values as the iter-
ative method of incremental adjustment proceeds. Before training is finished, many
epochs are generally needed (Svozil et al. 1997; Wu et al. 1992).
Back-propagation learning can be carried out in one of two fundamental ways
for a specified training set i.e. pattern mode and batch mode. Weight updating is
conducted after each training pattern is presented in the pattern mode of teaching
backpropagation. Weight updating is conducted in the back-propagation learning
batch mode after all the training examples are presented (i.e. after the whole epoch).
The pattern mode is preferred over the batch mode from an ‘on-line’ point of per-
spective, as it needs less local storage for each synaptic connection. In addition, given
that patterns are randomly presented to the network, the use of patterns by pattern
updating weights makes the search in weight space stochastic, making it less likely
that the backpropagation algorithm will be trapped in a local minimum. On the other
side, the use of training batch mode offers a more precise gradient vector estimate.
For example, pattern mode should be used in online process control, as not all train-
ing patterns are available in the given time. Ultimately, the comparative efficacy of
the two training methods relies on the issue solved (Russell et al. 1997).
Information flows through the network in forecast mode, from inputs to out-
puts. The network processes one example at a time, producing an estimate based
on the input values of the output value(s). The resulting error is used as an esti-
mate of the trained network’s predictive quality. One can begin with a practice set
in back-propagation learning and use the back-propagation algorithm to calculate
the network’s synaptic weights (Tonn 2002; Iturrarán and Parra 2014). The hope
is that the so-called neural network will become widespread. A network is said to
generalize well when the network-computed input-output relationship is right (or
almost right) for input/output models that have never been used in-network practice.
Generalization is not a mystical property of neural networks but can be contrasted
with the impact of excellent non-linear input information interpolation. When the
7.5 Neural Network Techniques 205
learning method repeats too many iterations (i.e. the neural network is over-trained or
over-fitted, there is no distinction between overtraining and overfitting), the network
can memorize the training data and thus be less willing to generalize comparable
patterns of input-output (Leite and Vidal 2011; Mahmood et al. 2017). The network
provides almost ideal outcomes from the practice set for examples but fails from the
test set for example. Overfitting can be contrasted with unsuitable polynomial degree
selection in polynomial regression. With noisy data, severe overfitting can happen
even when there are many more instances of training than weights. A sufficiently
big set of training instances is the fundamental condition for excellent generaliza-
tion. This training set must be representative of the set of all cases to which you
wish to generalize at the same time. The significance of this situation is due to two
kinds of generalization: interpolation and extrapolation (Masters 1995; Masri et al.
2000; Singh et al. 2016). Interpolation applies to instances encircled more or less
by neighboring instances of practice; all else is extrapolation. In specific, instances
outside the training information spectrum involve extrapolation. Interpolation can
often be performed reliably, but it is notoriously unreliable to extrapolate. There-
fore, adequate training information is essential to prevent the need for extrapolation.
Experimental design results in methods for choosing excellent training sets. Because
of a set quantity of training data, there are some efficient methods to avoid overfitting
and generalization (Wu et al. 1992).
An example from the Blackfoot field, Canada is presented here to increase the under-
standing of MLFN methods. As in the earlier case, the MLFN also applied to the
data, first is for composite trace and the second is for entire 3D datasets. For this
purpose, one composite trace extracted from near to well 08-08 and thereafter the
MLFN is applied on it to predict porosity trace. Figure 7.19 shows a comparison
Fig. 7.19 Comparison of well log porosity with actual porosity estimated using MLFN methods
206 7 Geostatistical Inversion
of predicted porosity and actual porosity. It is noticed from the Fig. 7.19 that the
predicted porosity is following the trend with actual porosity quite nicely. The aver-
age correlation coefficient is estimated to be 0.59 and RMS error is to be 4.7. These
values show that the MLFN technique is quite good enough to estimate subsurface
rock property although the correlation is low and error is high as compared with the
MLFN methods. This is not always true because the MLFN largely depends on input
parameters. If one can able to choose the perfect combination of input parameters,
then the methods can show better results as compared with multi-attribute regres-
sion. In this example, the analysis is presented to understand how MLFN works,
not to compare with other geostatistical methods discussed here. The analyses were
performed in different conditions and considering different input parameters hence
cannot compare explicitly with each other.
A cross plot between predicted porosity and actual porosity is also generated and
shown in Fig. 7.20. The distribution of scatter points demonstrates that the predicted
porosity is close to the actual porosity.
To analysis the 3D variation of porosity, the MLFN can be performed for the entire
seismic volume. It is time taking but can provide a detailed picture of the subsurface.
A cross-section at inline 27 is presented in Fig. 7.21. The inverted section shows
very high-resolution subsurface layers. Well, log porosity from well 08-08 is also
plotted on the predicted section. One can notice that the predicted porosity and well
log porosity matching very nicely with each other. The layer of the subsurface can be
easily identified with this section. The interpretation of the predicted porosity section
shows many patches of high porosity ranging from 0.12 to 0.15. One patch of high
porosity is matching with a well log anomaly zone which falls in between 1055 and
1070 ms time interval and predicted as a sand channel. This high porosity anomaly
Fig. 7.21 Cross-section of predicted porosity at inline 27 estimated using MLFN methods. The
anomalous zone is highlighted by the ellipse
is corroborating well with the low impedance anomaly and high porosity anomaly
estimated with other methods discussed so far. A 3D plot of predicted porosity is
provided in Fig. 7.22 to get the variation of the sand channel vertically as well as
horizontally. Further, to measure the horizontal variation of the sand channel, slice at
1065 ms time intervals are generated from predicted porosity volume and shown in
Fig. 7.23. The variation of the sand channel is highlighted by the arrow. It is noticed
from the figure that the sand channel varies from NE to SW direction which is already
shown for multi-attribute analysis cases.
Fig. 7.23 Slice at channel interval (1065 ms) highlighted channel trend
7.5 Neural Network Techniques 209
In Eq. 7.30, n represents the number of examples and D(x, xi ) is the distance
between the input point and each of the training points xi . The quantity D(x, xi ) can
be calculated as follows.
3
xj − xij 2
D(x, xi ) = (7.31)
j=1
σj
N
2
EV (σ1 , σ2 , σ3 ) = Li − Li (7.32)
i=1
This can be also noted that the prediction error depends largely on the choice
of the parameters σj . The validation error is minimized using a nonlinear conjugate
gradient algorithm with respect to the smoothing parameters.
Fig. 7.24 Comparison of predicted porosity with actual porosity estimated using PNN methods
is noticed that the predicted curves are following original curves very nicely for
all wells. Further, the cross-plot between predicted porosity and original porosity is
generated and shown in Fig. 7.25. The distribution of scattered points shows that
the inverted porosity is very close to the actual porosity which is also evident by
the correlation coefficient. The PNN shows superior results, i.e. lower prediction
error and lower validation error. The correlation coefficient is estimated to be 0.76
and RMS error is 3.93, indicating a good correlation and RMS. One can argue that
the probabilistic neural network predicts the logs with higher accuracy. The multi-
regression analysis predicted the logs with correlation 0.59 while the PNN predicted
them with correlation 0.76. A comparison of prediction error with validation error
for all wells are presented in Fig. 7.26. Both errors vary in a similar way with an
average value of 3.93.
Once the relation between the seismic attributes and the porosity records is deter-
mined, the PNN is applied to the data volumes. Figures 7.27 demonstrate a porosity
cube at inline 27. The well log porosity is also plotted over the predicted porosity
and it is found that both the porosities are matching very well. It is very possible to
distinguish the sand channel as a large anomaly of porosity. Note again the greater
resolution obtained through the use of the probabilistic neural network. The 3D view
of predicted porosity is shown in Fig. 7.28 and three-dimensional variation of the
sand channel are can be noticed with high porosity contrast. Figure 7.29 shows poros-
ity slice at the sand channel level and shows the horizontal distribution of the sand
7.5 Neural Network Techniques 211
Fig. 7.25 Cross plot of actual porosity and predicted porosity for well 08-07
Fig. 7.26 Variation of Prediction and validation error with well logs
channel. The high porosity zone in the seismic section is shown in Fig. 7.30. This
figure nicely highlighted the anomalous zone. From this analysis, one can easily
interpret the seismic section and can identify the visualization of the sand channel.
Apart from that, these sections can also be used to derive a relationship between
them. For example, the inverted impedance and predicted porosity section can be
cross plotted and can be used to estimate the relation between them by fitting a
straight line or higher-order polynomial. This relationship will be more meaningful
as compared with the relation derived from the well log curves at well location. The
well log derived equation provides a relationship at a particular location that cannot
212 7 Geostatistical Inversion
Fig. 7.27 Cross-section of predicted porosity at inline 27. The anomalous zone is highlighted by
an ellipse
Fig. 7.28 3D view of predicted porosity. The anomalous zone is shown in pink color
References 213
Fig. 7.30 Seismic cross-section highlighted the sand channel. The color bar indicates porosity
valid away from the well whereas the relationship derived from the cross-section is
valid anywhere in the region and give a better estimate (Maurya et al. 2017) of the
subsurface property.
References
Adeli H, Panakkat A (2009) A probabilistic neural network for earthquake magnitude prediction.
Neural Netw 22(7):1018–1024
Anderson JK (1996) Limitations of seismic inversion for porosity and pore fluid: lessons from
chalk reservoir characterization and exploration. In: SEG technical program expanded abstracts,
society of exploration geophysicists, pp 309–312
Barnes AE (1991) Instantaneous frequency and amplitude at the envelope peak of a constant-phase
wavelet. Geophysics 56(7):1058
214 7 Geostatistical Inversion
Barnes AE (1993) Instantaneous spectral bandwidth and dominant frequency with applications to
seismic reflection data. Geophysics 58(3):419–428
Barnes AE (1994) Theory of two-dimensional complex seismic trace analysis. In: SEG technical
program expanded abstracts, society of exploration geophysicists, pp 1580–1583
Bosch M, Mukerji T, Gonzalez EF (2010) Seismic inversion for reservoir properties combining
statistical rock physics and geostatistics: a review. Geophysics 75(5):75A165–75A176
Brown AR (1996) Seismic attributes and their classification. Lead Edge 15(10):1090
Calderón-Macías C, Sen MK, Stoffa PL (1997) Hopfield neural networks, and mean-field annealing
for seismic deconvolution and multiple attenuations. Geophysics 62(3):992–1002
Chambers RL, Yarus JM (2002) Quantitative use of seismic attributes for reservoir characterization.
CSEG Recorder 27(6):14–25
Chiles J (1988) Fractal and geostatistical methods for modeling of a fracture network. Math Geol
20(6):631–654
Chopra S, Marfurt KJ (2005) Seismic attributes—a historical perspective. Geophysics 70(5):3SO–
28SO
Chopra S, Marfurt KJ (2007) Seismic attributes for prospect identification and reservoir charac-
terization. Society of Exploration Geophysicists and European Association of Geoscientists and
Engineers
Cohen L, Lee C (1990) Instantaneous bandwidth for signals and spectrograms. In: International
conference on acoustics, speech, and signal processing, IEEE, pp 2451–2454
Dai H, MacBeth C (1994) Split shear-wave analysis using an artificial neural network. The first
Break 12(12):605–613
Dowla FU, Taylor SR, Anderson RW (1990) Seismic discrimination with artificial neural networks:
preliminary results with regional spectral data. Bull Seismol Soc Am 80(5):1346–1373
Doyen PM (1988) Porosity from seismic data: a geostatistical approach. Geophysics 53(10):1263–
1275
Draper N, Smith H (1966) Applied regression analysis. Wiley, New York
Dubrule O (2003) Geostatistics for seismic data integration in Earth models. Distinguished instructor
short course. Number 6. SEG Books
Eskandari H, Rezaee MR, Mohammadnia M (2004) Application of multiple regression and artificial
neural network techniques to predict shear wave velocity from well log data for a carbonate
reservoir, south-west Iran. CSEG Recorder 29(7):42–48
Haas A, Dubrule O (1994) Geostatistical inversion- a sequential method of stochastic reservoir
modeling constrained by seismic data. First Break 12(11):561–569
Hampson D, Todorov T, Russell B (2000) Using multi-attribute transforms to predict log properties
from seismic data. Explor Geophys 31(3):481–487
Hampson DP, Schuelke JS, Quirein JA (2001) Use of multiattribute transforms to predict log
properties from seismic data. Geophysics 66(1):220–236
Huang Z, Shimeld J, Williamson M, Katsube J (1996) Permeability prediction with artificial neural
network modeling in the Venture gas field, offshore eastern Canada. Geophysics 61(2):422–436
Iturrarán VU, Parra JO (2014) Artificial neural networks applied to estimate permeability, porosity
and intrinsic attenuation using seismic attributes and well-log data. J Appl Geophys 107:45–54
Jones TA, Helwick SJ (1998) Method of generating 3-d geologic models incorporating geologic
and geophysical constraints. US Patent 5:634–838
Langer H, Nunnari G, Occhipinti L (1996) Estimation of seismic waveform governing parameters
with neural networks. J Geophys Res Solid Earth 101(B9):20109–20118
Leiphart DJ, Hart BS (2001) Comparison of linear regression and a probabilistic neural network
to predict porosity from 3-D seismic attributes in Lower Brushy Canyon channeled sandstones,
southeast New Mexico. Geophysics 66(5):1349–1358
Leite EP, Vidal AC (2011) 3D porosity prediction from seismic inversion and neural networks.
Comput Geosci 37(8):1174–1180
Lindseth RO (1979) Synthetic sonic logs a process for stratigraphic interpretation. Geophysics
44(1):3–26
References 215
Luo X, Patton AD, Singh C (2000) Real power transfer capability calculations using multi-layer
feed-forward neural networks. IEEE Trans Power Syst 15(2):903–908
Macías D, Pérez-Pomares JM, García-Garrido L, Carmona R, Muñoz-Chápuli R (1998) Immunore-
activity of the ETS-1 transcription factor correlates with areas of epithelial-mesenchymal
transition in the developing avian heart. Anat Embryol 198(4):307–315
Mahmood MF, Shakir U, Abuzar MK, Khan MA, Khattak N, Hussain HS, Tahir AR (2017) A
probabilistic neural network approach for porosity prediction in the Balkassar area: a case study.
Journal of Himalayan Earth Science 50:1
Masri S, Smyth A, Chassiakos A, Caughey T, Hunter N (2000) Application of neural networks for
the detection of changes in nonlinear systems. J Eng Mech 126(7):666–676
Masters T (1995) Advanced algorithms for neural networks: a C++ sourcebook. Wiley, Hoboken
Maurya SP, Singh NP (2018) Comparing pre-and post-stack seismic inversion methods-a case study
from Scotian Shelf, Canada. J Ind Geophys Union 22(6):585–597
Maurya SP, Singh KH (2019) Predicting porosity by multivariate regression and probabilistic
neural network using model-based and colored inversion as external attributes: a quantitative
comparison. J Geol Soc India 93(2):207–212
Maurya SP, Singh KH, Singh NP (2019) Qualitative and quantitative comparison of geostatistical
techniques of porosity prediction from the seismic and logging data: a case study from the
Blackfoot Field, Alberta, Canada. Mar Geophys Res 40(1):51–71
Maurya SP, Singh KH, Kumar A, Singh NP (2017) Reservoir characterization using post-stack
seismic inversion techniques based on a real coded genetic algorithm. J Geophys 39(2):95–103
McCormack MD (1991) Neural computing in geophysics. Lead Edge 10(1):11–15
McCormack MD, Zaucha DE, Dushek DW (1993) First-break refraction event picking and seismic
data-trace editing using neural networks. Geophysics 58(1):67–78
McCulloch WS, Pitts W (1943) A logical calculus of the ideas immanent in nervous activity. Bull
Math Biophys 5(4):115–133
Murat ME, Rudman AJ (1992) Automated first arrival picking: a neural network approach 1.
Geophys Prospect 40(6):587–604
Poulton MM, Sternberg BK, Glass CE (1992) Location of subsurface targets in geophysical data
using neural networks. Geophysics 57(12):1534–1544
Pramanik AG et al (2004) Estimation of effective porosity using geostatistics and multi-attribute
transforms a case study. SEG 69(2):352–372
Romeo G (1994) Seismic signals detection and classification using artificial neural networks. Annals
of Geophysics 37:3
Röth G, Tarantola A (1994) Neural networks and inversion of seismic data. J Geophys Res Solid
Earth 99(B4):6753–6768
Russell B, Hampson D, Schuelke J, Quirein J (1997) Multiattribute seismic analysis. Lead Edge
16(10):1439–1444
Russell SA, Reasoner C, Lay T, Revenaugh J (2001) Coexisting shear- and compressional-
wave seismic velocity discontinuities beneath the central Pacific. Geophysical Research Letters
28(11):2281–2284
Schuelke JS, Quirein JA, Sag JF, Altany DA, Hunt PE (1997) Reservoir architecture and porosity
distribution, Pegasus field, West Texas—an integrated sequence stratigraphic-seismic attribute
study using neural networks. In SEG technical program expanded abstracts, society of exploration
geophysicists, pp 668–671
Schultz PS, Ronen S, Hattori M, Corbett C (1994) Seismic-guided estimation of log properties (Part
1: a data-driven interpretation methodology). Lead Edge 13(5):305–310
Singh S, Kanli AI, Sevgen S (2016) A general approach for porosity estimation using artificial
neural network method: a case study from Kansas gas field. Stud Geophys Geod 60(1):130–140
Soubotcheva N, Stewart RR (2004) Predicting porosity logs from seismic attributes using
geostatistics. CREWES Research Report -Volume 1:1–4
Specht DF (1990) Probabilistic neural networks. Neural Netw 3(1):109–118
Specht DF (1991) A general regression neural network. IEEE Trans Neural Netw 2(6):568–576
216 7 Geostatistical Inversion