Russell - Introduction to Seismic Inversion Methods

Introduction

to

Seismic Inversion Methods

Brian H. Russell Hampson-Russell SoftwareServices,Ltd. Calgary,Alberta

Course Notes Series, No. 2 S. N. Domenico, Series Editor

Societyof Exploration Geophysicists

Thesecoursenotesare publishedwithoutthe normalSEGpeerreviews. They havenot beenexaminedfor accuracyand clarity.Questionsor commentsby the readershouldbe referreddirectlyto the author.

ISBN 978-0-931830-48-8 ISBN 978-0-931830-65-5

(Series) (Volume)

Libraryof Congress CatalogCardNumber88-62743 Societyof Exploration Geophysicists P.O.

Box 702740

Tulsa, Oklahoma 74170-2740

¸ 1988 by the Societyof Exploration Geophysicists All rightsreserved. Thisbookor portionshereofmaynotbe reproduced in anyformwithoutpermission in writingfromthe publisher. Reprinted1990, 1992, 1999, 2000, 2004, 2006, 2008, 2009 Printed in the United States of America

]:nl;roduct1 on •o Selsmic I nversion •thods

Table

Bri an Russell

of Contents PAGE

Part

I

Introduction

1-2

Part

Z

The Convolution

Model

2-1

2.4 The Noise Component

2-2 2-6 2-12 2-18

Recursive

3-1

2.1 Tr•e Sei smic Model 2.2 The Reflection Coefficient 2.3

Part

Part

3

4

The Seismic

Wavelet

Inversion

- Theory

3.1

Discrete

3.2

Problems encountered

Inversion

3.3

Continuous

with

data

4-1

4. !

4-2 4-4 4-6 4-12 4-14

4.4 4.5

I ntroduc ti on resolution

Lateral resolution Noise attenuation

Recursive

Inversion

- Practice

5.3 Seismically derived porosity Sparse-spike Inversi on

6-1

6.1

6-2 6-4 6-22 6-30

The recursive

inversion

method

5.2 Information in the low frequency component 6

I ntroduc ti on

6.2 Maximum-likelihood 6.3 6.4 P art

Part

7

8

5-1 5-2 5-10 5-16

5.1

P art

3-2 3-4 3-8

Seismic Processing Considerati ons 4.2 Ampli rude recovery

5

real

Inversion

4.3 Improvementof vertical

Part

Series

aleconvolution and inversion

The L I norm method Reef Problem

I nversion applied to Thi n-beds

7-1

7.1 Thin bed analysis

7.Z Inversion comparison of thin beds

7-2 7-4

Model-based

8-1

Inversion

B. 1 I ntroducti

8.2

Generalized

on .

linear

inversion

8.3 Seismic1ithologic roodelling (SLIM) Appendix8-1 Matrix applications in geophysics

8-2 8-4 8-10 8-14

Introduction

Part

9

to Seismic

Travel-time

Inversion

Methods

Brian

9-1

Inversion

g. 1. I ntroducti

on

9.2 Numerical examplesof traveltime 9.3 Seismic Tomography

inversion

Part 10 Amplitude versus offset (AVO) Inversion 10.1 AVO theory 10.2 AVO inversion by GLI

9-2 9-4 9-10

10-1 10-2 10-8

Inversion

11-1

I ntroduc ti on

11-2 11-4

Part 11 Velocity

Theory and Examples Part 12 Summary

12-1

Russell

Introduction

to Seismic •nversion

Methods

Brian Russell

PART I - INTRODUCTION

Part

1 - Introduction

Page 1 -

1

Introduction

to Seismic

Inversion

Methods

Brian

Russell

I NTRODUCT ION TO SEI SMIC INVERSION METHODS ,

__

_•

i

Part

This

i

_

i

,

.

,

-

,

!

•

_,

l_

,

Introduction _

,

.

i.,.

_

.

course is intended as an overview of the current techniques used in

the inversion of seismic data. It would therefore seemappropriate to begin by defining what is meant by seismic inversion. The most general definition is

as fol 1 ows'

Geophysical inversion

involves

mapping the physical structure and

properties of the subsurface of the earth using measurementsmadeon the surface of the earth.

The above definition

work that is done in

is so broad that it encompassesvirtually

seismic analysis

course we shall primarily 'restrict

and interpretation.

all

the

Thus, in this

our discussion to those inversion

methods

which attempt to recover a broadband pseudo-acoustic impedance log from a band-1 imi ted

sei smic trace.

Another way to look at inversion is to consider it as the technique for creating a model of the earth using the seismic data as input. As such, it

can be considered as the opposite of the forwar• modelling technique, which involves creating

a synthetic seismic

section

based on a model of the earth

(or, in the simplest case, using a sonic log as a one-dimensional model). The relationship between forward and inverse modelling is shownin Figure 1.1.

To understandseismic inversion, we must first processes involved therefore

in

the

look at the basic

creation of seismic data. convolutional

model

understandthe physical Initially,

we will

of the seismic trace

in the

time and frequencydomains,consideringthe three components of this model: reflectivity,

seismic wavelet, and noise.

_

Part

I

- Introduction

m

i

Page 1 -

--.

2

Introduction to Seismic InverSion Methods

FORWARDMODELL I NG i

m

Brian Russell

INVERSEMODELLING(INVERSION) _

,

ß

ß

_

ß

EARTH MODEL

Input'

,

Process:

Output'

MODELLING

INVERSION

ALGORITHM

ALGORITHM

SEISMIC RESPONSE i

m

mlm

ii

EARTH MODEL i

ii

Figure1.1 Fo.•ard ' andsInverse Model,ling

Part

I - Introduction

Page I -

3

Introduction.

to Seismic Inversion Methods

Brian l•ussel 1

Once we have an understanding of these concepts and the problems which

can occur, we are in a position to

look at

the methodswhich are currently

ß

used to invert seismic data.

These methods are summarizedin Figure 1.2.

be on poststack seismic inversion where

primary emphasis of the course will the

ultimate

resul.t,

The

o

as was previously

Oiscussed, is a

pseudo-impeaance

section.

We will

start by looking at

inversion,

the

most contanonmethods of

which are based on single trace recursion.

these recurslye relationship

inversion

procedures,

it

the

unUerstand

important to look at the

between aleconvolution anU inversion, and how Uependent each

method is on the deconvolution scheme Chosen.

classical

is

To better

poststack

Specifically,

we will

consider

"whitening" aleconvolutionmethods, wavelet extraction methods, and

newer sparse-spike deconvolution methods such as Maximum-likelihood

deconvolution

and the

L-1

norm metboa.

Another important type of inversion methodwhich will be aiscussed is model-based inversion, where a geological moael is iteratively upUatedto finU the

best

fit

with the seismic data.

After this,

traveltime

tomography,will be discussedalong with several illustrative

inversion,

or

examples.

After the discussion on poststack inversion, we shall move into the realm of pretstack. These methoUs,still fairly new, allow us to extract parameters other than impedance, such as density and shear-wave velocity. Finally,

we will

aiscuss the geological aUvantages anU limitations

each seismic inversion roethoU,looking at examples of each.

Part

1 -

Introduction

Page i -

of

Introduction to SelsmicInversion Methods

Brian Russell

SEI SMI C I NV ERSI ON

.MET•OS,,,

POSTSTACK

PRESTACK

INVERSION

INVERSION

i

EF IEL D MODEL-BASED I RECURSIVEWAV TRAVELTIME

INVERSION • ,INVESION

INVERSION

,,

I METHODS ]

!TOMOGRAPHY) - "NARROWSPARSEBAND

Figure 1.2

Part

1 - Introuuction

--

LINEAR

NVERSIOUMETHODS i

i

SPIKE

A summaryof current inversion techniques.

Page 1 -

Introduction


Brtan Russell

PART 2 - THECONVOLUTIONAL MODEL

Part

2 - The Convolutional

Model

Page 2 -

Introduction


Part 2 -

Brian Russell

The Convolutional

Mooel

2.1 Th'e Sei smic Model

The mostbasic and commonly used one-Oimensionalmoael for the seismic trace is referreU

to as the convolutional

moOel, which states that the seismic

trace is simplythe convolutionof the earth's reflectivity with a seismic source function

with the adUltion of a noise component. In

equation

form,

where * implies convolution,

s(t) : w(t) * r(t) + n(t)s where

and

s (t)

= the sei smic trace,

w(t)

: a seismic wavelet,

r (t)

: earth refl ecti vi ty,

n(t)

: additive

noise.

An even simpler assumptionis to consiUerthe noise component to be zero, in which case the seismic tr•½e is simply the convolution of a seismic wavelet with t•e earth ' s refl ecti vi ty, s(t)

In

= w{t) * r(t).

seismic processingwe deal exclusively with digital data, that

data sampled at a constanttime interval.

is,

If weconsiUerthe relectivity to

consist of a reflection coefficient at each time sample(som• of which can be

zero), and the wavelet to be a smooth function in time, convolutioncan be thoughtof as "replacing"eachreflection. coefficient with a scaledversion of the waveletandsumming the result. The result of this processis illustrated

in Figures 2.1 and2.Z for botha "sparse"anda "dense"set of reflection coefficients.

Notice that convolution

with

the wavelet tends to "smear" the

reflection coefficients. That is, there is a total loss of resolution,which is the ability

to resolve closely spacedreflectors.

Part 2 - The Convolutional

Model

Page

Introduction

to Seismic Inversion

Nethods

Brian Russell

WAVELET:

(a) '*• •

:

-' ':'

REFLECTIVITY

TRACE:

Figure 2.1

(a)

Convolutionof a wavelet with a sparse"reflectivity. (a) •avelet. (b) Reflectivit.y. (c) Resu1ting Seismic Trace.

!

.

i

:

!

!

:

i

,

:

i

!

i

'?t

*

(b')

:

i

ß

i

c

o

Fi õure 2.2

o

Convolution of a wavelet with a sonic-derived

reflectivity. ,

Par•

i

o

, ß ....

!

,

m


(a) Wavelet. (b) Reflectivity. i

i

L_

Model

-

o

o

"dense"

(c) SeismicTrace

'

Page 2 -

3

Introduction

to Seismic

An alternate,

Inver'sion

Methods

Brian

but equivalent, way of

the frequency domain.

Russell

looking at the seismic trace is in

If we take the Fourier transform of

the previous

ß

equati on, we may write

S(f) where

= W(f) x R(f),

S(f) = Fouriertransform of s(t), W(f) = Fourier transform of w(t), R(f) = Fourier transform of r(t),

In the above equation we see that

ana f = frequency.

convolution becomesmultiplication

in

the frequency domain. However, the Fourier transform is a complex function, and it

is normal to consiUer the amplitude and phase spectra of the individual

components. The spectra of S(f) may then be simply expressed

esCf)= ew(f) + er(f), where

I •ndicates amplitude spectrum, and 0

In the

Figure 2.3

illustrates

the convolutional model

frequency domain. Notice that the time Oomainproblem of

resolution becomesone of loss of reOuceo by the effects


loss

of

frequency content in the frequency domain.

Both the high and low frequencies of the reflectivity

Part

.

other words, convolution involves multiplying the amplitude spectra

and adding the phase spectra.

in

indicates phase spectrum.

have been severely

of the seismic wavelet.

Mooel

Page ?. -

4

Introduction to Seismic Inversion Methods

AMPLITUDE

Brian Russell

SPECTRA

PHASE SPECTRA

w (f) I

I

-tR (f)

i i

, i.

I iit

!

loo

|11

s (f) i

I

i

i!

I

Figure 2.3

Part


Convolution in the frequency domain for the time series shown in Figure 2.1.

Model

Page 2 -

Introduction

2.g

to Seismic

The Reflection l_

_

,m

i

_

_

,

Inversion

Coefficient _

_

m_

_,•

,

Methods

Brian

Russell

Series _ _

ß

_

el

'The reflection coefficient series (or reflectivity,

as it is also called)

describedin theprevious sectionis oneof thefundamental physical concepts in the seismic method. Basically, each reflection coefficient maybe thought of as the res ponse of the seismic wavelet to an acoustic impeUance change within the ear th, where acoustic impedance is defined as the proUuct of compressi onal velocity and Uensity. Mathematically, converting from acoustic involves dividing the difference in the acoustic i ropedanceto re flectivity impedances by the sum of the acoustic impeaances. This gives t•e coefficient

at

the boundary between the two layers.

reflection

The equation is as

fo11 aws:

•i+lVi+l- iVi i where

•

Zi+l- Zi i+1

r = reflection

coefficient,

/o__density, V -- compressional velocity, Z -- acoustic impeUance, and

Layer i overlies Layer i+1.

Wemust also convert from depth to time by integrating the sonic log transit times. Figure •.4 showsa schematicsonic log, density log, anU resulting acoustic impedancefor a simplifieU

earth moael. Figure 2.$ shows

the resultof converting to thereflection coefficient seriesandintegrating to

time.

It should be pointed out that this formula is true only for the normal incidence case, that is, for a seismic wave striking the reflecting interface at right angles to the beds. Later in this course, we shall consider the case of

Part

nonnormal

inciaence.


Model

Page 2 -

6


STRATIGRAPHIC

Brian Russell OENSITY

SONICLOG

SECTION

30O

4OO

loo

200

2.0

3.0

,

ß

I

SHALE

LOG.

•T (•usec./mette)

.....

•

3600 m/s

DEPTH

_

SANOSTONE

ß

ß ß

ß

. .

'I

!

- .. ,

!_1

!

ß ß

!

v--I

!

UMESTONEI I I ! I ! I 1

V--3600 J

V= 6QO0

LIMESTONE

2000111

I

Fig. 2.4. BoreholeLogMeasurements. REFLECTWrrY

ACOUSTIC

VS TWO.WAY TIME

IMPED,M•CE (2•

(Y•ocrrv mm

,mm

mm

mm

rome

m

m

-----

V$ OEPTH

x OEaSn• 20K -.25 I

.am

O

Q.2S I

-.25 v

O '

+ .2S I

mm

SHALE .....

OEPTH

•--------'-[

SANDSTONE . . ... !

I

!11

,

I1

UMESTONE I I 1 I I I II i ! I 1 i I i SHALE •.--._--.---- •

1000m

-

1000 m

--

NO

•.'•

,•

LIMESTONE

- 20o0 m

2000 m

Fig.

2.5.

, ..

Creation of Reflectivity

Part g - The Convolutional Model

I SECOND

Sequence.

Page 2 -

7

IntroductJ

on 1:o Sei stoic Inversion

Our best method of

derlye

Herhods

observing

them from well log curves.

Bri an Russell

seJsm•c impedance and reflectivity

is •o

Thus, we maycreate an impedancecurve

by

multiplying together •he sonic and density logs from a well. Wemay•hen computethe reflectivlty by using •he formula shownearlier. Often, we do not have the density log available• to us and must makedo with only the sonJc. The

approxJmatJonof velocJty to •mpedance 1s a reasonable approxjmation, and

seemsto holdwell for clas;cics and carbonates(not evaporltes, however). Figure 2.6 showsthe sonic and reflectJv•ty traces from a typJcal Alberta well after they have been Jntegrated to two-way tlme. As we shall see later,

the type of

aleconvolution and inversion used is

dependent on the statistical assumptionswhich are made about the seismic reflectivity and wavelet. Therefore, howcan we describe the reflectivity seen in

a

well?

reflectivity

The

traditional

to be a perfectly

answer has always been that

we consider

random sequence and, from Figure •.6,

the

this

appears to be a good assumption. A ranUomsequencehas the property that

autocorrelation is a spike at zero-lag. autocorrelation are zero except the

its

That is, all the componentsof the

zero-lag value, as shownin the following

equati on-

t(Drt = ( 1 , 0 , 0 , .........

)

t zero-lag.

Let

2.7.

us test this idea on a theoretical

Notice that the

autocorrelation

of

random sequence, shownin

Figure

this sequence has a large spike at

ß

the zeroth lag, but that there is a significant noise component at nonzero lags. To have a truly random sequence, it must be infinite in extent. Also on this figure is shown the autocorrelation of a well log •erived

reflectivity. Wesee that it is even less "random"than the randomspike sequence. Wewill discuss this in more detail on the next page.

Part


Model

Page 2 -

8

IntroductJon

to Se•.s=•c Inversion

Methods

Br•an

Russell

RFC

F•g. 2.6. Reflectivitysequence derivedfromsonJc .log.

RANDOM

SPIKE SEQUENCE

AUTOCORRE•JATION OF RANDOMSEQUENCE

Fig.

2.7.

WELL LOG DERIVED REFLECT1vrrY

AUTOCORRELATION

OF REFLECTIVITY

Autocorrelat4ons of random and well log

der4vedspike sequences.

Part


Model

Page 2-

Introductlon

to Sei smic Inversion

Methods

Therefore, the true earth reflectivity

truly

Brian Russel 1

cannot be consideredas being

random. For a typical Alberta well we see a numberof large spikes

(co•responding to majorlithol ogic change)sticking up abovethe crowd.A good way to describethis statistically is as a Bernoulli-Gaussian sequence. The Bernoulli part of this term implies a sparsenessin the positions of the spikes and the Gaussianimplies a randomness in their amplitudes. Whenwe generatesuch a sequence,there is a term, lambda, which controls the sparsenessof the spikes. For a lambdaof 0 there are no spikes, and for a lambdaof 1, the sequence is perfectly Gaussian in distribution. Figure 2.8 shows a number of such series for different

typical Alberta well log reflectivity

values of lambda.

Notice that

a

wouldhavea lambdavalue in the 0.1 to

0.5 range.

Part


Model

Page 2 -

10

I ntroducti on to Sei smic I nversi on Methods

Brian Russell

It

tl

I

I

I

i

I

I

•11 I

LAMBD^•0.01

511 t

•tl

I

(VERY SPARSE)

11

311

I

4#

I

511 I

#1

TZIIE

LAMBDA--O.

I

(KS !

1

1,1

::.•"• •'•;'"' "";'•'l•' "••'r'•

-• "(11 I TX#E

(HS)

LAMBDAI0.5

LAMBDA--

1.0 (GAUSSIAN:]

EXAMPLESOF REFLECTIVITIES

Fig.

2.8.

Examplesof reflectivities factor

,

,

m

i

ß

to be discussed

using lambda

in Part

6.

i

Part 2 - The Convolutional Model

Page 2 -

11

Introduction

2.3

to Seismic Inversion ,Methods

The Seismic --

_

Brian Russell

Wavelet

ß

•

,

Zero Phase and Constant Phase Wavelets m _

m _

m

ß

m

u

,

L

m

_

J

The assumptiontha.t there is a single, well-defined wavelet which is convolved with the reflectivity to producethe seismic trace is overly simplistic. Morerealistically, the wavelet is both time-varying and complex in shape. However,the assumptionof a simple wavelet is reasonable, and in this

section

we shall

consider

several

types

of

wavelets

and

their

characteristics.

First,

let us consider the Ricker wavelet, which consists of a peak and

two troughs, or side lobes. The Ricker wavelet is dependentonly on its dominant frequency, that is, the peak frequencyof its a•litude spectrum or the inverse of the dominantperiod in the time domain(the dominantperiod is

found by measuringthe time from troughto trough). TwoRicker wave'lets are shownin Figures 2.9 and 2.10 of frequencies 20 and 40 Hz. Notice that as the anq•litude spectrumof a wavelet .is broadened,the wavelet gets narrower in the

timedomain, indicatingan increase of resolution.Ourultimatewaveletwould be a spike, with a flat amplitude spectrum. Sucha wavelet is an unrealistic goal in seismic processing, but one that is aimedfor. The Rtcker wavelets of

Figures 2.9

and 2.10

are also zero-phase, or

perfectly symmetrical. This is a desirable character.tstic of wavelets since the energy is then concentrated at a positive peak, and the convol'ution of the wavelet

with

a reflection

coefficient

will

better

resolve

that

reflection.

To

get an idea of non-zero-phase wavelets, consider Figure 2.11, where a Ricker wavelet

has been rotated by 90 degree increments, and Figure 2.12, where the

samewavelet has been shifted by 30 degree increments.

Notice that the 90

degree rotation

180 degree shift

displays perfect antis•nmnetry, whereas a

simply inverts the wavelet.

Part


The 30 degree rotations are asymetric.

Model

Page 2-

•2


Brian Russell

Fig.

2.9.

20 Hz Ricker

Wavelet'.

Fig.

•.10.

40 Hz Ricker

wavelet.

Fig.

2.11.

Ricker

wavelet

rotated

by 90 degree increments

Fig.

2.12.

Ricker

wavelet

rotated

by 30 degree increments

Part


Model

Page 2 -

13

Introduction

to Seismic

Of course,

Inversion

a typical

frequencies than that

fil•er

Methods

seismic wavelet

shownin Figure 2.13,

would be noticeable

a larger

Consider the

range of banapass

where we have passed a banaof frequencies has also had cosine tapers applied between 5

and 15 Hz, and between60 and 80 Hz. box-car.

contains

shownon the Ricker wavelet.

between15 and 60 Hz. The filter

that

Brian Russell

The taper reduces the "ringing" effect

if the wavelet amplitude spectrum was a

The wavelet of Figure 2.13 is

simple

zero-phase, and would be excellent as

a stratigraphic wavelet. It is often referred to as an Ormsbywavelet. Minimum Phase Wavelets

The concept of minimum-phaseis one that

is

vital

to aleconvolution, but

is also a concept that is poorly understood.

The reason for

understanding is that most discussions of

concept stress the mathematics

at

the

expense of

the

physical

the

interpretation.

use of minimum-phaseis adapted from Treitel

The

this lack of

definition

we

and Robinson (1966):

For a given set of wavelets, all with the sameamplitude spectrum,

the minimum-phase waveletis the onewhichhasthe sharpest leading edge. That is, only wavelets which have positive

The reason that

time values.

minimum-phase concept is important to us is

that

a

typical wavelet in dynamite work is close to minimum-phase. Also, the wavelet from the

seismic instruments

equivalent of the 5/15-60/80

in the aefinition

is

also

zero-phase wavelet is shownin Figure 2.14.

Part

As

used, notice that the minimum-phasewavelet has no component

prior to time zero and has its energy

possible.

minimum-phase. The minimum-phase

concentrated

as close to the origin

as

The phase spectrum of the minimum-waveletis also shown.


Model

Pa.qe 2 -

14

I•troduct•on to Seistoic!nversionNethods. ql

Re• R f1.38

Zero Phase I•auel•t

5/15-68Y88

-

0.6

- e.3e

{•

Trace

1

iii

...... ,

• .....

'

1

2be Trace

Reg 1)

Br•anRussell

min,l•

wavelet

Fig.

2.13. Zero-phase bandpass

Fig.

2.14. Minim•-phase equivalent of zero-phase wavelet shownin Fig. 2.13.

wavelet.

I

•/15-68/88 hz

18.00 p

Trace I

RegE

wayel

Speetnm

'188.88•

Trace1

0.8

188

_

!

m,m,

i

m

Part 2 -Th 'e Convolutional Model

i

Page 2-

15

Introduction


Brian Russell

Let us nowlook at the effect of different waveletson the reflectivity function itself.

Figure 2.15 a anU b shows a numberof different

conv6lved with the reflectivity in Figure Z.5.

(Trace 1) from the simple blocky model shown

The following wavelets have been used- high frequency

zero-phase (Trace •),

low frequencyzero-phase(Trace ½), high frequency

minimumphase (Trace 3), figure,

wavelets

low frequency minimum phase (Trace 5).

From the

we can make the fol 1owing observations:

(1) Low freq. zero-phase wavelet: - Resolution of reflections

- Identification

(Trace 4) is poor.

of onset of reflection

is good.

(Z) High freq. zero-phase wavelet: (Trace Z) - Resolution of reflections

- Identification

is good.


is good.

(3) Lowfreq. min. p•ase wavelet- (Trace 5) - Resolution

of reflections

i s poor.

- Identification of onset of reflection is poor. (4) High freq. min. phase wavelet: (Trace 3) - Resolution of refl ec tions is good.

- Identification


is poor.

Based on the aboveobservations, we wouldhave to consider the high

frequency,zero-phase waveletthe best, andthe low-frequency, minimum phase wavelet

Part

the

worst.


Model

Page 2 -

16

Introduction

!ql

RegR


Zer• PhaseUa•elet

•,'1G-•1•

Methods

14z

Russell

q2 RegC ZeroPhase 14aue16(' ' •'le-3•4B Hz e

F

- •.• ['

Brian

'

•,3 RecjB miniilium phue

'

q• Reg1) 'minimum phase "

'

•,leJ3e/4eh•

'

8

17 .•

e.e

(a)

•/••/'•-•"v--,._,, -r

e.•' ' "s•e'' ,m

,,

Tr'oce

[b) 700

Fig.

2.15.

Convolution of four different in (a) with trace I of (b). shown on traces

Part


Model

wavelets shown The results are

2 to 5 of (b).

Page 2 -

17

Introduction


Methods

Brian Russell

g.4 Th•N.oi se.Co.mp.o•ne ntThe situation

that has been discussed so far is the ideal case.

That is,

.

we have interpreted every reflection wavelet on a seismic trace as being an actual

reflection

from a

lithological

boundary.

Actually,

many of the

"wiggles"on a trace are not true reflections, but are actually the result seismic noise.

of

Seismic noise can be grouped under two categories-

(i) Random Noise - noise which is uncorrelated from trace to trace and is •ue mainly to environmental factors.

(ii)

CoherentNoise - noise which is predictable on the seismic trace but

is unwanted. An exampleis multiple reflection interference. Randomnoise can be thought of as the additive componentn(t) which was

seen in the equationon page 2-g. Correcting for this term is the primary reason for stackingour •ata. Stackingactually uoesan excellent job of removing ranUomnoise.

Multiples, one of the major sources of coherent noise, are causedby multiple "bounces" of the seismic signal within the earth, as shownin Figure 2.16. They may be straightforward, as in multiple seafloor bouncesor "ringing", or extremelycomplex,as typified by interbed multiples. Multiples cannot be thoughtof as additive noise andmustbe modeledas a convolution with the reflecti

vi ty.

Figure 2.17

shows the

theoretical

multiple

sequence which

generatedby the simple blocky modelshown on Figure •. 5. this

data,

it

is

Multiples maybe partially

powerful

the

multiples

be

removed by

stacking,

but

important that

elimination

technique.

aleconvolution, f-k filter.ing,

Such

techniques

and inverse velocity

If

we are

effectively

would to

be

invert

removed.

often require a more include

stacking.

predictive

These techniques

wil 1 be consi alered in Part 4.

Part


Model

Page 2 -

18


Fig.

2.16. Several multiple generating mechanisms.

TIME

TIME

[sec)

[sec)

0.7

REFLECTION COEFFICIENT

SERIES

Fig.

Brian Russell

2.17.

0.7

R.C.S. WITH

ALL

MULTIPLES

Reflectivi ty sequenceof Fig.

2.5.

with

and without multipl es.

.

Part 2 - The ConvolutionalModel

Page 2 -

19

PART 3 - RECURSIVE INVERSION - THEORY m•mmm•---'

.•

,-

-

-

'

•-

Part 3 - Recurstve Inversion - Theory

-

_

-

-

_-

_

Page 3 -

•ntroduct•on

to SeJsmic Znversion Methods

Brian Russell

PART 3 - RECURSIVE INVERSION - THEORY 3.1

Discrete ,

!

ß

Inversion ,

,

•

In section 2.2, we saw that reflectivity was defined in terms of acoustic impedancechanges. The formula was written:

Y•i+lV•+l ' •iV! 2i+1' Zi ri--yoi'+lVi+l+ Y•iVi---Zi..+l + Zi where

r -- refl ecti on coefficient,

/0-- density, V -- compressional velocity, Z -- acoustic impedance, and

Layer i overlies Layer i+1.

If we have the true reflectivity available to us, it is possible to recover the a.coustic impedanceby inverting the above formula. Normally, the inverse' formulation is simply written down,but here we will supply the missing steps for completness. First,

notice that:

Zi+l+Zi

Zi+1- Zt

2 Zi+1

Zi+l+Zi

Zi+1- Zi

2 Zf[

Zi+l+ Zi

Zi+l+Zi

Zi+l+Zi

I +ri- Zi+l +Zi + Zi+l +2i

Zi+l +Zi

Also

I-

ri--

Ther'efore

Zi+l Zi

l+r.

1

1

ill,

Part 3 - RecursiveInversion- Theory

ß

,

I

Page

Introduction to Seismic Invers-•onMethods

Brian Russell

pv-e-

TIME

(sec]

0.7

REFLECTION COEFFICIENT SERIES

Fig.

,

!

Part

3.1,

RECOVERED ACOUSTIC IMPEDANCE

Applying the recursiveinversionformulato a simple,and exact, reflectivity.

ß

3 - Recursive

Inversion - Theory

Page 3 -

!ntroductt •9r•

on to Se1 smJc ! nversi on Methods

;• • •;• • • •-••

9rgr•t-k'k9r9r• •-;• ;•

Or, the final

Brian

Russell

.................................................

•esult-

l+r i

Zi+[=Z

.

ß

This is called

the discrete

recursive

inversion

formula and is the basis

of many current inversion techniques. The formula tells us that if we know the acoustic impedanceof a particular layer and the reflection coefficient at the base of that layer, we may recover the acoustic impedance of the next layer. Of course we need an estimate of the first layer impedanceto start us off. Assumewe can estimate this value for layer one. Then

l+rl

Z2: Zli r1

,

Z3= Z2 11 +r2 -

To find the nth impedancefrom the first,

Figure reflection

3.1

shows the

coefficients

application

derived

in

and so on ...

r

we simply write the formula as

of the

section

recursive

2.2.

formula to

As expected,

the

the full

acoustic impedance was recovered. Problems •

ß

encountered

,

When the

m

with

i

i

•

real i

data !

m

recursive inversion formula is applied to real data,

that two serious problems are encountered.

we find

These problems are as follows-

(i) FrequencyBandlimi ti ng _

ß

Referring back to Figure 2.2 bandlimited

when

it

we see that

is convolved

with

the reflectivity

the seismic

wavelet.

is severely Both

the

low frequency componentsand the high frequency componentsare lost.

Part 3 - Recursive Inversion

- Theory

Page 3 -

4

"

Introduction to SeismicInversion Methods

0.2

0

Brian Russell

V•) 'V,•

•R

Vo:1000 m

Where: {ASSUME j•: l)

--• V,•= 1000 i-o.t

R = +0.2

-- 1500 •ec'. m (a) - 0.1

'•0.2

R• R=

Vo=1000m

= 818ii•.m R•=-0.1 -'+ ¾1 R =+0.2 R: -0.1

Figure 3.2 Effect of banUlimitingon reflectivity, single reflection coefficient,

where(a) shows

anU (b) showsbandlimited

refl ecti on coefficient. I

__

___

i

_

i

Part 3 - Recursire Inversion - Theory

i

m

i

m

I

Page 3 -

Introduction

to Seismic

(ii)

Noise

The

inclusion

Inversion

of coherent

makethe estimate• reflectivity

Methods

or random noise

let

us first

into

the seismic

Russell

'trace

will

deviate from the true reflectivity.

To get a feeling for the severity inversion,

Brian

of

the above limitations

use simple models. To illustrate

on recursire

the

effect

of

bandlimiting, consider Figure 3.Z. It shows the inversion of a single spike (Figure 3.2 (a)) anUthe inversion of this spike convolved with a Ricker wavelet (Figure 3.2 (b)). Even with this very high frequency banUwidth

wavelet, we have totally lost our abil.ity to recover the low frequency componentof the acoustic impedance. In Figure 3.3 the model derived in

minimum-phase wavelet.

section Z.2 has been convolved with a

Notice that the inversion of the data again shows a

loss of the low frequency component. The loss of the low frequency component is the most severe problem facing us in the inversion of seismic data, for is extremely

Oifficult

to

directly

recover

spectrum, we may recover muchof

the

it.

original

At

the

it

high end of the ß

frequency content using

deconvolution techniques. In part 5 we will address the problem of recovering the low frequency component. Next,

consider

sources, but will

reflectivity. reflection train

the

problem of noise.

always tend

to

interfere

This noise with

our

may be from many recovery of the true

Figure 3.4 showsthe effect of adding the full (including

multiple

transmission losses) to the model reflectivity.

As we can see on the diagram, the recovered acoustic impedancehas the basic shape as the true acoustic impedance, but becomesincreasingly

same

incorrect

with depth. This problemof accumulatingerror is compoundeU by the amplitude problemnsintroduced by the transmission losses.

Part 3 - Recurslye Inversion - Theory

Page 3 -

6

Introduction

to Seismic Invers,ion Methods

Brian Russell

pv-•,

TIME

0.?


Fig.

3.3.


The effect

SYNTHETIC

INVERSION

(MWNUM-PHASE WAVELET)

OF SYNTHETIC

of bandlimiting on recurslye inversion.

TIME

TIME

(see)

(re.c)

0.7


Fig.

Part 3 - Recursive

3.4.


The effect

Inversion

- Theory

R.C.S. WITH ALL MULTIPLES

of noise on recursive


inversion.

Page 3 -

Introduction


3.3 Continuous

Brian Russell

Inversion

A logarithmic

relationship

is

often used to

approximate the above

formulas. This is derived by noting that we can write r(t) function in the following way:

as a continuous

r(t) -- Z(t+dt) Z(t+dt)+- Z{t) _ 1 d Z(t) Z(•) - •' z'(t)

ß Or

! d In Z(t)

r(t) = •

The inverse

formula

dt is

thust

Z(t)=Z(O) exp 2y r(t)dt. 0

The precedingapproximation is valid if case.

r(t) <10.3• which is usually the

A paper by Berteussen and Ursin (1983), goes into muchmore detail

the continuous versus discrete approximation.

Figures 3.5 and 3.6 from their

paper showthat the accuracy of the continuous inversion algorithm is 4% of the correct value between reflection

If

on

coefficients

within

of -0.5 and +0.3.

our reflection coefficients are in the order of + or - 0.1, an even

simplerapproximation maybe made by dropp'ing the logarithmicrelationship: t

1d Z(t)•_==• Z(t)--2'Z(O) fr(t) dt r(t)--• -dr VO

Part 3 - Recursive Inversion - Theory

Page 3 -

8


m

i

,,

,m

I +gt -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3

Fig. 3.5

Brian Russell

I

I

½xp (26•)

IIIII

Difference

0.0 0.05 0.11 0.18 0.25 0.33 0.43 0.• 0.667 0.8182 1.0 1.222 1.500 1.86

0.14 0. I? 0.20 0.25 0.30 0.37 0.45 0.•5 0.670 0.8187 1.0 1.221 1.492 1.82

-0.14 -0.12 -0.09 -0.07 -0.05 -0.04 ' --0.02 --0.01 -0.003 --0.0005 0.0 0.001 0.008 0.04

0.4 0.5

2.33 3.0

2.23 2.7

o.1 0.3

0.6

4.0

3.3

0.7

0.7 0.8 0.9 1.0

5.7 9.0 19.0 co

4.1 5.0 6.0 7.4

1.6 4.0 13.0 •o

Numericalc•pari son of discrete and continuous i nversi

on.

(Berteussen and Ursin, 1983)

$000 }m MPEDANCE (O ISCR. ) r-niL O

${300 -•O IFFERENCE o

SO0

OI FFERENCE( SCALEDUP)

T •'•E

Fig. 3.6

t SECONOS

C•parisonbetween impedance c•putatins basedona discrete and a continuousseismic •del.

(Berteussenand Ursin, 1983)

Part 3 - Recursire .Inversion - Theory

Page 3 -

Introduction'to

Seismic Inversion

Methods

Brian Russell

PART 4 - SEISMIC PROCESSINGCONSIDERATIONS

Part 4 - Seismic Processing Considerations

Page 4 -

1

•ntroduction

4.1

to Seismic •nvers•on Methods

B.r.ian Russell

Introduction

Havinglookedat a simple model'of the seismic trace, anu at the recursire inversion alogorithmin theory, we will nowlook at the problem of

processingreal seismiceata in order to get the best results fromseismic inversion. We may group the key processingproblemsinto the following categories:

( i ) Amp 1i tu de recovery. (i i) Vertical

resolution

improvement.

(i i i ) Horizontal resol uti on improvement. (iv)

Noise elimination.

Amplitudeproblemsare a majorconsiderationat the early processing stagesandwewill look at both deterministicamplitude recoveryandsurface consistent

residual static time corrections.

Vertical resolution improvement

will involve a discussion of aleconvolution and wavelet processing techniques.

In

our discussion of horizontal resolution we will look at

the resolution

improvement obtained in migration, using a 3-D example.Finally, wewill consider several approaches to noiseelimination, especially the elimination of multi pl es.

Simply stateu, to invert our seismic data we usually assume the And to arrive at an one-dimensional model given in the previous section. approximationof this model (that is, that each trace is a vertical, band-limited reflectivity function) we must carefully process our data with these

considerations in minU. Figure 4.1 showsa processing flow which could

be useUto do preinversion processing.


Page 4 -

2

Introduction


Methods

Brian Russell

INPUT RAW DATA

DETERMINISTIC AMPLITUDE

CORRECTIONS

,.

SURFACE-CONSI STENT AMPt:ITUDE ANAL'YSIS

_•m

mlm

SURFACE-CONS ISTENT

DECONVOLUTIO, NFOLLOWED

BY HI GH RESOIJUTI.ON DECON i

i

SURFACE-CONSISTENT STATI CS ANAIJYSIS

VELOCITY ANAUYSIS

APPbY STATICS AND VEUOCITY

MULTIPLE

ATTENUATION

STACK ß

•

MI GRATI ON ,

Fig.

ll

,

ß

'

ß

Simpl i fiedi nversi onprocessing flow.

4.1.

I

,

_

i

11

,

m


-

--

m

_

•

•

,11

Page 4 -

3

Inl;roducl:ion 1:oSeJsmlc Invers1on Nethods

4.2

BrJan Russell

Am.p'l i tu.de..P,.ecovery

The most dJffJcult job in the p•ocessing of any seismic line

is

ß

•econst•uctingthe amplJtudes of the selsmJc t•aces as they wouldhavebeenJf the•e

were no dJs[urbJnginf'luences present.

We normally make the

simplJficationthat the distortionof the seJsmic amplJtudes maybe put into three main categories'sphe•Jcaldivergence, absorptJon,and t•ansmJssion loss.

Based on a consideration of

these three factors, we maywrJte aownan

approximatefunctJonfor the total earth attenuation-

At: AO* (b

/ t) * exp(-at), where

anU

Thus, if data, formula.

the

t

= time,

At = recorded amplitude, A0 = true ampl i tude, a,b

= constants.

we estimate the constants in the above equation from the seismic

true amplitudesof the data coulUbe recoveredby using the inverse The deterministic

amplitude correction and trace

to

trace

mean

scaling will accountfor the overall gross changesin amplitude. However, there may still be subtle (or even not-so-subtle) amplitude problems associatedwith poorsurface conditions or other factors. To compensate for these effects, it is often advisable to computeand apply surface-consistent

gain corrections. This correction involves computing a total gain value for each trace and then decomposing this single value in the four components

Aij= SixRjxG kxMkX •j, where A = Total amplitude factor,

M = Offset

component,

S = Shot component,

X = Offset

distance,

R:

i,j

Receiver

component,

G = CDP component, and


= shot,receiver

pos.,

k = CDP position.

Page 4 -

Introduction to Seismic .Inversion Methods

Brian Russell

SURFACE SUEF'A•

AND

CONSIb'TEh[O{

T |tV•E :

,RiL-rER ß

Fig.

4.2.

,

Surface and sub-surface geometryand surface-consistent decomposition. (Mike Graul). ,

Part 4 - SeismicProcessing Considerations

Page 4 -

5

Introduction


Figure 4.g

Brian Russell

(from Mike Graul's unpublished course notes) shows the

geometryusedfor this analysis. Notice that the surface-consistent statics anti aleconvolutionproblemare similar. For the statics problem,the averaging

can be •1oneby straight summation.For the amplitude problemwemust transform the above equation into additive form using the logarithm:

InAij= InSi +InRj+InGk+lnkMijX•. The

problem can then be treated exactly the same way as in the

statics

case. Figure 4.3, fromTaner anti Koehler (1981), showsthe effect of doing surface consistent amplitude and statics corrections.

4.3 I•mp. rov. ement_ o.[_Ver. t.i.ca.1..Resoluti on Deconvolution

is

a

process by which an attempt is made to

seismic wavelet from the seismic trace,

remove the

leaving an estimate of reflectivity.

Let us first discussthe "convolution"part of "deconvolution" starting with the equation for the convolutional model

st--wt* r t

st= the seismictrace,

where

wt= the seismicwavelet, rt= reflection coefficient series, * = convolution operation. In

the

frequency domain S(f)

The

procedure reflection

deconvol ution

and consists coefficients.

• W(f) x R(f)

process

of

is

.

simply the

reverse

of the convolution

"removing" the wavelet shape to

reveal

the

We must design an operator to do this, as in the

fol 1owlng equati on-

rt: st* o

whereOr--operator-- inverseof wt .

Part 4 - Seismic Processing Considerations ,

Page 4 -

6

Introduction

to Seismic

Inversion

Methods

ii

Brian

Russell

11

i

ß

1'

ii

'..,•' •,

,"

"

"

ß d.

Preliminary stack bet'ore surface consistent static andomplilude corrections.

Fig.

4.3.

Stacks with and without

ß Stockwith surfaceconsistent staticand amplitudecorrections.

surface-consi stent

corrections. (TaneranuKoehler,1981). ß


,

Page4 - 7

Introduction


Brian Russell

In the frequency domain,this becomes R(f) = W(f) x 1/W(f) .

After this extremely simple introduction, it

may appear that the

deconvolutionproblemshouldbe easy to solve. This is not the case, and the continuing research into the problem testifies to this. There are two main problems. Is our convolutionalmodelcorrect, and, if the modelis correct, can we derive the

true

wavelet

from the

data?

The

answer to the first

questionis that the convolutionalmodelappearsto be the best modelwe have comeup with so far. The main problem is in assuming that the wavelet does not vary with time. In our discussionwe will assumethat the time varying problemis negligible within the zoneof interest.

The secondproblemis much more severe, since it requires solving the ambiguousproblemof separatinga wavelet and reflectivity sequencewhenonly the seismic trace is known. To get around this problem, all deconvolution or

wavelet estimation programsmakecertain restrictive assumptions,either about the wavelet or

the

reflectivity.

There are

methods: those which make restrictive

two classes of deconvolution

phase assumptions and can be considered ,

true wavelet processingtechniquesonly whenthese phaseassumptionsare met, and those

which do not

make restrictive

phase assumptions and can be

consideredas true wavelet processingmethods. In the first category are (1) Spiking deconvolution, (2) Predictive

deconvolution,

(3) Zero phase deconvolution, and (4) Surface-consi stent deconvoluti on.


Page 4 -

Introduction


Brian Russell

(a)

Fig.

A comparison of non surface-consistent and surface-consistent decon on pre-stack data. {a) Zero-phase deconvolution. {b) Surface-consistent soikinB d•convolution.

4.4

(b), Fig. 4.5

Surface-consistent decon comparisonafter stack. (a) Zero-phase aleconvolution. (b) Surface-consistent deconvol ution.

'--'-

,

ß

,

,•

,t

ß

ß

_

,

,

_ _

,,

,

Part 4 - .Seismic Processing Consioerations

,_

,

,

Page4 -

Introduction

to Seismic Invers.ion Methods

Brian Russell

In the second category are found (1) Wavelet estimation using a well 1og (Strat Decon).

(Hampson and Galbraith 1981) (2) Maximum-1ikel ihood aleconvolution.

(Chi et al, Let

us

illustrate

the

lg84) effectiveness

surface-consi

stent

aleconvolution.

surface-consi

stent

scheme involves

components. di recti

ons-

must therefore

We common

con,non offset

Referring

the

of

one

of.

to Figure 4.•,

convolutional

average

source, commonreceiver,

over

four

notice

proauct

different

ways to perform it.

that

a

of four geometry

common depth point (CDP), and

(COS). The averaging must be performed iteratively

are several different

the methods,

and there

The example in Figures 4.4 ana 4.5

shows an actual surface-consi stent case study which was aone in the following way'

(a) Computethe autocorrelations

of each trace,

(b) average the autocorrelations in each geometry eirection to get four average autocorrel ati OhS,

(c) derive and apply the minimum-phase inverse of each waveform, and (•) iterate through this procedure to get an optimumresult. Two points to note

when you are

looking

at

the

case study are the

consistent definition of the waveform in the surface-consistentapproachan• the subsequent improvementof the stratigraphic

interpretability

We can compareall of the above techniques

page.

of the stack.

using Table 4-1 on the next

The two major facets of the techniques which will

be comparedare

the

wavelet estimation procedure and the wavelet shaping procedure.


Page 4 -

10


Table

4-1

Comparison of Deconvolution MethoUs m

METHOD

Spiking

Brian Russell

m

ß

ß

m

WAVELET ESTIMATION

WAVELETSHAPING

Min.imum phase assumption Randomrefl ecti vi ty

Ideally shapedto spike. In practice, shapedto minimum

assumptions.

phase,higherfrequencyoutput.

Predi cti ve

No assumptions about

Does not whiten data well.

Deconvol uti on

wavelet•

Removes short andlong period

Deconvol ution

multiples.

Does not affect

phaseof wayelet for long lags. ..1_,

Zero

Phase

Deconvol utton

Zero phaseassumption. Randomrefl ectt vi ty

Amplitude spectrumi$

assumption.

whi tened.

m

Phase is not altered.

Canshapeto desired output. Surface-cons.

Deconvolution

Minimumor zero phase. Randomreflecti vi ty

Ampli rude spectrum i s

assumption.

whitenedless than in single

Phasecharacteri s improved.

trace

methods.

Stratigraphic

No phase assumption.

Phase of wavelet is zeroed.

Deconvol ution

However, well must match

Amplitude spectrum not

sei smi c.

whi tened.

Maximum-

No phase assumption.

Phase of wavelet is zeroed•

L ik el i hood

Sparse-spike assumption.

Amp 1i rude spectrumi s

deconvol ution


whi tened.

Page 4

11'

Introduction

4.4

to Seismic

Lateral

Inversion

Methods

Brian

Russell

Resol uti on

The complete three-dimensional (3-D) diffraction problemis shownin Figure 4.6 for a modelstudy taken fromHerman,et al (1982). Wewill look'at line 108, which cuts obliquely across a fault and also cuts across a reef-like structure.

Note that

it misses the second reef

structure.

Figure 4.7 shows the result of processing the line. In the stacked section we maydistinguish two types of diffractions, or lateral events which do not represent true geology.

The first

type are due to point reflectors

the plane of the section, and include the corners at the base of the reef structure

secondtype are out-of-t•e-plane

in

sides of the fault and the sharp

which was crossed by the line.

diffractions,

The

often called "side-swipe". This

is most noticeable by the appearance of energy from the second reef booy which was not crossed. In the two-dimensional (2-D) migration, we have correctly removed

the

2-D

out-of-the-plane problems.

diffraction

diffractions.

The final

patterns,

but

are

The full

3-D

migration

migrated section has also

still

bothere•

the

corrects for these

accounted for

positioned evehts such as the obliquely dipping fault.

by

incorrectly

This brief summaryhas

not been intended as a complete summaryof the migration procedure, but rather

as a warning that structural

migration {preferably 3-D) mustbe performedon complex

lines for the fol 1 owing reasons:

(a)

To correctly position dipping events on the seismic section, and

(b)

To remove

diffracted

Although migration

events.

can compensatefor someof the

lateral

resolution

problems, we must rememberthat this is analogous to the aleconvolution problem in that not all of the interfering

be aware that the

interference,

true

effects may be removed. Therefore, we must

one-dimensional

seismic trace,

free of any lateral

is impossible to achieve.


Page 4 -

12

Introduction


Brian Russell

71

lol

I 131

101

131

(a] 3-

108

D MODEL

ß

ß ß

ß ß ß

..................................

LINE .............................

.........................................

....................................

{hi 8•8•0 Fig. i

4.6.

LAYOU•

3-D model experiment. mm

_

ml

j

(Herman et al, 1982).

mm

Part 4 • Seismic Processing Considerations

Page 4 -

13

Introduction

4.5

Notse


Attenuation

As we' discussed in an earlier

either

Brian Russell

•andom 'or coherent.

section, seismic noise can be classified

Random noise

is

as

reduced by the stacking process

quite well unlessthe signal-to-noise ratio dropsclose to one. In this case, a coherencyenhancement programcan be used, which usually involves sometype of trace mixing or FK filtering. However,the interpreter mustbe aware that any mixing of the data will "smear" trace amplitudes, makingthe inversion result on a particular trace less reliable. Coherent

noise

is much more difficult

sources of coherent noise is multiple

Two of

to eliminate.

One of

major

interference, explained in section 2.4.

the major methodsused in the elimination of multiples

filtering

the

are

the FK

method,and the newerInverse Velocity Stacking method. The Inverse

VeiocityStackingmethod involvesthe followingsteps: (1) Correct the data using the proper NMOvelocity, (2) Model the data as a linear

sumof parabolic shapes,

(This involves transforming to the Velocity domain),

(3) Filter outthe parabolic components with a moveout greaterthansome pre-determined limit (in the order of 30 msec),and (4) Perform the inverse transform.

Figure 4.8, taken from Hampson (1986), showsa comparisonbetweenthe two methodsfor a typical multiple problemin northernAlberta. The displays are all' co•on offset stacks.

Notice

that

although both methods have performed

well on the outside traces, the Inverse Velocity Stacking methodworks best on

the inside traces. Figure4.9, also fromHampson (1986), showsa comparison of final

stacks with and without multiple attenuation.

It is obvious 'from this

comparison that the result of inverting the section whichhas not had multiple attenuation

would be to introduce spurious velocities

into the solution.

The

importanceof multiple elimination to the preprocessingflow cannot therefore be overemphasized.

m

i

i

m

,

i

.

i

Part 4 - Seismic Processing Consideration•

m

_

i

i

_

L

,=•m__

_

i

m

Page ½ -

ß

•

14


Brian Russell.

!lilt tiiti ll!1111iitt i)tt iltli ii/lit tttl• ill

(b] LINEld8 - 2-D MIGRATION

IIIIIIll!!1111111111111it I!1111111 I!11111111111illl illIiIIIIIIIIIil!111111tllilil!illlllll!111illlll [1111111111111111111111111 III!!1111 I!111111111111111 IIIIIIIilllllllll!111111111111111111111 ?•111[•i•• IIIIIIIII !1111111111111111 IIII!IIIiill•illlllillllllllllliillllllllllllh

•.,

}!l!iilll •lllllilllllll i!iiJ :illllllllllllilitiilillit!illllll ,o

111lllllllllllllllllllll1111llllll Iilllllll!ll!llll I111 illllllllilllllllllllllllllllllllllii{lilllllllllllll

"•fillllllllll!1111illi!111 IIIIIIIII IIIIIII1111111111 IIIIIlillilllllll!1111!1!111111111 '• Col LINE 108 - 3-D MIGR•ATION

F•g. 4.7. Migrationof modeldata shown in F•g. 4.6.

- -(Hermanet al, 1982). Part 4 - Seismic ProcessingConsiderations ß

Page 4 -

15

Introduction

to Seismic

Inversion

Methods

Brian

AFTER

AFTER

F-K MULTIPLE

INVERSE VELOCITY STACK MULTIPLE ATTENUATION

INPUT

Russell

ATTENUATION J. '

"

' ')'%':!•!t!'!11!1'1

';.•m,:'!:',./-•-•l- •r'm--

"';;:.m;: .... ,;lliml;..•

all

.

m#l

Fig, 4.8.

Commonoffset

stacks

calculated

from data before multiple

attenuation, after inverse velocity and after F-K multiple attenuation.

stack multiple attenuation, (Hampson, 1986)

888

Zone d Interest

1698 -4 Secondreal-dataset conventional stackwithout multiple attenuation.

'•",• ......

•,,t.•/:,.•t.,.

•

lee

,ii%' .t

---';•-•' "" • ""

).

;•,•<,:u(•:'J,.•J L,.•.,!I',,', ,'; •, • •, , ,

.% '.

• • ' "'

"•'"t"

1•%';J• ß' ß'

,

,,,, • •

I• •,,,,,ß..... '1"' ',''.•..;•t(•' •,)"•,'.m,,•""•. •,,,•• '•;••

•t•,

..•'•t,..'•"•'i•'

ß.... -....

•-

; -'".' ,•..'.... '. 2•>.': '..'•, •;,%"'•1

'" "'

Zone of

,,, .t•iill••)•.•);•l',"P,'•)'•"•'".•r'"mm"•""•P"• "•)r'" ,,..,.,,..,,,_. •,,.,....•.,...,..,...,..,.•..,....,,,.,.,.. gt••' •..,,,. ' '" •- ..... ,•

Interest

.. ,,, ••,•,•p}•h•?.•,.•,, r•.•,•. •.}•.U.•,.., ,nm, ";' ........ ,,,, ../•.• ,•.•'•, .'•-•%.... ,'......•.•_ ,,,,.,.•,,,,,.., .l•,• ), •.•, ' ,,{. • ,•,•m,l,•, .'•' ',•'•...r'• .... '.""::•"'•'•""""="'""•" Fig. 4.9.

Second real data stack after inverse velocity multiple attenuation. (Hampson, 1986)


stack

Page 4 -

16

Introduction

to Seismic Inverslon

Methods

Brian Russell

PART 5 - RECURSIVE INVERSION - PRACTICE _

Part

5 - Recursive

Inversion

_

_

_

_

- Practice

..

.

.•

,•

_

_

Page 5 -

i

Introduction

5.1

to Seismic

The Recurslye

Inversion

Inversion

Methods

Brian

Russell

Method

We have nowreached a point

where we may start aiscussing the various

algorithms currently usedto invert seismicdata. Wemustremember that all these techniquesare baseUon the assumptionof a one-aimensionalseismic trace model. T•at is, we assumethat all the corrections which were aiscussed in section 4 have been correctly applied, leaving us with a seismic section in

whic• each trace represents a vertical, band-limiteU reflectivity this section we will

series.

In

look at someof the problems inherent in this assumption.

The most popular techniquecurrently used to invert seismic Uata is

referred

.

to as recursire

inversion

and goes under such trade

names as SEISLOGana

VERILOG. The basic equations used are given in part 2, anU can be written

Zi+ 1 Zi

ri--Zi+l+ Zi where

and

The seismic data are

<===__===> Zi+l =Zi

LIJ

r i = ith reflection

,

coefficient,

Zi --/•Vi= density x velocity. simply assumeato fit

the forward model and is

inverted usingthe inverse relationship. However,as wasshownin section 3, one of t•e key problemsin the recursire inversion of seismic data is the loss of the low-frequency component. Figure 5.1 shows an exampleof an input seismic

section

aria the

resulting pseuao-acoustic impeaance without

the

incorporationof low frequency information. Notice that it resemblesa phase-shifteU version of the seismic•ata. The questionof introUuclng the low frequencycomponent involvestwo separateissues. First, wheredo we get the low-frequency component from, ana, second,howao we incorporateit?

Part

5 - Recurslye

Inversion

- Practice

Page 5 -

2.

Introduction


Brian Russell

1171121e9leS1ol 92 93

i••'••11• I I Ittltl =:::•:::::::-•--lll[l•1111t• •'• •1tlllttllltl Ill•l1t1 l!IIit! 'ti! llltfltll!!l• I I !1!t•n•'i•l, , •" • --'' ' .

ßß ,• '•

,

•.•- -.

•]

.•,.,• •4•

,•, .... _•.

• •

• f •• .

•• .•

............ .........

-::-=

• •-•,• ., • •,l•,,..,,, •,•-,, • ••..•,, ,,•,,.,,• ,•,•,I,.•.•..,• ....... •.•,,• . •.•.•.. .,%•.• ' .,•, ,.•,.... •i• • ....... ....... ••,•;.• • .,• .., .....• ....•...,• •.,.,,,,..•,.., ,,,.•,•,•.•.• ....•.,• • .,............. . •'2 •q•, •,,• ,.. J.I .,,, • .•,• ....... : ..•..... :,.. . ,• . ,'-.. ....

• ,•'•,:•'•'

....2•Y•'

___ • .... ..• ,.• -

• .......

• ... • ß•-ß •]• •,••'•,•,..,• ,•v•,• :•,.•.•l.;•.•.'..•l•l;ql .n.................... ...;•'. •;....: ••'•,J'•, .. •,• ................... '•'',,•:'•l•,fm; •,•' ,,,' ',',•, ",,','",',' •":•'•'•"•m• i•q•'t•'•'•a ....:•., •,•],' .•' ' - '""W',-•

•. ß•,'•.

,

:•,• m•,•'.

•l••J• •":•!•

•- (• ,, •

-

- ,,•.-..•.•.,'.;.',-,.. .........

7•*

. :c•i•

••'

•'•:; •, /....................

. ....

-(•-•( •.•,••(•'••'•"•:•"•'•7

'• . , .

x•{

• ....

.................. •m:'•'

,_

• ............

• ...................

•'•:.,•

......

•....- ......... •" ß7•'

"•''•

.... ; '_ ],;,'• ; '-•-•,••-----m'•l•

r'u'" •$• .... ,......r ... •<• • ß • - ' •'•' /.•"• ,,,_ . ; .... •,.:•.-.............. • ...... •%--=: . .•..

=". .... ' . .'••'•q•

7' • '• • '. '

.----

"•. •q• ..... ........... •.... ,........... •• .....

• ';•u

,••.,,•-.•,,•?'..%•.,

•*•'•d•ti',i

l•l•l'i'/lt' i•"'; •:•;•t•l,•i•21.•.l•'*.'•.'l•,•-•ii•.'•'..•,•:b-''?

•?•'•'• ;••q• •.('•'•'"•",•h/•'•'} ••'•'•"' c'((•'•'".......... ß"

• .....

"•

'1 '•

'

' '

.... •, ß , -'

'

•

.....

--.- -••_ •

'

••"'•I'i•I• ,,.•_.'.';'".

........

:: ::

-

...... ß

•'.•-•-• '. ,,,'• ,,'•-<•., ,. ,,• ,

(a) OriœinalSeismic Data. Heavylines indicatemajorreflectors. N

N

N

'"

"

0.7 0.7 0.8

0.8 0.9

'I

10

i

!l

I

1.0 I 1 I

1.2

12

.I .!

1.3

!.3

i

!

1 4

1.4

1.5

1.5 1.6

1.7 1.7

I

I

18 i

I

(b)

19

Recursiveinversion of data in (a). ß

Figure Part 5 - Recursive Inversion - Practice

5.1

(Galbraith and Millington, 1979) Page 5 -

3


Brian Russell

Thelow frequencycomponent can be foundin oneof three ways' (1) From a filtered

sonic log

The sonic log is the best wayof derivinglow-frequency information in the vicinity of the well. However, it suffers fromtwomainproblems'it is usually stretched with respectto the seismicdata andit lacks.a lateral

component. Theseproblems, discussed in GalbraithandMillington(1979), are solved by using a stretching algorithm which stretches the sonic log information to fit

the seismic data at selected control points.

(2) From seismic velocity

analysis

In this case, interval velocities are derived from the stacking velocity functions along a seismicline usingDix' formula. The resulting function

will be quite noisyandit is advisable to do someformof two-dimensional filtering on them. In Figure 5.2(a), a 2-D polynomialfit has beendone to smoothout the function.

This final

set of

traces represents the filtered

interval velocity in the 0-10 Hz rangefor eachtrace and may be added directly to the invertedseismictraces. Refer to rindseth (1979), for more de ta i 1 s.

(3) Froma geological model Using all incorporated.

available sources, a

blocky geological model can be built

and

This is a time-consuming method.

Part 5 - Recursire Inversion - Practice

Page 5 -

4.

Introduction

to Seismic InversiOn Methods .

Brian Russell

.

70000

GOOO0

$0000

(pvl 4oooo

'/sgc

( b)

$oooo ZOOO0

,

I0000

/ --V..308 (PV)* 3460 ,

(a) i

Figure 5.2

.

ß

VELocrrY SURFACE 2rid ORDERPOLYN• Frr s •mTZ•eH CUTFtT• tRussell and Lindseth,

Part 5 - RecursiveInversion- Practice

Page5 -

5

1982).

Introduction


Brian Russell

Second, the low-frequency component can be addedto the high frequency

component by either adding reflectivity stage or the impedance stage. In section 2.3, it wasshownthat the continuousapproximationto the forwardand inverse equations was given by Inverse Equation

Forward Equati on

n Z(t) r(t) =•-1 d 1dt -

<::==> Z(t)

t

=Z(O) exp 2•0r(t)dt.

Sincethe previoustransforms are nonlinear(because of the logarithm), Galbraith and Millington (1979) suggestthat the addition of the low-frequency

component shouldbe madeat the reflectivity stage. In the SEISLOG technique they are addedat the velocity stage. However, dueto other considerations, this should not affect the result too much.

Of course, we are really interested in the seismic velocity rather than

the acoustic impedance. Figure5.2(b), fromLindseth(lg79), shows that an approximatelinear relationship exists between velocity and acoustic impedance, given by V = 0.308 Z + 3460 ft/sec.

Notice

that

this

relationship

is good for carbonates and clastics and

poor for evaporitesand shouldtherefore be usedwith caution. A moreexact relationship may be found by doing crossplots from a well close to the prospect. However, usinga similar relationshipwemayapproximately extract velocity information from the recoveredacoustic impedance.

Figure 5.3 showslow frequencyinformationderived from filtered

sonic

logs. The final pseudo-acoustic impedancelog is shownin Figure5.4 including the low-frequency component.Notice that the geologicalmarkersare moreclearly visible on the final inverted section.

Part 5 - Recurslye Inversion

- Practice

Page 5 -

6

Introduction


Figure 5.3

Methods

Brian Russell

LowFrequency comDonent derived from"st.reched:'sonic loœ.

0.7 0.8 0.9

0.9

l.O

1.0

1.1

1.1

1.2

1.2

1.3

1.:)

1.4

1.4

1.5

I$

1.6

1.6

1.7

1.7

1.8 19

19

Figure 5.4

Final inversioncombinin•Figures 5.1(b) and 5.3. Lines indicate major reflectors.

(Galbraith and Millington, 1979) Part

5 - Recursive

Inversion

- Practice

Page 5 -

7

Introduction


In sugary,

Methods

Brian Russell

the recursive methodof seismic inversion may be given by the

fol 1owing flowchart'

I

i

i

INTRODUCE LOW FREQUENCIES •) I•.v••o••DO-•CO••c • ,

ICORRECT TO PSEUDO ' VELOCITIES ß

,

CONVERT TO DEPTH I

Recursi ve Inversion ,

.

_

Procedure ß

.,

.

i

A commonmethod of display used for inverted sections is to convert to actual

interval

transit

times.

coloured according to a lithological

These transit

times are then contoured and

colour scheme. This is an effective way

of presentingthe information• especially to those not totally familiar'with normal

seismic

sections.

Part 5 - Recursive Inversion

- Practice

Page 5 -

8

Introduction

to Seismic

Inversion

(a)

Methods

Brian Russell

Frequency (e)

1

(b)

Fig.

(a) Frequency response of a theoretical

differentiator.

(b) Frequencyresponseof a theoretical integrator. (Russell and Lindseth, !982 ) ,m

Part

5 -Recursire

Inversion

- Practice

,i

m

ml

Page 5 -

,

,

9

Introduction

5.2

to Seismic

Inver.si.on

Methods

Brian

Russell

I nfor.marl o.nI ?•_Th.e. Lo..w .F.r.equ.e. ncycompo..ne. nt The key factor which sets inverted data apart from normal seismic data is

the inclusion of the low frequency component,regardless of howthis component

is introduced. In this section we will look at the interpretational advantages of introducing this component. The information in this section is taken from a paper by Russell and Lindseth (1982). We

start

by

assuming

reflectivity-impedance

the

extremely

simple

moael

for

the

relationship which was introduced in part 5.1. However,

we will neglect the logarithmic relationship of the more complete theory (this is justifiea for reflection coefficients less that 0.1), so t•at t

_1dZ(t)<=__==> Z(t) =2Z(O)j• 0r(t)at'

r(t) - • dt-

If we consider a single harmonic component, we may derive the frequency response of this tel ationship, which is

dejwt -dt "-- jwejwt<===> . jwtdt= -jw

eJWt

where

w--

21Tf,

In words.,differentiation introducesa -6 riB/octaveslope from.the high end of the spectrumto the low, and a +90 degree phase shift. Integration introduces a -6 dB/octave slope from the low end to the high end,

and a -90

degree phase shift. Simpler still, differentiation removeslow frequencies and integration puts them in. Figure 5.5 illustrates these relationships. But how aoes all this effect our geology? In illustrated three basic geological models' (1) Abrupt 1i thol ogi c change,

Figure 5,6 we have

ß

(2) Transitional

lithologic

change, an•

(3) Cyclical change.

Part

5 - Recursire

Inversion

- Practice

Page 5 -

10


(A)MAJORLITHOLOGIC CHANGE

Brian Russell

(B)TRANSITIONAL LITHOLOGIC CHANGE

Vl

I i

I. I

I I

V1

I i I

i

V:V•+KZ

i

(C)CYCLICALCHANGE

v•

Fig.

Part

5.6.

5 - Recursire

_

!

Threetypesof lithological models' (a) Major change, (b) Transitional, (c) Cyclical. (Russell andLindseth, 1982).

Inversion

- Practice

Page 5-

11


Brian Russell

Wemayillustrate the effect of inversiononthesethreecasesbylooking at both seismic anUsonic log Uata. To showthe loss of high frequencyon the

sonic log, a simple filter

is used,andthe associated phaseshift is not

introUuced.

To start with, considera major1ithologic boundary as exempl i lieu by the Paleozoicunconformity of Western Canada,a changefroma clastic sequence to

a carbonate sequence.Figure5.7 shows that mostof the information aboutthe

largestepin velocityis containeU in theD-10liz component of the soniclog. In Figure5.8, the seismicdataandfinal Uepthinversion are shown.On the seismic data, a major boundaryshowsup as simply a large reflection coefficient, whereas,on the inversion, the large velocity step is shown. RAWSONIC

FILTERED SONICLOGS

VELOCITY FT/SEC

0

10000

10-90HZ

O-IOHZ

O-CJOHZ

TIME

0.3-

0.5-

Fig. 5.7. Frequency components of a sonic log.

(Russell and Lindset•, 1982).

!

L

, , ,

I

I ß

[ I

L

Part 5 - Recursire Inversion - Practice

Page 5 -

12

Introduction

to Seismic

Inversion

Methods

Brian

Russell

o'- .

ß

(a)

.%;

DEPTH

SEISLOG DEPTH ß

o

..... ßlOP OF "' . ß""I:'ALEOZOIC -425'

(b)

Fig.

..... Part

5.8.

_.........

5 - Recursive

Major litholgical'change, Saskatchewan example. (a) Sesimic s_ection, (b) Inverted section.

_(R_q•sell ....and Li,pqse_th,_•!98_2)___ Inversion

- Practice

Page 5 -

13


To illustrate

transitional

and cyclic

Brian Russell

change, a single examplewill be

used. Figure. 5.9 showsa soniclog from an offshore Tertiary basin, illustrating the rampswhichshow a transitional velocity increase, and the rapidly varyingcyclic sequences. Noticethat the 0-10 Hzcomponent contains all the information about the ramps, but the cyclic sequenceis containedin

the 10-50 Hz component.Onlythe Oc component is lost from the cyclic component uponremovalof the low frequencies. Figure 5.10 illustrates the samepoint usingthe original seismicdata andthe final depthinversion. In summary, the informationcontainedin the low frequencycomponentof the soniclog is .lost in the seismic data. This includessuchgeological information as the dc velocity component,large jumps in velocity, and linear

velocity ramps. If this informationcould be recoveredandincluUeaduring the inversion process, it wouldintroducethis lost geological information.

Fig. 5.9. Soniclog showing cyclic andtransitionalstrata. (Russell

Part

5 - Recurslye

Inversion

- Practice

and LinOseth, 1982) Page 5 -

14

Introduction

to Seismic

Inversion

Methods

Brian

Russell

SEISMIC SECTION-CYCUC & TRANSITIONAL STRATA

(a)

(b) i

1-3500

ß

Part

5 - Recursive

Inversion

- Practice

Page 5 -

15

Introduction

5.3


Brian Russell

Sei smical ly Derived Porosi ty --

ILI

,

ß

We have shownthat

I

seismic data may be quite

adequately inverted to

pseudo-velocity (and hence pseudo-sonic)information i f our corrections and assumptions are reasonable.

true

sonic log

Thus,

information

we may try to treat

and extract

the inverted

petrophysical

data

data as

from

it,

specifically porosity values. Angeleri and Carpi (1982) have tried just this, with mixed results. The flow chart for their procedure is shown in Figure 5.11. In their chart, the Wyllie formula and shale correction are given by:

At

where

--transit

time for fluid saturated rock,

Zstf= porefluid transittime, btma:rockmatrixtransittime, Vsh = fractionalvolume of shale,and

btsh: shaletransittime. The derivation

of

porosity

control.

Figure 5.12 shows the

porosity

for

was tried on a line which had good well plot

each of three wells. from

the

seismic

well

log

Notice that the fit

clean sands and very poor in the dirty information

of

sands.

section

only

porosity versus seismic is reasonable in

the

Thus, we mayextract porosity under

the

most

favourable

conditions, notably excellent well control and clean sand content.

Part 5 - Recurslye Inversion - Practice

Page 5 -

16


Brian Russell

!•ILI61C .AT&'$[IS'MI• .AT&'

F ']w[tt 'ill ]

I-"'• ''' m.,,•, _,ml.111,l lit

-[ ,gnu mill i' •ill. Utl..I •%lOtOG

I IIITEIPllETATII i

l!

WlltK

:

t

Fig.

' .

5.11. Porosity eval uati on flow diagram.

(AngeleriandCarpi, 1982). ,

WELL

__

ClII

ß"

2

PNIIVI

o..- OPt

, e

e

WELL

, ß I

e

.

e

e

.

.

ß

e

.....

' e

e

WELL

poeoItrv

ß e

,

3

e

I

i

I e

e

,e

ß

i

CPI

-e

i

ß

ß

ß

e

1.4

1.7

1.8,

1.9

.

Fig.

5.12. Porosity profiles from seismic data and borehole data.

Shalepercentage is al so displayed. (Angel eri andCarpi, 1982). i

Part

5 - Recursire

Inversion

- Practice

Page 5 -

,

17

Introduction

to Sei stoic Inversion

Methods

Brian

Russel 1

PART 6 - SPARSE-SPIKE INVERSION • {

Part 6 - Sparse-spike Inversion

•

......

•

I

]

m

•

m

6-

1

Introduction to Seismic Inversion Me.thods

6.1

Brian Russell

Introduction

Thebasic theoryof maximum-1 ikeli hooddeconvol ution (MLD)wasdeveloped by Dr. Jerry Mendeland his associatesat USCanUhasbeenwell publicised (KormyloandMendel,1983;Chiet el, 1984). A paperby Hampson and Russell (1985) outlined a modification of maximum-likelihood Ueconvolutionmelthod ,

which allowed the methodto be moreeasily applied to real seismic •ata. One of the conclusionsof that paper wasthat the methodcould be extenoedto use

the sparsereflectivity as the first step of a broadbandseismic inversion technique.This technique,which will be termedmaximum-likelihood seismic inversion, is discussed later

You will

in these notes.

recall that our basic model of the seismic trace is

s(t) = w(t) * r(t) + n(t), where

s(t)

: the seismic trace,

w(t) : a seismic wayelet,

Notice that the solution to the

there

r(t)

: earth reflectivity,

n(t)

= addi tire

above equation

are three unknownsto solve for.

and

noise.

is indeterminate,

However, using certain

since

assumptions,

the aleconvolution problem can be solved. As we have seen, the recursire method of seismic inversion is basedon classical aleconvolution techniques, which assumea randomreflectivity

and a minimumor zero-phase wavelet.

They

produce a higherfrequency waveleton output,but neverrecoverthe reflection coefficient series completely. More recent aleconvolution techniquesmaybe

groupedunder the categoryof sparse-spike meth•s. That is, they assumea certain modelof the reflectivity

and make a wavelet estimate based on this

assumption.

Part 6 - Sparse-spike

Inversion

6-

2

Introduction


Methods

Brian Russell

ACTUAL REFLECTIVITY

POISSON-GAUSSIAN I,:,

I

..

SERIES OF LARGE EVENTS

--F

GAUSSIAN

BACKGROUND

OF SMALL

SONIC-LOG

EVENTS

REFLECTIVITY

EXAMPLE

Figure 6.1 The fundamental assumption of the maximum-likelihood method.

Part 6-

Sparse-spike Inversion

6-

3

Intr6duction to Seismic Tnvetsion Methods

Brian Russell

These techniques include(1) btaximum-Likel ihood deconvolutton

and inversion.

(2) L1 norm deconvolution and inversion.

(3) Minimumentropy deconvol ution (MEO). From the point of view of seismic inversion, sparse-spike methods have an

advantage over classical methodsof deconvolutionbecause the sparse-spike estimate, with extra constraints, the

reflectivity.

We

deconvolution, and will

will

Maximum-Likelihood i

i

m

focus

ß

m

m

m

ß

ß

ß

m

_

_

l!

Figure 6.1 illustrates

.

.

•

on

maximum-likelihood

the L1 normmethod of Dr.

Doug

and Inversion

m

I

Maximum-Li kel i hood Deconvoluti I

bandwidth estimate of

not be discussed in these notes.

Deconvolution

!

initially

then move on to

O1denburg. The MEDmethod will 6.2

can be used as a full

am

..

the

_

ß

on I

_

fundamental assumption of Maximum-Likelihood

deconvolution, which is that the earth' s reflectivity

is composed of a series

of large events superimposedon a Gaussian backgroundof smaller events.

This

contrasts with spiking decon, which assumesa perfectly randomdistribution reflection

Figure

coefficients.

The real

sonic-log

reflectivity

at the bottom of

6.1 shows that in fact this type of model is not at all

Geologically,

the

large

events

correspond to

of

unreasonable.

unconformities

and major

ß

1i thol ogic boundaries.

From our assumptions about the model, we can derive an objective

function

whichmaybe minimizedto yield the "optimum" or mostlikely reflectivity. and wavelet

combination consistent with the statistical

assumption. Notice

this method gives us estimates of both the sparse reflectivity ,,


that

and wavelet.

m


Brian Russell

INPUT

WAVELET

SPIKE SIZE'

SPl• ••:

9.19

50.00

REFLECTIVITY

NOISE' 39.00 NOISE

OB,.ECTIVE'

98.19

Figure6.2(a) Objectivefunctionfor onePoSsible solutionto inputtrace.

INPUT

WAVELET

SPIKE S!7_F: 6.38

SPIKE DENSIq'•, 70.85

REFLECTIVITY

NOISE: 81.• 5

NOISE

OBJECTIVE:158.98

Objective function fora second possible solution toinput

Figure 6.2(b)

trace. Thisvalueis higher than6.2(a),.indicating a less 1ikely solution.

!

,

,,


6-

5

Introduction


Brian Russell

The objective function j is given by

ß

k=l

- R2

N2

k=l

r(k)

where

- 2mln(X)-

= reflection

2(L-re)In(i-A)

coeff.

at kth

sample,

m = numberof refl ecti OhS, ß

L : total numberof samples,

N : sqare root of noise variance, n : noise at kth sample, and

• = likelihood that a given sample has a reflection.

Mathematically,

the

expected behavior

of

the

objective function is

expressed in terms of the parametersshownabove. No assumptions are made

aboutthe wavelet. The reflectivity sequenceis postulatedto be "sparse", meaning that the expected numberof spi•es is governed by the parameter lambda, the ratio of the expected numberof nonzer. o spikes to the total number of trace samples. Normally, lambda is a numbermuchsmaller than one. The

other parametersneededto describe the expected behavior are R, the RMS•size

of the large spi•es, andN, the RMSsize of t•e noise. With these parameters specified,

any glven deconvolution sol ution

can be examined to see.whether it

is likely to be the result of a statistical processwith thoseparameters.For example,if the reflectivity estimate has a number of spikesmuchlarger than the expectednumber,then it is an unlikely result. In

simpler terms, we are looking for the solution

with

the minimum

number of spikesin its reflectivity and t•e lowestnoisecomponent.Figures 6.2(a) and 6.2(b) showtwopossiblesolutionsfor the sameinput synthetic

trace. Noticethat the obje6tive functionfor the onewith the minimum spike structure

is indeed the lowest


value.

6-

6

Introduction to Sei smic I nversi.on Methods

Bri an Russel1 Synthetic

Reflectivity

Original

I,

Model

I

I terati

,1.2.

ill. -.I

,i.

on I

I terati

on 2

Iteration

3

I teration

4

Iteration

S

Iteration

6

Iterati

on 7

Figure 6.3.

The Sinl•le Most Likely Addition (SMLA)algorithm illustrated for a simple reflectivity

Part 6 - Sparse-spi ke Inversion

model.

6-

7

Introduction


Brian Russel1

Of course, there maybe an infinite numberof possible solutions, and it

would take too muchcomputertime to look at eachone. mTherefore, a simpler method is used to arrive at the answer. Essentially,

we start with an initial

wavelet estimate,es'timatethe sparsereflectivity, ' improve the wavelet and iterate through this sequence of steps until an acceptably low objective function is reached. This is shownin block form in Figure 6.4. Thus, there

is a two step procedure-havingthe waveletestimate, updatethe reflectivity, and then,

In

having the

reflectivity

estimate,

update the

wavelet.

These procedures are illustrated on model data in Figures 6.3 an• 6.5. Figure 6.3, the proceUure for upUating the reflectivity is shown. It

consists of adding reflection

coefficients oneby one until an optimumset of

"sparse" coefficients hasbeenfound. Thealgorithmusedfor updatingthe reflectivity

is callee the single-most-likely-addition

algorithm (SMLA)since

after each step it tries to find the optimumspike to add. Figure 6.5 shows the procedure for updating the wavelet phase. The input model is shownat the

top of the figure, and the up•ated reflectivity and phaseis shownafter one, two, five, and ten iterations. Notice that the final result compares favourably

with the model wavelet.

WAVELET

ESTIMATE

IMPROVE ES•TE

WAVELET

ESTIMATE

REFLECTIVITY

Fiõure 6.4.

Theblockcomponent method of solvingfor both reflectivity andwavelet. Iterate around the loop unti 1 converRence.


6-

8

Introduction

to Seismic Invers.ion Methods

Wayel et

Brian Russell

ReflectiVity '

Synthetic

INPUT MOD INITIAL

Ill

Fi õure 6.5.

,I ,

CUESS

TEN ITERATIONS

The procedure for updatinõ the wavelet in the maximum-likelihood

Between each iteration

iter. ation on reflectivity

method.

above, a separate

(see Fiõure 6.3)

has been done.


6-

9

Introductionto SeismicInversionMethods

Figure 6.6

is

an exampleof the

Brian Russell

algorithm

applied

to a synthetic

seismogram. Notice that the major reflectors have been recovered fairly

and that the resultant trace matchesthe original trace quite accurately. course, the reflection

smaller

reflection

coefficient

coefficients

are

missing

well

Of

in the recovered

series.

Let us nowlook Cretaceous gas play

at in

comparison between the

some real data. Southern Alberta. input

anU

The first example is a' basal Figure 6.7(a) and (b) shows the

output

stack

from the

aleconvolution

procedure. Also shownare the extracted and final wavelet shapes. The main things to note are the major increase in detail (frequency content) seen in the final stack, and the improvementin stratigraphic content.

Figure 6.8 is a comparisonof input and output stacks for a typical Western

Canada basin seismic line.

The area is an event of interest

between

0.7 anU 0.8 seconds, representing a channel scour within the lower Cretaceous.

Although the scour is visible on both sections, a dramatic improvementis seen in the resolution

of the infill

Within the central

of

this

channel

on the deconvolved section.

portion of the channel, a .positive reflection

lateral extent of five traces is clearly

with a

visible and is superimposed on the

Uominant negative trough.

INPUT:

V.

,.: --

ESTIMATED:

ttl J':ll'j' ß

"'" Figure 6.6

" Synthetic


seismogram test.

6-

10

Introduction


Methods

Brian Russell

'SONIC

SYNTHETIC 0.5

LOG iZ.i

0.6

0.7

0.8

EXTRACTED

(a) Initial

WAVELET

seismicwith extracted wavelet.

0.5

0.6

.

0.8

(b)

Final deconvolved seismic with zero-please wavelet. ....

Part 6-

Figure


-

-_

6.7 __

._

11

Introduction

This

to Seismic

is

quite

Inversion

possibly

Methods

Brian

a clean channel sand and may or

Russell

may not be

prospective. However, this feature is entirely absent on the input stack. Overlying the channel is a linear anomalywhich could represent the 'base of a gas sand, and is muchmore sharply lateral

and vertical

Finally

reflectivity. reflections

This

is

the

shown in

deconvolved output

Figure 6.9.

both in a

the base of

the

are

present.

and estimated

the

Although some of the subtle

missing from this estimated reflectivity,

that all the main reflectors

clearly

on the output section,

sense.

we have taken

are

defined

It

is

there is no doubt

interesting

channel (at 0.7;- seconds)and the

to note how

base of

the

6-

12

postulated gas sand on top of the channel have been delineated.



Brian Russell

INPUT

DECONVOLVED

STACK

STACK 0.6

0.7

0.8

0.9 Figure 6.8

An input stack over a channelscour and the resul ting deconvol ved sei smic.

ESTIMATED

DECONVOLVED

REFLECTIVITY

STACK

0.6

0.7

0.8

0.9 Figure 6.9

The deconvolved result from Figure 6.8 and its


estimated reflectivity.

m

13

Introduction

to Seismic

Maximum-Likel

An

ihood

obvious

reflectivity

Inversion

Methods

Brian

Russell

Inversion

extension

of

the

theory

to

invert

the

"blocky" impedancefrom

to Uevise a broad-band or

data (Hampson and Russell, 1985). impedanceZ(i) maybe written

is

Given the reflectivity,

r(i),

es ti mated

the

seismic

the resul ting

Z(i)=Z(i_l )[1+r(i)] 1 - r(i)

Unfortunately, from MLD produces

additive

noise.

'

of thi s formula to the reflectivity

application

res ults,

unsatisfactory

especially

this

estimate

spectrum.

the

presence of

Although the MLD algor it•m'extrapol ares outsi de the bandwidth

of the wavelet to produce a broad-band reflectivity

of

in

estimates

The result

degraaed by noi se at the

is is

that

the

while

estimate, the reliability low frequency

short

wavelength

end of

features

the

of the

impedancemay be properly reconstructed, the overall trenu is poorly resolvea.

This is equivalent to saying that the times of the spires on the reflectivity estimate are better

of

resolved than their

amplituaes.

In order to stabilize the reflectivity estimate, independent knowleUge the impedancetrenU may be input as a constraint. Since r(i) < l, we can

derive

a convolutional

reflectivity,

type

equation

between

acoustic

impeUance anU

written

In Z(i)

= 2H(i) * r(i)

+ n(i),

where Z(i) = the knownimpedancetrend,

• i <0 H(i)

:

• i >0 and

n(i) : "errors" in the input trend.

_

Part 6 - Sparse-spike

Inversion

6•

14

Introduction

to Seismic

Inversion

Figure 6.10

Figure 6.11

Methods

Brian

Russell

Input Model parameters.

Maximu•m-L i kel i hood i nversi on result from Figure 6.10.

ß


.m

__

6-

lb

Introduction

The

to Seismic

Inversion

error series n(i)

'Methods

reflects

the fact

Brian

that the trend

Russell

information

is

approximate. We now have two measured time-series:

the seismic trace, T(i),

and the

own wavelet

log

of

impedanceIn Z(i),

each with

its

and noise

parameters. The objective function is modified to contain two terms weighted by their relative noise variances. Minimizing this function gives a solution

for r(i) whichattemptsa compromiseby simultaneouslymoUellingthe seismic trace while

conformingto the knownimpedancetrend.

If

both the

seismic

noise and the impedance trend noise are modelledas Gaussiansequences,their respective variances become"tuning" parameterswhich the user can modify to shift the point at which the compromise occurs. That is, at one extreme only

the seismicinformationis usedand at the ot•er extremeonly the impedance trend.

In our first example,the methodis tested on a simple synthetic. Figure 6.10 showsthe sonic log, the derived reflectivity, the zero-phase wavelet used to generate

the

synthetic,

and finally

the

synthetic

itself.

This

example was usedinitially becauseit truly represents a "blocky" impedance (and therefor.e a "sparse" reflectivity) and therefore satisfies the basic assumptions of the method.

In

Figure

6.11 the maximum-likelihood inversion result

is

shown.

In

this casewe haveuseda smoothedversion of the sonic velocities to provide the constraint.

A visual

comparisonwoulU indicate

that

the extracteU

velocity profile corresponds very well to the input. A moredetailed comparisonof the two figures showsthat the original and extracted logs do not matchperfectly. T•ese small. shifts are due to slight amplitude problems on the extracted reflectivity.

It is doubtful that a perfect match could ever

be obtai neU.


6-

16

Introduction

to Seismic

Figure 6.12

Inversion

Creation of a seismic model from a sonic-log.

Figure 6.13 •-

_

!

......

ii__

Bri an Russel 1

Methods

Inversion result from Figure 6.12. -

-

Part 6 - Sparse-spi•e Inversion

i

-

•_!

mm

i

i

i

ß

i

i

!

It_l

I

17


Let us nowturn our attention to

Brian Russell

a slightly

more realistic

synthetic

example. Figure 6.12 showsthe applicationof this algorithmto a sonic-log derivedsynthetic. At the' top of the figure we see a sonic log with 'its reflectivity sequencebelow. (In this example,we have assumedthat the density is constant, but this is not a necessary restriction.) The

reflectivity wascbnvolved with a zero-phase wavelet,bandlimitedfrom10 to 60 Hz, andthe final syntheticis shown at the bottomof the figure. The results

Figure 6.13.

of

The initial

the middle panel,

the maximum-likelihood inversion method are

log is shownat

sbown in

the top, the constraint is shownin

and the extracted resull• is shownat the

bottom of

the

diagram. In this calculation, the waveletwasassumed known. Notethe blocky nature of the estimatedvelocity profile compared with the actual sonic log profile. Again, the input and output logs do not matchperfectly. The fact that the two do not perfectly matchis due to slight errors in the reflectivity sizes whichare amplified by the integration process,and is partially the effect of the constaintused. Theconstraintshownin Figure

6.13 wascalculatedby applying a 200 ms smootherto the actual log. In practice, nearby well

this information could be derived from stacking velocities or from control.


6-

18

Introduction

to Seismic

Inversion

Methods

Brian

Russell

* !

Figure 6.14

An input seismic 1ine to be inverted.

:

ß

'. eel'?

e4dl

Figure 6.15

Maximum-Liklihood reflectivity

estimate

from

seismic in Figure 6.14.


6-

19

Introauction


Brian Russell

Finally, we showthe results of the algorithm applied to real seismic data. Figure 6.14 showsa portion of t•e input stack. Figure 6.15 showsthe •D extracted reflectivity. Figure 6.16 shows the recoveredacoustic impedance,wherea linear ramphas beenusedas the constraint. Notice that the inverted section •isplays a "blocky" character, indicating that the major features of the impedance log havebeen successfullyrecovered. This blocky

impedance canbe contrastedwith the more traditional narrow-band .inversion procedures, whichestimate a "smoothed" or frequencylimited version of the impedance.Finally, Figure 6.17 shows a comparison betweenthe well itself and the

inverted

section.

In summary, maximum-likelihoodinversion is a procedurewhich extracts a broad-band estimate of the seismic reflectivity

and, by the introduction of

1inear constraints, al lows us to invert to an acoustic impedancesection which retains the major geological features of boreholelog data.


6-

20

Introduction to Seismic Inver.sion Methods

Figure 6.16

Brian Russell

Inversion of reflectivity shownin Figure 6.15.

SEISMIC INVERSION WELL

+

Figure 6.17

SONIC LOG

A comparison of the inverted seismic data and the sonic log at well location.


..

21


6.3

Brian Russell

The L 1 Norm Method --

__LI

_

_

_

i

.

Another method of- recursive,

single trace inversion which uses a

"sparse-spike" assumption is the L1 normmethod, developed primarilyby Dr. DougOldenburgof UBC. and Inverse Theory andApplications (ITA). This method is also often referred to as the linear programming method,and this can lead

to confusion. Actually, the two namesrefer to separateaspectsof the method. Themathematical modelusedin the construction of the algorithm is the minimizationof the L1 norm. However,the methodusedto solve the problem is linear programming. The basic theory of this methodis found in a paper by Oldenburg, et el (1983). The first part of the paper discussesthe noi se-free convolutional model,

x(t) --w(t) * r(t),

where x(t) = the seismictrace, w(t) --the wavelet, an•

r(t) -- the reflectivity. The authors point out that if

a

high-resolution

aleconvolution

is

performedon the seismictrace, the resulting estimateof the reflectivity can be thought of as an averagedversion of the original reflectivity,

as shownat

the top of Figure6.18. This averaged reflectivity is missingbotht•e high andlow frequency range,andis accurateonly in a band-limiteacentral range of frequencies. Although there are an infinite numberof ways in which the missing frequencycomponentscan be supplied, Oldenburg, et al (1983) show that

we can reduce this nonuniqueness by supplying more information

to

the

problem, such as the layered geological model

•= 0 if t •l• r(t)--•,rj6(t-l•), where j--!

Part 6-

Sparse-spike

Inversion

=1 ift:•

, an• .

6-

22

Introduction


Brian Russell

b

ß I

ß m

ß m

•

e

f

o

0.0

T.IJdE•(•J

1m

.50

joo

j25

FR F.,O [HZJ

I I

!

Figure 6.18

Part 6-

Synthetic test of L1 NormInversion, moUified fro•.q Oldenburget al (1983). (a) Input impedance, (b) Input reflectivity, (c) Spectrumof (b), (d) Low frequencymodeltrace, (e) Deconvolutionof (•), (f) Spectrumof (U), (g) Estimatedimpedancefrom L1 Norm method, (•) Estimated reflectivity, (i) Spectrumof (•).


6-

23

Introduction

to Seismic

Inversion

•.le.thods

Brian

Russell

Mathematically, the previous equation is considered as the constraint to

the inversion problem. Now, the layered earth model equates to a "blocky" impedancefunction, which in turn equates to a "sparse-spiKe" reflectivity function.

The

above constraint

will

thus

restrict

"sparse" structure so that

extremely fine

reflection

not be fully

coefficients,

will

our inverted

structure,

result

to

a

such as very small

inverted.

The other key difference in the linear programmingmethod is that the L1 norm is minimized

rather

than the L2 norm.

The L1 norm is defined

of the absolute values of the seismic trace.

is

as the

sum

TrueL2 norm, on the other hand,

defined as the square root of the sumof t•e squares of the seismic

trace

values. The two norms are shownbelow, applied to the trace x:

x1 :

xi

and

x2:

xi

i--1

i:1

The fact that the L1 norm favours a "sparse" structure is shown in

following simple example. (Taken from the notes to Dr. convention course'

"Inverse

seismograminversion").

theory with application

the

Oldenburg's1085 CSEG to

aleconvolution and

Let f and g be two portions of seismic traces, where'

f: (1,-1,0) and g : (0,%•,0) . The L2 norms are

therefore'

The L1 normsare given by'

fl -- 1 + 1 : 2

and gl =

Notice that the L1 norm of wavelet

g

is

' smaller than the L1 norm of f,

whereas the L2 norms are both the same. Hence, minimizing the L1 norm would

reveal that g is a "preferred" seismic trace based on it's


sparseness.

6-

24

Introduction

to Seismic

Inversion

Methods

Brian

Input sei smic data

(a)

(b)

Estimated refl ec ti vi ty

(c)

Figure 6.19

Russell

Final

impedance

L1 14ormmetboOapplied to real seismic data,

(Walker and Ulrych, 1983) Part 6 - Sparse-spike Inversion

6-

25

Introduction to Seismic Inversion MethoUs

Brian Russell

Several other authors had previously consideredthe L1 normsolution in

deconvolution (Claerboutand Muir, 1973, andTayloretal., 1979), however, they considered the problem in the timedomain.Oldenburg et al.w suggested solving the problemusing frequencydomainconstraints. That is, the reliable

frequencyband is honored whileat the same timea sparsereflectivity is created.

The results of their. algorithm on synthetic data are shownat the

bottom of Figure6.18. Theactual implementation of the L1 algorithmto real seismicdata has been done by Inverse Theory andApplications(ITA). The processingflow •or the linear programming inversion methodis shownbelow.

InterPreter'= CMPStacl=

r(t)©w(t)

t

ß

,i

i

I,,ico,ect,', ,o,' Residu Pm'm,e o,w(t) I ß

i

i

i i

i

I

i

i

I Fourier Trans•• of<•r (t)> I i

Scale Data

Const. mints. From $tackins•_V'elocitles I Con,straints From 'Well Logs I ii &

i

Unear Programing Invemion

Assume r(t)ß• n;)(t-•q),isaspame, reflection series. Minimize the sum of absolute reflection strengU•. FulFBand Reflectivity Series r (t) Signal to Noise Enhancement and Display Preparation Integration to Obtain Impedance Sections

Figure 6.19(b) The L1 Norm(Linear Programming.) Method. (Oldenburg,1985). Part 6-


6-

26

Introduction

TSN

tO0

to Seismic Inver. s,ion Methods

90

80

70

60

Brian Russell

50

40

30

20

tO

1,2

1,3

1,4

1,5

1,6

1,7

1,:8

.2,0

2ø2

Figure 6.20

Inputseismicdatasectionto L1Norm inversion. (O1denburg, 1985'


6-

27

Introduction


Figure 6.19 showsthe

Methods

application

Brian Russell

of

the above technique to an actual

seismic line from Alberta. The data consist of 49 traces with a sample rate of 4 msecand a 10-50 Hz bandwidth. The figure showsthe linear programming reflectivity

and impedanceestimates below the input

should be pointed out that a three the

final

results

in

both

seismic section.

It

trace spatial smootherhas been applied to

cases.

Finally, let us consider a dataset fromAlberta which has been processeU through the LP inversion method. The input seismic is shownin Figure 6.2D and the final inversion in Figure 6.21. The constraints useU here were from well

log

data.

In the final

superimposed on the final

inversion notice that the impedance has been

reflectivity


estimate using a grey level

scale.

6-

28


Brian Russell

1.6

1.7

1.8

1.9

2.0

2.1

2.2

Figure 6.21

Part 6 -Sparse-spike

Reflectivity and grey-level plot of impedancefor the L1 Norminversion of data in Figure 6.20.

Inversion

(O1denburg, 1985

6-

2-9

Introduction

6.4


Reef P roblee ß

Br•an Russell

_

Onthe next fewpages 'is a comparisonbetweena recursiveinversion procedure (Verilog) anda sparse-spikeinversionmethod(MLD). The sequence !

of pages includes the following:

- a sonic log and its derived reflecti vt ty, - a synthetic seismogramat both polarities, - the original seismic line, showingthe well location, - the Verilog inversion, - the MLD inversi

and

on.

BaseUon the these data handouts, do the following interpretation exerc i se:

([) Tie the synthetic to the seismicline at SP76. (Hint- use reverse polari ty syntheti c).

(g) Identify and color the following events in the reef zone- the Calmar shale (which overlies the Nisku shaly carbonate), - the 1retort shale, and

- .the porous Leduc reef.

(3) Comparethe reefal events on the seismic and the two inversions. Use a blocked off version of the sonic log.

(4) Determine for parallelism which section tells you the most about the reef

Part 6-

zone?


6-

30


Rickel,

g

Phas•

3g

26

Hz

Ns,

REFL. COEF.

Brian Russell

DEPTH lib

VELOCI •¾ Eft,/sec.

...,--

...,--

...m

$11qPLE I HTI3tViILAliPLI •IIi)E

Figure 6.22

Part 6-

I

Sonic Log and synthetic


2 Ns.

tiC. Ilql•. - Sonic

Pei.•ri es onlg

at the reef well.

6-

31

Introduction

to Seismic

Inversion

Methods

Brian

Russell

.47

'49

5! 55 57

?! '?•

;'5 ??

B5

$5

99

•41

ß

i

Part 6-

Sparse-spike

Inversion

32

Introduction

to Seismic

Inversion

Methods

Brian Russell

***__********************************************************

0.7

Part 6-

0.8

.la.El


1 .:0

1.1

1..2

1. E)

6-

33

Introduc%ion

0.?

[o Seismic

IB.G

Inversion

O.G

Meltotis

1.0

Brian

1.1

1.2

ß

i i

i i

I •

1.•

eet I i e

I i

eI

Russell

1.4 ei ßl e

ß

e

I

1

25 27 31

•x

$1

$?

67 ?.1 ?$ ß

i i

'111

113

119

123 127 I

131

I .

I

137 ß

IIIII

i•i

iii ß

.

Part 6 - Sparse-spike Inversio•

ß

i iiiiiiiiii

I1•11•1111 ii

lee-

,

IIII

ß

.

..,

34


Brian Russell

PART 7 - INVERSION APPLIED TO THIN BEDS

Part 7 - Inversion applied to Thin Beds

Page 7-

I

Intro4uction

7.1


Methods

Brian

Russell

Thin Bed Analysis

One of the problemsthat we have identified in the inversion of seismic traces is the loss of resolution caused by the convolution of the seismic wavelet with the earth's reflectivity. As the time separation between reflection

coefficients

becomessmaller, the interference between overlapping

wavelets becomesmore severe.

effect

Indeed,

of reflection coefficients

in Figure 6.19 it was shownthat the

one sampleapart and of opposite sign is to

simply apply a phaseshift of 90 degrees to the wavelet. In fact, the effect is

more of

a differentiation

of the wavelet,

which alters

the

amplitude

spectrumas wel1 as the phase spectrum. In this section we will look closer at the effect thin

of wavelets on thin

beds and how .effectively

we can invert

these

bed s.

The first (1973).

comprehensivel'ook at thin bed effects was done by Widess

In this paper he used a model which has becomethe standard for

discussing thin beds, the wedgemodel. That is, consider a high velocity laye6 encasedin a low velocity layer (or vice versa) and allow the thickness of the layer to pinch out to zero. Next create the reflectivity responsefrom the impedance,and convolvewith a wavelet. The thickness of the layer is given in terms of two-waytime through the layer and is then related to the dominantperiod of the wavelet. The usual wavelet used is a Ricker becauseof the simpl i city of its shape.

Figure 7.1 is taken from Widess' paper and showsthe synthetic section as the thickness of the layer

decreases

from twice

wavelet to 1/ZOth of the dominant period.

the dominant period of the

(Note that what is refertea to as a

wavelengthin his plot i s actually twice the dominantperiod). A few important points can be noted from Figure 7.1.

First,

the wavelets start

interfering

with eackotherat a thicknessjust below two dominant periods,but remain Clistinguishable down to about one period.

Part 7 - Inversion

applied

to Thin Beds

Page 7-

g

Introduction

to Seismic

Inversion

Methods

Brian

Figure 7.1

Effect

Russell

of bed thickness

on

reflection waveshape,where (a) Thin-bed model, (b) Wavelet shapesat top

PI•OPAGA! ION I NdC ACnOSS TK arO) .

and bottom re fl ectors,

•'-----•).z _1

(c) Synthetic seismic

I

model, anU (d) Tuning parameters

as measured

from

resul ting waveshape.

--t

(C)

5O

(D)

,

2.0

,.

THINBEDREGIME J PEAK-TO-TROUGH/ AMPLITUDE

0.8

1.0

< 0.4

/ -0.4

,• i

-40

\ .

0

.

.

. .

20

40

MS

TWO-WAY TRUE THICKNESS (MILLISECONDS)

Figure 7.2

A typical detection and resolution cha•t used to interpret bed thickness from zero phase seismic data.

('Hardage, 1986) .

..

_

i

i

,,

,

i

_

-

-


-

-_-

-

_

-

_ .....

Page 7-

l.

3

_

Introduction

to Seismic

Inversion

Methods

Brian

Russell

Below a thicknessvalue of oneperiodthe waveletsStart merginginto a single wavelet, increase

and an amplitude increase

is

observe•.

This

amplitude

is a maximumat 1/4 period, and decreases from this point down... The

amplitude is appraochingzero at

1/•0

period, but note that the resulting

waveform is a gO degree phase shifted version of the original wavelet.

A more quantitative way to measurethis information is to plot the peak to trough amplitude difference and i sochron across the thin bed. This is done in Figure 7.•, taken from Hardage (1986). This diagram quantifies what has already been

seen qualitatively

the

seimsic

section.

That

is

that the

amplitude is a maximumat a thickness of 1/4 the wavelet dominant period,

also that this is the lower isochron limit. considered

to

obtain fully

be

the thin bed threshhold,

resolved reflection

and

Thus, 1/4 the dominant period is below which it

is

difficult

to

coefficients.

7.2 In.versionCamparison of T.hinBees ß

To test out this theory, a

using

both

recursire

thin

bed model was set up and was inverted

inversion and maximum-likelihood

impedance model is shown in Figure

7.3,

aleconvolution.

and displays a velocity

The

decrease in

the thin bed rather than an increase. This simply inverts the polarity of Widess' diagram. Notice that the wedge starts at trace 1 with a time thickness

sample.

of

100 msec and thins down to a thickness of 2 msec,.or

The resulting

synthetic

seismogram is shownin Figure 7.4.

.one

time

A 20 Hz

'Ricker wavelet was usedto create the synthetic. Since the dominant period

(T) of a 20 Hz Ricker is 50 msec,the wedge hasa thicknessof 2T at trace 1, T at trace 25, T/2 at trace 37, etc.

Parl•'7 - 'inverslYn'ap'plled 1•oThin'-Beds .....

Page 7 --'4

'•-

Introduction

to Seismic

4

8

12

Inversion

16

Methods

20

24

Brian

28

32

36

40

44

Russell

48

lOO

200

3OO

400

ß

500

Figure7.3 True impedance fromwedge model.

o

lOO

200 .

300 ß

400

500

Figure 7.4

Wedgemodel reflectivity 20 HZ Ricker

convol ved with

wavelet.

Part 7 - Inversion applied to Thin BeUs

Page 7-

5

Introduction

to Seismic

First,

the

Inversion

Methods

let us consider the effect

Brian

of

performing a recursire inversion on

wedgemodel. The inversion result is shownin Figure 7.5.

low frequency component was not added into better

show the effects

was also felt

little

of the initial

that the

information

with recursire

addition

inversion..

Note that the

the solution of Figure 7.5, to

recursire phase of the inversion.

of

to this test.

the

low frequency

the

thickness of

It

componentwould ado

Notice that there'are two major

First,

Russell

problems

the beU has only been

resolved downto about 25 msec, which is 1/2 of the dominantperiod. Remember, that this is a two-way time, has

been resolved

therefore

we say that the bed thickness itself

down to 12.5 msec,

or

1/4

period.

This

theoretical

resolution limit is the sameas that of Widess. Also, the top of the weUge appears "pulled-up" at the right side of the plot as the inversion has trouble with the interfering

wavelets.

that there are actually recursire

inversion

A second problem is that, although we know

only three distinct

has estimate•

at

velocity

least

units in the section, the

seven

in the vertical

=irection.

ß

This

result

is Uue to the banu-limited

Uescriptively,

nature

of the Ricker

wavelet.

More

every wiggle on the section has been interpreted as a velocity.

ß

Next,

consider

a maximum-likelihood inversion

constraint used was simply a

linear

ramp.

In

of

the

weOge.

The

this case, the shape of the

ß

wedge has been much better inversion.

That

is,

defined,

However, notice that

the maximum-likelihood

due to the broad-band

the resolution

inversion

inverted

has still

method also failed

bed thickness below 1/4 dominant period. recursively

limit

nature

of

the

been observeU.

to

resolve

the

The "pUll-up" observedon the

section is also in evidence here.

In summary,even though sparse-spike methods give an output section that

is visually

more appealing than recursively

appear to be a way to break the

inverted sections, there does not

low resolution

limit

of 1/4 of the dominant

sei smic peri od.

_


i

Page 7-

_

i

mk

6

Introduction to Seismic Inversi.on Methods

4

8

12

16

20

24

28

Brian Russell

32

36

40

44

48

o

300-

400.

Figure 7.5 Recursive inversion of wedgemodelshownin Figure 7.4.

4

8

'

12 16 ' • i '

20 ' I i

24

28

32

36

40

44

48

100-.................. 300

400

500

,,

Figure 7.6

.

Maximum-likelihoodderived impedanceof wedgemoUel shown i n Figure 7.4.


Page 7-

7

[ntroductJon


Methods

Br•an Russel•

PART 8 - MODEL-BASED INVERSION _

Part

8 - Model-based

-

Inversion

_

-

m

m

L

ß

....

•

Page 8 -

Introduction

8.1

to Seismic

Introducti

In the

information

Inversion

Methods

Brian

Russell

on

past

sections

directly

of

the

course,

we have derived

reflectivi-ty

from the seismic section and used recursire

produce a final velocity versus depth model.

inversion to

We have also seen that these

methods can be severely affected by noise, poor amplitude recovery, and the band-limited nature of seismic data. will

be included

In first

in

the

final

inversion

and comparing the model to

our

approac• is

intuitively data

itself.

seismic data.

We shall then use the

comparison between real and modeled data to iteratively

the model in such a way as to better

this

result.

this chapter, we shall consider the case of builaing a geologic moUel

results of t•is

of

That is, any problems in the data itsel f

match the seismic data.

shown in Figure 8.1.

Notice

that

update

The basic idea

this

method is

very appealing since it avoids the airect inversion of the seismic On the other hand, it may be possible to come up with a model

that matchesthe data'very well, but is incorrect. (This can be seen easily by noting that there are infinitely manyvelocity/depth pairs that will result in the sametime value.) This is referred to as the problem of nonuniqueness. To implement the approach shown in fundamental questions.

First,

Figure 8.1,

we need to answer two

what is the mathematical relationship

the model data and the seismic data?

Second,

between

how do'we update the' model? We

shall consider two approaches to theseproblems,the generalizedlinear inversion (GLI) approach outlined in CooRe and Schneider (1983}, and the Seismic Lithologic (SLIM) method which was developed in Gelland and Larner (1983).

Part

8 - Model-based

Inversion

Page 8 -

2

Introduction

to Seismic

Inversion

Methods'

Brian

Russell

CALCULATE UPDATE IMPEDANCE

ERROR

ERROR SMALL ENOUGH

NO

YES

SOLUTION = ESTIMATE

Model Based Invemion

Figure 8.1

Flowchart for the model based inversion

Part

8 - Model-based

Inversion

technique.

Page 8 -

Introduction

8.2

%o Seismic

Generalized

Linear

Inversion

Methods

Brian

Russell

Inversion

The generalized linear inversion(GLI) method is a methodw•ich can be. applied to virtually any set of geophysicalmeasurementsto determine the geologicalsituation whichproduced these results. That is, given a set of geophysicalobservations, the GLI method will derive the geological model which best fits

t•ese observations in a

least squares sense.

Mathematically,

if we express the model and observations as vectors

M:(m 1,m 2, ..... , mk) T=vector ofkmodel parameters, and T: (t1,t2, ..... , tn)T: functional

of n observations.

between the model and observations can be expressed

Then the relationship in the

vector

form

t i = F(ml,m2, ...... , mk) Once

the

functional

relationship

ß i : 1, ...

, n.

has been derived

observations and the model, any set of

between

model parameters will

the

produce an

ß

output. But what model?GLI eliminatesthe needfor trial

and error by

analyzing the error betweenthe model output and the observations, and then perturbing the model parameters in such a way as to produce an output which will produceless error. In this way, we may iterate towards a solution. Mathematically' F(M)

)F(MO)

= F(Mo)+ aT

•M,

MO-Initial•odel,

where

M:

true earth model,

AM:

changein model parameters,

F(M) : observations,

F(Mo):calculated valuesfrominitial

model, and

•)F(M O)

.2 •

Part

8 - Model-based

Inversion

= change in calculated values.

Page 8 -

4

Introduction to Seismic Inversion Flethods

Brian Russell

IMPEDANCE

(GM/CM3)(FT/SEC) X1000

IMPEDANCE

4.6

41.5AMPLITUDE

41.5 4.6

ß

ml

i

I

41.5 4.6

41.5

i

ß

,

-

ii

,i

i,

i

i

ii

,

ß

ß

,

,

i

:. __

b

Figure 8.2

A synthetic

c

d

e

test of the GLI approach to model based

inversion.

(a) Input impedance. (b) Reflectivity derived from (a) with added multiples. (c) Recurslye inversion of (b). (d) Recurslye

inversion

of (b)convolved

(e) GLI inversion of (b).

Part

8 - Model-based

I nversi on

with wavelet.

(Cooke and Schneider, 1983)

Page 8 -

Introduction

to Seismic

Inversion

Methods

Brian

Russell

But note that the error between the observations and the computed values i s simply

•F = F(M)- F(MO). Therefore,

the above equation can be re expressed as a matrix equation •F

where

= A AM,

A:

matrix with

of deri vatives n rows anU k columns.

The soluti on to the above equation would appear to be

•M

-1

= A

where A-l: matrixinverseof A.

•F,

However, since there are usually more observations than parameters (that

is, n is usually greater therefore

does

not

than k)

have

a

the matrix A is

true

inverse.

This

usually not square and is

overdeterminedcase. To solve the equation in that

referred

to

as

an

case, we use a least

squares solution often referred to as the Marquart-Levenburg method(see Lines and Treitel

(1984)).

The solution

is given by

•M: (AT'A)-IA TZ•F. Figure 8.1 can be thought of as a flowchart of the GLI methodif we make the impedanceupdate using the methodjust described. However, we still must derive

the

observations. convol utional

functional

relationship

The simplest

necessary to relate

solution

the

which presents itself

model to

the

is the standarO

model

s(t)

= w(t) * r(t),

where r(t)

= primaries only.

Cooke and Schneider (1983) use a modilied version of the previous formula

in which multiples

and transmission losses are modelled. Figure 8.2 is a

composite from their paper showingthe results of an inversion applied to

a

single synthetic impedance trace.

Part

8 - Model-based

Inversion

Page 8 -

6

Introduction

to Seismic

Inversion

•

Methods

Brian

'

Russell

• IMP.EDANCE x1OOO (GM/CM3)(FT/$EC)

ß ._ ,.o

. .:-:

. .:.

.

ß ß ,, ,,

. .• .....

. .:.........•:: :...., ...... .. :...lO ? "e'.

: :........:..:.-.-_-

........

,,

, ß ....-.

-.

4': -':::.•/-.:.!i!i..::..':.. :.:......:.':i•i.'-'-:.. :...•.•. '......'..'.. :.' }::! - ..'.:"

•

,

.

'

300M$

.

,

Figure 8.3

2-D model to test GLI algorithm.

The well on the right encountersa gas sand while the well on the left does not. (Cooke and Schneider, 1983) IMPEDANCE

AMPLITUDE

(GM/CM3! (FT/SEC)X1000 10

,,,.l

Figure 8.4

Model

traces

derived

from

m)del in Figure 8.3. {Cooke and Sc)•neider, 1983)

Part 8 - Model-based •nversion

Figure 8.5

38 10

A

38

B

GLI inversion of model traces. Compare with sonic log on right side of Fi•iure 8.3. (Cooke and Schneider, 1983) Page 8 -

7

Introduction

to Seismic

Inver. sJon Methods

Brian

Russell

In Figure 8.2, notice that the advantage of incorporating multiples in the solution is that, although they are modelled in co•uting the seismic

response,they are not included in the model parameters. This is a big advantage over recursire methods, since those methods incorporate multiples into the solution if they are not removedfrom the section. Another

important

feature

of

this

particular

method is

the

the

parameterizationused. Instead of assigninga different value of velocity at each time sample, large

geological blocks were defined.

assigneda starting impedancevalue, impedancegradient, time.

This

Each block was

and a thickness in

reduceU the numberof parameters and therefore

simplified

the

computation.However,there is enough flexibility in this modellingapproach to

derive a fairly

detailed geological inversion.

Wewill

now look at both a

syntheticandreal examplefromCookeandSchneider(1983). A 2-0 synthetic example was next considered by Cooke and Schneider (1983). Figure 8.3 showsthe model, which consisted of two gas sands encased in shale.

One well encountered the

sand and the other missed.

The impedance

profile of the discovery well is shownon the right. Figure 8.4 shows synthetic traces over the two wells, in which a noise component has been added.

Finally,

Figure 8.5 showsthe initial

guess and the

final

solution,

for which the gradients have been set to zero. Notice that although the solution is not perfect, the gas sand has been delineated.

Part

8 - Model-based

Inversion

Page 8 -

8


Brian Russell

I

I

YES

___•__•J ' '

- _

FINALMObEL

._ x•,

• .... r -• •;•,•

-.-'%•..

-cx-r. .,'•_;'•.:. ß-,• . . . .. .-. t .•..

Figure 8.6

I11 ustrated

flow chart

for

the SLIM method.

(Western GeophysicalBrochure)

Part 8 - Model-based Inversion

Page 8 -

9

! ntroducti

on to Sei stoic !nver si on Methods

Brian

Russel 1

8.3 Sei_smic L_ithologicModelling(,SLIM) Although the n•thod outlined in Cookeand Schneider (1983) promise, it has not, commercially.

as

far

However,

fully

released,

than

the

direct

of

a seismic

flowchart

the

layers

of

variable

points along the line.

similar

and

is

Modeling (SLIM) method of

perturbation

of

the

SLIM method taken from a Western

as in the GLI method, an initial velocity,

Also, the

of a model rather

section.

created and comparedwith a seismic section. of

very

Although the details of the algorithm have not been

inversion

Notice that,

aware, been implemented

appears

Seismic Lithologic

the method does involve

Figure 8.6 shows a

brochure.

author is

one method that

commercially available is the

Western Geophysical.

as this

showed much

density,

geological model is

The model is defined as a series and thickness at

seismic wavelet is either

various

control

supplied (from a

previous wavelet extraction procedure) or is estimated from the data. synthetic model is then comparedwith

the

The

seismic data and the least-squared

error sum is computed. The model is perturbed in such a way as to reduce the error,

and the process is repeated until

convergence.

The user has total geological information

control over the constraints and may incorporate from any source. The major advantage of this method

over classical

methods is

recurslye

that

noise

in the seismic

section

is not

incorporated. However,as in the GLI method,t•hesolution is nonunique. The best examples of

Gelland and Larner (1983). initial

applying

this

method to

Figure 8.7 is taken from their paper and shows an

Denver basin model which has 73 flat

boundaries

of

a

real data are given in

sonic log.

layers derived from the major

Beside this is the actual stacked

data

to

be

Page 8 -

10

inverted.

Part

8 - Model-based

I nversi

on

Introduction

to Seismic

Inversion

Methods

Brian

Russell

1kit

,1.4

1.4

1.6

2.0

2.0

Stack

Initial Figure 8.7

Left' Init)al Denver Basinmodelseismic. Right: Stacked section from DenverBasin. (Gelfand andLarner, 1983).

lkft

.4

1.6

1.8

2.0 •

1.8

Field data

Synthetic

Figure 8.8

Left:

Reflectivity

F•na• SLIM JnversJon of Figure 8.7 spl iceU into Right- Final reflectivity from ' -- _• -- __--__ii m ' -' (Gelfand and Larner, Part 8 - Model-based Inversion

data shown 1n field data. inversion. 1983). • .......... Page 8 -

2.0 .m: 11

Introduction

to Seismic

Inversion

Methods

Brian

Russell

In Figure 8.8 the stack is again shown in its most complexregion, with the final synthetic data is shownafter 7 iterations through the program. Notice the excellent agreement. On the right hand si•e of Figure 8-.9 is the final reflectivity section from which the pseudoimpedanceis derived. Since this reflectivity is "spi•y", or broad band, it already contains the low

frequencycomponent necessary for full inversion. Finally, Figure 8.10 shows the final inversion compared with a traditional recursire inversion.

Note the

'blocky' nature of the parameterbasedinversion when comparedwith the recurs i ve i nvers i on.

In

summary,

which can reflectivity

parameter

be thought is

of

extracted.

based inversion i s an iterative

as a geology-baseddeconvolution since the full I•

propagatedthroughthe final result

Part

8 - Model-based

Inversion

model1ing scheme

has the

advantage

as in recursire

that

errors

are

not

inversion.

Page 8 -

12

Introduction

w

to Seismic

Methods

500-ft N •

lkft

1

Inversion

Brian

114 mile S

Russell

114 mile SE

-

.5

.7

l m ß

1.9

Figure 8.9

Impedance section

derived

from SLIM inversion

of

Denver Basln 1 ine showni n Figure 8.7.

{GelfanUand Larner, 1983) W

lkft

50011 N

114mileS

ß

ß

114mileS

.•

E

1.7

19

Fi gure 8.10

Traditional

recursire

inversion

of Denver Basin line

from Fi gur.e 8.7.

(Gelfana anU Larner, 1983) Part

8 - Model-based

Inversion

Page 8 -

13

Introduction

to Seismic

Inversion

Appendix 8-!

Methods

Brian

Russell

Mat_r_ix .appljc.•at. ions_inGe•ophy.s.ics

Matrix theory showsup in every aspect of geophysicalproocessing.Before looking at generalized matrix theory, let us considerthe application of matrices to

the solution of a linear equation, probably the most important

application.

For example, let

3x1+2x2 : 1, and x1- x2 = 2. By inspection, we see that the solution is

However,we Could .haveexpressed the equations in the matrix form A

X

2

x1

1 -1

x2

3

=

y,

or

The sol ution

is,

ß

therefore x

=

-1

A

x1 x2

8 - Model-based

Inversion

1 . -2

1

1

1

-1/5

or

Part

y,

3

Page 8 -

14

Introduction


Me.thods

In the above equations we had the and the

problem t•erefore

of dimension N.

same numberof equations as unknowns

had a unique answer. In matrix terms, this

that the problem can be set up as a

vector

Brian Russell

square matrix of dimension N x N times a

However, in geophysical problemswe are Uealing with

the real earth anU the equations are never as nice.

fewer

means

equations than

unknowns (in

Generally,

which case the

we either

situation

is

have

called

underdetermined) or more equations than unknowns(in which case the situation

is is

calleU overdetermined). In geophysicalproblems,the underUetermined case of

little

interest

overde termi neU

to

case is

us since

there

is

no unique

of muchinterest since it occurs in

solution.

the

The

following

problems: (!)

Surface

consi stent

resi dual statics,

(2) Lithological modelling,and (3) Refracti on model 1i ng. The overdetermined system of categories-

consistent

extending our earlier

equations • can be split into two separate inconsistent. These are best described by

an•

example.

(a) Con•s.i s••t Overd..etermined Lin.earEqua.t. ion.s

In

this

equations

case

are

we

simply

reUunUant equations may

have more equations

than

unknowns, but

scaled versions of t•e others.

In

this

the extra

case,

the

simply be eliminated, reducingthe prø•lemto the

square matrix case. For examp.le, consider adding a third equation to our earl ier example,

so that

3x1+2x2 : 1, x1- x2 : 2, anU

Part

8 - Model-based

Inversion

5x1- 5x2 : 10.

Page 8 -

15

Introduction

to Seismic

This may be written

Inversion

Methods

Brian

Russell

in matrix form as

2

x1 x

o

$

But notice that the third equation is simply five times the second, and therefore conveys no new information.

equations back to the original

We may thus

reduce the system of

form.

(b) Inco, ns, is,ten•tO•verd. e•ermine. d L.i.near Equa•i.on?

In this case the

extra

equations are not scalea

versions

of other

equations-in the set, but conveyconflicting information. In this case, there is no solution to the problem which will

solve all the equations.

This is

usually the case in our seismic wor• and indicates the presence of measurement noise and errors. As an example, consider a modification to the preceding equations, so that

3x1+2xZ -- 1, x1- x2 -- Z,

ana 5x1- $x2 = 8. This may be written

3

in matrix form as

2 -I -5

Part

8 - Model-based

Inversion

xI

1

x2

2

8

Page 8 -

16

Introduction

to Seismic

Inversion

He.thods

Brian

Russell

'Now the third equationis not reducible to either of the other two, ana an alternate

of

least

solution

solution must be found.

The most popular aproach is the method

squares, which minimizes the sumof the squared error and the observed results.

That is,

if

we set the error

between the to

e=Ax-y,

then we si reply mini mize n

eTe-(eI , ez, .......

, en)

=

ei ß

2

Le. Re expressing the 'preceding equation in terms of the values x, y, and A, we

have

ß

E = eTe= (y - Ax)T(y - Ax)

= yTy_ xTATy _ yTAx + xTATAx. We then solve the equation bE_

bxi The final

solution to

the

least-squares

problem is given by the normal

equati OhS

AT Ax = AT y or

Part

8 - Model-based

Inversion

x = (ATA)-lATy .

Page 8 -

17

Introduction


Br•an Russell

PART g - TRAVELTIME INVERSION

ml ß

Part

g - Traveltime

Inversion

i

ii

Page 9 -

I

Introduction

to Seismic

Inversion

Methods

Sei smi c Travel

9,1

Brian Russell

Inversion

time

Introduction _

In

_•

L_

,

_.

this

_

section we will

look at a type of inversion

several names, incluUing traveltime tomography.

The last

term

inversion,

that

raypath inversion,

goes under ana seismic

tenUs to be overuseU at the moment, so it

is

important to use the term correctly. In section 9.3 we shall showan example whic• may be considerea as seismic tomography. As all of t•e other names suggest,

however,

seismic

measurements to infer

the structure

extracteU are velocities

structure

traveltime

and

can be derived.

of the earth.

depths,

Initial)y,

but Jr'has becomeobvious that

inversion

this

measurements can be used effectively

uses a set of traveltime The parameters

aria [herefore

which

a gross model of earth

this was considered the ultimate accurate

set

are

of velocity

goal,

versus depth

to constrain other types of

inversion.

For example, the'velocities could be used as the low frequency componentin recursire inversion, or as the velocity control for a depth migration. The way in which traveltime set of times from a dataset.

inversion is carried out is to first

These picks

pick

a

m•y come from any of three basic

types of seismic datasetsSurface

seismic measurements

- shots and geophones on the surface, VSP measurements

- shots. on surface,

geop•ones in well,

and

Cross-hol e measurements

- s•ots anU geophonesboth in well. Once the times have been picked, they must be made to fit subsurface.

In the next section,

we will

a model of the

look at some straighforwara

examples

of using traveltime picks in order to resolve the earth's velocity and depth structure.

Part

9 - Travel

time

I nversi on

Page 9 -

2

Introduction

to Seismic

Inversion

(a)

Methods

Brian

Russell

(b)

Figure 9.1 Travel paths through a single,

constant velocity

block.

(a.) Surface recording, (b) VSP recording, (c)Cross-hole recording

$

R

(b)

s•

(c)

Figure 9.2 Travel paths through two blocks of slightly

differing

veloc-ity.

(a) Surface recording, (b) VSP recording, (c) Cross-hole recording

Part

9 - Travel

time

Inversion

Page 9 -

Introductfonto SefsmJ c Inversfon Methods

Brjan Russell

9.• Numerical_Exa•mples of .Travelti• In•v. ersion

Consider the simplestpossiblecase, a constant velocityearth. Figure 9.1 showsthe travel paths that wouldresult from the three geometry configurationsgiven a square area of dimension L by L. Note that the traveltimes in Figure 9.1 would simply be: (1) Surface sei smic'

t--Z

L p

or

p-- t/Z

L,

(z) vsP-

t --•L p

or p -- t /i•L, ana

t=Lp

or

(3) Cross-hol e'

p=t/L,

where p -- ! / V.

Obviously, all three sets of measurementscontain the sameinformation.

However,if the velocity(or slowness p) andthe depth are both unknown, neitheronecanbe determined froma singletimemeasurement. Anevengreater ambiguity comes into play if we havea single measurement but morethan one

box. In Figureg.g this situation is shown. Notice that the equationsnow would involve three unknownsand only one measurement.

A moregeneralmodelis proposed in Bishopet al (lg85) an• Bor•inget al (1986). The earth is representedas a number of boxesof constantsize and

velocity. Although thevelocity of eachboxis a constant, thevelocity may vary from box to box. This is shownin Figure 9.3. The objective is thus to computethe seismic travel path through each box using the traveltime

measurements. A keyproblem hereis howto allowthe rays to travel through the boxes. The first order approximationwould be straight rays with no bending. However,i f Snell's law is use4, the problembecomes moredifficult to sol ve.

Part-g

Travel i'i me'-'in'ver'si on .....

Pag• g' '-

4'

Introduction

to Seismic

Inversion

Source

Figure 9.3

Methods

Brian

Russell

Receiver

Separation of the earth into small constant vel oci ty blocks for sei stoic travel time inversion.

(Bording et al,

Part

g - Traveltime

Inversion

1986)

Page 9 -

5


Brian Russell

Let us apply the straight ray approximationin the simplecaseof having simply two blocks of different velocity. In this case, we have coupled together both surface and VSP measurements. Two possible recorUing arrangementsare shownin Figures 9.4 and 9.5.

The situations illustrated

are

obviously oversimplified since we have assumeda straight ray approximation in

both boxes. That is, there is no refraction at the velocity discontinuity, and the reflection point is directly at the center of the two boxes. However, if we assumethat the velocities vary only slightly, this approximation is reasonable.

Let us start with the situation illustrated by Figure 9.4. In this case, t•ere is a single shot with geophonesboth on the surface and in a borehole at the base of the layer. If we assume that the sides of the boxes are unity in length (1 cm or m or km. m), the travel time equations are (1) For the. raypath from S to R

wherePl: 1/velocityin box1 P2: 1/velocityin box2

(Z) Fortheraypath fromS to R2:

t2= q•Pl+• P2. 2

2

Thus, the total problem can be expressed in matrix form as:

•

•

Pl

•r• •]• •

2

P2

tl :

t2

or Ap: t .

The solution to the previous equation is then

p = A-lt. Unfortunately,

a quick try at solving the above equation will show that

the Ueterminant of A is

Physically,

this

is

O,

which

telling

means that

the inverse is nonrealizable.

us that the two travel

paths spene equal

proportions of their paths in eac• box.

Part

9 - Travel

time

Inversion

Page 9 -

6


Brian Russell

P,,

S

x

P•"% P,'v,

Figure 9.4

Surface

for a single shot.

$•

R!

Figure 9.5

and VSP raypaths

St

Surface and VSP raypaths for two separate

Part 9 - Travel time Inversion

shots.

Page 9 -

7

Introduction


Methods

Brian Russell

A simple way to remedy this situation

raypath.

This is shownin Figure 9.5.

In

shot one-half a box length to t•e left

case, the traveltime

is to move the shot for the second.

this situation,

we have moveUthe

for the recorUing in t•e hole.

In this

equations are

(1) Forthe raypathfromS1 to R1:

tl: 1•Pl+•P2 (2) Forthe raypathfromS2 to R2In this case notice from the diagramthat

tan 0 : 1/1 $ : 2/3 = 0 6667, or B : 33 69o Thus

cos 0 = 0.8320

and (see figure)

x = 1/(2 x 0.832) = 0.6 y = 3/(•

x 0.832) = 1.8

y-x=l.2 Therefore

t2:1.2 Pl + 0.6 P2 ' Thus, the total

1.2

with

problem can be expressed in matrix form as'

•[•

Pl

tl

0.6

P2

t2

sol ution

Pl

0.6

- 2

t1

-1.2

2

t2

1

o.85

P2

Problem' Try to solve the above equation

and 1.1 kin/sec.

T•at is,

when the two velocities

work out the traveltimes

are 1.0

and plug them into the last

matrix equati on.

Part

9 - Travel time

Inversion

Page 9 -

8


Brian Russell

Pick seismic

Estimatevelocityat well usingsoniclogandVSP

reflectiontimes,t

Estimate V(x,z)byusing V(xo,z),

Initial Model

the reflection traveltimes and the

theassumption of verticalrays Start with

top layer

Layer

Computer forward

Stripping .Inversion PerturbV(x,z) byleastsquares or manually

Add another

layer

Final Seismic Model

Figure 9.6

modeltraveltimes, f, by normalraytracing

It- fll'

layersbeen

ii

Model iscomplete,, I

A possibleflowchart for seismictraveltim inversion. (Lines et al,

Part 9 - Travel time Inversion

1988)

Page 9 -


Brian Russell

9.3 Seism•i •c.T..omo. graphy The term tomography was first used in the medicalfield for the imaging of humantissue using Nuclear MagneticResonance (NMR) and other physical measurements. In the seismic field

the

velocity

field

it

has come to mean the reconstruction

of the earth by the analysis of traveltime

of

measurements.

Excellent overviews of tomographyare g.iven in Bording et al (1986), and Lines et al (1988}. You will find t•at the latter paper introduces the term "cooperative inversion" since both seismic and gravity measurementsare used in

the inversion, but that much of the technique used by the authors

can be

consi alered sei smic tomography.

Figure 9.6, taken from the paper'by Lines et al (1988), shows the chart that they propose for performing be considered quite general,

traveltime

inversion.

even though many variations

of it

flow

This method can are used in the

industry. Basically, the process starts with an estimate of the model which, in the flowchart shownin Figure 9.6, is deriveU from the sonic log and VSP measurements. Next, traveltime

picks are made from the seismic data.

In this

case, stacked CDPdata is usecl, but the shot profiles (or CDPprofiles)

could

,

also be used. refraction

As well,

arrivals.

travel time In

raytraced, and an error

traveltimes.

the is

next

picks

can be made from VSP data and

stage of

the

computed between the

process,

the

model is

computed and observea

Based on the error computed, a new model is computeU. This

done using the GLI

technique

described

procedure shownin Figure 9.6,

in Chapter 7 of these notes.

is

In the

the inversion is done layer by layer until

the

model is complete.

Although any traveltime

Stewart

inversion can be considered tomography, Dr.

(personal communication) points out that to be analagous with

medical field,

where physical

measurements are

Rob

the

taken completely around t'he

imaged object, a true seismic tomography experiment would involve aata on more than one side of the portion seismic

Part

of

the

earth

to

be

imaged,

such as surface

and VSP.

9 - Travel time

Inversion

Page 9 -

10


Brian Russell

WlC

•0 _

m

mß

i

roll i

-1Z

..2

0

x I• vsP SOURCE

• 3-D SOURCE • GEOPHONE

Surface geometry for tomographicimaging example.

Figure 9.7

(Chiu and Stewart,

Well C VSP

Une8901 898921 89D•B8960 898980 89CDP Une

CDP

0.0

.......... •::•=•'•:"•.::"--'::':.-:'::.i•r.:iE)•".Z•!;.".•h. •.

0.1 ---':......

-" '•'•":'":

1987)

Depth (m) 185 9O7

205 0

fi'•L .o.• :•(;:•

460

730

895

•1o

?6o

895

ß.• mo•w,• .'•.' •

•

.......

•

.• --'-..

oJ 0.4 o.5 ß

.

..

(b)

(a) Fi gue 9.8

(a) Picked events on 3-D seismic.. (b) Picked events

on VSP.

(Stewart and Chui, .....

Part

9 - Traveltime

Inversion

Page 9 -

11

1986)

Introduction


Brian Russell

An example of using multiple datasets for seismic tomography is found in

Stewart and Chiu (1986), and Chiu and Stewart (1987).

The objective was a

Glauconitic

Since

channel

sand which

containeU heavy oil.

this

was a

development survey, a lot of measurementswere available to imagethe subsurface, including well log data, VSP, and 3-D seismic. Figure 9.7 shows

the •ensity of information along a portion of one seismic line. Figure 9.8 shows the various datasets used in the tomographicimaging. Figure 9.8(a) showsthe stacked seismic data with the key events indicated and Figure 9.8{b) showsthe picked VSP from well C. Finally, Figure 9.9 showsthe well l'ogs and synthetic

from a different

well,

clearly

indicating

the Glauconitic

channel.

The tomographic technique involved picking events from both the VSP first

arrivals and the prestack 3-D seismic data. Traveltime inversion wasdone by the

technique described

in Chiu and Stewart (1987).

The method involves

starting with a simple modelof the subsurfaceand perturbing this modelusing the errors between the picked traveltimes and the raypath times through the model. This raytracing

method differs

from the

method shown in

is done a nonzero source to receiver

offset,

Figure

9.6 since

and also the VSP data.

To test the method, Chiu and Stewart created a synthetic model. Figures g.10(a) and (b) show raytrace plots for the VSP and surface Uata, respectively,

through this model.

ZERO PHASE BANDPASS 10/15 - 80/110 Hz RFC

NORMAL

DENSITY (kg/m 3 )

VE-OCITY(m/sec) SOIl

lime

030O

(sec) till

Figure g. g

ß

Wel1 log curves and synthetic

IO#

ß

showing Glauconitic

channel.

(Stewart and Chiu 1987) Part

9 - Travel

time

I nversi on

Page 9 -

12

Introduction


Methods

Brian Russell

offset (krn) 1.o

0.0

#

2.0

$Oul•CE

A •OPHONE'

i

i

i

i

i

i

(a) Figure 9.10

{•)

(a) Surface raypaths through model used to test inversion. (b) VSP raypaths through model.

1.0 ,

ii

!

2.0 ! i 1! -

1987)

Voity

Offset (km) 0.0

(Chiu and Stewart,

0.0

2.0

4.O

---

a

Figure 9.11

Results of tomographicinversion of model data

using VSPand surface data. Part

9 - Travel

time

Inversion

(Chiu and Stewart,

1987)

Page 9 -

13


Brian Russell

Figure9.11 shows the results of the inversionprocessusingboththe VSP and surface seismicdata. To makethe test morerealistic,

randomnoise was

added to travel'time picks. Notice that thecorrectresulthasbeenobtained in

four

iterations.

Let us now return to the case study described initially. The final velocity/depth model is shownin Figure 9.12. Notice that the velocities fit

quite well with the averagedsonic log velocities. This velocity model was usedto produceboth a depthmigrated seismic section, shownin Figure 9.13, and a full seismicinversionbasedon the maximum-likelihood technique. The final inversion is not showndue to colour reproduction limitations. As can be seen in Figure 9.13,

the Glauconitic

delineated. The depth tie is also excellent.

channels have been well

The conclusion that the authors

makeis that if severaltypesof geophysical measurements can be intergrated, the result is an improvedproduct. Each set of data acts as a constraint the

Part

on

others.

9 - Travel

time

Inversion

Page 9 -

14

Introduction

to Seismic Inversi on Methods

Offset(km)

-1.0

0.0

Brian Russell

Velocity Oan/s)

LO

0.0

3.0

6.0

TC)iliO6RAPmC (C) INVERSION

- SONCL06 (1:2)

Figure 9.12

Results of tomographic inversion of G1auconitic channel.

(ChiUandStewart,1987) .ß ß: .•m,,, ...

.......... ' 't ''•"',' l..,.,.,;,,•.

J•

ß -.--:'

._:_.4sl•l •,_, i!',i,? ,a•.. ,:.I.,,t.:, ?

800 .. : 900

:'":

""' ""' .......

'

Depth

(m) 1000 11oo 1200

1300 1400

Fi gure g. 13

Depth migration of seismic aatashown in-Figure 9.8(a). Tomographic velocities of Figure9.12 havebeenused.

(Stewart an• Chiu, Part

9 - Travel time

Inversion

Page 9 -

15

1986)


Brian Russell

PART 10 - AMPLITUDE VERSUS OFFSET INVERSION

Part 10 - Amplitude versus Offset

Inversion

Page 10 -

1

Introduction

to Seismic

Inversion

Methods

Brian

Russell

10.1 AV.O Theo.•y.._ Until

traces.

result

now, we have discusseO only the inversion of zero-incidence seismic

That is, we have considered

of

a seismic ray striking

each reflection

coefficient

to be the

the interface between two layers

at

zero

degrees. In this case, the 'reflection coefficient is a simple function of the P-wave velocity and density in each of the layers.

The formula, which we have

seen many times, is simply

i+lvi+ - ivi

zi+- zi

ri= Yoi+iVi+l+ yO iVi where

r:

•Zi+l + Zi reflection

coefficient

yo: density, V = P-wave vel oci ty,

Z: and

Layer i overlies

When we allow the seismic ray to strike

angles,

as

results.

coefficient

in

a commonshot recording,

In this case, there is P- to

becomes a function

density of each of the layers.

derived:

reflected

acoustic impedance, Layer i+1.

the boundary at nonzero incidence

a much .more complicated S-wave conversion

of the P-wave velocity,

situation

and the reflection

S-wave velocity,

and

Indeed, there are now four curves that can be

P-wave amplitude, transmitteU P-wave amplitude,

reflected

S-waveamplitude,and transmittedS-wave amplitude. The variation of ß

amplitude

with

Poisson's ratio,

offset

also

involves

another

physical

parameter

called

which is related to P-and S-wave velocity by the formula

•' =-

(Vp / VS• 2- Z .

Poisson's ratio can theoretically

vary between 0 and 0.5.

Part 10 - Amplitude versus Offset Inversion

Page 10 -

2

Introduction


$i

Brian Russell

Sr

at, •t BOUNDARY

(X2'•2

t

$•

Figure 10.1

Reflectedandtransmittedrays createdwhena P-wave strikes the boundarybetweentwo layers.

(Waters, 1981).

•o,2+, -•sin2•, - ' 'cos2•,- •x,n-•:/ •D,/•-cos2+,/ Figure 10.2

Zoeppritzequationswhichdescribethe amplitudes of the rays shownin Figure 10.1.

(Waters, 1981).

Part 10 - Ampli rude versus Offset

Inversion

Page 10-

3

Introduction


Methods

Brian Russell

The equations from which the ampl'itude variations can be derived are callea

the Zoeppritz

equations. They are derived from the

continuity

of

displacement and stress in both the normal and tangential directions across an interface

between two layers.

Figure 10.1 showsthe seismic rays

across a

boundary, and Figure 10.2 gives the final form of the equations. They are taken from • textbook by Waters (however, someof the signs were wrong, and

they are fixed

in

the diagram).

Since we have four equations with four

unknowns,they can be rearranged in the form of a ½ x 4 matrix equation Ax--y with

soluti

on

x = A-ly . Over the years, several effects.

authors

have discussed amplitude versus offset

However, these authors concluded that the effect would be negligible

on seismic data.

In

a landmark paper, Ostrander (1984) showedthat for a

significant changein Poisson's ratio, a major changein the P-waveamplitude coefficient can be seenas a function of offset. This Poisson's ratio change is

most noticeable

in a gas sand, where the ratio .

encasing shales to as low as 0.1 that,

in

the

gas sand itself.

such extreme cases, the P-wave reflection

positive to negative for increase

in

can change from 0.4 in

Ostrander showed

coefficient

a decrease in Poisson's ratio

in P-wave velocity,

or from negative to positive

the

can go from

coupled with an

for an increase

in

Poisson's ratio coupled with a decrease in P-wave velocity.

Figure 10.3(a) showsthe gas sand model that Ostrander used and Figure 10.3(b)

shows the result of

reflection shows that

coefficients. this

effect

can

amplitude

versus

offset modelling of the P-wave

Figures 10.5(a) and (b), also taken from Ostrander, inUeed

be

observeU

on

a

common offset

stack.

Figure 10.5(a) showsa stackeOseismic section witl• three apparent "bright spot" anomalies. Unfortunately, only wells A anU B were productive. The three common offset stacks, shownin Figure 10.5(b), indicate that only locations A ,

and B actually Uisplay an AVOeffect.


Page 10 -

4


Brian Russell

GAS*':*"*•* t•t

Vl•Z=8.000 /32 -"2.14

;..•.•::., SHALE

:o.,

ß

$=10.000

•---'

/4)3 =2.40 (•'3 =0.4

Figure 10.3 (a)

Synthetic

gas sand model.

(Ostrander,

1984)

0.41 IN SAND 0.3

0.2

t..,

0.1

ANGLE OF INCIDENCE - •,-'-e'

0

I0 o

20 o

30 ø

40 ø

0

..,

NO

GAS

.,,,,,,,.,.ooo.o.o..... ..ooo., ßo.,.,o.,-'. oø,.,,,o .,,o*oo...... -0.2

-0.3

-0.,4

Figure 10.3 (b) Computedreflection coefficients as a function of offset for reflections

from top and bottom interfaces

of model

s•own in Figure 10.3 (a). • ,

,

,

(Ostrander, IlL

_

--

Part 10 - Amplitudeversus Offset Inversion

_,

11

i

,

i

m m

im

,

1984)

ß

Page 10 -

Introduction

The

to Seismic

Inversion

Methods

method useU to identify

this effect

Brian

is only partly

Russell

qualitative,

and

can be diagrammedas shownbelowINPUT

BUI L D MODEL

SEISMIC

SHOT PROFILES

COMPUTE

COFFSTACK

SYNTHETIC

im

...

VISUAL

.

COMPARISON

m

MANUALLY

•

CHANGE PARAMETERS

-

,m

NO

MODEL MATCHES REALITY .

.

Figure 10.4

Flow

Chart

for

Manual

AVO Inversion

Obviously, this visual meth'od of comparison leavesmuchto be Oesired. We will therefore look at several methods for the qualitative

inversion of AVO

data,

the

both

of

normal-incidence

which have been looked at

previously

in

context

of

inversion. ß


Inversion

Page 10-

6

Introduction

1•0

0.0


170

160

1•50

140

Methods

130

Brian Russell

120

110

100

90

•0

70

• , -- • , • • " ! :• "•• :. ' '' -' ' : •.... •- ' -• I; " • • -•'I'*O..'ß,•m * .... •. re-!:ß ß' l, 't ß .". '.'..' ' t'-' ß I , I i I, mm I.. i ' ,. '

'

..,.,,.,..•..,.,,.... •,",',,-. ..... --,,,,'.,,.-,',,,,, ..... •,,.,..... ,,,' ......... ,-•~., ....... •,,.'""'....'!.'_'•"' ..... •,,., ............. '.... ß

'

i

i i.

ii I

I

ß IIII

II

i

ß

0.0

'i: _--.,..•.•...'_'_-•..• ,:•:;.-': ', • ..=•..•:.' -':."•i "=-" '.•.•.'-'?'?'?'?'?'•L--•,•.•'•.•.•..,':: .... ß•-..c• .._ .o..•:?;i•. ' --?.-•" ":'.. : ' • ...... il •;,.•:.?:•-•=;.•'.'•..='.:•:-1: . . ,•.•. .-...,_•....:,•.•.:..•,.•;,..-.

•

•

.*... .... .

ß

.- -...--- . ; ...... ..• .•:..;,-.;:.•..

0.5 •'.l_'.-•. .: : •_•_•..._. .... _

_•:_- .•...4. .....
. .:.;. ,""-•.•r..•_-•".::

-:.:...:..._...• ........ _:.....;._.•.= •,'*' ....... ':'.:-r-'_.•.•; . . 0.5

1.0

I.,...,:,•....,;.•. .........:.,...-..

ß

...

....,..,

t,,,

"_,,,d,.•, I•l.leeile*e,,I,'t ! :litI•ll""' IIt•'•l';I;,d

..........12.0 .

..

,•e.•. Illl•lll.1111el'•ßß..le

•............ ;.; .... :.....i:-; - .. • '."-'•1, ø.'•1•-.

ß'.:--;•=....=:;;:.1• .... ...,;•.".....

2.5

'......... "•"":": ....................... 2..5 ß ß -..... ,.....-.,.-.,.-..-,.,...........,............

,..,- .1,• ..----?'"'1

"' '

Figure 10.5(a)

ß........;:.-.'..

.

"1;;::= ....:":" .....;""':"' ""1'.'-'-'-' ' "•':,';;;;":',:: ....:"

ß

'}i:.;.=i.•;-."':;::.'.:•: '.' ß

ß

.......;.::.;:..

S•acked seismic line showing"brigh• spot" anomalies. Loca%ions

A an•

B are

known gas.

(Osl:rander,

!984)

SP 80

":l:1111il•

.... ,e*'11:,l:,

ol,, titIll, ....

ß

.

eellie

;;;;;;;;i;il ,,

' ' ß

.

I•

..111tl•

' ........

6•$•'

1012'

.

,

6952'

1012 ß

Fiõure 10.5(D)

Common offset sl;acksover locations A, B, andC from stackedsection in Figure 10.5/a). Notice the AVO increase

on A and B.

(Os!;ranOer,

Par[ 10 - A,•pli[ude

versus Offset

Inversion

1984)

Page 10 -

7

Introduction

to Sei stoic Inversi

10.2 AVO Inversion

Recall

on Methods

Bri an Russell

by GLI

that

in the theory of generalized linear inversion

were three important components'a geological relationship

(.GLI)

there

modelof the earth, a physical

between the earth and a set of geophysical measurements,anU a ß

set of geophysicalobservations. This methodWasdiscusseain both chapters8 anU

9,

applied

to

stacked data

respectively.

Now, let us apply the

wil 1 be cal led

AVO inversi

In

secti on 10.1,

inversion method to

an• traveltime unstacked data.

inversion, The result

on.

the three components needed to perform GLI inversion

on

AVOUatawe,redescribed. Ourmodel of the earth is a series of layers with t•e el astic 'parameters of P-and S-wave vel oci ty, density, and Poi sson'S ratio. Our physical relationship

between this

derived using the Zoeppritz equations.

model ana seismic CDPprofiles And, finally,

picked amplitudes and times of events on a CDPprofile

was

the observations are the

or common offset

stack.

By computingderivatives from the Zoeppritz matrix, it is possible to set up a GLI solution data.

This

to t•e AVO problem similar solution

to the solution

found for zero-offset

is

aF(Mo•)

FIM): F(M D)+ •)M bM, where

Mo:initial earthmodel, M:

AM:

true earth model,

change in model parameters,

F(M) -- AVOobservations,

F(MO):Zoeppritz valuesfrominitial model,and

•)F(M O)

i)--••: change in calculated values.

The implementation is simply a variation

of

the manualmethod, anO is

sinownon the next page.


Page 10 -

8


Brian Russell

INPUT SEISMIC SHOT PROFILES

COFFSTACK J

COMPUTE

SYNTHETIC

i

PICK

STORE COMPUTED

AMPLITUDES

AMPLI TUDES

i

t •

COMPUTE

COMPUTE MODEL

ERROR

PARAMETER CHANGE USI NG GLI

NO

ERROR

YES

MODEL MATCHES

REAbITY ,

Figure 10.6

AVO inversion by the GLI method


Inversion

Page 10 -

Introduction


Brian Russell

Wewill nowlook at an example of GLI inversion of offset

data.

amplitude versus

First, consider the integrated well logs shownon the left hand

panel of Figure 10.7. Actually., only the sonic log or P-wavelog wasrecorded in the field. The density log was derived from the sonic using Gardner's equation, the Potssonratio wasfixed at 0.25, and the S-wavewasderived from the

P-wave and Potsson ratio logs.

On the logs, three

layers

have been

blocked at depth and a significant Poisson's ratio changehas been introduced in the middle block. On the right hand side of Figure 10.7, notice that the amplitudeversusoffset curves have been displayed for the third layer. As predicted earlier, the P-wave reflection coefficient displays a strong increase of amplitude with offset.

Figure 10.8 showsthe sameset of blockedlogs on the left, but showsthe seismic responseof the amplitude change on the right. This synthetic was produced by simply'replacingthe zero-incidenceamplitudeswith the amplitudes derived from the Zoeppritz calculations. The events between600 and 700 msec display a pronounced amplitude changewit h offset.

ZOEPPRII'Z $IHI:LE INTERFFJCE

TESTLO TESTLO TESi'DE lE•-S t•;$TPO

Eq, m,nt: 3 Ti,•:

$7• Depth: 795

589

1998

Of*f's•:c' ..........................

.

i Reflected

.....

•

2•

2,5

4;8

P-Wave

Transmitted P-Wave(-9.8)

...... Ref'l ected S-Wave ........... Transmitted S-Wave ,,,

ß

,

mm

i

m

Figure 10.7 Blockedwell logs on left,

with computedZoeppritz curves

for layer 3 on right. _

Part 10 -^mplJtude versus Offset Inve•sJon

_

Page 10-

10

Introducti

on to

Sei smi c Inversi

on Methods

Brian

Russel 1

20EPPRITZ •EFLECTIUITY tlODEL

TESTLOG TESTLOG IEH$ITY1 u•Ym u•/m 9/½½

......

MODEL1 (meters.)

S-I,IFIUE! POISSOH1 EU usam

,L

909

727

545

363

.

181

-

ß.•..o...-.. ............... ,

....,....

•ee-. •......................... -.•" ......... ......... .............. •.'.........,.

:: ......................... •...• ........... " 268

268

Figure 10.8

2.5

468

LeftRight:

.$

A "blowup"of the blockedlogs shownin Figure 10.7. A synthetic

commonoffset

stack and the AVO curves

shown on the right of Figure 10.7.


Page 10-

11

Introduction


Methods

However, does the change seen in

Brian Russell

Figure 10.8 reflect

the reality

of the

situation? Figure 10.9 showsa set of CDPgathers which correspond to the model. The gathers are a realistic modelled dataset and were generatedwith no changein Poisson's ratio. Since the gathers are noisy and contain fewer

tracesthanthe synthetic CDP profile shownin Figure10.8, theywereusedto create

a commonoffset stack.

The geometry of this st.ack

is

described

in

Ostrander'spaper, and the resulting gathersare often referredto' as Ostrander gathers. Traces within a CDP/offset window were gathered and stacked, resulting in increased signal to noise. Figure 10.10 showsa display of the logs, synthetic model, and common offset stack. The mismatch in amplitudes is now obvious. ß

ß Next,

the

amplitudes of

the

event on the

contanon offset

stack

correspondingto the event displayed in Figure 10.7 were picked. The event above the

anomalous layer was also picked.

The picks were then

used along

with the computedamplitude versus offset curve to invert the data by the GLI method. In the inversion, two parameters were allowed to vary- the Poisson's ratio

in the layer of interest,

and a

scalar

which relates the magnitude of

the seismic picksto themagnitude of theactual'amplitudes.

Figure 10.9

CDPgathers from a seismic dataset corresponding to synthetic shownin Figure 10.8.

Part 10 -Amplituae

versus Offset Inversion

Page 10 -

12

I ntroducti onto SeistoicI nversi onMethods

Bri an Russel 1

I NVER$]ON FULL HOIIEL

TESTL TESTL I)EN$I $-WI:IUPOI$$ _ us/m

u•/•

g/cc

us/•

rIO]]EL1 (meters) EU

909

727

545

353

181

0

COFFSTK1( n,elers )

838

6•4

498 :)32

li•E;

ItßßIII 50

I 2•0

2•0

2.5

4•0

Figure 10.10

,5

A comparison of the synthetic coneon offset stack from Figure 10.8 {middle panel) with a con,nonoffset stack created from the CDPgathers of Figure 10.9 (right panel). T•e left panel showsthe blocked well logs from which the synthetic

was created.


Page l(J-

13

0


Brian Russell

The results of this inversion are shown in Figure 10.11. The figure

shows the change in Poisson's ratio beforeandafter inversion(dashed line before, solid line = after) on the left hand side. Onthe upperright is shownthe matchbetween the observed picksin the upperlayer (shown as small

squares) andthefinal theoreticalcurve{s•own asa solidline). Thelower right showsthe samething for the lower layer.

Finally, Figure10.1Zshowsthe comparison betweenthe coanon offset stackandthe syntheticmodelafter the modelhasbeenrecomputed with the new amplitude changes fromtheupdated Poisson's ratio. Noticethe improvement in the

match.

.

II•'RSIOH

SIN•E

LI•ER:

!

ß

I101•ELI

Ewnt (2) P.bove Laver

e .e•

0.042

6.666

.

.

..

O•'•'r•

Event (3) Belo4aLaver

0.048

OO

0

6.624

70O

i

6,8

Figure 10.11

i

Poi•s•s

i

Ratio

i

e.•

6.606

e

S5e

O('•set

( m)

Theresultsof a GLI inversionbetween the computed, amplitudes

of Figure10.7andthepickedamplitudes fromtheconmon

offset stackof Figure10.10. Thedashed line onthe plot onthe left is the Poisson'sratio beforeinversion,andthe solid line is after inversion. Theplots on the right show the newcomputed curveswiththe picks(squares) superimposed.


Page 10-

14

Introduction

to Seismic

Inversion

Methods

Brian

Russell

IHUER$IOH FULL MODœL

MOIlEL2(me•er$)

EU

Figure 10.12

909

727

545

363

COFFSTKI ( meters )

lB1

B

838

664

498

332

166

A replot of Figure 10.10, where the synthetic has been

recomputedusing the newPoisson's ratio value.


Inversion

Page 10 -

15

[nt•oduct•

on t• $e• sm•c [nve•sJ on Methods

B• an Russe• ]

PART 11 - VEI:OCITY INVERSION

Part 11 - Velocity

Inversion

Page 11 -


Brian Russell

Part 11 - Veloci..tyI.n.v. ersi on __

11.1

! ntroduc ti on The

last

inversion.

topic to be discussedin these notes is the topic of velocity

Alth

acutally

fit

ough this technique is

referred

to as inversion, it does not

in to the narrower category of inversion techniques that we have

been dtscussing in this course. These techniqueshaveall involved inputting

a stacked, or unstacked, seismic dataset and inverting to a velocity depth section.

seismic

The output of

section

amplitudes,

the

velocity

inversion

properly positioneU in depth, but still

and still

band-limiteU.

versus

described here is the

plotted

as seismic

As such this technique is closer

to that

of depth migration.

In

this

section,

we will

inversion, and then look at a

this

at

few examples.

the theory of

summary. Our discussion here will

11

- Velocity

Inversion

and

there

follow

velocity

An excellent review article

subject is given in Bleistein and Cohen(1982).

theory of the method is reviewed

Part

look briefly

is

In this

article,

also an extensive

on

the

literature

that article.

Page 11 -


Brian Russell

KII.OFEœT

-2

-1

0

KILOFEET

I

2

-2

-I

0

I

2

,

(a) Figure 11.1

(b)

The effect of the velocity inversion methodon synthetic

data. (a) A "buried focus"effect, (b) The output from the velocity

inversion

method.

(Bleistein

KILOFEET o

1

-1

C)

and Cohen 1982 )

KILOFEET o

1

ß

ii'1 uJ LL m

o ....

(a) Figure 11.2

(b)

A second example of the effect

synthetic data.

of velocity

inversion on

(a) Input section with diffraction,

(b) Output from velocity

inversion. (Bleistein

and Cohen

1982 )

......

Part 11 - Velocity

Inversion

Page 11 -

Introduction

to Seismic

Inversion

Methods

Brian

Russell

11.2 Theory and .Examples

The velocity inversion procedure is referred to as an inverse

scattering

problem, in which the interior of the earth is mapped by inver.ting the observations from multiple acoustic sources. (This is a long way of saying that the seismic

section

is

inverted!)

method is the acoustic wave equation.

and classical migration is

that

Thus,

the

point for this

The difference between this

perturbation

techniques

transforms are used rather than downwardcontinuation The initial

starting

technique

and integral

of the wave equation.

work in this area was done by NormanBleistein

and Jack Cohen

at the University of Denver. In their initial paper, Cohenand Bleistein (197g), they employed only a perturbation technique in the inversion of seismic data. In simple terms, this technique involves using a constant velocity in the wave equation, perturbing this constant velocity by a small amount, and then, by observing the backscattered wavefield, solving for the perturbed velocity.

This methodsolves for only the reflection

strength

of

the mappedinterfaces.

In their morerecent paper, Bleistein solution

was proposed which al.so solves for

refraction. inversion

and Cohen(1982), a moreaccurate

Clayton and Stolt of

seismic

data.

(1981) Their

transmission losses

and

have applied a similar methodto the

method is referred

to as

the

Born-WKBJ

method, and thus this approach to inversion is often cal led Born inversion.

Despite the differences in the mathematicsbetweenthe velocity inversion methodsand migration methods, the results look very similar to those of

migration. For example,Figure 11.1, fromBleistein andCohen(1982), shows the input an• inverted

result

for

a g-D buried focus.

Note that, as in

migration, the "bow-tie" has been imagedto a synclinal feature.

Part 11 - Velocity

Inversion

Page 11 -

4

Introduction

•o Seismic Inversion

q ! ()f!. []

I,;,(11). iJ

ll;11]ll. ()

ß

ß

Methods

,

I I 3l)!]. l,

Brian Russell

I ;i'llroll. IJ.

I I,•I O0. II

I I11.)[11J. fl

•

ß

½J

ß

(a)

qlOO .-.,,

6500

8900

!

!

I 1 300 -

1 37.r.,P

16' ,'!..[]

t

18b•G

,

c) ß

(b)

Figure 11.3

The effect

of velocity

inversion on real data.

(a) Input section (MarathonOil),

(•) Output section. (Bleistein

Part 11 - Velocity Inversion

and Cohen 1982)

Page11 -

5


Brian Russell

Figure 11.2, also from Bleistein and Cohen (1982), showstl•e velocity inversion of a diffraction tail from a geological discontinuity. Notice that the diffraction tail has been "collapsea", again as in migration.

Finally, Figure 11.3 showsan example of applying the velocity inversion techniqueto a real dataset. Again, note the similarity with classical depth migration. The fact that this section is plotted as wiggle trace only makes the plot di fficul t to evaluate.

In summary,this technique cannotbe classed with the other methodswhich have been discussed in this'course due to its similarity with depth migration. However, research in'this area is continuing at a steady pace, and the technique promises much for the future.

Part

11 - Velocity

I nversion

Page11 -

6

Introduction to Seismic Inversion Nethods

Brian Russell

PART 12 - SUMMARY

Part 12 - Summary

Page 12 -

i

Introduction


Methods

Brian Russell

lZ.1 Sgmmary In inversion

these

notes,

we have

of seismic data.

reviewed the current

The basic

model

methods used

in

the

used in most of these methods is

the one-dimensionalmodel, which states that the seismic trace is simply the convolution of a zero phase wavelet

the

earth's

acoustic

shownin Figures

with a reflectivity

impedanceprofile.

12.1,

12.2,

sequence derived from

Flowcharts for these methoUs are

and 12.3.

Let

us initially

summarizethe

advantages and disadvantages of the three methodsof single trace which

inversion

have been discussed:

(1) Recursire _

!

Inversion

,•

_

-

•

Advantage s:

(i)

Utilizes

the complete seismic trace in its calculation.

(l i ) A robust procedure when used on clean seismic data. (iii)

Output is in wiggle trace format similar

to seismic data.

Di sadvantages:

(i)

Errors

are

propagated through the recurslye solution

if there are

phase, amplitude, or noise problems. (i i) The low frequency componentmust be derived from a separate source.

(2) Spar. se-SP.i kg_.Invers. ion Advantages-

(i)

The data

itself

is

used in

the

calculation,

as

in

recursi

ve

i nver si on.

(ii) (iii)

A geological looking inversion is produced. The low frequency information

is included mathematically

in

solution.

Part 12 - Summary

Page 12

the

Introduction


Methods

Brian Russell

BAND-LIMITED SEISMIC TRACE INTRODUCE LOW FREQUENCY COMPONENT

REFL

COEFF.

I

INVERT

I

TOIMPEDANCE

IMPEDANCE

SCALE TO VELOCITY AND DEPTH

DISPLAY

Fiõure

12.1

Band-Limited

Inversion

(Recursive)

ß

Part

12 - Summary

ß

Page 12-

,

3

Introduction

to Seismic

Inversion

Methods

Brian Russell

,

Dôsadvantages'

Statistical

(i)

nature of the sparse-spike methodsused are subject to

probl eros i n noisy Uata.

Final output lacks muchof the fine detail seen on recursively inverted data. Only the "blocky" componentis inverted. (3) Model-Base• I nver si on

Advantages'

(i) A complete solution, including low frequencyinformation, is possible to

(ii)

ob rain.

Errors are distributed

(iii)

through the solution.

Multiple and attenuation effects can be modelled.

Di sadvantages'

(i) A completesolution is arrived at iteratively andmayneverbe reached ( i.e. the solution maynot converge). It is possible that more than one forward modelcorrectly fits data (nonuniqueness).

(ii)

The

other

methods which

were considered

were

traveltime

the

inversion,

velocity inversion, and amplitude versus offset inversion. All are important methods, but cannot be compareddirectly with the three previous methods (comparingapples with oranges?). ,

The traveltime

accurate velocity constraint

for

either

inversion

method was

an excellent

versus depth model. These velocities one of

the

classical

method for finding an make

an

excellent

inversion methodsor for a depth

migration.

Part

12 - Summary

Page 12 -

4


BrianRussell

EXTRACT

SPARSE REFLECTIVITY INTRODUCE LINEAR CONSTRAINTS

INVERT TO IMPEDANCE

I

vELøcmTY!

DL_•.••_•. AND m

Fiõure

Part 12 - Sugary

12.2

i

i

m

i

Broad-BandInversion (Sparse-Spike)

Page 12-

Introduction


Methods

Brian Russell

The velocity inversion methodwas shownto be very similar to depth migration. The output from this method could therefore be used as input to one of the other

Finally,

three

classical

amplitude versus offset inversion adds an extra dimension to the

inversion problem since it velocity future,

methods of inversion.

inversion but still

is

method.

truly This

has a number of

a lithologic

inversion

method is definitely hurdles

the

rather than a method of

the

to overcome. This author's humble

opinion is that once the interpreter is able to do a complete lithological inversion on their seismic datasets, the other methodswill be replaced. The other conclusion from this course is that the more separate datasets (surface seismic, VSP, well log, gravity, etc..) the interpreter can use in an inversion,

the better

Part 12 - Sumnary

the final

product will

be.

Page 12-

6

Introduction


Brian Russell ß

MODEL TRACE

CALCULATE ERROR

ERROR SMALL ENOUGH

IMPEDANCE ESTIMATE

UPDATE IMPEDANCE

NO

YES

ON

=ESTIMATE

Fiõure

Part 1'•- •'" .... Summary

12.3

Mode 1-Based

Inversion

Page 12-

Introduction

to Sei stoic Inversion

Herhods

Brian Russell

REFERENCES

Angeleri, G.P., and Carpi, R., 1982, Porosity prediction from

seismic data'

Geophys.Prosp., v.30, p.$80-607.

Berteussen, K.A., and Ursin, B., 1983, Approximate computation of

the acoustic impedancefrom seismic datap.

Geophysics, v. 48,

1351-1358.

Bishop, T.N.,

Bube, K.P., Cutler,

Resnick, J.R.,

Shuey, R.T.,

R.T., Langan, RT., Love, P.L.,

SpinUler,

Tomographic determination of velocity

D.A., and Wyld, H.W., 1985,

and depth in laterally

varying media- Geophysics, v. 50, p. 903Bleistein,

N., and Cohen, J.K.,

Present status,

Bording, R.P.,

198•, The velocity inversion problem-

new directions:

Lines, L.R.,

923.

Geophysics, v.47. p.1497-1511.

Scales, J.A.,

ana Treitel,

S., 1986,

Principles of travel time tomography' SEGContinuing EUucation notes, Geophysical inversion and applications.

Chi, C., Mendel, J.M., and Hampson,D., 1984, A computationally fast approach to maximum-likelihood aleconvolution: p.

Geophysics, v. 49,

550-565.

Chiu, S.K.,

and Stewart,

R.R.,

dimensional seismic velocity

seismic profiles,

1987, Tomographic determination

structure

of three-

using well logs, vertical

and surface seismic data:

Geophysics, v.52,

p. 1085-1098. Claerbout,

J.F.,

and Muir, F.,

Geophysics, v. 38, p.

Part 12 - Summary

1973, Robust Modeling with erratic

data:

8Z6-844.

Page 12 -

Introduction

Clayton,


R.W., and Stolt,

acoustic reflection

Cohen, J.K.,

R.H.,

data:

and Bleistein,

Brian Russell

1981, A born WKBJinversion

Geophysics, v. 46,

method for

1559-1568.

N, 1979, Velocity inversion procedure for

acoustic waves: Geophysics, v. 44, p. 1077-1087. Cooke, D.A.,

and Schneider, W.A., 1983, Generalized linear

of reflection Galbraith,

J.M.,

seismic data: and Millington,

Geophysics, v. 48, p. G.F.,

inversion

665-676.

1979, Low frequency recovery in

the inversion of seismograms: Journal of the CSEG, V. 15, p. 30-39. Gelland, V., and Larner, K., 1983, Seismic litholic

modeling:

presented at the 1983 convention of the CSEG,Las Vegas. Graul, M., Deconvolution and wavelet processing:

Unpubished SEG course

notes.

Hardage, R., 1986, Seismic Stratigraphy: London

Geophysical Press,

- Amsterdam.

.Hampson, D., and Galbraith, M., 1981, Wavelet extraction by sonic-log correltation:

Journal of the CSEG, v. 17, p. 24-

Hampson, D., 1986, Inverse velocity Journal of the CSEG, V. 22, p. Hampson, D., and Russell,

42.

stacking for multiple

elimination:

44-55.

B., 1985, Maximum-Likelihood seismic

inversion (abstract no. SP-16)- National CanauianCSEGmeeting, Ca.!gary, Alberta.

Part 12 - Summary

Page 12 -

9

Introducti on to Sei stoic Inversi on Methods

Bri an Russell .

Herman, A.J.,

Anania, R.M., Chun, J.H.,

Jacewitz,

C.A.,

and

Pepper, R.E.F., 1982, A fast three-dimensionalmodeling technique and fundamentals of three-dimensional frequency-Uomainmigration: Geophysics, v. 47, p.

Jones, I.F.,

1627-1644.

and Levy, S., 1987, Signal=to-noise ratio enhancementin

multi channel seismic data via the Karhunen-Loeve transform,

Geophysical Propecting, v. 35, p. Kormyl o, J.,

anu Mendel., J.M.,

deconvolution-

IEEE Trans.

12-32.

1983, Maximum-likelihood seismic

on Geoscience and Remote Sensing,

v. IT - 28, p. 482 - 488.

Lines, L.R.,

Schultz, A.K., and Treitel,

of geophysical data:

Lines, L.R.,

Geophysics, v. 53, p. 8-

and Tritel,

anU its application

S., 1988, Cooperative inversion

S.,

1984, A review of least-squares

to geophysical problems'

Prospecting, v. 32, p.

20. inversion

Geophysical

159-186.

Lindseth, R.O., 1979, Synthetic sonic logs - a process for stratigraphic interpretation:

Geophysics, v. 44, p. 3-

26.

Oldenburg, D.W., 1985, Inverse theory with applica.tion to aleconvolution and seismogram inversion.

Oldenburg, D.W., Scheuer, T.,

impedancefrom reflection Ostrander,

W.J.,

Unpublished course notes.

and Levy, S., 1983, Recovery of the acoustic

seismograms:Geophysics, v. 48, p. 1318-1337.

1984, Plane wave reflection

at non-normal angles of incidence:

Part 12 - Summary

coefficients

for gas sands

Geophysics, v. 49, p. 1637-1648.

Page 12 -

10

Introduction %oSeismic Inversion Methods

Brian Russell

Russell, B.H., and Lindseth, R.O., 1982, The information content of synthetic sonic logs - A frequencydomainapproach- presentedat the 1982 convention of the EAEG, Cannes, France.

Shuey,R.T., 1985, A simplification of the Zoeppritz equations: Geophysics,v. 50, p. 609-614. Stewart, R.R., andChiu, S.K.L., 1986, Tomography-based imagingof a heavyoil reservoir usingwell-logs, VSPand3-D SeismicdataJournal of the CSEG, v. 22, p.

73-86.

Taner, M.T., an• Koehler, F., 1981, Geophysics,v. 46, p. 17-22.

Surface

consistent

corrections-

Taylor, H.L., Banks,S.C., and McCoy,J.F., 1979, Deconvolutionwith the L1 norm: Geophysics,v. 44, p.

39-52.

Trei tel, S., and Robinson,E.A., 1966, The design of hi gh-resolution

digital filters ß IEEETransactions on Geo•cience Electronics, v. GE-4, No. 1, p.

25-38.

Walker, C., and Ulrych, T.J., 1983, Autoregressive recovery of the acoustic impedanceGeophysics,v. 48, p. 1338- 1350. 1981, Reflecti on seismology, a tool for energy resource

Waters, K.H.,

(second edition)-

exploration

Wiley, NewYork.

Western Geophysical Co., 1983, Sei smic Li thol ogic MoUeli ng:

Technical

Brochure.

Widess, M.B., 1973, How thin is a thin bed?- Geophysics, v. 38,

p. 1176-

Part

1180.

12 - Summary

Page 12 -

11

ISBN: 978-0-931830-65-5

90000

306

ISBN 978-0-931830-48-8 ISBN 978-0-931830-65-5

(Series) (Volume)

Russell - Introduction to Seismic Inversion Methods

Recommend Documents