Determination of macro subsurface models by generalised inversion

MACRO SUBSURFACE MODELS

BY GENERALIZED INVERSION

P.M. van der Made

TR diss

LIST OF CONTENTS

GLOSSARY

1 THE PROBLEM 1

1.1 Interpretation of seismic data and time to depth conversion 1

1.2 Detailed inversion of seismic data 1 1.3 Target oriented processing 2

1.4 Macro models 5 1.5 Problem definition: macro model estimation 5

1.6 Outline 5

2 POSSIBLE SOLUTIONS 7

2.1 Dix'and Hubral's methods 7 2.2 Depth reconstruction with image rays 8

2.3 Macro model estimation using generalized inversion theory 8 2.3.1 Velocity estimation algorithms using generalized inversion

theory as described in literature: are they suitable for macro

model estimation? 9

2.4 Conclusion 10

3 GENERALIZED INVERSION THEORY AND ITS APPLICATION

TO MACRO MODEL ESTIMATION 13

3.1 Generalized inversion 13

3.1.1 Formulation of the inverse problem 13

3.1.2 Optimization 14 3.1.3 Features of generalized inversion 17

3.2 Application to macro model estimation 19 3.3 Features of the objective function in macro model estimation for

two types of data vectors 20 3.3.1 The objective functions 20

3.3.2 Discussion 25 3.4 Conclusion 25

4.1.2 A priori information on the model description 34

4.2 A priori information on the data 36 4.3 Discussion of the various forms of a priori information 37

4.4 Resolution analysis 38 4.5 Conclusion 41

PARAMETERIZATION OF MACRO MODELS 43

5.1 Geometry parameterization 43 5.2 Velocity parameterization 47

5.3 Constraints 48 5.3.1 Connection constraints 49

5.3.2 Constraints to exclude crossing layer interfaces 52

5.3.3 Constraints on the velocity 53

5.4 Conclusion 53

RAYTRACING 55 6.1 Finding the ray 57

6.1.1 Raytracing by shooting 59 6.1.2 Minimum-time raytracing 59 6.1.3 Comparison of the two methods 61

6.1.4 A hybrid solution 61 6.2 Raytracing in structurally complex media 62

6.2.1 Finding the interface sequence 62

6.2.2 Multiple arrivals 63 6.3 Application to the traveltime inversion scheme: the traveltime as

MACRO S I M M MODELS

BY GENERALIZED INVERSION

PROEFSCHRIFT

ter verkrijging van de graad van doctor aan de Technische Universiteit Delft,

op gezag van de Rector Magnificus, prof.dr. J.M. Dirken,

in het openbaar te verdedigen ten overstaan van een commissie, aangewezen door het College van Dekanen op dinsdag 26 januari 1988 te 16.00 uur door

PIETER MAARTEN VAN DER MADE

geboren te Jakarta natuurkundig ingenieur

Gebotekst Zoetermeer/ 1988

TR diss

1608

Copyright © 1987. by Jason Gcosysiems, Delft, The Netherlands.

All rights reserved. No pan of this publication may be reproduced, stored in a retrieval system or transmitted in any form or by any means, electronic, mechanical, photocopying, recording or otherwise, without the prior written permission of the author, P.M. v.d. Made. Jason Gcosysiems B.V., P.O. Box 596, 2600 AN Delft, The Netherlands.

OP-DATA KONINKLUKE BIBLIOTHEEK, DEN HAAG v.d. Made. Pietcr Maarten

Determination of Macro Subsurface Models by Generalized inversion / Pieter Maarten van der Madc-IS.I. : s.n.| (Zoctcrmecr: Gcbolcksi). - 111.

Thesis Delft. - With rcf. - With summary in Dutch ISBN 90-9001995-2

SISO 567.2 UDC 550.34.01(043.3)

Subject headings: generalized inversion, seismic interpretation, seismic processing, travelling inversion. macro model, raytracing

7 APPLICATION OF TRAVELTIME INVERSION IN

PRACTICAL PROCESSING SCHEMES: OVERVIEW 67

7.1 The tracking problem 67 7.1.1 Automated pre-stack tracking 67

7.1.2 Inversion of stacking information 69 7.1.3 Recursive processing using downward continuation 69

7.2 Recursive processing versus simultaneous inversion 71

7.3 The statics problem 76 7.4 Integrated processing schemes 80

8 APPLICATION OF TRAVELTIME INVERSION IN

PRACTICAL PROCESSING SCHEMES:

INVERSION OF STACKING INFORMATION 87

8.1 Reconstruction of pre-stack traveltime from stacking information . . . . 87

8.2 Direct inversion of stacking information 88 8.3 Hyperbolic move-out assumption: implications 89

9 PRACTICAL SYNTHETIC AND FIELD DATA EXAMPLES 97 9.1 Synthetic example of the estimation of a salt dome structure 97

9.1.1 Inputdata 97 9.1.2 Initialmodel 99 9.1.3 The inversion loop 101 9.1.4 Final model and quality control 101

9.1.5 Discussion 103 9.2 North Sea field data example 105

APPENDIX A, B, C I l l

REFERENCES 120 SUMMARY 124 SAMENVATTING 125 CURRICULUM VITAE 127

GLOSSARY

The following symbols and abbreviations are used in this thesis. Unless otherwise defined in the text, they have the following meaning:

I ] : denotes column vector

I |T : denotes row vector

AT transpose of matrix A

(A) jj : denotes element i,j of matrix A, i is the row index, j is the column index

estimate of parameter (vector) p denote integers, used for indices i, j , k, 1, m, n I, J, K, L, M, N 1-D 2-D 3-D s c G O x z a CMP GSG GRG SNR To Vs t sd(x} var{x| covar(x,y) pdf

denote integers, used for ranges One dimensional Two dimensional Three dimensional slowness interval velocity geometry objective function horizontal position depth tangent

Common Midpoint (gather) Common Shot Gather Common Receiver Gather Signal to Noise Ratio post-stack traveltime stacking velocity traveltime standard deviation of x variance of x covariance of x and y probability density function

CHAPTER 1

THE PROBLEM

1.1 INTERPRETATION OF SEISMIC DATA AND TIME TO DEPTH CONVERSION

One of the goals of the seismic method is to aid the exploration and production of oil-, gas- and other mineral reservoirs in the subsurface. The idea is to acquire an image of the subsurface without (extensive) drilling.

Unprocessed seismic data do not give a useful image of the subsurface. The most common processing performed on seismic data is CMP-stacking. It reduces the data volume and enhances the signal-to-noise ratio. Moreover, the data is organized in such a way that an (almost) 1:1 image of the subsurface is obtained: the structure of the subsurface is visualized. However, this image is still in the time domain. Since the goal of seismic exploration is eventually to drill for oil or gas, an image in the depth domain is required. In practice, the complexity of this problem is underestimated: after sophisticated, and detailed interpretation in the time domain is done, very crude time to depth inversion (e.g. using Dix' formula (1955)) is carried out. Also the problem of time to depth conversion has been neglected in literature on seismic processing for a long time. Only in the past few years has it become fully recognized as an essential step that deserves closer attention.

1.2 DETAILED INVERSION OF SEISMIC DATA

Detailed time to depth conversion can be carried out in two ways: 1. Seismic depth migration

2. Inverse scattering

The goal of both processes is to invert the seismic data to estimate the detailed parameters of the subsurface. The techniques are mainly different with respect to the parameters estimated. In inverse scattering, the elastic properties of the medium (P- and S- wave velocity and mass density) is estimated, whereas in seismic migration the reflectivity distribution of the subsurface is estimated. The relation between the two approaches has been described for the linearized case by Berkhout (1984).

However, neither procedure solves the subsurface estimation problem completely. For depth migration a model of the velocity in the subsurface is needed as input. The reason is that for optimal focussing at a certain depth level, the propagation effects to that depth level have to be known. Obviously propagation depends on velocity (of propagation), so we need a velocity model beforehand.

Moreover, since seismic data is bandlimited, with both procedures a bandlimited depth model is obtained. The main problem is that the low frequency information, and consequently the trend information is missing. As a result, both procedures need an input subsurface model containing the low frequency information on the velocity and other seismic parameters (density, absorption) of the subsurface.

In state of the art seismic processing, where depth migration is common practice, the input velocity model is often obtained with very crude (and inefficient) methods. In practice, these models are not sufficiently able to provide good depth migration results, especially in complex subsurface structures. However, the required velocity model does not need to contain detailed information on the velocity field, since only the propagation effects have to be explained. As is well known, the detail in the subsurface can not be estimated from propagation information. Vice versa, it is not necessary to know the detail in the subsurface (velocity field) to explain the propagation effects of seismic waves.

1.3 TARGET ORIENTED PROCESSING

A new approach to seismic processing is the so-called target oriented processing. A first example is target oriented, detailed lithology estimation, as carried out in the Delft consortium project PRINCEPS (1982). The procedure is visualized in figure 1.1. The objective is to estimate the lithology parameters of a target zone (e.g. an oil reservoir) in detail. By stripping the seismic surface data from the propagation effects of the overburden (i.e. the zone between the surface and the target zone), undistorted data for the estimation of the target zone is obtained. The procedure to remove the propagation effects from the data is called redatuming. A second example of target oriented processing, is target migration, as carried out in the Delft consortium project TRITON (1986). In this procedure seismic data is redatumed to the target

zone, and after that migrated (see figure 1.2). This procedure is more efficient than traditional migration, because the data is only migrated over the target zone.

Redatuming, which is related to depth migration, needs a velocity model as input, for the same reasons it is needed in depth migration: propagation in the subsurface is primarily dependent on the velocity. However, as in migration, this model does not have to contain detailed information on the velocity field of the subsurface.

seismic data redatuming to target zone / macro model / of overburden detailed lithology estimation in target zone

H

Figure I.l: Target oriented, detailed lithology estimation.

seismic data

I

1

Figure I J a : The macro layers of the subsurface.

detail macro

macro layer 1 macro layer 2 macro layer 3 macro layer 4

Figure 1.3b: A vertical cross section of (I) the macro model (thick line), and (2) the band-limiled. detailed model, obtained with detailed inversion of seismic data, superposed on the macro model. The macro model contains the low frequency information, and is needed as input to the detailed inversion.

1.4 MACRO MODELS

Obviously, inversion from seismic data as described in the previous section is not sufficient to solve the total estimation problem. Both inverse scattering and seismic migration need a background medium as input. This model should contain the low frequency information that is not present in the seismic data, but does not have to contain detailed information (see figure 1.3). The model should contain only velocity information that is relevant for the propagation of seismic waves, not the reflection of seismic waves. This is what we call a macro model of the subsurface. The question is now, where do we get the low frequency information from? The only way to obtain it is, obviously, by including other information. The most practical way is to use geological knowledge on the area under investigation, combined with the interpretation of the seismic data by a geologist. However, this interpretation is done on the stacked seismic data, and consequently the information needed is not available in the depth domain. The geological model should be mapped to the depth domain to obtain the macro model of the subsurface.

1.5 PROBLEM DEFINITION: MACRO MODEL ESTIMATION

As discussed in the previous sections, there are a large number of reasons for requiring a macro model of the seismic medium parameters in the subsurface as a function of depth and lateral co-ordinate(s). Hence the problem I attempt to address in this thesis is the following: How can we get such a model?

In this dissertation, I propose to concentrate on techniques which estimate the seismic P-wave velocity only. Estimation of the S-wave velocity, anisotropy, mass density, and attenuation are not considered.

A number of different solutions to the problem will be put forward and compared. I will proceed to discuss one of these solutions, traveltime inversion, in detail. Traveltime information is sufficient to estimate the macro velocity model. The amplitude information in the data is needed for the estimation of detailed velocity models (or for density and absorption estimation). Synthetic and field data examples of macro model estimation by traveltime inversion will be presented.

1.6 OUTLINE

In chapter 1 the problem I propose to discuss is defined as follows: How can we find a macro model of the depth interval velocity structure of the earth. In the same chapter, the need of such a model is also justified: it is needed as input in seismic processing and interpretation.

In chapter 2 the possible solutions to the problem stated in chapter 1 are discussed. Both conventional and new techniques to find a macro subsurface model are mentioned.

An especially interesting class of estimation schemes is based on generalized inversion theory. In the rest of this dissertation the focus will be on these techniques. In chapters 3 to 6 the theory

on which these schemes are based is treated. To apply general inversion theory, first the problem to be solved has to be mathematically defined as a set of equations, relating the data to the parameters to estimate. In chapter 3 it is discussed how these equations can be solved. Often these solutions are non-unique. To resolve this, additional data have to be introduced. This is discussed in chapter 4. There are 2 prerequisites for defining the set of equations:

• A parametric description of the subsurface model has to be defined. In chapter 5 the parametric macro subsurface model is discussed.

• A mathematical relation between the data and the model parameters has to be defined. In the case of traveltime inversion, on which the focus of the remainder of this thesis will be, this can be done with raytracing. In chapter 6, raytracing is discussed.

In chapter 7 and 8 is shown how traveltime inversion can be applied in practical processing. The main problem is how to measure the traveltimes from seismic data. There are a number of possible ways, which are discussed in chapter 7. A procedure which has proven to be very practical, is the inversion of stacking information. It is discussed in detail in chapter 8. In chapter 9, synthetic and field data examples of the inversion of stacking information are given.

CHAPTER 2

POSSIBLE SOLUTIONS

In this chapter I intend to make an inventory of techniques that could possibly solve the problem defined in section 1, the estimation of the macro model. I do not pretend that this inventory is complete. The various techniques mentioned will be discussed with respect to:

• accuracy; • efficiency; • stability;

• how they deal with the non-uniqueness of the seismic problem; • the data (domain) they work on.

2.1 DIX' AND HUBRAL'S METHODS

In state of the art processing, the macro subsurface model is still often obtained with Dix's formula (1955). The input data to this technique are traveltimes picked from the stacked section and the stacking velocity with which they are obtained (see also appendix A and Wapenaar and Berkhout, 1985):

T * V2 - T * V2 p _ / O.i s.i O.i-1 s.i-1

T0 , i ~ O.i-1 i Z

i = S(

0,,-V-.)*

.c

estimated depth of reflector i,

estimated interval velocity of the layer above reflector i, stacked traveltime of reflector i (T0 0 = 0),

stacking velocity of reflector i.

This technique is based on a very simple model of the subsurface. The main assumption is that the earth consists of plane horizontal reflectors. A second assumption is that the pre-stack

To.i V..i

reflection traveltime is a hyperbolic function of the offset (hyperbolic move-out). This is only true for plane reflectors in a medium with a constant velocity. If the move-out curves are not hyperbolic, Dix's method works well only for small offset data. However, from small offset data little resolution for the stacking velocity can be obtained.

The technique is, however, very efficient. Only some simple formulas have to be applied to the data.

An extension to this method is the procedure proposed by Shah (1973) and Hubral (1976). This method can handle dipping reflectors. However, hyperbolic move-out is still assumed, and for non-constant velocity media only small offset data may be used. Furthermore, it is assumed that the layer interfaces are not curved. As input, in addition to the stacked traveltimes and stacking velocities, the time dip of the stacked time horizons are needed. The procedure works recursively downward into the subsurface.

In practical situations the use of time dips is a problem, as they are difficult to measure from the seismic data, and are liable to be very sensitive to noise. Moreover, since the procedure is recursive, the estimation errors may grow accumulatively as deeper reflectors are estimated. In standard seismic processing, this technique is not very commonly used.

2.2 DEPTH RECONSTRUCTION WITH IMAGE RAYS

A procedure that is often used in the industry for the construction of depth models, is the following:

1. Time migration of the stacked section;

2. Traveltime picking from the time migrated section;

3. Recursive reconstruction of the depth model using image rays.

In this procedure the interval velocities are not estimated. They have to be known a priori, and are often obtained by either using Dix's formula, or from the well-log. In this thesis I will not discuss this procedure further, since I intend to focus on procedures that estimate both depth and velocity.

2 . 3 MACRO MODEL ESTIMATION USING GENERALIZED INVERSION

THEORY

A new class of schemes for the estimation of depth and velocity (and other seismic parameters) is becoming more and more accepted in seismic industry and research. They involve generalized inversion theory or parameter estimation theory. It is assumed that a model of seismic wave propagation is known, i.e. the seismic data can be described as functions of the depth- and seismic parameters. These functions are generally non-linear. The objective is to solve the set of equations, which will result in a model of the earth in terms of depth and seismic parameters (e.g. velocity) that optimally fits the data. The procedure to do this is iterative, since the

equations are non-linear. Non-linear inversion theory is discussed in chapter 3.

As mentioned in chapter 1, it is not possible to obtain unique subsurface parameter estimates from the seismic data alone There are many solutions that fit the data. To resolve this problem extra information has to be added. This information can often be obtained from well-logs or from geological knowledge of the area. In generalized inversion theory this information can be treated explicitly as a priori information on the subsurface parameters to be estimated. In dedicated estimation schemes such as Dix's and Hubral's method, a priori information on the parameters to be estimated can not be handled. They operate on the seismic data only.

The earth is often assumed to be a sparse system of (homogeneous) layers, which is a form of a priori information as well. Usually, the sparseness assumption is based on geological knowledge. Dix's and Hubral's schemes use sparse subsurface models, since they invert for a set of horizons, picked by an interpreter, only. In generalized inversion schemes, this information can also be handled by choosing a proper model description (sparse model). This is discussed in chapter 4. The sparseness assumption is of vital importance in macro model estimation. The reason is, that it compensates for the lack of low frequencies in the seismic data. The various forms of a priori information will be discussed in detail in chapter 4.

Summarizing, generalized inversion theory is very well suited for macro model estimation. The geological a priori information can be incorporated, which means that the low frequencies missing in the seismic data can be included in the model.

2.3.1 Velocity estimation algorithms using generalized inversion theory as described in literature: are they suitable for macro model estimation?

There have already been a lot of techniques published which incorporate (non-linear) inversion theory for velocity estimation. Examples are for instance: Tarantola (1984), Bishop et al. (1985), Landa et al. (1986), Gjoystdal and Ursin (1981), Gray and Golden (1983), Chiu et al. (1986), Sword (1986). The techniques may be classified on basis of three criteria (see also table 2.1): Author Tarantola Sword Landa Bishop Gray Gjoystdal Chiu Subsurface model gridded gridded sparse gridded/sparse sparse sparse sparse Forward model wave eq. raytracing raytracing raytracing raytracing raytracing raytracing Data seismic traces seismic traces seismic traces picked trav.t. picked trav.t. picked trav.t. picked trav.t. Table 2.1: Classification of velocity estimation techniques based on generalized inversion, with respect to the parameterization of the subsurface model used, the forward model used, and the input data.

1. Kind of models estimated: detailed gridded models or sparse macro models. 2. Forward modeling used: raytracing or wave equation.

3. Data domain used: picked traveltimes or seismic traces.

The term gridded model refers to a model that consists of velocity cells, each cell with an independent velocity parameter. The velocity cells are positioned on a grid, with a small grid distance. As a result, the gridded models contain many more parameters than sparse macro models, which are parameterized on basis of geological knowledge of the area. In this section it is investigated whether the techniques described in the literature are suitable for our purpose, i.e. macro model estimation.

Tarantola's technique (1984) is meant primarily to estimate the detail in the velocity field, and works directly on the seismic traces. This technique is not suitable for our purpose, since we are interested in macro model estimation. As a matter of fact, this technique needs a macro model as input, as is discussed in chapter 1.

Techniques such as Bishop et al. (1985), Gray and Golden (1985), are based on raytracing and need picked (pre-stack) traveltimes as input. The difference between the two techniques is the definition of the subsurface model. Gray and Golden use a sparse macro model, whereas Bishop et al. use a hybrid form of a sparse model and a detailed gridded subsurface model. In the latter case the velocity field is parameterized by (square) velocity cells on an equidistant grid, but only a sparse set of reflectors is used.

Landa's technique is also based on raytracing, but it works directly on the seismic traces. The idea of this technique is to find a sparse subsurface model that gives an optimal coherency of the data along the traveltime curves generated (with raytracing) in this subsurface model (model driven coherency optimization). The disadvantage of the procedure is that it is only accurate if a reasonable estimate of the subsurface is already available. Only small refinements to this estimate can be found. This is discussed in more detail in chapter 3.

Sword's method (1986), uses traveltimes and ray parameters (T-p data) as input, and the model estimated consists of velocity cells on a grid. The procedure to pick the T-p data from the seismic traces is automatic, and is an integral part of the estimation method.

Finally, Gjoystdal and Ursin propose a hybrid technique in which iteratively (1) the velocity field within a layer is estimated with inversion theory, and (2) the layer boundaries are reconstructed (recursively) in a way similar to Hubral (1976). This means that this technique is liable to yield substantial inaccuracies for deeper layers (see section 2.1).

2.4 CONCLUSION

As stated in chapter 1, we are interested in schemes which estimate global models of the subsurface (macro models). This means that schemes which assume a sparse system of layers in the subsurface are very well suited for our purpose: A macro description of the major layers

is estimated, and the detail is neglected. The methods of Dix and Hubral estimate such models, but have the disadvantages listed above.

In recent literature, techniques that use generalized inversion for macro model estimation are described by Chiu et al., Landa et al., and Gray and Golden. They seem very promising, because they do not have the problems of Dix' and Hubral's methods. These techniques will be discussed in more detail in the next chapters. It will be shown that they still have limitations, or that they are incomplete. The main problem is that the models estimated do not include steep dips, faults and other structural irregularities.

All these techniques have in common that raytracing is used to describe the forward model, i.e. the relation between subsurface parameters and the data. With raytracing, complex structures can be handled. Moreover, raytracing is fast compared to other seismic modeling techniques. Because we are dealing with macro models (no thin layers) the assumptions made in raytracing are not severely violated.

The input data to the techniques described by Gray and Golden, and by Chiu et al. are picked traveltimes. The technique described by Landa uses the seismic traces as input. In chapter 3 the advantages of using traveltimes as input to a macro model inversion scheme rather than the seismic traces themselves are discussed. The disadvantage of using picked traveltimes as input is, obviously, that die traveltimes have to be picked from the seismic traces. How this problem can be solved is discussed in section 7.

In the next chapters, a total strategy for obtaining the macro subsurface model from seismic data is proposed. The problems addressed are:

• inversion theory (chapter 3),

• incorporation of a priori information (chapter 4), • macro model parameterization (chapter 5), • raytracing (chapter 6),

CHAPTER 3

GENERALIZED INVERSION THEORY

AND ITS APPLICATION TO

MACRO MODEL ESTIMATION

In this chapter, first the concept and main features of generalized inversion are discussed. Second, the application to macro model estimation is discussed. Two types of input data to a macro model estimation procedure are compared: seismic traces and seismic traveltimes. The seismic applications are discussed from an inversion viewpoint only.

3.1 GENERALIZED INVERSION 3.1.1 Formulation of the inverse problem

An estimation problem can be formulated as a set of equations which have to be solved (Jackson, 1972). The equations (forward model) describe the data as functions of the unknowns (model parameters). The data can be of any kind. In the seismic case, these functions are based on the seismic wave propagation in the subsurface model:

d, = d^pj + n, d. = d.(p) + n. 2 2 2 (3.1) dn = d (p) + n with: P n; data point i,

unknowns (model parameters), noise or uncertainty of data point i.

An estimation problem can be formulated as the minimization of a measure of the data mismatch (objective function). If a least squares (12) norm of the data mismatch vector is used as an objective function, it can be written as follows (Jackson, 1979):

where:

O : objective function, p : model parameter vector, dn, : measured data vector, d(p) modeled data vector

Ou : covariance matrix of the data noise (uncertainties).

Note that minimization of the 12-norm implies, in a stochastic formulation, that the noise or uncertainties of the data have a gaussian probability density function (pdf), and that the maximum likelihood solution is searched for. If another pdf of the uncertainties is assumed, a different norm will be the result. However, in practical situations, no knowledge is available on the distribution of the data noise, and the gaussian assumption is as good as any other. Moreover, the stochastic formulation is no necessity. From a practical viewpoint, the only thing one is interested in inversion is the minimization of the data mismatch, and that is precisely expressed in formula (3.2).

For a discussion of stochastic inversion, the reader is referred to Tarantola and Valette (1982a and 1982b), and Duijndam (1987).

3.1.2 Optimization

To minimize an objective function we need to search for the minimum. Algorithms which do that are called optimization or search algorithms. These algorithms can be subdivided in two categories:

1. Local search algorithms, 2. Global search algorithms.

For a detailed discussion of local optimimization the reader is referred to Scales (1985). Here only the relevant main features are mentioned.

Local optimization algorithms search iteralively for a minimum of the objective function starting from an initial estimate. Often, the local gradient and Hessian of the objective function is used to find the search direction (Newton's method). If the equations to solve are linear, and the 12-norm is used, this gives the information to find the minimum directly. In all other cases, this gives only an approximation to ihe correct search direction. A line search along the search direction has to be carried out, until the minimum along this direction is found. After that, the procedure is repeated iteratively, until the minimum of the objective function is found.

The gradient of the objective function can be calculated as follows (Scales. 1985):

g = 2 JTCd d' ( d ( p ) - d J (3.3)

where:

g : gradient of objective function,

J : Jacobian of the scaled data vector. Each row contains the derivatives of a data point to the model parameters.

The Hessian of the objective function is (Scales, 1985):

H = 2 JTC ^ j + 2 ( d ( p ) - dm) C ^ T (3.4)

where:

H : Hessian of the objective function,

T : vector containing the Hessians of the data vector elements. Tj is the Hessian of d,

If the equations d = d(p) are linear, the Hessian is (Scales, 1985):

H = 2 JTCd d' j (3.5)

For non-linear problems, equation (3.5) gives a good approximation to the true Hessian if the data mismatch is not too high (Scales, 1985). The Gauss-Newton minimization procedure, which is an adaption of the Newton method, uses equation (3.5) to approximate the Hessian. The advantage is that only the first derivatives of the data have to be calculated.

Other adaptions of the Newton method are the so-called Quasi Newton methods. In these methods the Hessian does not have to be calculated. In this procedure an initial approximation of the Hessian (often the identity matrix) is updated dynamically during the optimization process. Of course, the derivatives of d(p) have to be calculated to construct the gradient and the Jacobian. The way in which these are calculated depends, of course, on the kind of forward problem involved. In chapter 6. is described how these can be calculated for traveltime inversion.

Alternative procedures are the conjugate gradient method and the steepest descent method. These procedures only use the gradient of the objective function. The advantage of the procedures is that the Hessian matrix does not have to be stored and manipulated. This is important for systems with a large number of parameters. Tarantola (1984) uses for instance the steepest descent method in his inversion scheme, because he estimates a huge number of parameters. In macro model estimation the number of parameters is limited, and Newton's method is feasible.

The disadvantage of local search algorithms is that they only look for the minimum locally. Hence, these algorithms are liable to get stuck in a local minimum. If local minima are present in the objective function within the feasible parameter space, a global search algorithm should be used. Global search algorithms step gridwise or randomly through the (feasible) parameter space, in such a way that all local minima are detected, and the global minimum can be selected. An example of a global search algorithm is the "simulated annealing" approach, amongst others described by Kirkpatrick et al. (1983). These procedures can overcome the local minima problem. However, they are very expensive in terms of computing time. A further discussion is beyond the scope of this thesis.

Figure 3.1: The iterative behaviour of an optimization algorithm. 2 0 0 -x 400 600-800 200 400 600 800 1000 1200 1400 1600 1800 lateral distance (m)

Figure 3.1a: The true subsurface model with the data rays.

2 0 0 - 400- o-v 600-800 r 200 400 600 800 1000 1200 1400 1600 1800 lateral distonce (m)

Figure 3.1b: The initial model.

200 400 600 800 1000 1200 1400 1600 1800 lateral distance (m)

In this thesis, local search algorithms are used because the goal is to find a macro model estimation procedure that is not only feasible on super computers.

In figure 3.1 an example of the iterative behaviour of the optimization algorithms is given. It involves the estimation of a subsurface layer using seismic traveltime data. The optimization algorithm used is a Newton algorithm, i.e. a local search algorithm which uses the gradient and the Hessian of the objective function.

In figure 3.1a, the rays to generate the synthetic data are displayed. The traveltimes along these rays are used as input data (dm) to the inversion.

The iterative procedure is started with the initial model displayed in figure 3.1b. In figures c, d and e the major iterations are displayed. The final iteration yields a model that is almost equal to the "true" model.

In figure 3.1 f the value of the objective function (calculated with eq. 3.2) for the initial estimate (index 0), and for the 3 iterations (indices 1,2 and 3) is displayed.

3.1.3 Features of generalized inversion

In this section a short overview is given of the main features of generalized inversion, that are relevant for macro model estimation:

• Simultaneous inversion of different data sets: The formulation of the inverse problem in equations (3.1) and (3.2), makes it possible to invert any kind of data. Also, different types of data can be inverted simultaneously:

di i = dn<P> + nii du = du(p) + n,2 diN ■ V P ) + niN d2, = d2,(p) + n2, (3.1a) d22 = d22(p) + n22 d2M = d2M(P> + n 2 M with: dn d2i P nii " 2 i

data point i of data set 1, data point i of data set 2, unknowns (model parameters),

noise or uncertainty of data point i of data set 1, noise or uncertainty of data point i of data set 2.

Data set 1 may for instance contain seismic data, and data set 2 Vertical Seismic Profiling data. Of course, more than 2 data sets can be inverted simultaneously in the same manner.

Q. V O-i 2 0 0 - 4006 0 0 -8 0 0 - '*''" f ** i r r 1 1 1 T 200 400 600 BOO 1000 1200 1400 1600 1800 lateral distance (m)

Figure 3.Id: The second iteration.

2 0 0 c 4 0 0 -a V 6 0 0 -800 _{1 1 1 1 1 r} 200 400 600 800 1000 1200 1400 1600 1800 lateral distance (m)

Figure 3.1e: The third and final iteration. The final result is almost equal to the true model of figure 3. la.

0 1 2 3

• A priori information: Often, on an estimation problem some a priori information on the model parameters is available (e.g. from well-logs or geological knowledge). The simultaneous- inversion feature may also be used to include a priori information on the subsurface to be estimated. This is discussed in chapter 4.

• Resolution analysis: When more than one data set is inverted simultaneously, the influence of the data sets may be identified and quantified with the so called resolution matrix. This is discussed in chapter 4.

• Accuracy: In principle, any accuracy can be obtained with generalized inversion. The obtained accuracy depends on the number of iterations allowed to the optimization algorithm. • Automatic procedure: Using generalized inversion, completely automatic algorithms can be

built. Once all the data is entered, the optimization procedure is completely automatic. • Efficiency: A drawback is that the optimization procedure, certainly when the equations

(3.1) are very non linear, may take many iterations and a large amount of computation time. • Non-recursive: With generalized inversion, non-recursive estimation procedures can be

built. As discussed in chapter 2, recursive algorithms are inaccurate for deeper layers.

3.2 APPLICATION TO MACRO MODEL ESTIMATION

To apply generalized inversion in macro model estimation, a sparse representation of the subsurface should be chosen (see chapter 2, section 2.3). In chapter 4 this fact is further demonstrated. In chapter 5 is discussed how a suitable parametric representation of the macro model can be defined. Here, it is sufficient to mention that the (sparse) parametric macro model consists of 2 types of parameters:

• Geometry parameters (pG), which describe the geometry of the interfaces between the macro

layers.

• Velocity parameters (pc), which describe the velocity field within the layers. The velocity

field within one layer should be parameterized such that the velocity is either constant or very smoothly varying, which means that only low frequency variations are allowed.

Note, that this implies quite a different parameterization than usually implemented in seismic-inversion. Usually, grid parameterization is used, where the parameters are velocity (and density or other seismic parameters) on a equidistant grid, and the coordinates of the grid points are not parameters to be estimated.

If we apply 12-norm generalized inversion to macro model estimation, we simply have to put the macro model parameters in the parameter vector of equation (3.2).

P «■ [PG Pel

where:

pc : velocity parameter vector,

p : total parameter vector.

What remains to be defined, is the data vector d in equation (3.2), and the relation between data and parameters, the "forward model" d(p). First the data vector, and after that the forward model must be chosen. For the data vector, we restrict ourselves to two possibilities:

1. Seismic traveltimes, picked from the seismic traces by an interpreter. 2. The seismic traces themselves.

In the next section, the most suitable data vector is investigated. It will be concluded that seismic traveltimes are optimally suited for macro model estimation.

There are various choices open for the forward model. A thorough investigation of the various forward modeling algorithms known from literature is out of the scope of this thesis. In this thesis, raytracing is used as forward modeling scheme. Raytracing is discussed in chapter 6.

3.3 FEATURES OF THE OBJECTIVE FUNCTION IN MACRO MODEL ESTIMATION FOR TWO TYPES OF DATA VECTORS

To investigate the most suitable data vector, the behaviour of the objective function for the two types of data vectors mentioned in the previous section are discussed. In chapter 2, macro model inversion schemes from literature for both types of data vectors are mentioned:

1. Macro model estimation with seismic traveltimes as data vector elements: Chiu et al. (1986), Gray and Golden (1983),

2. Macro model estimation with seismic traces as data vector elements: Landa (1986).

In the following discussion, it is shown that using seismic traces the behaviour of the objective function as a function of the parameters is such that problems will arise in minimization with a local search algorithm. To clarify the point, inversion schemes are also discussed which use seismic traces to estimate a detailed gridded model (e.g. Tarantola, 1984). Note that these procedures are actually outside the scope of this thesis, because they are not meant for macro model estimation.

3.3.1 The objective functions

In this subsection, the objective functions for three types of inversion schemes are defined. In the next subsection, their behaviour as a function of the parameters will be discussed.

• The objective function for traveltime inversion using a macro- model parameterization:

0 = ( t ( p ) - t JTC ;i ,( t ( p ) - tm) (3.6)

where:

O : objective function,

p : macro model parameter vector (see section 3.2), i.e. a vector containing a set of parameters describing the geometry of the layers (pG), and the velocity field

tn, : measured data vector, i.e. a vector containing the picked traveltimes from a seismic section,

t(p) : modeled data vector, i.e. a vector containing the modeled traveltimes corresponding to the measured traveltimes,

C„ : covariance matrix of the data noise: The picking errors in the traveltimes. • The objective function for inversion of seismic traces using a macro-model parameterization

("model driven coherency optimization"):

o = -£(i/r),

(tO

I<

macro model parameter vector, as in formula (3.6) (see also section 3.2), amplitude recorded at receiver j from source i at time t,

modeled pre-stack traveltime from source i to receiver j , source index,

j : receiver index,

T : time gate around modeled traveltime.

It can be seen that the objective function is of a different nature. It is not the 12-norm of the data mismatch vector (dm - d(p)). Also, the data vector d(t) is not used completely in the

optimization process. Only the data samples in a time gate (-T/2.T/2) around the modeled traveltimes are used.

• The objective function for inversion of seismic traces using a detailed, gridded parametric model ("inverse scattering"):

0 = ( d ( p ) - d j C j ( d ( p ) ^ j (3 .8 )

where:

O : objective function,

p : model parameter vector. The model parameters form a gridded subsurface model. Every parameter is the velocity al a grid point of an equidistant grid (see also section 3.2). The grid distance meets the Nyquist criterion in the wavenumber domain (see chapter 4);

Figure 3.2: Local minima for 3 types of inversion procedures. The one dimensional seismic data lo invert (figure 3.2b) generated in the model depicted in figure 3.2a.

Depth — » ■

Figure 3.2a: True model.

Time — ^

-Figure 3.2b: Seismic data.

c ;

B

Figure 3.2c: Initial griddcd model fof the inverse scattering procedure. The parameters lo estimate are the reflector amplitudes at every grid point (14 parameiers). The initial values are 0.

E

h

1—£ f—I ♦ 1 1 1

5 6 f 8 % 10 11 12 13 V

3 4

Figure 3.2d: The inverse scatlering estimation result.

.t

II

Parameter

t

Depth

Figure 3.2f: The initial model for model driven coherency optimization.

1

Depth —►

Figure 3.2g: The estimation result from model driven coherency optimization.

Depth parameter

Figure 3.2h: The objective function for model driven coherency optimization as function of the depth parameter.

t

Time — ►

Figure 3.2i: The travelume picked from the seismic data in figure 3.2b.

1

Figure 3.2j: The initial model for travelume inversion.

1

Depth —>•

Figure 3.2k: The estimation result from travcltimc inversion. The initial model of figure 3.2j was used.

«♦ \ /

-—. f— \ e

Depth parameter — ►

: measured data vector. The elements are the samples of the (pre-stack) traces; : modeled data vector. The elements correspond to the samples of the (pre-stack)

traces. The data is modeled using the wave equation; : covariance matrix of the data noise (uncertainties). 3.3.2 Discussion

In model driven coherency optimization, the objective function as function of the parameters has many local minima. The cause is the so-called cycle skipping. When two traces are cross correlated, peaks of wavelets can easily give a local maximum in the correlation with side lobes of wavelets (or side lobes with side lobes, see also figure 3.2f, g, and h). Consequently, minimization procedures are liable to get stuck in one of the local minima, and thus may never find the global minimum corresponding to the best attainable subsurface estimate.

In inverse scattering, this problem is less likely to be encountered. The reason is that no macro model parameterization is used in the inversion. The model used has a sufficient number of degrees of freedom to bypass the local minima problem (see figure 3.2c, d, and e). Note, that a minimal data mismatch does not guarantee a realistic estimation result. Because the data is bandlimited, the inverse scattering procedure lacks resolution to get unique estimation results. Traveltime inversion does not have the local-minima problem either, in spite of the use of a macro model. The ambiguities caused by cycle skipping have been removed from the data by traveltime picking (see also figure 3.2i, j , k, and 1). The cycle skipping problem has now, however, been moved to the traveltime picking process. How to solve that is discussed in chapter 7.

Of course some local minima may still be present in the objective function for traveltime inversion (or inverse scattering). It is almost impossible to make any predictions about that. Our experience shows, however, that local minima are no major problem in traveltime inversion.

3.4 CONCLUSION

The 12-norm of the data mismatch seems to be a good objective function in practical situations. Local, gradient (and Hessian) based search algorithms are to be preferred over global search algorithms, if we want to design an estimation scheme that is feasible on small computers. Since we propose to use sparse macro-model parameterization, the data vector elements should be the picked seismic traveltimes when a local, gradient based search algorithm is used.

CHAPTER 4

HARD AND SOFT A PRIORI INFORMATION

As is well known, the seismic problem is non-unique in its solution. An example is given in figure 4.1. In this problem many subsurface models fit the same seismic data set. To resolve this, extra (a priori) information has to be added. This information can come from well-logs, or from geological knowledge of the area. Also, a hypothesis based on the interpretation of the seismic section of the interpreter may serve as a priori information. Often, the assumption is made that the earth is a sparse system of reflecting (homogeneous) layers. This is also a form of a priori information.

Resolving the non-uniqueness of the seismic problem is not the only reason to use a priori information. Another reason for using it, is simply that it is there: to achieve an optimal estimation result, all available information should be used. In this section, the various ways to use a priori information are discussed in general, and how they can be used in macro model estimation in particular.

As discussed in chapter 3, an estimation problem can be formulated as a set of equations which have to be solved. This is the basis of general inverse theory. The equations (forward model) describe the data as functions of the unknowns (model parameters). In seismics the functions are based on the seismic wave propagation in the subsurface model.

s, = s,(p) + n, s2 = s2(p) + n2 ( 4 1 ) Sf, = sN(p) + nN where: Si P Si(p)

seismic data point i,

unknowns (model parameters), relation describing the forward model,

noise or uncertainty of data point i. It may also include the error in the "theory", i.e. the error in the forward model Sj(p).

Figure 4.1: An example of a non-unique seismic inversion problem. o.ooo j= o.soo-J3 i .ooo " ^ 1 S O O -l f -l _ -L2_ 1OOO 2 0 0 0 3 0 0 0 • o [ « r a d l l t a n o * ( r n )

Figure 4.1a: The seismic data. In this hypothetical case only zero offset data from one horizontal reflector is used. 5 0 0 . 1 O O O -- 1 S O O 2 0 0 O H 2 S O O O SOO-. 1 O O O H 1 S O O -2 0 0 0 2 S O O 1500. 3000. m/s m/s 2000 2000 2500. 3000. m/s m/s ..Ï' ".«'■ ! (b) (d) l O O O 2 O 0 O 3 0 0 0 4-000 l a t e r a l d i i t o n c t ( r n )

In principle, a priori information can be imposed on 2 entities: 1. The subsurface model.

2. The data.

4.1 A PRIORI INFORMATION ON THE SUBSURFACE MODEL

A priori information on the subsurface model can be implemented in two ways:

1. By defining constraints on the model parameters themselves. This is the classical way of using a priori information in generalized inversion.

2. By the choice of parameterization. Although a priori information is often implicitly included this way, it is not always recognized as such.

In the next two subsections, both types of a priori information will be discussed.

4.1.1 A priori information on the model parameters

If for instance a well-log is available, this gives us direct information on the model parameters we want to estimate from the seismic measurements. When general inversion theory is used, this information can be utilized by adding it in the form of extra equations to the set of equations describing the information in the seismic data.

s, = s,(p) + nt, s2 = s2(p) + ns 2 sN = y p ) + ns N a, = a,(p) + na, (4.2) h = ^ P ) + na,2 aM= aM<P>+ ". ,M where: Si a, P n

seismic data point i, a priori data point i,

unknowns (model parameters),

noise or uncertainty of seismic data point i, noise or uncertainty of a priori data point i.

The functions a^p) describe the relation between the a priori data and the model parameters. Usually, the a priori information is specified on the parameters directly:

a. = p. + n .

i * i a,i S.I

n..i

a = I p + n where:

a : a priori parameter vector, p : parameter vector, I identity matrix,

na : vector containing the uncertainties of the a priori parameter vector.

Note that the a priori parameter vector is, in principle, not the same as the initial parameter vector mentioned in chapter 3. The initial parameter vector is merely a staning point for the estimation procedure, and does not hold information on the final estimate. However, often the initial parameter vector is chosen equal to the a priori parameter vector (if available).

As described in section 3 the total set of equations may be solved by minimizing the weighted 12-norm of the data mismatch vector.

0=(d(p)-d

)

ci(d(p)-d

)

.

where: O P

objective function, model parameter vector,

measured data vector, containing both seismic and a priori data: dm = |s, s2... sNa, a2.... aM],

d(p) : modeled data vector, containing both seismic and a priori modeled data: d(p) = (s,(p)... sN(p) a,(p).... aM(p)|,

Q j : covariance matrix of the uncertainties in both seismic and a priori data.

In most cases, it is reasonable to assume that the a priori data is stochastically independent from the seismic data.

In this case, the objective function can also be formulated with the contributions of the seismic data and a priori data separated:

O = ( s ( p ) - sm) V ( s ( p ) - sm) + T (4.4)

(a(P)-O

« ( ^ P ) - 0

measured seismic data vector: s,,, = [sj s2 .... s^],

modeled seismic data vector: s(p) = [s^p) S2(p)... sN(p)],

covariance matrix of the uncertainties in the seismic data, input a priori data vector: a,,, = |a] a2 .... a^jj,

a(p) : modeled a priori data vector: a(p) = |a,(p) a2(p)... aN(p)],

CM : covariance matrix of the uncertainties in the a priori data.

Note that minimization of the 12 norm of the a priori data implies, in the stochastic formulation, that a gaussian probability density function (pdf) of the uncertainties in the a priori data is assumed (see also chapter 3). However, this may not be an satisfactory assumption. In many situations the uncertainties in the a priori data may be with absolute certainty within a given range. If the probability within this range is uniform, this would result in a box-car pdf. Of course any hybrid form of the two types of pdf are possible. In most practical cases, the box-car pdf and the gaussian pdf are the only two types used, since no information is available to make further refinements. A priori data with boxcar pdf's are often called hard constraints. A priori data with gaussian pdf's are often called soft constraints.

Hard constraints on the model parameters can not be implemented as in equation (4.4), because this equation supports only a gaussian pdf of the a priori information. The implementation of hard constraints is usually done by using "constrained optimization". Constrained optimization is described in detail in many text books (e.g: Scales, 1985). Further discussion is outside the scope of this thesis.

In figure 4.2, an example of inversion using soft a priori information is given. The data in figure 4.1a is inverted. The solution of this problem, based on the seismic data only is non-unique. We shall assume that a priori information is available on the subsurface model to be estimated, as depicted in figure 4.2a. The a priori information is "soft" (Gaussian pdf), and the standard deviation of the a priori depth parameters are:

sd{z)=200m ,

and the standard deviation of the velocity parameter is: sd(c} =200m/s .

The standard deviations of the zero offset traveltimes are: sd{t0) = 10ms .

In figure 4.2b, the data mismatch of the a priori model is depicted. Obviously, this model does not fit the seismic data. In the inversion procedure, this model is adapted in such a way that it does fit the seismic data (see figure 4.2d), and stays as close as possible to the a priori model (see figure 4.2c), or, to be more precise, the objective function of eq. (4.4) is minimum. Using the a priori information, the subsurface model that fits the seismic data optimally and is closest to the a priori model is chosen.

A typical example of hard constraints is given in figure 4.3. The parameters of this scheme are the depths of the layers Zj. Obviously they may never cross each other. The constraint equations are:

Figure 4.2: Inversion using soft a priori information. The model to be estimated is parameterized by 2 depth parameters (at the left and right boundaries of the plot), and 1 velocity parameter (the interval velocity in the upper layer).

soo-f

1 0 0 0 2 0 0 0 3 0 0 0 l a t e r a l d l e t a n e e ( m )

Figure 4.2a: The a priori model.

£ i.c

2 . 0 0 0

-ö i o o o a o o o 3 0 0 0 lateral d l e t o n a e ( m )

Figure 4.2b: The data mismatch of the a priori model. The lower dots (connected by the horizontal line) are the input data. The upper doLs are the modeled data. The vertical lines denote the data mismatches for each zero-offset data point.

1800. 3000. m / s m / s 9 O 0 -1 OOO1 5 0 0 -2 0 0 0 1OOO 2 0 0 0 3 O 0 0 * O 0 0 l a t e r a l d l i t a n c * ( r n )

Figure 4.2c: The final estimate, after simultaneous inversion of the seismic data and a priori data.

B x& , ' n | o . s o , £ 1 . O O O Sfi 1 . O O O -2 0 0 0 J O O O A O O O l a t e r a l d i s t a n c e ( m )

Figure 4.2d: The data mismatch of the final estimate. Now the modeled data arc equal to the input data, so no data mismatches can be seen anymore.

Figure 4.3a: The parametric model. The parameters are:

Z]: the depth of the first interface; 7^. the depth of the second interface.

probability Figure 4.3b: The 2 dimensional pdf. The shaded area has the probability 0, and is excluded as inversion result. The white area has uniform, larger than 0, probability.

0 < z . , - z . <°°, for i = 1,2

The probability that Zj+j - z, is within this range, is of course uniform, resulting in a boxcar

pdf.

4.1.2 A priori information on the model description

There are many ways to parameterize the subsurface model. As already mentioned in previous chapters, there are two major categories:

1. Gridded, detailed (micro) models. 2. Sparse, macro models.

In detailed models, the subsurface is usually parameterized by velocity cells, each cell with an independent velocity parameter (see figure 4.4). In the 2-D case, often square cells are taken. The dimensions of the cells are usually small compared to the wavelength of the seismic data:

Ax<c . / ( 4 * f ) (4 5%

min v m a x ' 1'♦•-»/

where: Ax

'max

width of velocity cell,

maximum frequency in seismic data, minimum velocity in subsurface.

This is the well known anti aliasing criterion for depth migration operators for primary reflection data.

Formula (4.5) implies that detailed models have the degrees of freedom to Tit any (seismic) data set with a maximum frequency fmax. This means that no extra information is added to the

seismic data by using such a model description. If the grid distance is,

A x > cm„ / ( 4 * fm i n) ( 4.6 )

where:

fmjn : minimum frequency in seismic data, cmax : maximum velocity in subsurface,

the model parameterization itself adds low frequency information to the solution outside the data bandwidth, if normal incidence reflection data is used. For non-normal incidence data a more complex criterion for the grid distance applies. The grid criterion can be determined for these cases by investigating kx- k, (x- and z- component of the wave vector, see e.g. Berkhout, 1987)

representation of the data for the model under investigation. Using such a representation, the lowest wave numbers kx and kz can be identified, and the corresponding maximum Ax and Az

(horizontal and vertical grid distance) can be derived, which add no a priori information. Note that inversion of seismic reflection data, using an equidistant grid with a grid-distance for which inequality (4.6) holds, will certainly yield aliasing, and the inversion results will be unreliable.

c, C2 C3 C4

Cl - 1 Cl CC1 Cl . 2

Figure 4.4: A gridded velocity model. Each cell has an independent velocity parameter, ci is the velocity in cell i.

In macro model estimation, it is assumed that the subsurface consists of a sparse set of layers, and that the layers have a constant or smoothly varying velocity. Usually these assumptions are based on geological knowledge. We can parameterize a model such that deviations from that assumption are not allowed (see chapter 5). Here, it is sufficient to describe the main features of such a parameterization. As already mentioned in chapter 3, the parameteric model consists of 2 types of parameters (see figure 4.5):

• Geometry parameters (pG), which describe the geometry of the interfaces between the macro

layers.

• Velocity parameters (pc), which describe the velocity field Cj(x,z) within the layers. The

velocity field within one layer should be parameterized such that the velocity is either constant or very smoothly varying, which means that only low frequency variations are allowed.

Usually, the layers have dimensions for which inequality (4.6) certainly holds. Consequently. low frequency information is added to the data outside the data bandwidth (if normal incidence reflection data is used). However, the information is introduced in a sensible way, because the parameterization is based on geological a priori knowledge, instead of using rectangular velocity cells on an equidistant grid.

By the choice of parameterization, hard constraints are imposed on the solution, to the extent that the solution is sparse. The "sparsity constraint" is essential in macro model estimation, since it gives the low frequency information on the subsurface model that is not present in the seismic data. Consequently, using sparse macro models in inversion schemes, more resolution for the model parameters is obtained than with detailed models. For the 1-D case, this is also shown in theory and with examples by van Riel and Berkhout (1985).

4.2 A PRIORI INFORMATION ON THE DATA

The seismic data are usually contaminated with a considerable amount of random and correlated noise. Obviously, the quality of the subsurface estimate depends on the signal to noise ratio of the data. The influence of the noise on the estimates can be reduced by interpretation of the data, as is shown for two types of inversion schemes:

1. Inversion of seismic traces, 2. Inversion of seismic traveltimes.

Consider an inversion scheme, utilizing an objective function of the type found in equations (4.4), and in which the data vector consists of the samples of the seismic traces. As can be seen in equation (4.4), the objective function is inversely weighted by the (co)variances of the data. For every sample, the interpreter may define the degree to which it is corrupted by noise, and assign it an assorted (co)variance. Samples wich are interpreted as having a low signal to noise ratio should be assigned a high (co)variance and vice versa. Since the samples are inversely

weighted with the covariance matrix, samples with a low covariance will be of greater influence on the subsurface estimate than samples with a high covariance. Consequently, the subsurface estimate is less influenced by the data noise. Obviously, assigning (co)variances to all samples of seismic traces is a enormous task, rendering this approach to a large extent impractical. Picking seismic traveltimes is an extreme form of the classification of seismic data as described above. Only the peaks of the major reflections are selected from the seismic traces and their traveltimes are used as data in a traveltime inversion scheme. All other samples of the seismic traces remain unused or are assigned, as it were, an infinite variance.

Note that by picking the peaks of the seismic events, the sidelobes are ignored. Consequently, the frequency content of the data is enhanced by the interpretation of the data, or, in other words, the interpretation adds a priori information to the data.

4.3 DISCUSSION OF THE VARIOUS FORMS OF A PRIORI INFORMATION

Summarizing, a priori information can both be imposed on the model parameters and on the data. A priori information on the model parameters can be introduced in two ways:

• By imposing constraints on the model parameters,

• By the choice of parameterization. There are two important classes of parameterization: Macro-models and detailed models.

A priori information on the data can be introduced by picking the seismic events. Since this results in a sparse representation of the data, it is the equivalent of the macro model for the subsurface parameters in the data domain. It seems logical to use sparse data when a sparse model is used. So, picked seismic data seem to be the appropriate data for macro model estimation.

In table 4.1, the table presented in chapter 3 (table 3.1) is rearranged to show how a priori information is used in techniques described in literature.

sparse data (picked traveltimes) sampled data (seismic traces)

sparse model sampled model

(macro model) (detailed model)

Gray, Chiu, (Bishop) Gjoystdal, (Bishop)

Landa Tarantola, Sword

Table 4.1: Classification of velocity estimation techniques based on generalized inversion theory. The classification is made with respect to the a priori sparseness assumption imposed on the data and the subsurface model.

It can be seen that most techniques for macro model estimation indeed use picked traveltimes as input. There are some exceptions: The technique described by Landa et al. (1986), uses sampled

data, although a sparse macro model description for the subsurface is used. The reason is that, due to the low signal to noise ratio of seismic data, the interpretation and event picking is often a difficult process (see chapter 7). However, to construct an initial macro subsurface model, an interpretation of the data is unavoidable. Of course, this interpretation does not have to be accurate, since it involves only an initial model. Landa's technique can best be used for a final update of a depth model obtained with traveltime inversion, when the picked traveltimes are inaccurate.

Bishop's technique, which utilizes picked traveltimes as input, uses a hybrid model description. The velocity field consists of velocity cells on a grid (gridded model approach), but only a sparse set of reflectors is used (sparse macro model approach).

4.4 RESOLUTION ANALYSIS

If a priori information is added to an estimation problem, it may be desirable to be able to afterwards identify which features of the estimated model are determined from the a priori data, and which from the seismic data. This is called resolution analysis. If the a priori information is added in the form described in section 4.1.1, i.e. by adding extra equations, the resolution analysis may be performed using the resolution matrix (Jackson, 1979):

where:

R : resolution matrix,

Js : Jacobian of the seismic data vector,

Ja : Jacobian of the a priori data vector.

This matrix follows from an linearized model of the data as function of the parameters. The linearization is performed at a given realization of the parameter vector p0, e.g. the final

estimate. The linearized forward equations are, disregarding the additive noise: s - s(p0) = Js (p - p,j) a-a(p0) = Ja( p - p0) or: As = Js Ap (4.7a) Aa = Ja Ap

The least squares estimator operating on the seismic data only (Jackson, 1979) is:

A P e - P X V ' K V

$<£*•> (4.7b)

where Ape is the estimated parameter perturbation vector.

A

p

= (

I < J

J J c ^ J,)-

(jj c;

J

) A

or:

Ap =RAp

(4.7c)

where R is the resolution matrix of equation (4.7).

The meaning of the resolution matrix can be understood as follows: This matrix maps, in a linear sense, the unknown "true" subsurface parameters on the estimated subsurface parameters. Ideally, this matrix should be equal to the identity matrix, which means that unknown "true" subsurface parameters are mapped on the estimated parameters by the seismic data in an undistorted and unique sense. If the matrix is not equal to the identity matrix, the parameters are not determined uniquely by the seismic data only. For instance, if row i of R contains zero elements only, parameter i is completely unresolved by the seismic data. Consequently, this parameter is determined completely by the a priori data. A more detailed discussion on the interpretation of resolution matrices is given by van Riel (1982), and van der Made (1984a). An example is given in figure 4.6. The resolution matrix is calculated for the inversion problem

Figure 4.6: The resolution matrix for the estimation problem of figure 4.2. None of the diagonal elements are equal to 1, so none of the parameters are fully resolved. The values of the off diagonal elements are a measure for the cross-correlation of the parameters.