#### Improving seismic interpretation: a high-contrast approximation to the reflection coefficient of a plane longitudinal wave

Pet.Sci.
Improving seismic interpretation: a high-contrast
Yin Xingyao 0
Zong Zhaoyun 0
Wu Guochen 0
0 China University of Petroleum (Beijing) and Springer-Verlag Berlin Heidelberg 2013
amplitude versus offset (AVO) analysis and inversion in exploration geophysics. However, the weak properties contrast hypothesis of those linearized approximate equations leads to big errors when the two media across the interface vary dramatically. To extend the application of AVO analysis and inversion to high contrast between the properties of the two layers, we derive a novel nonlinearized high-contrast (A PP wave is a reflected compressional wave from an incident compressional wave (P-wave).) This novel approximation is derived from the exact reflection coefficient equation with Taylor expansion for the incident angle. Model tests demonstrate that, compared with the reflection coefficients of the linearized approximations, the reflection coefficients of the novel nonlinearized approximate equation agree with those of the exact PP equation better for a high contrast interface with a moderate incident angle. Furthermore, we introduce a nonlinear direct inversion method utilizing the novel reflection P-wave velocities, S-wave velocities, and densities in the upper and lower layers across the interface. This nonlinear inversion algorithm is able to estimate the inverse of the nonlinear function in terms of model and suitability of this novel approximation for a high contrast interface, and we still could estimate the six parameters across the interface reasonably when the parameters in both media across the interface vary
High-contrast interface; reflection coefficient; amplitude variation with angle; multi-
1 Introduction
The Zoeppritz equation (Zoeppritz and Erdbebnenwellen,
1 9 1 9 ) a n d i t s a p p r o x i m a t i o n s a s t h e f u n d a m e n t a l
mathematical formulae for describing the amplitudes of
PP reflected waves from P-wave incident plane waves in
exploration geophysics under plane wave approximation
play an important role in AVO analysis/inversion (Smith
2012b). The Zoeppritz equation gives the precise values of
intrinsic nonlinearity makes it less appropriate in practical
applications. Therefore, linearized approximations with
different parameterization of the Zoeppritz equations are
different types of linearized approximations see Russell et
al (2011). The linearized approximations are derived under
the hypothesis of weak property contrasts between layers or
limited incident angle. However, these assumptions do not
hold especially at unconformities or at interfaces between
al, 2011). Therefore, in this paper, we attempt to derive an
approximation of the PP reflection coefficient to adjust to
high contrast situations.
We utilize the forward modeling and nonlinear inversion
method to test the feasibility and suitability of this novel
approximation. Two forward modeling models with different
degrees of property contrast are established and we compare
the reflection coefficients with the novel approximation,
exact Zeoppritz equation and linearized approximation,
respectively. As for the inversion method, an artificial
neural network nonlinear direct inversion is introduced to
estimate the six parameters with the novel approximation
as a forward solver. The artificial neural network nonlinear
a kind of nonlinear direct inversion approach rather than an
optimization approach. It has been proved that inversion is
inverse of G( ) , which is the forward solver. It can provide
several solutions like the multiple realizations in stochastic
inversion by Bayesian inference (Buland and Omre, 2003).
(2009). In appendix A hereunder, we will give the necessary
description of this method for the nonlinear inversion problem
with the novel approximation equation.
2 Modeling
The general theory of the P-wave reflection has been
widely discussed in the literature, so we shall reproduce
here only that required for an understanding of the notation
and terminology that we will use in this paper. For the cases
of incident longitudinal waves polarized in the plane of
be expressed as
(Aki and Richards, 1980)
,
RPP
the transmission angles of the longitudinal wave and shear
RPP is the reflection coefficient of the
VP1, VS1 and 1 are the P-wave velocity,
S-wave velocity and density in medium 1, and VP2, VS2 and 2
is the ray parameter.
The linearized approximation of Eq. (1) is given by Aki
and Richards (1980) as,
(2)
(3)
(4)
where
V
VP
VP
VS
VS
RPP
1 4VS2 p2
black), linearized approximation equation (2) (dashed blue)
and our novel nonlinearized approximation equation (6) (red
dots), and we can see that reflection coefficients with these
three equations show good similarity in the weak contrast
case. Fig. 2 displays the result of the model two, we can see
show high errors compared to that from the exact equation,
and the reflection coefficient from our novel approximation
still shows high similarity to that from the exact equation at a
moderate incident angle.
3 Nonlinear inversion
To test the possibility of estimating parameters with
our novel approximation, we introduce an artificial neural
2009). The inversion is formulated in a direct inversion
scheme utilizing Eq. (6) as the forward solver. We restrict our
computational domain to various two-layer models to test the
effectiveness and potential of our novel approximation in high
contrast media. Similar to Rabben et al (2008), we attempt to
Taking the Taylor expansion for incident angle of Eq. (1),
we can express RPP in a closed form as,
VP1VP22VS1 1 2
2
VP1VP22VS2 2
3
8VS41 13 16VS21VS22 12 2
8VS42 1 2
2
A1VP1VP2 2
2VP1VP22VS1VS2 1 2
2
8VP2VS21VS2 12 2
4VP22VS1VS2 1 2
2
8VP2VS1VS22 1 2
2
2VP2 2
A2
1
(6)
the wavelet estimation and the convolution in the modeling.
We utilized the exact equation (1) as the synthetic model.
The model vector m comprises of P-wave velocities, S-wave
velocities and densities in the upper and lower layers. The
observed data d
incident angles of Eq. (1).
The model parameters and the data parameters can be
related through the nonlinear forward mapping,
d
G m
Model One
Here, we attempt to search for all possible solutions
of model parameters to satisfy the observed reflection
method, we suppose the inverse to G m is G d .
Although the inverse mapping may not exist in entire spaces,
is so smooth that the inverse
inside which the mapping G
of G does exist,
G d
m
The artificial neural network direct inversion method is
an inverse (not optimizing) algorithm, utilizing numerical
approximation of Eq. (8) in empirically constrained
subspaces. Supposing there exists an inverse mapping
and its numerical approximation inside these subspaces, it
works simultaneously with a population of several so-called
individuals. Each individual contains a parameter vector, a
data vector and the model error. The model error is used for
and for their sorting from the best to the worst model. The
computation records already evaluated and tested models so
that these models can be reused later. Repeated usage of some
models generates the possibility of efficient inversion with
minimum number of forward evaluations. Besides, several
solutions can be expected with this algorithm. Details of
reproduce the algorithm in Appendix A but only to a level
required for an understanding of the notation and terminology
that we use in our examples and discussion section.
4 Examples
Various two-layer models are established in the inversion
test. The first one is a gas sand/shale model. The P-wave
velocity (VP1), S-wave velocity (VS1) and density (Density
3,
respectively, while the P-wave velocity (VP2), S-wave velocity
(VS2) and density (Density 2) in the shale sand is 3,048 m/s,
1,244 m/s and 2,400 kg/m3, respectively, and we refer to this
model, Eq. (1) is used to generate RPP at different incident
angles to simulate observed data. The introduced nonlinear
inversion method is then used to generate ten solutions for
each of the six parameters. With the introduced inversion
method, the results of parameters estimation are displayed
in Fig. 4. Taking VP1 for example, the left figure shows the
comparison between the estimated solutions and the true
value. The sequence numbers from 1 to 10 indexes the ten
solutions for VP1, the value at sequence number 11 gives the
average value of the ten estimated values, and the value at
sequence number 12 gives the true value. The right figure
displays the relative error between the estimated solutions and
the true value. The sequence numbers from 1 to 10 indexes
the relative error between each estimated solution for VP1 and
the true value, sequence number 11 gives the relative error
between the average value of all ten values and the true value.
4, we can see that all six parameters can be inverted well and
the relative error is around 3% for each parameter.
from medium 1 to medium 2 in the second model. It shows
respect to the parameters in medium 1. Fig. 6 displays the
result of estimating six parameters with the second model
30
25
%
i,o 20
t
a
r
e
ng 15
a
h
c
itve 10
a
l
e
R 5
0
1 2 3 4 5 6 7 8 9 10 11 12
Sequence number
1 2 3 4 5 6 7 8 9 10 11
Sequence number
10
%
,ro 6
r
r
e
e 5
v
it
a
le 4
R
9
8
7
3
2
1
0
10
9
8
7
%
,ro 6
r
r
e
e 5
v
it
lea 4
R
3
2
1
0
3500
3000
2500
2000
s
/
m
,1
VP1500
1000
500
0
3500
3000
2500
/s 2000
m
,2P
V1500
1000
500
0
1 2 3 4 5 6 7 8 9 10 11 12
Sequence number
1 2 3 4 5 6 7 8 9 10 11
Sequence number
1 2 3 4 5 6 7 8 9 10 11 12
Sequence number
1 2 3 4 5 6 7 8 9 10 11
Sequence number
(Continued)
1 2 3 4 5 6 7 8 9 10 11 12
Sequence number
1 2 3 4 5 6 7 8 9 10 11
Sequence number
1 2 3 4 5 6 7 8 9 10 11 12
Sequence number
1 2 3 4 5 6 7 8 9 10 11
Sequence number
when the maximum incident angle is 31 degree. We can
see that even in high contrast media, the inversion with the
exact reflection coefficient equation still estimates the six
parameters reasonably.
Fig. 7 to Fig. 12 display the inversion results with the 2-D
surface model. Taking the P-wave velocity in upper medium
for example, Fig. 7(a) - 7(c) display P-wave velocity in upper
medium of the true model, inverted result and the relative
error between the true model and inverted result, respectively.
in Fig. 7 to Fig. 12. From the inverted results, we can see
that, with the high-contrast approximation and the nonlinear
inversion algorithm, we can obtain reasonable inversion
results, and the relative error is around 4% for each parameter.
3500
3000
2500
3
1 2 3 4 5 6 7 8 9 10 11 12
Sequence number
1 2 3 4 5 6 7 8 9 10 11
Sequence number
5500
5000
4500
4000
3500
/s 3000
m
,2P2500
V
2000
1500
1000
500
0
5500
5000
4500
4000
3500
/s 3000
m
,1S 2500
V
2000
1500
1000
500
0
5500
5000
4500
4000
3500
/s 3000
m2S2500
V
2000
1500
1000
500
0
3
2
1
0
10
9
8
7
,r%6
o
r
re 5
e
v
it
la 4
e
R
3
2
1
0
3
2
1
0
1 2 3 4 5 6 7 8 9 10 11 12
Sequence number
1 2 3 4 5 6 7 8 9 10 11 12
Sequence number
1 2 3 4 5 6 7 8 9 10 11
Sequence number
1 2 3 4 5 6 7 8 9 10 11 12
Sequence number
1 2 3 4 5 6 7 8 9 10 11
Sequence number
10
9
8
7
%
,r 6
o
r
r
ee 5
v
it
la 4
e
R
3
2
1
0
10
9
8
7
3
2
1
0
,r%6
o
r
r
e 5
e
v
it
la 4
e
R
(Continued)
5500
5000
4500
4000
3 /m3500
g
,k 3000
1
y
its 2500
n
e
D2000
1500
1000
500
0
5500
5000
4500
4000
3 m3500
/
g
,k 3000
2
itsy 2500
n
eD2000
1500
1000
500
0
1 2 3 4 5 6 7 8 9 10 11 12
Sequence number
1 2 3 4 5 6 7 8 9 10 11
Sequence number
(a)
(b)
(c)
3.5
3.0
2.5
2.0
1.5
1.0
0.5
20 40 60 80 100 20 40 60 80 100
In line In line
Fig. 8 P-wave velocities in the lower medium (a) True model, (b) Inverted result, and (c) Relative error
20 40 60 80 100
Inline
20 40 60 80 100 20 40 60 80 100
Inline Inline
Fig. 9 S-wave velocities in the upper medium (a) True model, (b) Inverted result, and (c) Relative error
20 40 60 80 100
In line
20 40 60 80 100
In line
(c)
20 40 60 80 100
In line
(a)
20 40 60 80 100
In line
20 40 60 80 100 20 40 60 80 100
In line In line
Fig. 11 Densities in the upper medium (a) True model, (b) Inverted result, and (c) Relative error
2180
2160
2140
2120
2100
10
20
30
ine 40
l
ss 50
roC 60
70
80
90
2220
2200
5 Conclusions
In this paper, we derived a high-contrast approximation
of the exact PP reflection coefficient in terms of six
parameters including P-wave velocities, S-wave velocities
and densities in upper and lower layers around a reflector.
We utilized the forward modeling and inversion method to
test the validity and feasibility of this novel approximation.
Forward modeling tests demonstrated the priority of the novel
approximation to the linearized approximation in reflection
coefficient modeling. A nonlinear direct inversion method
was introduced to estimate the six layer parameters around
multi-parameters with our novel nonlinearized approximation
of exact reflection coefficient equation could still get
reasonable inversion result even when the parameters in both
Acknowledgements
We would like to acknowledge the sponsorship of the
National 973 Program of China (2013CB228604) and
the National Grand Project for Science and Technology
acknowledge the support of the Australian and Western
Ve n t u r e P a r t n e r s , a s w e l l a s t h e We s t e r n A u s t r a l i a n
Energy Research Alliance (WA:ERA), and Foundation
WX2013-04-01).
Buland A and Omre H. Bayesian linearized AVO inversion. Geophysics.
prediction from seismic prestack data.
Geophysics. 2008
. 73(3):
C13-C21
Technical Program Expanded Abstracts. 2004
continental sandstone interbeds and deep volcanic rocks. Applied
Karimi O, Omre H and Mohammadzadeh M. Bayesian closed-skew
Gaussian inversion of seismic AVO data for elastic material
amplitude polynomial methods for characterisation of hydrocarbon
seismic data from unconsolidated sediments containing gas hydrate
and free gas. Geophysics. 2004. 69(1): 164-179
Rabben T E, Tjelmeland H and Ursin B. Nonlinear Bayesian joint
Rimstad K, Avseth P and Omre H. Hierarchical Bayesian lithology/
poroelasticity. Geophysics. 2011. 76(3): C19-C29
nonlinear equations. Technical Computing Prague. ISBN
978-807080-733-0. Praha. 2009. p.90
Shuey R T. A simplification of the Zoeppritz equations. Geophysics.
inversion of PP reflections from plane interfaces using effective
Smith G C and Gidlow P M. Weighted stacking for rock property
estimation and detection of gas. Geophysical Prospecting. 1987.
Ulvmoen M and Omre H. Improved resolution in Bayesian lithology/
fluid inversion from prestack seismic data and well observations:
Ayzenberg M, Tsvankin I, Aizenberg A, et al. Effective reflection
coefficients for curved interfaces in transversely isotropic media.
Ulvmoen M, Omre H and Buland A. Improved resolution in Bayesian
lithology/fluid inversion from prestack seismic data and well
Bortfeld R. Approximations to the reflection and transmission
B73-B82
Ursin B, Bauer C, Zhao H S, et al. Combined seismic inversion and
gravity modeling of a shallow anomaly in the southern Barents Sea.
The artificial neural network inversion was initially
direct inversion approach rather than optimization approach.
It mainly contains the following steps.
Problem initialization
starting model population, even if doing so is easy. We just
mimin < mi < mimax
(A-1)
The starting population of models is generated from
the defined range with uniform probability. The number of
models in the starting population q is not very important
because it will change during iterations. A suitable choice
can be 30 q 60 . At each iteration, the current population
of models M B mB , dB , err B is sorted according to the
individual errors between the models and the candidate
solution. The diameter of the population R defines the size
of a subspace, inside which the next population of models
will be generated, and the index of the prediction function
ip specifies the prediction method used for predicting the
candidate solution, including linear regression (ip=1),
radial basis function network (RBFN) (ip=2), and Kriging
prediction (ip=3). Both of these parameters can be tuned
during the inversion, but in the beginning they are both set
to 1.
Prediction of population and candidate solution
There is a geometrical criterion in distinct iteration cycles
before population predicting, one model is selected as the
center of the population (mC), and the other surrounding
models are located randomly in the distance R measured from
the center of the population.
The prediction population is generated in such a way
that the center is located close to the expected solution, and
surrounding models are located randomly along the surface of
a hypersphere with diameter R and center mC. Provided both
the diameter R and the center mC are known, new models of
the predicting population can be obtained as follows. Firstly,
The matrix tensor Cm whose value on the principal diagonal
line is defined as the square of the difference between the
maximum and minimum of each model parameter can be
decomposed using a Choleski decomposition as Cm=L LT.
Secondly, a random six dimensional unit vector g is
generated. Then, the proposed candidate model can be
expressed as,
mC
RL g
(A-2)
In any case when the candidate model mC is outside the
parametric hypercube, it is projected along the direction
(mg mC) to the closest face of the parametric hypercube.
The archive of already evaluated models is checked and the
model {mk, dk, errk} is selected, and this model is the closest
archive model to mg and still not connected to the predicting
population. We need to compute the distance sg between this
model and the candidate model. If the distance sg is smaller
than the distance between neighboring surrounding models,
the model {mk, dk, errk} is connected to the population,
otherwise, the candidate model mg is evaluated and is
connected to the predicting population and copied to the
archive for future use. Finally, the second step made 1
times to obtain a population of total size q.
The prediction population above can be used for
estimating the solution (8) with the prediction algorithm.
The prediction algorithm is implemented for the local
approximation of the inverse mapping. Different prediction
algorithms with different index ip can be selected to estimate
the solutions, however the best solution is often to use
different prediction algorithms even inside each individual
inverse problem. Therefore, in our case, we use the prediction
algorithm in a cyclic manner according to the variable ip.
Freeman and Co. 1980 Alemie W and Sacchi M D. High-resolution three-term AVO inversion by means of a Trivariate Cauchy probability distribution . Geophysics.
Geophysics . 2003 . 68 ( 4 ): 1140 -1149 Vedanti N and Sen M K. Seismic inversion tracks in situ combustion: A case study from Balol oil field , India. Geophysics . 2009 . 74 ( 4 ): B103 -112 Wang H Y , Sun Z D , Wang D , et al. Frequency-dependent velocity prediction theory with implication for better reservoir fluid its application . SEG Technical Program Expanded Abstracts . 2008 elastic impedance for hydrocarbon-sand discrimination . Geophysics . Zhang R, Sen M K and Srinivasan S . A prestack basis pursuit seismic inversion . Geophysics . 2013 . 78 ( 1 ): R1 - R11 lithology and fluid discrimination . SEG Technical Program Expanded Abstracts . 2010 Zhu X and McMechan G A. Elastic inversion of near- and postcritical reflections using phase variation with angle . Geophysics . 2012 .
1919 . 1 : 66 - 84 Zong Z Y , Yin X Y and Wu G C. Elastic impedance variation with Engineering . 2012a . 9 ( 3 ): 247 Zong Z Y , Yin X Y and Wu G C. AVO inversion and poroelasticity with P- and S-wave moduli . Geophysics . 2012b . 77 ( 6 ): N29 - N36