# Data modeling by inversion

Series | Investigations in Geophysics |
---|---|

Author | Öz Yilmaz |

DOI | http://dx.doi.org/10.1190/1.9781560801580 |

ISBN | ISBN 978-1-56080-094-1 |

Store | SEG Online Store |

## J.1 The generalized linear inversion

We shall start by making the following formal statement:

*Given an observed data set* **d**, *estimate a set of parameters* **p** *which are used to construct a model* **d′** *of the observed data set* **d**, *such that the difference between the observed data set* **d** *and the modeled data set* **d′** *is minimum based on a specific norm*.

For most geophysical applications, the norm for minimization usually is chosen to be *L*_{2}, in which case the error energy in data modeling is made minimum in the least-squares sense. Another norm that is appropriate in some applications is *L*_{1}, in which case the magnitude of the error in data modeling is made minimum.

To form a mathematical framework for the objective stated above, we need a *model equation* that relates the modeled data with the model parameters to be estimated as

**(**)

where **d′** is the modeled data vector, **p** is the model parameter vector, and **L** is the matrix that relates the modeled data vector to the model parameter vector.

The error vector **e** is defined as the difference between the modeled data vector and the observed data vector

**(**)

Substitute equation (**J-1**) into equation (**J-2**) to obtain

**(**)

Following Lines and Treitel ^{[1]}, the least-squares solution for equation (**J-3**) can be determined. First, the cumulative squared error *S* is expressed as

**(**)

where *T* is for transpose. By substituting for **e** from equation (**J-3**), we get

**(**)

Expand the right-hand side to get

**(**)

Now differentiate both sides of equation (**J-4c**) with respect to **p** and observe the requirement for least-squares minimization that *∂S*/*∂***p** = 0

**(**)

Apply matrix transpose and rearrange the terms

**(**)

which yields the desired least-squares solution

**(**)

where **L ^{T}L** is the covariance matrix and (

**L**)

^{T}L**is the least-squares (also called generalized linear) inverse of**

^{−1}L^{T}**L**.

Equation (**J-6a**) represents the generalized linear inverse (GLI) solution to the parameter vector **p**. This solution is widely used in many stages of seismic data analysis. Examples include deconvolution (optimum wiener filters), residual statics corrections (residual statics corrections), refraction statics corrections (refraction statics corrections), and the discrete Radon transform (the radon transform). In most applications, the solution needs to be constrained since the covariance matrix **L ^{T}L** needs to be inverted. The constrained solution is given by

**(**)

where *β* is called the damping factor and **I** is the identity matrix.

In some applications, the generalized linear inverse problem is formulated in the frequency domain. Then, the unconstrained solution is given by

**(**)

and the constrained solution is given by

**(**)

where the asterisk denotes complex conjugate, and **p**, **d**, and **L** are complex.

The computationally suitable technique to solve for the parameter vector **p** is dictated by the properties of the coefficient matrix **L**; that in turn is influenced by how the parameter vector **p** is defined. The coefficient matrix **L** may be sparse with many of its elements being zero as for the residual and refraction statics. It may be near singular with the elements of two rows being very similar in size as for the discrete Radon transform. The covariance matrix **L ^{T}L** may have the Toeplitz structure with its elements being diagonally symmetric as for deconvolution and the discrete Radon transform. It may also be diagonally dominant as for deconvolution. In geophysical applications, techniques to solve for the parameter vector

**p**in equations (

**J-6a**,

**J-6b**) or (

**J-7a**,

**J-7b**) include Levinson recursion, conjugate gradient, Gauss-Seidel and singular-value decomposition.

## J.2 The GLI formalism of deconvolution

Deconvolution is fundamentally a data modeling technique. Specifically, we model a 1-D seismogram for a minimum-phase estimate of the source wavelet, to predict multiples, and ultimately obtain an estimate of white reflectivity series (optimum wiener filters). The mathematical foundation of deconvolution is provided in Appendix B. Here, we shall apply the theory of the generalized linear inversion outlined in the previous section to describe deconvolution within the framework of 1-D seismic inversion. Inversion techniques at higher dimensions will be summarized in the next section.

We shall consider designing a least-squares inverse filter *f*(*t*) that converts a wavelet *w*(*t*) to a desired form *d*(*t*) such that the difference *e*(*t*) between the actual output *y*(*t*) and the desired output *d*(*t*) is minimum in the least-squares sense. The zero-delay unit spike is a special case of the desired output *d*(*t*). Other forms of *d*(*t*) can also be considered, such as a zero-phase band-limited wavelet. The model equation for deconvolution is given by

**(**)

Consider the discrete form of equation (**J-8**), with *w*(*t*) represented by the *m*–length time series (*w*_{0}, *w*_{1}, *w*_{2}, …, *w*_{m−1}), and *f*(*t*) represented by the *n*–length time series (*f*_{0}, *f*_{1}, *f*_{2}, …, *f*_{n−1}).

Equation (**J-8**) can then be expressed in matrix form

**(**)

This is the equation of complete transient convolution. Define the output vector on the left-hand side by **y**, the coefficient matrix on the right-hand side by **L**, and the filter vector by **f**. Equation (**J-9**) takes the compact form

**(**)

Follow the steps from equations (**J-2**) through (**J-6**) to obtain the least-squares solution for the filter vector **f**:

**(**)

Consider the special case of a three-point wavelet (*w*_{0}, *w*_{1}, *w*_{2}). Set up the **L** matrix of equation (**J-9**) for this special case

**(**)

with its transpose **L ^{T}** given by

**(**)

Now multiply the two matrices to get

**(**)

Compute the first three lags of the autocorrelation (*r*_{0}, *r*_{1}, *r*_{2}) of the wavelet (*w*_{0}, *w*_{1}, *w*_{2}):

**(**)

and note that the elements of the covariance matrix **L ^{T}L** given by equation (

**J-12c**) are the first three autocorrelation lags of the wavelet (

*w*

_{0},

*w*

_{1},

*w*

_{2}) given by equations (

**J-13**). Thus, for the general case, we note that

**(**)

where (*r*_{0}, *r*_{1}, *r*_{2}, …, *r*_{n−1}) are the first *n* autocorrelation lags of the input wavelet series (*w*_{0}, *w*_{1}, *w*_{2}, *w*_{m−1}).

For the special case of a desired output vector **d**

**(**)

obtain the matrix product **L ^{T}d**

**(**)

Now compute the first three lags of the crosscorrelation (*g*_{0}, *g*_{1}, *g*_{2}) of the desired output (*d*_{0}, *d*_{1}, *d*_{2}, *d*_{3}, *d*_{4}) with the input wavelet (*w*_{0}, *w*_{1}, *w*_{2})

**(**)

**(**)

**(**)

and note that the elements of the matrix **L ^{T}d** given by equation (

**J-15b**) are the first three lags of the crosscorrelation of the desired output (

*d*

_{0},

*d*

_{1},

*d*

_{2},

*d*

_{3},

*d*

_{4}) with the input wavelet (

*w*

_{0},

*w*

_{1},

*w*

_{2}) given by equations (

**J-16a**,

**J-16b**,

**J-16c**). Thus, for the general case, we note that

**(**)

where (*g*_{0}, *g*_{1}, *g*_{2}, …, *g*_{n−1}) are the first *n* lags of the crosscorrelation of the desired output **d** with the input wavelet **w**.

By substituting equations (**J-14**) and (**J-17**) into equation (**J-11**), we get

**(**)

This is the discrete form of the classic Wiener-Hopf integral equation to estimate the least-squares shaping filter **f** that converts an input wavelet **w** into a desired form **d**.

In general, the autocorrelation matrix **L ^{T}L** in equation (

**J-11**) is diagonally dominant since

*r*

_{0}is the largest value of the autocorrelation lags (

*r*

_{0},

*r*

_{1},

*r*

_{2}, …,

*r*

_{n−1}). Nevertheless, practical implementations of equation (

**J-11**) often require adding a small fraction

*ε*of the zero-lag of the autocorrelation to the diagonal elements of the matrix

**L**:

^{T}L

**(**)

where *ε* is called the prewhitening factor and **I** is the identity matrix.

Note that the autocorrelation matrix **L ^{T}L** given by equation (

**J-14**) is of Toeplitz form. A typical length for the deconvolution filter

**f**in equation (

**J-18**) is between 40 and 80 samples. This makes the size of the autocorrelation matrix

**L**to be between 40 × 40 and 80 × 80. If the autocorrelogram is computed from an input seismogram represented by a single trace, then the length of the input data vector is typically 1000 samples. The ratio of the length of the input seismogram used in computing the autocorrelation lags to the filter length should be no less than 8.

^{T}LEquation (**J-11**) may be solved for the filter vector **f** using Levinson recursion (Section B.6). It may also be solved by the conjugate gradient method ^{[2]} ^{[3]} ^{[4]}. The conjugate gradient method by Hestenes ^{[5]} described and exemplified by Wang and Treitel ^{[2]} is outlined below.

Refer to equation (**J-11**) and define

**(**)

and

**(**)

so that

**(**)

Our goal is to estimate the filter vector **f** recursively starting with an initial estimate **X _{0}** which may be defined as a null vector. The recursive estimate is terminated when the residual vector

**R**defined by

_{i}

**(**)

becomes a null vector itself. Here, **X _{i}** is the parameter vector estimate after the

*i*th iteration. If the autocorrelation matrix

**A**has dimensions

*n*×

*n*, the conjugate gradient method yields the solution

**f**after

*m*<

*n*iterations.

The recursive scheme requires the following initial values

**(**)

**(**)

and

**(**)

The recursive scheme computes the energy of the residual after the *i*th iteration

**(**)

To begin the recursion, set *i* = 0 and compute the residual error *c*_{0} using equation (**J-24a**). Then, insert the values for *c*_{−1} from equation (**J-23a**) and *c*_{0} from equation (**J-24a**) into ^{[2]}

**(**)

to get the value for *b*_{−1}. Next insert the values for **P**_{−1} from equation (**J-23b**), **R**_{0} from equation (**J-23c**), and *b*_{−1} from equation (**J-24b**) into

**(**)

to get a value for **P**_{0}. Then, apply the recursions given by

**(**)

**(**)

and

**(**)

to get an estimate **X**_{1} of the filter vector **f** given by

**(**)

and the associated new residual vector **R**_{1} given by

**(**)

Return to the beginning of the recursion and repeat the computations given by equations (**J-24a**) through (**J-24h**). Stop the iteration after a specified value of *m* < *n*.

I have tested this algorithm and compared the results with those obtained from Levinson recursion. For a filter length *n*, conjugate gradient gives the exact solution derived from Levinson recursion after *m* = *n* iterations. Nevertheless, an approximate solution that is acceptable without much compromise can be achieved by *m* < *n*. This means that for a large filter length *n*, conjugate gradient can be more efficient than Levinson recursion.

## J.3 Applications of the GLI technique

In this section, we shall review selected applications of generalized linear inversion (GLI) for data modeling as part of a seismic data processing sequence. Specifically, we shall compile a summary of the theoretical discussions in previous appendixes on deconvolution (Sections B.5 and B.8), residual statics corrections (Section C.4), refraction statics corrections (Section C.8), and the discrete Radon transform (Section F.3) within the context of the GLI solution given by equations (**J-6a**,**J-6b**) or (**J-7a**, **J-7b**).

For each application, we shall make reference to the model equation, input data vector **d**, the coefficient matrix **L**, the parameter vector **p** and their specific sizes, and the technique to solve for the parameter vector.

The model equation for *surface-consistent deconvolution* expressed in the frequency domain is given by (Section B.8):

**(**)

where is the logarithm of the amplitude spectrum of the modeled seismogram are the log-amplitude spectra of the individual components in the time domain — *s _{j}*(

*t*),

*g*(

_{i}*t*),

*h*(

_{l}*t*), and

*e*(

_{k}*t*), respectively. The term

*s*(

_{j}*t*) is the waveform component associated with source location

*j*, the term

*g*(

_{i}*t*) is the component associated with receiver location

*i*, and the term

*h*(

_{l}*t*) is the component associated with offset dependency of the waveform defined for each offset index

*l*= |

*i − j*|. The fourth component

*e*(

_{k}*t*) represents the earth’s impulse response at the source-receiver midpoint location,

*k*= (

*i*+

*j*)/2.

Consider a data set with *n _{s}* shot locations and

*n*channels, so that the total number of traces is

_{c}*n*×

_{s}*n*. For

_{c}*n*×

_{s}*n*values of the actual spectral components at frequency

_{c}*ω*, and

*n*shot locations,

_{s}*n*receiver locations,

_{r}*n*midpoint locations, and

_{e}*n*offsets, we have the following set of model equations:

_{h}

**(**)

Write equation (**J-26**) in matrix notation, as in equation (**J-1**). Here, is the column vector of (*n _{s}* ×

*n*)-length on the left-hand side in equation (

_{c}**J-26**),

**L**is the sparse matrix with dimensions (

*n*×

_{s}*n*) × (

_{c}*n*+

_{s}*n*+

_{h}*n*+

_{e}*n*), and

_{r}**p**is the column vector of (

*n*+

_{s}*n*+

_{h}*n*+

_{e}*n*)-length on the right-hand side of the same equation. Except the four elements in each row, the

_{r}**L**matrix contains zeros. Consider an example of

*n*= 1000,

_{s}*n*= 240,

_{c}*n*= 60,

_{h}*n*= 2000, and

_{e}*n*= 1000. You would then have 240 000 equations to estimate 4060 parameters — more equations than unknowns.

_{r}We want to estimate for each frequency *ω* the model parameters **p** such that the difference between the actual spectral component and the modeled spectral component is minimum in the least-squares sense. Follow the steps from equation (**J-2**) to (**J-7a**) for the complex case and derive the GLI solution for surface-consistent deconvolution, The parameter vector **p** that contains the spectral components which are associated with the source and receiver locations, offset dependency, and earth’s impulse response, is solved for each frequency component *ω*. Results from all frequency components are then combined to obtain the terms in equation (**J-25**). The surface-consistent spiking deconvolution operator to be applied to each trace in the data set is then the minimum-phase inverse of *s _{j}*(

*t*) *

*g*(

_{i}*t*) *

*h*(

_{l}*t*).

The size of the matrix **L ^{T}L** is (

*n*+

_{s}*n*+

_{h}*n*+

_{e}*n*) × (

_{r}*n*+

_{s}*n*+

_{h}*n*+

_{e}*n*). Consider 500 frequency components for which equation (

_{r}**J-28**) is to be solved for. For the example specified above, you would have to solve a matrix equation of the size 4060 × 4060 for each of the 500 frequency components. A practical scheme for solving this equation is based on the Gauss-Seidel iterative scheme described in Section B.8.

To estimate *surface-consistent shot and receiver residual statics*, we model traveltime deviations on moveout-corrected CMP gathers. The model equation for residual statics estimation is given by (Section C.4)

**(**)

where are the modeled traveltime deviations associated with a reflection event on moveout-corrected CMP gathers, *s _{j}* is the residual statics shift at the

*j*th source location,

*r*is the residual statics shift at the

_{i}*i*th receiver location,

*G*is the structure term at the

_{k}*k*th midpoint location, and is the residual moveout at the

*k*th midpoint location.

Consider a data set with *n _{s}* shot locations and

*n*channels, so that the maximum number of picks is

_{c}*n*×

_{s}*n*. For

_{c}*n*×

_{s}*n*picks of

_{c}*t*, and

_{ij}*n*shot locations,

_{s}*n*receiver locations,

_{r}*n*midpoint locations, we have the following set of equations:

_{G}

**(**)

which in compact matrix notation is given by **t′** = **Lp** as in equation (**J-1**). Here **t′** is the column vector of (*n _{s}* ×

*n*)-length in equation (

_{c}**J-28**),

**L**is the sparse matrix in the same equation with dimensions (

*n*×

_{s}*n*) × (

_{c}*n*+

_{s}*n*+

_{r}*n*+

_{G}*n*), and

_{G}**p**is the column vector of (

*n*+

_{s}*n*+

_{r}*n*+

_{G}*n*)-length on the right-hand side of the same equation. Except for the three elements in each row, the

_{G}**L**matrix contains zeros. Consider an example of

*n*= 1000,

_{s}*n*= 240,

_{c}*n*= 1000, and

_{r}*n*= 2000. You would then have 240 000 equations to estimate 6000 parameters.

_{G}Follow the steps from equation (**J-2**) to (**J-6**) and derive the GLI solution for surface-consistent residual statics, **p** = (**L ^{T}L**)

**. Here**

^{−1}L^{T}t**t**denotes the column vector of (

*n*×

_{s}*n*)-length that represents the traveltime deviations picked from moveout-corrected CMP gathers. The size of the matrix

_{c}**L**is (

^{T}L*n*+

_{s}*n*+

_{r}*n*+

_{G}*n*) × (

_{G}*n*+

_{s}*n*+

_{r}*n*+

_{G}*n*). For the example specified above, you would have to solve a matrix equation

_{G}**p**= (

**L**)

^{T}L**of the size 6000 × 6000. A practical scheme for solving this equation is based on the Gauss-Seidel method (Section C.4).**

^{−1}L^{T}tTo estimate *surface-consistent shot and receiver intercept time anomalies*, and thus obtain shot and receiver refraction statics, we model refracted arrival times derived from picking first breaks (refraction statics corrections). The model equation for the refracted arrivals is given by (Section C.8)

**(**)

where *T _{j}* and

*T*are the intercept time anomalies at shot and receiver locations, respectively, and

_{i}*s*is the bedrock slowness. Here, we shall consider the special case of surface shots. Hence, for

_{b}*n*shot/receiver stations the parameter vector is

**p**: (

*T*

_{1},

*T*

_{2}, …,

*T*;

_{n}*s*).

_{b}Consider a data set with *n* shot locations and *n _{c}* channels, so that the total number of traces is

*n*×

*n*. Assume that you will pick the first breaks from half the number of the traces. For

_{c}*n*×

*n*/2 picks of

_{c}*t*, and

_{ij}*n*+ 1 parameters

**p**: (

*T*

_{1},

*T*

_{2}, …,

*T*;

_{n}*s*), we have the following set of equations:

_{b}

**(**)

which in matrix notation is **t′** = **Lp** as in equation (**J-1**). Here **t′** is the column vector of (*n* × *n _{c}*/2)-length in equation (

**J-30**),

**L**is the sparse matrix in the same equation with dimensions (

*n*×

*n*/2) × (

_{c}*n*+ 1), and

**p**is the column vector of (

*n*+ 1)-length on the right-hand side of the same equation. Except for the three elements in each row, the

**L**matrix contains zeros. Consider an example of

*n*= 1000 and

*n*= 240. You would then have 120 000 equations to estimate 1001 parameters.

_{c}Follow the steps from equation (**J-2**) to (**J-6**) and derive the GLI solution for surface-consistent refraction statics, **p** = (**L ^{T}L**)

**. Here**

^{−1}L^{T}t**t**denotes the column vector of (

*n*×

*n*/2)-length that represents the observed (picked) refracted arrival times in equation (

_{c}**J-30**). The size of the matrix

**L**is (

^{T}L*n*+ 1) × (

*n*+ 1). For the example specified above, you would have to solve a matrix equation

**p**= (

**L**)

^{T}L**of the size 1001 × 1001. A practical scheme for solving this equation is based on the Gauss-Seidel method (Section C.4).**

^{−1}L^{T}tIn some applications, instead of directly modeling the observed data, such as refracted arrivals, it may be preferable to perturb the parameter vector **p** by a small amount **Δp** and compute the change in the model vector **Δt′**. Perturbation of the parameter vector associated with the earth is the basis for tomographic inversion. Since it is the change in the parameter vector and not the data vector that we estimate, tomographic application of the GLI technique may be appropriately considered an example of earth modeling rather than data modeling. In Section J.6, we review reflection traveltime tomography for updating earth model parameters. Here, we will refer to the model equation for refraction tomography (Section C.9)

**(**)

where is the amount of change in the difference between the observed traveltimes *t _{ij}* and the initial estimate of the modeled traveltimes as a result of the change in the parameter Δ

*p*; Δ

*s*and Δ

_{wj}*s*are the slowness perturbations within the weathering layer at shot and receiver locations, respectively; and Δ

_{wi}*s*is the slowness perturbation within the bedrock layer. The coefficient terms

_{b}*Z*and

_{j}*Z*and

_{i}*X*are functions of the weathering layer thickness, critical angle of refraction at shot and receiver locations and shot-receiver separation (Section C.9).

_{ij}Examine the structure of equation (**J-31**) and note that, instead of modeling refraction traveltimes by way of equation (**J-29**), we model the change in traveltimes by way of equation (**J-31**) and thus estimate the near-surface parameters. Hence, for *n* shot/receiver stations, the parameter vector is **Δp** : (Δ*s*_{w1}, Δ*s*_{w2}, …, Δ*s _{wn}*, Δ

*s*). We refer to the scheme based on equation (

_{b}**J-31**) as the

*iterative*GLI solution.

Consider a data set with *n* shot locations and *n _{c}* channels, so that the total number of traces is

*n*×

*n*. Again, assume that you will pick the first breaks from half the number of the traces. For

_{c}*n*×

*n*/2 picks of

_{c}*t*and

_{ij}*n*+ 1 parameters

**Δp**: (Δ

*s*

_{w1}, Δ

*s*

_{w2}, …, Δ

*s*, Δ

_{wn}*s*), we have the following set of equations:

_{b}

**(**)

We write these equations in matrix notation as in equation (**J-1**) as

**(**)

where **Δt′** is the column vector of (*n* + *n _{c}*/2)-length on the left-hand side of equation (

**J-32**),

**L**is the sparse matrix in the same equation with dimensions (

*n*+

*n*/2) × (

_{c}*n*+ 1), and

**p**is the column vector of (

*n*+ 1)-length on the right-hand side of equation (

**J-32**). Except for the three elements in each row, the

**L**matrix contains zeros. For the example of

*n*= 1000 and

*n*= 240, you would have 120 000 equations to estimate 1001 parameters.

_{c}The GLI solution to equation (**J-34**) satisfies the requirement that the energy of the error vector **e** = **Δt** − **Δt′** is minimum and is given by

**(**)

where **Δt** denotes the column vector of *n* + *n _{c}*/2-length that represents the difference between the observed (picked) refracted arrival times and the initial estimate of the modeled times in equation (

**J-32**). The size of the matrix

**L**is (

^{T}L*n*+ 1) × (

*n*+ 1). For the example specified above, you would have to solve a matrix equation

**p**= (

**L**)

^{T}L**of the size 1001 × 1001. A practical scheme for solving this equation is based on the Gauss-Seidel method (Section C.4).**

^{−1}L^{T}tWe now review *data modeling by the discrete Radon transform*. Let *d*(*h, t*) represent a CMP gather in the offset space (the plane of offset versus two-way traveltime) and *u*(*v, τ*) be the transformed data in the velocity space (the plane of stacking velocity versus two-way zero-offset time). The model equation for the hyperbolic Radon transform in the Fourier transform domain is given by (Section F.3)

**(**)

where *ω′* is the Fourier dual of *t′*, with *t′* = *t*^{2}, *t* being the two-way traveltime. Also, *v* is the hyperbolic moveout velocity and *h* is the (half) offset.

Equation (**J-33**) can be written in the matrix form **d′** = **Lu** as in equation (**J-1**) for each *ω′* component of *d′*(*h, ω′*) and *u*(*v, ω′*), where **L** now is a complex matrix of dimensions *m* × *p*:

**(**)

and **d′** and *u* are complex vectors of lengths *m* and *p*, respectively, where *m* and *p* are the number of offsets and the number of velocities used in computing the discrete Radon transform. Note that the elements of the **L** matrix only depend on the geometry of the input CMP gather and the range of velocities used in constructing the velocity-stack gather *u*(*v, τ*).

Follow the steps described by equations (**J-2**) through (**J-7a**) to obtain the GLI solution for each Fourier component of *u*(*v, ω′*) given by **u** = (**L ^{T*}L**)

**, where the asterisk denotes complex conjugate. The size of the matrix**

^{−1}L^{T*}d**L**is

^{T}L*m*×

*p*. For

*m*= 120,

*p*= 120, and 500 Fourier components, you would have to solve a matrix equation

**p**= (

**L**)

^{T}L**of the size 120 × 120 for each of the 500 Fourier components. The elements of the matrix**

^{−1}L^{T}t**L**given by equation (

**J-36**) may have similar magnitudes in two different rows because of the small separation of offsets associated with the input CMP gather. As such, the matrix can be near singular and requires a solution based on singular-value decomposition (Section F.3).

Finally, slant-stack transform is a special case of the generalized discrete Radon transform based on linear moveout equation. For details of the generalized linear inversion formulation of the slant-stack transform, see Multichannel filtering techniques for noise and multiple attenuation#F.3 Mathematical foundation of the discrete radon transform|Section F.3]].

## J.4 Dix conversion

We shall develop the traveltime equation for a CMP recording geometry and an earth model that comprises horizontal isovelocity layers (Figure J-1). Then, we shall derive the Dix equation for interval velocities computed from rms velocities associated with the horizontal layers. From the geometry of Figure J-1, note that the half-offset *h* is given by

or

**(**)

From Snell’s law, the ray parameter associated with the raypath in Figure J-1 is

**(**)

where *v _{k}* is the interval velocity for the

*k*th layer. It follows that

**(**)

Also note that

**(**)

where Δ*τ _{k}* is the two-way zero-offset time along the vertical raypath within the

*k*th layer. Substitute equations (

**J-37c**) and (

**J-37d**) into equation (

**J-37a**), and for convenience, define the full offset

*x*= 2

*h*

**(**)

From the geometry of Figure J-1, note that the two-way traveltime *t _{n}*(

*x*) along the raypath from the source

*S*to the reflection point

*R*on the

*n*th interface back to the reflector

*G*is given by the sum of the traveltimes Δ

*t*along the raypath segment within each layer

_{k}or

**(**)

From Snell’s law given by equation (**J-37b**), it follows that

**(**)

Substitute equations (**J-37d**) and (**J-39b**) into equation (**J-39a**)

**(**)

Equations (**J-38**) and (**J-40**) are the raypath equations for an earth model with horizontal isovelocity layers. Based on source-receiver reciprocity, the traveltime can be represented by a curve that is symmetric with respect to offset *x* (Taner and Koehler, 1969)

**(**)

which also is referred to in Section C.1. We need to expand the square-root term in equations (**J-38**) and (**J-40**) into the Taylor series to derive the expressions for *t*^{2}(*x*) and the even powers of *x* displayed in equation (**J-41**) ^{[6]}. The series expansions are

**(**)

and

**(**)

Now, write the infinite series in equations (**J-42a**) and (**J-42b**) as a summation

**(**)

and

**(**)

In equation (**J-43a**), the terms (*a*_{1}, *a*_{2}, *a*_{3}, …) are the fractional coefficients (1, 1/2, 3/8, …) displayed in equation (**J-42a**). Similarly, in equation (**J-43b**), the terms (*b*_{0}, *b*_{1}, *b*_{2}, …) are the fractional coefficients (1, 1/2, 3/8, …) displayed in equation (**J-42b**). Interchange the order of the summations in equations (**J-43a**) and (**J-43b**) to obtain

**(**)

and

**(**)

Define two terms, *A _{j}* and

*B*

_{j}

**(**)

**(**)

and rewrite equations (**J-44a**) and (**J-44b**) as

**(**)

and

**(**)

The even powers of *x* that we need to substitute into equation (**J-41**) can also be expressed as a double summation of the form given by equation (**J-46a**)

**(**)

where *m* = 1, 2, 3, …, and *A*_{j, 2m} is a recursive combination of *A _{j}* of equation (

**J-45a**)

^{[6]}.

Similarly, that we need to substitute into equation (**J-41**) can be expressed as a double summation of the form given by equation (**J-46b**)

**(**)

where *B*_{j, 2} is a recursive combination of *B _{j}* of equation (

**J-45b**)

^{[6]}.

By using equations (**J-47a**) and (**J-47b**) in equation (**J-41**), and setting the terms on both sides of the equation with like powers of *p* equal to another, the coefficients, *C*_{0}, *C*_{1}, *C*_{2}, …, can be determined in terms of the model parameters — interval velocities *v _{k}* and layer thicknesses defined in terms of two-way vertical times Δ

*τ*. While details can be found in Hubral and Krey

_{k}^{[6]}, here, we shall consider a truncated form of equation (

**J-41**) to derive the Dix equation.

Specifically, we shall consider the following form of equation (**J-41**)

**(**)

Expand the series in equations (**J-44a**) and (**J-44b**) up to the second power of the ray parameter *p* and set *a*_{1} = 1, *b*_{0} = 1, and *b*_{1} = 1/2 to obtain

**(**)

and

**(**)

Refer to equation (**J-37b**) and note that by retaining only the terms up to the second power of *p*, we are making a small-spread approximation as in equation (**J-48**).

Define

**(**)

and

**(**)

where *τ _{n}* is the total two-way zero-offset elapsed time along the vertical raypath from the surface down to the

*n*th interface, and

*V*is the rms velocity at the

_{n}*n*th interface measured along the vertical raypath (Figure J-1).

Use the definitions given by equations (**J-50a**) and (**J-50b**) to rewrite equations (**J-49a**) and (**J-49b**) to yield

**(**)

and

**(**)

Now, square both sides of equations (**J-51a**) and (**J-51b**), and again, retain only the terms up to the second power of *p*

**(**)

and

**(**)

Finally, substitute equations (**J-52a**) and (**J-52b**) into equation (**J-48**)

**(**)

and set the terms with like powers of *p* on both sides of the equation equal to one another to determine the coefficients *C*_{0} and *C*_{1}

**(**)

and

**(**)

Back substitution into equation (**J-48**) yields the familiar hyperbolic normal moveout equation

**(**)

We now return to equation (**J-50b**) and split it in the following manner:

**(**)

By the definition of the rms velocity given by equation (**J-50b**), note that

**(**)

Substitute equation (**J-56b**) into equation (**J-56a**)

**(**)

Finally, rearrange the terms and note that Δ*τ _{n}* =

*τ*

_{n}− τ_{n−1}to obtain the Dix equation

**(**)

which is equation (**1**) of the main text. Equation (**J-57**) is the basis for Dix conversion of rms velocities to interval velocities. Given the rms velocities measured from the surface down to the (*n* − 1) and *n*th layer boundaries along the vertical raypaths and the associated two-way zero-offset times *t*_{n−1} and *t _{n}*, the layer velocity within the interval between the two layer boundaries can be computed using equation (

**J-57**). Keep in mind the underlying assumption that the earth model comprises horizontal isovelocity layers. Although reflector dip can be incorporated into the Dix equation

^{[6]}, in practice, it is most desirable to work with rms velocities estimated from gathers generated by prestack time migration. Additionally, keep in mind that the rms velocities used in equation (

**J-57**) are based on straight-ray assumption; thus, ray bending at layer boundaries is not accounted for in Dix conversion. Finally, since equation (

**J-48**) is an approximation to equation (

**J-41**), the offset range used in estimating the rms velocities

*V*and

_{n}*V*

_{n−1}corresponds to a small spread.

## J.5 Map processing

A map is defined as a 2-D surface *g*(*x, y*). Depending on the quantity being mapped, *g*(*x, y*) may have many types of units; for example, gravitational attraction (mGal), magnetic intensity (gamma), elevation, or even times picked along marker horizons from seismic data. Here, we set the convention that the positive *x*-axis points eastward and the positive *y*-axis points northward. A discrete map function is represented by a grid of mesh points over the *x − y* plane. These mesh points are spaced commonly at equal intervals in the *x* and *y* directions. For many types of mapping, *g*(*x, y*) is a smooth function and such a map can be analyzed in the Fourier transform domain. However, there are situations (as for isochron and structure maps) in which the map function has discontinuities that represent faulting.

Maps usually are created from irregularly spaced observation values. Thus, the map function at a particular grid point must be computed by some fitting procedure. A local plane surface fitting procedure is described later in this section.

Figure J-2 shows the contour map of the test surface used in this study. Void grid points have been filled with the arithmetic mean of the map function. Before map creation, some correction may be applied to observed data, followed by various types of editing. The main topic of this section is the next processing stage, which involves various 2-D processing techniques, each of which is aimed at achieving a particular interpretive goal. To clarify the objectives in digital processing, consider the properties of the quantity represented by a map function. In particular, consider the gravity problem. Ultra-long wavelength anomalies generally are associated with variations in crustal thickness. Moderately long wavelengths usually are caused by variations in the basement topography. Medium to short wavelength anomalies are related closely to local tectonic disturbances such as folding. Finally, very short wavelength anomalies have a variety of sources; some are because of small, near-surface features, while others are just noise of various types. Note that gradation in wavelength variations is related somewhat to the depth of the source that causes the anomaly. The longer the wavelength, the deeper the anomaly; the shorter the wavelength, the shallower the anomaly.

Similar physical interpretations are possible with other map functions, such as the seismic (time) horizon map in Figure J-2. Some of the very short wavelengths here result from near-surface effects that are manifested as residual statics. Also note in Figure J-2 that some moderately long wavelengths correspond to structural undulations that exist in the area. The most important observation made from this map is that most of the features with different wavelengths are not spatially isolated, but are superimposed. This characteristic is common for all types of map functions, whether gravity, magnetic, elevation, or time. To separate the effects of different features from each other, we must analyze them in terms of wavelength. This wavelength analysis also is a way to discriminate the depth of the source for some types of data such as potential-field measurements.

The basic motivation behind digital processing of maps is to separate anomalies. Simple 2-D smoothing and wavelength filtering are techniques for separating anomalies. Vertical derivatives and analytic continuation also are useful for enhancing certain anomalies so that they appear more pronounced on the map. Some techniques for separating anomalies by their wavelengths are described next.

Before applying digital processing techniques, it is best to study the map in terms of wavelength composition. The 2-D amplitude spectrum of a map is an excellent tool for recognizing not only the wavelength content, but also the orientation of various components. The most useful display is the color contour plot of the amplitude spectrum (Figure J-3) from which various bands of wavelengths are distinguished clearly. Pink represents long, beige represents moderate, and yellow represents short wavelength anomalies. For a 2-D real function, such as a map, the amplitude spectrum is antisymmetric. Thus, only a pair of quadrants (the first and second) of the amplitude spectrum needs to be displayed.

A simple 2-D smoothing operation is the easiest way to obtain a map that represents the regional anomaly. Figure J-4 shows the map in Figure J-2 after 2-D smoothing. Most of the very short wavelength features present in the center of the original map have been eliminated. Depending on the interpreter’s goal, this output might be a satisfactory regional map. (However, a smoother regional map often is desired.)

Smoothing basically is done by computing the average value of the grid points that fall onto a ring of some desired radius. The center of the ring coincides with the output point. There is no reason why more than one ring should not be considered. For *n* concentric rings with *m _{i}* points over the

*i*th ring, the average value

*g*of the quantity

_{i}*g*that is being mapped is

_{ij}

**(**)

Thus, the cumulative average *g* over *n* rings is

**(**)

Weighting factors, which depend on the distance from the center of the rings, often are used in smoothing algorithms. In these cases, equation (**J-58b**) becomes

**(**)

where *w _{i}* are the weights. The residual anomaly is defined by

**(**)

where *g*_{0} is the grid value at the center of the rings. Figure J-5 shows the regional anomaly obtained by using 15 rings. In general, the more rings, the more smoothing of the data. Note that the output of the ring method is relative. Depending on the nature of the objectives of map processing, Figure J-5 may or may not be a good estimate of the regional anomaly. In most cases, determining what is regional and what is residual is a matter of interpretive judgment.

The hypothetical map in Figure J-6a has ten different anomalies of various shapes and orientations. Figure J-6b is the 2-D amplitude spectrum of this map with (*k _{x}*,

*k*) as the Fourier duals of (

_{y}*x, y*). Anomaly 3, which has infinite wavelength in the north-south direction (

*y*-axis) and a wavelength of 2

*AA′*in the east-west direction (

*x*-axis), lies on the

*k*-axis in the transform domain. Anomaly 7 also lies on the

_{x}*k*-axis, but is farther from the origin since it has a shorter wavelength component (2

_{x}*BB′*) in the

*x*direction. Anomaly 8 tilted counterclockwise by an angle

*θ*degrees from the

*y*direction has an apparent wavelength component 2

*CC′*in the

*y*direction and 2

*DD′*in the

*x*direction.

Note the following relations:

**(**)

where *λ _{x}* =

*DD′*and

*λ*=

_{y}*CC′*are the wavelength components. By definition

and

which, after substitution into equation (**J-60a**), become

**(**)

Equations (**J-60a**) and (**J-60b**) imply that map trend *MM′* is orthogonal to transform trend *TT′*.

Unlike the simple, linear anomalies 3, 7, and 8, rounded anomalies 1, 2, and 10 map onto the entire transform plane. Moreover, because of their equal size, anomalies 1 and 2 cannot be distinguished on the plane of 2-D amplitude spectrum, even though they are isolated in the space domain. However, the spatial extent of these two anomalies is much wider than that of anomaly 10. Thus, they are confined to lower wavenumbers. Finally, elongated features like anomalies 4, 5, 6, and 9 map onto the transform plane at certain dominant directions that are orthogonal to their respective spatial trends. The dots on the *k _{x} − k_{y}* plane in Figure J-6b depict the maximum energy associated with the anomalies.

Gridding involves fitting a locally plane surface to a set of control points around each grid output point. Consider a plane-surface fit in the least-squares sense:

**(**)

The least-squares error is

**(**)

where *g* is the observed value at the grid point (*x, y*), and *M* is the number of observations at and around that grid point. For local plane fitting, *M* usually is set to 9 points. We want to find a set of (*a*_{0}, *a*_{1}, *a*_{2}) for which *L* is minimum;

**(**)

By substituting equation (**J-61**) into equation (**J-62**), we have

**(**)

Then, by doing the differentiations expressed by equation (**J-63**), we get the following set of simultaneous equations

When put into matrix form, we obtain

**(**)

Equation (**J-65**) is solved for the set of coefficients (*a*_{0}, *a*_{1}, *a*_{2}).

Unlike the hypothetical anomalies in Figure J-6a, real data consist of anomalies of various shapes and orientations that are superimposed on each other. However, in the transform domain, the real anomalies can be separated in terms of their wavelength contents and orientations; something that cannot be achieved in the space domain. Thus, the transform domain provides a way to apply various filtering operations to a map. A band of wavelengths can be passed-regardless of orientation, as shown by the radial filter transfer function in Figure J-7a. Anomalies that are oriented in the range (*θ*_{1}, *θ*_{2}) from the northerly direction can be passed regardless of their size, using a directional filter whose transfer function is given by Figure J-7b. Finally, a band of wavelengths with a directional orientation can be passed. This is achieved by a filter with a transfer function as shown in Figure J-7c. In practice, as in other filter design techniques, the transfer functions shown must be tapered in the neighborhood of cutoff wavelengths for optimal performance.

Once the transfer function is designed to suit the purpose, the filter itself can be applied to the map in the transform or the space domain. In the transform domain, the transfer function is multiplied with the 2-D Fourier transform of the map. Subsequent inverse transformation yields the filtered map. To apply the filtering in the space domain, first inverse Fourier transform the filter’s transfer function to obtain its 2-D impulse response. Two-dimensional convolution of this impulse response with the map yields the filtered map. Note that filters whose transfer functions are depicted in Figure J-7 do not cause any phase shift; i.e., anomalies are not displaced in the space domain after filtering.

Before applying the 2-D filters, the possible bands of the wavelengths present on the amplitude spectrum of the map to be filtered must be determined. Four regions were distinguished from the color plot (Figure J-3) of the amplitude spectrum of the test surface. We have the regional component down to 21-km wavelengths. A subregional part of the spectrum between the 21- and 12-km wavelengths is defined. A moderate range of wavelengths is between 12 and 8 km, and the residual components are beyond the 8-km cutoff wavelength. This residual region can be subdivided; anything less than 4 km corresponds to very short wavelength anomalies. These anomalies may be caused by various types of noise.

The bands that were just determined now are used as cutoff wavelengths of the radial filter transfer function. A given map can be scanned with a suite of low-pass filters and several filtered maps can be produced, each with a potentially unique interpretational value. A result from such a filtering operation is shown in Figure J-8. Note the resemblance between Figure J-8 and the regional map (Figure J-5) that was obtained from the ring method described in this section. Although not demonstrated here, note that as the bandwidth of the filter is increased, more and more of the short wavelength anomalies are included, making the output less and less regional in character.

A filter scan is not limited to low-pass filtering only. High-pass, band-pass, and band-reject filters can be applied to maps to achieve a particular interpretational goal. For example, we may want to obtain a residual map that is free of the very short wavelength anomalies that usually are considered noise. This is accomplished by properly selecting the band-pass filter. A high-pass filter is the proper choice to retain all of the short wavelength region of the spectrum.

Directional filters scan the map of interest at various angles to emphasize a particular trend that may exist in the data. Figure J-9 shows the result of one of the several directional filters applied to the test surface in Figure J-2. A low-pass radial filter was cascaded with each directional filter so that any regional trend could be delineated in the area. The results of directional filtering indicate the existence of a prominent regional trend approximately N30°W (Figure J-9). In some cases, a certain band of wavelengths may have one dominant trend that is different from that of another band of wavelengths. This situation may imply a change in the tectonic setting over the geologic history in the area.

## J.6 Reflection traveltime tomography

We want to update the parameters of an earth model that has already been estimated by a suitable combination of inversion methods discussed in this chapter. Reflection traveltime tomography is an inversion method for estimating the earth model parameters from the reflection traveltimes associated with the observed seismic data. The reflection traveltime from a source at the surface to a reflection point at the subsurface and back to a receiver at the surface is represented by an integral of the traveltime segments along the raypath that depend on the earth model parameters themselves. This makes the direct inversion of the traveltimes to estimate the earth model parameters a nonlinear problem. Nevertheless, based on a formal linearization of the eikonal equation, which is a ray-theoretical approximation to the scalar wave equation (Section H.2), it can be shown that small changes in reflection traveltimes are *linearly* related to small changes in earth model parameters ^{[7]}. We are thus motivated to use reflection traveltime tomography not to estimate an initial earth model, but to *update* an already estimated model.

To develop a theory for the model update, we shall set the following framework for our strategy:

- Define the earth model parameters in terms of slowness functions and depths at the boundaries of the layers included in the initial model.
- Assume that the initial model is made up of horizontal layers with laterally invariant model parameters.
- Use the reflection times from CMP data associated with the layer boundaries included in the model and update the earth model such that the discrepancy between the modeled reflection times and the actual reflection times is minimum in the least-squares sense.
- Estimate the
*changes*in the model parameters, rather than the parameters themselves. - Perturb the initial model in two parts — slowness perturbation and depth perturbation.
- Perturb the model parameters while preserving the offset values of the seismic data.

Consider the raypath geometry in a horizontally layered earth model shown in Figure J-10. We want to derive the reflection traveltime equation for the raypath from the reflection point *R _{i}* at the interface

*z*in the subsurface to the receiver location

_{n}*G*at the surface

_{j}*z*= 0. This one-way raypath is associated with a CMP gather at midpoint location

*y*and half-offset

*h*. Our earth model consists of

*n*layers above the reflection point

*R*.

_{i}The modeled traveltime segment from *A* to *C* along the raypath *R _{i}G_{j}* within the

*k*th layer is

**(**)

or

**(**)

since *AB* = *z _{k} − z_{k−1}* and

*AC*=

*AB*sec

*θ*. In equation (

_{k}**J-66a**),

*v*is the interval velocity for the

_{k}*k*th layer and

*x*is the lateral distance from the midpoint location

*y*. Also, in equation (

**J-66b**),

*s*= 1/

_{k}*v*is the slowness of the

_{k}*k*th layer and

*θ*is the angle of incidence at the

_{k}*k*th layer boundary.

The modeled traveltime from the reflection point *R _{i}* on the

*n*th interface to the receiver

*G*is then the sum of the traveltime segments given by equation (

_{j}**J-66b**) within

*n*layers between the points

*R*and

_{i}*G*:

_{j}

**(**)

We now derive the expression for the half-offset *h _{n}* between the midpoint location

*y*and the receiver location

*G*. Refer to Figure J-10 once again and note that the offset segment

_{j}*BC*is given by

**(**)

or

**(**)

The half-offset *h _{n}* is the sum of the lateral segments within the

*n*layers above the reflection point

*R*given by equation (

_{i}**J-68b**):

**(**)

The raypath used in deriving the traveltime and offset equations (**J-67**) and (**J-69**), respectively, is described by the ray parameter *s* as

**(**)

which is the horizontal component of the slowness.

Consider an initial estimate of the parameter vector **p** : (…, *s _{m}*, …,

*z*, …) along the raypath from

_{m}*R*to

_{i}*G*in Figure J-10, where 1 ≤

_{j}*m*≤

*n*. We want to minimize the difference between the observed times

*t*and the modeled times by iteratively perturbing the initial estimate of the parameter vector. A change Δ

_{n}*p*in the parameter vector will change the modeled times as follows:

**(**)

The error in modeling the traveltimes is given by

**(**)

Substitute equation (**J-71**) into equation (**J-72**)

**(**)

Now define the difference Δ*t _{n}* between the observed traveltimes

*t*and the initial estimate of the modeled traveltimes and rewrite equation (

_{n}**J-73**) as

**(**)

The second term on the right is the amount of change in Δ*t _{n}* as a result of the change in the parameter Δ

*p*. Define this term as and rewrite equation (

**J-74a**) once more as

**(**)

where

**(**)

The parameter vector **p** : (…, *s _{m}*, …,

*z*, …), with 1 ≤

_{m}*m*≤

*n*, is composed of the slowness variable

*s*and the depth variable in equation (

_{m}**J-75a**) actually has two parts — one part associated with the slowness perturbation and the other associated with the depth perturbation, given by

**(**)

where are the contributions of the slowness perturbation and depth perturbation, respectively.

To compute the traveltime perturbations caused by the slowness and depth perturbations, we shall follow the elegant formulation by Sherwood ^{[8]}. The slowness perturbation Δ*s _{m}* in layer

*m*, where 1 ≤

*m*≤

*n*, yields a traveltime perturbation given by

**(**)

Since the offset of the input data is preserved in the perturbation, note from Figure J-10 and equation (**J-70**) that, while perturbing the slowness, the propagation angle also is perturbed. That is why we include the second partial differential term on the right-hand side of equation (**J-76a**). In fact, slowness perturbation of a single layer *m* causes perturbation of propagation angles within *all* of the *n* layers along the traveltime path. And that is why we have the summation in the second term on the right-hand side of equation (**J-76a**).

Compute the partial derivatives in equation (**J-76a**) using equation (**J-67**) to yield

**(**)

Substitute equation (**J-70**) into the second term on the right-hand side of this equation:

**(**)

Apply the same perturbation given by equation (**J-76a**) as for the traveltime equation (**J-67**) to the offset equation (**J-69**) to obtain

**(**)

Compute the partial derivatives in equation (**J-77a**) using equation (**J-69**) to yield

**(**)

Note that this summation is the same as the summation in the second term on the right-hand side of equation (**J-76c**). A direct consequence of the rule for preserving the offset value during perturbation leads to setting Δ*h _{n}*(

*s*) = 0 in equation (

_{m}**J-77b**). This means that the second term on the right-hand side of equation (

**J-76c**) vanishes, giving us the final expression for the slowness perturbation:

**(**)

We now evaluate the traveltime perturbation caused by a depth perturbation Δ*z _{m}* in layer

*m*, where 1 ≤

*m*≤

*n*. This perturbation is given by

**(**)

Since the offset of the input data is preserved in the perturbation, note from Figure J-1 that, while perturbing the depth, the propagation angle also is perturbed. Again, that is why we include the second partial differential term on the right-hand side of equation (**J-79a**). Since depth perturbation of a single layer *m* causes perturbation of propagation angles within *all* of the *n* layers along the traveltime path, we have the summation in the second term on the right-hand side of equation (**J-79a**).

Compute the partial derivatives in equation (**J-79a**) using equation (**J-67**) to yield

**(**)

A perturbation in the thickness defined by (*z _{m} − z*

_{m−1}) also causes a change in the thickness (

*z*

_{m+1}−

*z*), but in the opposite direction. That is why we have the two parts in the first term of the right-hand side of equation (

_{m}**J-79b**) — (

*s*sec

_{m}*θ*) and (

_{m}*s*

_{m+1}sec

*θ*

_{m+1}). Substitute equation (

**J-70**) into the second term on the right-hand side of this equation to obtain

**(**)

Apply the same perturbation given by equation (**J-79a**) as for the traveltime equation (**J-67**) to the offset equation (**J-69**):

**(**)

Compute the partial derivatives in equation (**J-80a**) using equation (**J-69**) to yield

**(**)

The two parts of the first term of the right-hand side of this equation (tan *θ _{m}*) and (tan

*θ*

_{m+1}) are justified as for equation (

**J-79b**).

A direct consequence of the rule for preserving the offset value during perturbation leads to setting Δ*h _{n}*(

*z*) = 0 in equation (

_{m}**J-80b**). This then leads to

**(**)

which, upon substitution into equation (**J-79c**) yields

**(**)

A further substitution of equation (**J-70**) then is made to get

**(**)

Simplification of the terms gives us the final expression for the depth perturbation as

**(**)

We now combine, using equation (**J-75b**), the slowness perturbation given by equation (**J-78**) and the depth perturbation given by equation (**J-82**) and obtain

**(**)

Examine the structure of equation (**J-83**) and note that, instead of modeling reflection traveltimes, we model the change in traveltimes, and thus estimate the change in the model parameters. The parameter vector is **Δp** : (Δ*s _{m}*, …, Δ

*z*, …), where 1 ≤

_{m}*m*≤

*n*.

Finally, write equation (**J-83**) in matrix form for all the traveltime perturbations to obtain

**(**)

where

**(**)

and

**(**)

Write equation (**J-84**) in compact matrix notation:

**(**)

where **Δt′** is the column vector on the left-hand side of equation (**J-84**), **L** is the sparse matrix with nonzero elements *Z _{m}* and

*S*given by equations (

_{m}**J-85a**,

**J-85b**), and

**Δp**is the column vector on the right-hand side of equation (

**J-84**). Except for the two elements in each row, the

**L**matrix contains zeros.

The generalized linear inversion (GLI) solution to equation (**J-86**) satisfies the requirement that the energy of the error vector

**(**)

is minimum and is given by (Section J.1)

**(**)

where **Δt** denotes the column vector that represents the difference between the observed reflection traveltimes and the initial estimate of the modeled times, and *T* denotes matrix transposition.

We now outline the procedure for a tomographic model update based on the theory described in this section:

- Perform prestack depth migration using the initial model and generate a set of image gathers.
- Compute the residual moveout for all offsets along events on image gathers that correspond to the layer boundaries included in the model, and thus build the traveltime error vector
**Δt**. For instance, you may have 10 layers, 1000 CMPs with a fold of 30. This means that the length of the traveltime error vector**Δt**is 300 000. - Define the initial model by a set of slowness and depth parameters, and construct the coefficient matrix
**L**in equation (**J-88**) by computing the nonzero matrix elements*Z*and_{m}*S*given by equations (_{m}**J-85a**) and (**J-85b**). For instance, you may have 10 layers and each layer may be defined by 50 pairs of slowness and depth values in the lateral direction. This means you would have 1000 parameters in your model space. It also means that for the example given in step (b), you have 300 000 equations to solve for 1000 parameters. The solution for this overdetermined system is given by equation (**J-88**). - Estimate the change in parameters vector
**Δp**, by way of the GLI solution given by equation (**J-88**). - Update the parameter vector
**p**+**Δp**. - Iterate steps (a) through (e) as necessary to minimize the discrepancy between the modeled and actual traveltimes.

## J.7 Threshold for velocity-depth ambiguity

Velocity-depth ambiguity states that an error in depth is indistinguishable from an error in velocity. Specifically, we shall investigate under what circumstances two earth models defined by layer velocities and reflector geometries in depth are indistinguishable. We shall adopt the derivation by Lines ^{[9]}. Albeit the underlying theory to be presented is somewhat simplistic compared to the comprehensive analysis by Bickel ^{[10]}, it does provide a practical insight to important aspects of velocity-depth ambiguity.

Consider a single, horizontal reflector at depth *z* in a medium with constant velocity *v*. The reflection traveltime associated with a CMP recording geometry of a fixed offset 2*h* is given by the hyperbolic moveout equation

**(**)

Perturbation of velocity *v* by Δ*v* and depth *z* by Δ*z* causes a change in traveltime *t*(*v, z*) by Δ*t* expressed as

**(**)

First, differentiate both sides of equation (**J-89**) with respect to velocity *v*

**(**)

and simplify by way of equation (**J-89**)

**(**)

Next, differentiate both sides of equation (**J-89**) with respect to depth *z*

**(**)

and simplify, again, by way of equation (**J-89**)

**(**)

Substitute the partial differentials given by equations (**J-91b**) and (**J-92b**) into equation (**J-90**), then normalize with respect to traveltime *t*

**(**)

For the case of zero offset, *t* = 2*z*/*v*; hence, equation (**J-93**) takes the special form

**(**)

Equation (**J-94**) states a rule of fundamental importance: *for a zero-offset case, when the perturbation in velocity* Δ*v*/*v is the same as the perturbation in depth* Δ*z*/*z, no change occurs in traveltime; thus, velocity-depth ambiguity is infinite*.

Now consider two earth models — an unperturbed model defined by the variables (*z, v*) and a perturbed model defined by the variables (*z* + Δ*z*, *v* + Δ*v*), such that the two models are indistinguishable at zero offset. This means that the zero-offset times *t*_{0} associated with the two models are identical. Rewrite equation (**J-89**) in terms of the constant zero-offset time *t*_{0} = 2*z*/*v*

**(**)

Perturb the velocity *v* by Δ*v* and evaluate the change in time Δ*t* to get

**(**)

Evaluate the partial differential by way of equation (**J-95**) to obtain

**(**)

Expand the square-root term by the Taylor series up to the second order to get

Finally, retain only the term with power of two in *h*, normalize with respect to traveltime *t*_{0} and rearrange the variables to obtain the expression for the absolute value of the fractional change in velocity Δ*v*/*v*

**(**)

where *z* = *vt*_{0}/2. Equation (**J-97**) states another rule for velocity-depth ambiguity: *for a reflector at depth z, velocity ambiguity defined by* Δ*v*/*v can be made smaller by increasing the offset h and decreasing the error in the time pick defined by* Δ*t*/*t*_{0} *at offset h*. Note that for the zero-offset case, *h* = 0, the velocity ambiguity is infinite. This conclusion also was reached earlier by way of equation (**J-94**).^{[11]} ^{[12]}

## References

- ↑ Lines and Treitel (1984), Lines, L. R. and Treitel, S., 1984, Tutorial: A review of least-squares inversion and its application to geophysical problems: Geophys. Prosp., 32, 159–186.
- ↑
^{2.0}^{2.1}^{2.2}Wang and Treitel, 1973, Wang, R. J. and Treitel, S., 1973, The determination of digital Wiener filters by means of gradient methods: Geophysics, 38, 310–326. - ↑ Koehler and Taner, 1985, Koehler, F. and Taner, M. T., 1985, The use of the conjugate-gradient algorithm in the computation of predictive deconvolution operators: Geophysics, 50, 2752–2758.
- ↑ Claerbout, 1992, Claerbout, J. F., 1992, Earth soundings analysis — processing versus inversion: Blackwell Scientific Publications.
- ↑ Hestenes (1956), Hestenes, M. R., 1956, The conjugate gradient method for solving linear systems: Proc. Sympo. in Applied Math., VI, McGraw-Hill Book Co.
- ↑
^{6.0}^{6.1}^{6.2}^{6.3}^{6.4}Hubral and Krey (1980), Hubral, P. and Krey, T., 1980, Interval velocities from seismic reflection time measurements: Soc. Expl. Geophys. - ↑ Aldridge, 1994, Aldridge, D. F., 1994, Linearization of the eikonal equation: Geophysics, 59, 1631–1632.
- ↑ Sherwood et al., 1986, Sherwood, J. W. C., Chen, K. C., and Wood, M., 1986, Depths and interval velocities from seismic reflection data for low-relief structures: Proc. Offshore Tech. Conf., 103–110.
- ↑ Lines, 1993, Lines, L., 1993, Ambiguity in analysis of velocity and depth: Geophysics, 58. 596–597.
- ↑ Bickel, 1990, Bickel, S. H., 1990, Velocity-depth ambiguity of reflection traveltimes: Geophysics, 55, 266–276.
- ↑ Landa, E., Thore, P., Sorin, V., and Koren, Z., 1991, Interpretation of velocity estimates from coherency inversion: Geophysics, 56, 1377–1383.
- ↑ Thorson, J. R., Gever, D. H., Swanger, H. J., Hadley, D. M., and Apsel, R. J., 1987, A model-based approach to interval velocity analysis: 57th Ann. internat. Mtg., Expanded Abstracts, 458–460.

## See also

- Introduction to earth modeling in depth
- Models with horizontal layers
- Model with low-relief structure
- Model with complex overburden structure
- Model building
- Model updating
- Exercises