# Data modeling by inversion

Series Investigations in Geophysics Öz Yilmaz http://dx.doi.org/10.1190/1.9781560801580 ISBN 978-1-56080-094-1 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 L2, 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 L1, 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

 ${\displaystyle \mathbf {d} ^{\prime }=\mathbf {Lp} ,}$ (J-1)

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

 ${\displaystyle \mathbf {e} =\mathbf {d} -\mathbf {d^{\prime }} .}$ (J-2)

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

 ${\displaystyle \mathbf {e} =\mathbf {d} -\mathbf {Lp} .}$ (J-3)

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

 ${\displaystyle S=\mathbf {e^{T}e} ,}$ (J-4a)

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

 ${\displaystyle S=(\mathbf {d} -\mathbf {Lp} )^{\mathbf {T} }(\mathbf {d} -\mathbf {Lp} ).}$ (J-4b)

Expand the right-hand side to get

 ${\displaystyle S=\mathbf {d^{T}d} -\mathbf {p^{T}L^{T}d} -\mathbf {d^{T}Lp} +\mathbf {p^{T}L^{T}Lp} .}$ (J-4c)

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

 ${\displaystyle -\mathbf {d^{T}L} +\mathbf {p^{T}L^{T}L} =\mathbf {0} .}$ (J-5a)

Apply matrix transpose and rearrange the terms

 ${\displaystyle (\mathbf {L^{T}L} )\mathbf {p} =\mathbf {L^{T}d} ,}$ (J-5b)

which yields the desired least-squares solution

 ${\displaystyle \mathbf {p} =(\mathbf {L^{T}L} )^{\mathbf {-} 1}\mathbf {L^{T}d} ,}$ (J-6a)

where LTL is the covariance matrix and (LTL)−1LT is the least-squares (also called generalized linear) inverse of 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 LTL needs to be inverted. The constrained solution is given by

 ${\displaystyle \mathbf {p} =(\mathbf {L^{T}L} +\beta \mathbf {I} )^{\mathbf {-} 1}\mathbf {L^{T}d} ,}$ (J-6b)

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

 ${\displaystyle \mathbf {p} =(\mathbf {L^{T\ast }L} )^{\mathbf {-} 1}\mathbf {L^{T\ast }d} ,}$ (J-7a)

and the constrained solution is given by

 ${\displaystyle \mathbf {p} =(\mathbf {L^{T\ast }L} +\beta \mathbf {I} )^{\mathbf {-} 1}\mathbf {L^{T\ast }d} ,}$ (J-7b)

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 LTL 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

 ${\displaystyle y(t)=w(t)\ast f(t)}$ (J-8)

Consider the discrete form of equation (J-8), with w(t) represented by the m–length time series (w0, w1, w2, …, wm−1), and f(t) represented by the n–length time series (f0, f1, f2, …, fn−1).

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

 ${\displaystyle {\begin{pmatrix}y_{0}\\y_{1}\\y_{2}\\\vdots \\y_{m+n-1}\end{pmatrix}}={\begin{pmatrix}w_{0}&0&\ldots &0\\w_{1}&w_{0}&&\\w_{2}&w_{1}&w_{0}&\\\vdots &w_{2}&w_{1}&\ddots \\w_{m-1}&\vdots &w_{2}&\ddots \\0&w_{m-1}&\vdots &\\\vdots &&w_{m-1}&\\0&&&w_{m-1}\end{pmatrix}}{\begin{pmatrix}f_{0}\\f_{1}\\f_{2}\\\vdots \\f_{n-1}\end{pmatrix}}.}$ (J-9)

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

 ${\displaystyle \mathbf {y} =\mathbf {Lf} .}$ (J-10)

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

 ${\displaystyle (\mathbf {L^{T}L} )\mathbf {f} =\mathbf {L^{T}d} .}$ (J-11)

Consider the special case of a three-point wavelet (w0, w1, w2). Set up the L matrix of equation (J-9) for this special case

 ${\displaystyle \mathbf {L} ={\begin{pmatrix}w_{0}&0&0\\w_{1}&w_{0}&0\\w_{2}&w_{1}&w_{0}\\0&w_{2}&w_{1}\\0&0&w_{2}\end{pmatrix}},}$ (J-12a)

with its transpose LT given by

 ${\displaystyle \mathbf {L^{T}} ={\begin{pmatrix}w_{0}&w_{1}&w_{2}&0&0\\0&w_{0}&w_{1}&w_{2}&0\\0&0&w_{0}&w_{1}&w_{2}\\\end{pmatrix}}.}$ (J-12b)

Now multiply the two matrices to get

 ${\displaystyle \mathbf {L^{T}L} ={\begin{pmatrix}w_{0}^{2}+w_{1}^{2}+w_{2}^{2}&w_{1}w_{0}+w_{2}w_{1}&w_{2}w_{0}\\w_{0}w_{1}+w_{1}w_{2}&w_{0}^{2}+w_{1}^{2}+w_{2}^{2}&w_{1}w_{0}+w_{2}w_{1}\\w_{0}w_{2}&w_{0}w_{1}+w_{1}w_{2}&w_{0}^{2}+w_{1}^{2}+w_{2}^{2}\end{pmatrix}}.}$ (J-12c)

Compute the first three lags of the autocorrelation (r0, r1, r2) of the wavelet (w0, w1, w2):

 ${\displaystyle {\begin{array}{l}r_{0}=w_{0}^{2}+w_{1}^{2}+w_{2}^{2}\\r_{1}=w_{0}w_{1}+w_{1}w_{2}\\r_{2}=w_{0}w_{2},\end{array}}}$ (J-13)

and note that the elements of the covariance matrix LTL given by equation (J-12c) are the first three autocorrelation lags of the wavelet (w0, w1, w2) given by equations (J-13). Thus, for the general case, we note that

 ${\displaystyle \mathbf {L^{T}L} ={\begin{pmatrix}r_{0}&r_{1}&r_{2}&\ldots &r_{n-1}\\r_{1}&r_{0}&r_{1}&\ldots &r_{n-2}\\r_{2}&r_{1}&r_{0}&\ldots &r_{n-3}\\\vdots &\vdots &\vdots &\ddots &\vdots \\r_{n-1}&r_{n-2}&r_{n-3}&\ldots &r_{0}\end{pmatrix}},}$ (J-14)

where (r0, r1, r2, …, rn−1) are the first n autocorrelation lags of the input wavelet series (w0, w1, w2, wm−1).

For the special case of a desired output vector d

 ${\displaystyle \mathbf {d} ={\begin{pmatrix}d_{0}\\d_{1}\\d_{2}\\d_{3}\\d_{4}\\\end{pmatrix}},}$ (J-15a)

obtain the matrix product LTd

 ${\displaystyle \mathbf {L^{T}d} ={\begin{pmatrix}w_{0}d_{0}+w_{1}d_{1}+w_{2}d_{2}\\w_{0}d_{1}+w_{1}d_{2}+w_{2}d_{3}\\w_{0}d_{2}+w_{1}d_{3}+w_{2}d_{4}\\\end{pmatrix}}.}$ (J-15b)

Now compute the first three lags of the crosscorrelation (g0, g1, g2) of the desired output (d0, d1, d2, d3, d4) with the input wavelet (w0, w1, w2)

 ${\displaystyle g_{0}=w_{0}d_{0}+w_{1}d_{1}+w_{2}d_{2},}$ (J-16a)

 ${\displaystyle g_{1}=w_{0}d_{1}+w_{1}d_{2}+w_{2}d_{3},}$ (J-16b)

 ${\displaystyle g_{2}=w_{0}d_{2}+w_{1}d_{3}+w_{2}d_{4},}$ (J-16c)

and note that the elements of the matrix LTd given by equation (J-15b) are the first three lags of the crosscorrelation of the desired output (d0, d1, d2, d3, d4) with the input wavelet (w0, w1, w2) given by equations (J-16a, J-16b, J-16c). Thus, for the general case, we note that

 ${\displaystyle \mathbf {L^{T}d} ={\begin{pmatrix}g_{0}\\g_{1}\\g_{2}\\\vdots \\g_{n-1}\end{pmatrix}},}$ (J-17)

where (g0, g1, g2, …, gn−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

 ${\displaystyle {\begin{pmatrix}r_{0}&r_{1}&r_{2}&\ldots &r_{n-1}\\r_{1}&r_{0}&r_{1}&\ldots &r_{n-2}\\r_{2}&r_{1}&r_{0}&\ldots &r_{n-3}\\\vdots &\vdots &\vdots &\ddots &\vdots \\r_{n-1}&r_{n-2}&r_{n-3}&\ldots &r_{0}\end{pmatrix}}{\begin{pmatrix}f_{0}\\f_{1}\\f_{2}\\\vdots \\f_{n-1}\end{pmatrix}}={\begin{pmatrix}g_{0}\\g_{1}\\g_{2}\\\vdots \\g_{n-1}\end{pmatrix}}.}$ (J-18)

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 LTL in equation (J-11) is diagonally dominant since r0 is the largest value of the autocorrelation lags (r0, r1, r2, …, rn−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 LTL:

 ${\displaystyle (\mathbf {L^{T}L} +\epsilon \mathbf {I} )\mathbf {f} =\mathbf {L^{T}d} ,}$ (J-19)

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

Note that the autocorrelation matrix LTL 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 LTL 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.

Equation (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

 ${\displaystyle \mathbf {A} =\mathbf {L^{T}L} }$ (J-20a)

and

 ${\displaystyle \mathbf {B} =\mathbf {L^{T}d} ,}$ (J-20b)

so that

 ${\displaystyle \mathbf {Af} =\mathbf {B} .}$ (J-21)

Our goal is to estimate the filter vector f recursively starting with an initial estimate X0 which may be defined as a null vector. The recursive estimate is terminated when the residual vector Ri defined by

 ${\displaystyle \mathbf {R} _{i}=\mathbf {B} -\mathbf {AX} _{i}}$ (J-22)

becomes a null vector itself. Here, Xi is the parameter vector estimate after the ith 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

 ${\displaystyle c_{-1}=1,}$ (J-23a)

 ${\displaystyle \mathbf {P} _{-1}=0,}$ (J-23b)

and

 ${\displaystyle \mathbf {R} _{0}=\mathbf {B} -\mathbf {A\ X} _{0}.}$ (J-23c)

The recursive scheme computes the energy of the residual after the ith iteration

 ${\displaystyle c_{i}=\mathbf {R} _{i}^{T}\mathbf {R} _{i},\quad i=0,1,2,\cdots ,m.}$ (J-24a)

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

 ${\displaystyle b_{i-1}=c_{i}/c_{i-1}}$ (J-24b)

to get the value for b−1. Next insert the values for P−1 from equation (J-23b), R0 from equation (J-23c), and b−1 from equation (J-24b) into

 ${\displaystyle \mathbf {P} _{i}=\mathbf {R} _{i}+b_{i-1}\mathbf {P} _{i-1}}$ (J-24c)

to get a value for P0. Then, apply the recursions given by

 ${\displaystyle \mathbf {Q} _{i}=\mathbf {AP} _{i},\quad m\leq n,}$ (J-24d)

 ${\displaystyle d_{i}=\mathbf {P} _{i}^{T}\mathbf {Q} _{i},}$ (J-24e)

and

 ${\displaystyle a_{i}=c_{i}/d_{i}}$ (J-24f)

to get an estimate X1 of the filter vector f given by

 ${\displaystyle \mathbf {X} _{i+1}=\mathbf {X} _{i}+a_{i}\mathbf {P} _{i}}$ (J-24g)

and the associated new residual vector R1 given by

 ${\displaystyle \mathbf {R} _{i+1}=\mathbf {R} _{i}-a_{i}\mathbf {Q} _{i}.}$ (J-24h)

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):

 ${\displaystyle {\tilde {X}}_{ij}^{\prime }(\omega )={\tilde {S}}_{j}(\omega )+{\tilde {H}}_{l}(\omega )+{\tilde {E}}_{k}(\omega )+{\tilde {G}}_{i}(\omega ),}$ (J-25)

where ${\displaystyle {\tilde {X}}_{ij}^{\prime }(\omega )}$ is the logarithm of the amplitude spectrum of the modeled seismogram ${\displaystyle x_{ij}^{\prime }(\omega ),\ {\text{and}}\ {\tilde {S}}_{j}(\omega ),{\tilde {H}}_{l}(\omega ),{\tilde {E}}_{k}(\omega )\ {\text{and}}\ {\tilde {G}}_{i}(\omega )}$ are the log-amplitude spectra of the individual components in the time domain — sj(t), gi(t), hl(t), and ek(t), respectively. The term sj(t) is the waveform component associated with source location j, the term gi(t) is the component associated with receiver location i, and the term hl(t) is the component associated with offset dependency of the waveform defined for each offset index l = |i − j|. The fourth component ek(t) represents the earth’s impulse response at the source-receiver midpoint location, k = (i + j)/2.

Consider a data set with ns shot locations and nc channels, so that the total number of traces is ns × nc. For ns × nc values of the actual spectral components ${\displaystyle {\tilde {X}}_{ij}}$ at frequency ω, and ns shot locations, nr receiver locations, ne midpoint locations, and nh offsets, we have the following set of model equations:

 ${\displaystyle {\begin{pmatrix}\ \\\ \\\ \\\vdots &\\{\tilde {X}}_{ij}^{\prime }&\\\vdots \\\ \\\ \\\ \\\end{pmatrix}}={\begin{pmatrix}\ \\\ \\\ \\\cdots 1\cdots \ \ 1\cdots \ \ 1\cdots \ \ 1\cdots \\\ \\\ \\\ \\\end{pmatrix}}{\begin{pmatrix}\vdots \\{\tilde {S}}_{j}\\\vdots \\\vdots \\{\tilde {H}}_{l}\\\vdots \\\vdots \\{\tilde {E}}_{k}\\\vdots \\\vdots \\{\tilde {G}}_{i}\\\vdots \end{pmatrix}}.}$ (J-26)

Write equation (J-26) in matrix notation, ${\displaystyle \mathbf {{\tilde {X}}^{\prime }} =\mathbf {Lp} }$ as in equation (J-1). Here, ${\displaystyle \mathbf {{\tilde {X}}^{\prime }} }$ is the column vector of (ns × nc)-length on the left-hand side in equation (J-26), L is the sparse matrix with dimensions (ns × nc) × (ns + nh + ne + nr), and p is the column vector of (ns + nh + ne + nr)-length on the right-hand side of the same equation. Except the four elements in each row, the L matrix contains zeros. Consider an example of ns = 1000, nc = 240, nh = 60, ne = 2000, and nr = 1000. You would then have 240 000 equations to estimate 4060 parameters — more equations than unknowns.

We want to estimate for each frequency ω the model parameters p such that the difference between the actual spectral component ${\displaystyle \mathbf {\tilde {X}} }$ and the modeled spectral component ${\displaystyle \mathbf {{\tilde {X}}^{\prime }} }$ 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, ${\displaystyle \mathbf {p} =(\mathbf {L^{T\ast }L} )^{\mathbf {-1} }\mathbf {L^{T\ast }{\tilde {X}}} .}$ The parameter vector p that contains the spectral components ${\displaystyle {\tilde {S}}_{j},{\tilde {G}}_{i},{\tilde {H}}_{l}\ {\text{and}}\ {\tilde {E}}_{k},}$ 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 sj(t) * gi(t) * hl(t).

The size of the matrix LTL is (ns + nh + ne + nr) × (ns + nh + ne + nr). Consider 500 frequency components for which equation (J-28) is to be solved for. For the example specified above, you would have to solve a matrix equation ${\displaystyle \mathbf {p} =(\mathbf {L^{T\ast }L} )^{\mathbf {-1} }\mathbf {L^{T\ast }{\tilde {X}}} }$ 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)

 ${\displaystyle t_{ij}^{\prime }=s_{j}+r_{i}+G_{k}+M_{k}x_{ij}^{2},}$ (J-27)

where ${\displaystyle t_{ij}^{\prime }}$ are the modeled traveltime deviations associated with a reflection event on moveout-corrected CMP gathers, sj is the residual statics shift at the jth source location, ri is the residual statics shift at the ith receiver location, Gk is the structure term at the kth midpoint location, and ${\displaystyle M_{k}x_{ij}^{2}}$ is the residual moveout at the kth midpoint location.

Consider a data set with ns shot locations and nc channels, so that the maximum number of picks is ns × nc. For ns × nc picks of tij, and ns shot locations, nr receiver locations, nG midpoint locations, we have the following set of equations:

 ${\displaystyle {\begin{pmatrix}\vdots \\\vdots \\\vdots \\t_{ij}^{\prime }\\\vdots \\\vdots \\\vdots \\\end{pmatrix}}={\begin{pmatrix}\ \\\ \\\ \\\ \\\cdots 1\quad \cdots 1\quad \cdots 1\quad \cdots x_{ij}^{2}\cdots \ \\\ \\\ \\\ \\\ \\\ \\\end{pmatrix}}{\begin{pmatrix}\vdots \\s_{j}\\\vdots \\r_{i}\\\vdots \\G_{k}\\\vdots \\M_{k}\\\vdots \end{pmatrix}},}$ (J-28)

which in compact matrix notation is given by t′ = Lp as in equation (J-1). Here t′ is the column vector of (ns × nc)-length in equation (J-28), L is the sparse matrix in the same equation with dimensions (ns × nc) × (ns + nr + nG + nG), and p is the column vector of (ns + nr + nG + nG)-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 ns = 1000, nc = 240, nr = 1000, and nG = 2000. You would then have 240 000 equations to estimate 6000 parameters.

Follow the steps from equation (J-2) to (J-6) and derive the GLI solution for surface-consistent residual statics, p = (LTL)−1LTt. Here t denotes the column vector of (ns × nc)-length that represents the traveltime deviations picked from moveout-corrected CMP gathers. The size of the matrix LTL is (ns + nr + nG + nG) × (ns + nr + nG + nG). For the example specified above, you would have to solve a matrix equation p = (LTL)−1LTt of the size 6000 × 6000. A practical scheme for solving this equation is based on the Gauss-Seidel method (Section C.4).

To 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)

 ${\displaystyle t_{ij}^{\prime }=T_{j}+T_{i}+s_{b}x_{ij},}$ (J-29)

where Tj and Ti are the intercept time anomalies at shot and receiver locations, respectively, and sb is the bedrock slowness. Here, we shall consider the special case of surface shots. Hence, for n shot/receiver stations the parameter vector is p : (T1, T2, …, Tn; sb).

Consider a data set with n shot locations and nc channels, so that the total number of traces is n × nc. Assume that you will pick the first breaks from half the number of the traces. For n × nc/2 picks of tij, and n + 1 parameters p : (T1, T2, …, Tn; sb), we have the following set of equations:

 ${\displaystyle {\begin{pmatrix}\vdots \\\vdots \\\vdots \\t_{ij}^{\prime }\\\vdots \\\vdots \\\vdots \\\end{pmatrix}}={\begin{pmatrix}\ \\\ \\\ \\\cdots \quad 1\quad \cdots \quad 1\quad \cdots \quad x_{ij}\\\ \\\ \\\ \\\end{pmatrix}}{\begin{pmatrix}\vdots \\T_{j}\\\vdots \\T_{i}\\\vdots \\s_{b}\end{pmatrix}},}$ (J-30)

which in matrix notation is t′ = Lp as in equation (J-1). Here t′ is the column vector of (n × nc/2)-length in equation (J-30), L is the sparse matrix in the same equation with dimensions (n × nc/2) × (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 nc = 240. You would then have 120 000 equations to estimate 1001 parameters.

Follow the steps from equation (J-2) to (J-6) and derive the GLI solution for surface-consistent refraction statics, p = (LTL)−1LTt. Here t denotes the column vector of (n × nc/2)-length that represents the observed (picked) refracted arrival times in equation (J-30). The size of the matrix LTL is (n + 1) × (n + 1). For the example specified above, you would have to solve a matrix equation p = (LTL)−1LTt of the size 1001 × 1001. A practical scheme for solving this equation is based on the Gauss-Seidel method (Section C.4).

In 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)

 ${\displaystyle \Delta t_{ij}^{\prime }=Z_{j}\Delta s_{wj}+Z_{i}\Delta s_{wi}+X_{ij}\Delta s_{b},}$ (J-31)

where ${\displaystyle \Delta t_{ij}^{\prime }}$ is the amount of change in the difference between the observed traveltimes tij and the initial estimate of the modeled traveltimes ${\displaystyle t_{ij}^{\prime }}$ as a result of the change in the parameter Δp; Δswj and Δswi are the slowness perturbations within the weathering layer at shot and receiver locations, respectively; and Δsb is the slowness perturbation within the bedrock layer. The coefficient terms Zj and Zi and Xij are functions of the weathering layer thickness, critical angle of refraction at shot and receiver locations and shot-receiver separation (Section C.9).

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 : (Δsw1, Δsw2, …, Δswn, Δsb). We refer to the scheme based on equation (J-31) as the iterative GLI solution.

Consider a data set with n shot locations and nc channels, so that the total number of traces is n × nc. Again, assume that you will pick the first breaks from half the number of the traces. For n × nc/2 picks of tij and n + 1 parameters Δp : (Δsw1, Δsw2, …, Δswn, Δsb), we have the following set of equations:

 ${\displaystyle {\begin{pmatrix}\vdots \\\vdots \\\vdots \\\Delta t_{ij}^{\prime }\\\vdots \\\vdots \\\vdots \\\end{pmatrix}}={\begin{pmatrix}\ \\\ \\\ \\\cdots \quad Z_{j}\quad \cdots \quad Z_{i}\quad \cdots \quad X_{ij}\\\ \\\ \\\ \\\end{pmatrix}}{\begin{pmatrix}\vdots \\\Delta s_{wj}\\\vdots \\\Delta s_{wi}\\\vdots \\\Delta s_{b}\end{pmatrix}}.}$ (J-32)

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

 ${\displaystyle {\boldsymbol {\Delta }}\mathbf {t^{\prime }} =\mathbf {L} {\boldsymbol {\Delta }}\mathbf {p} ,}$ (J-33)

where Δt′ is the column vector of (n + nc/2)-length on the left-hand side of equation (J-32), L is the sparse matrix in the same equation with dimensions (n + nc/2) × (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 nc = 240, you would have 120 000 equations to estimate 1001 parameters.

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

 ${\displaystyle {\boldsymbol {\Delta }}\mathbf {p} =(\mathbf {L^{T}L} )^{\mathbf {-1} }\mathbf {L^{T}} {\boldsymbol {\Delta }}\mathbf {t} ,}$ (J-34)

where Δt denotes the column vector of n + nc/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 LTL is (n + 1) × (n + 1). For the example specified above, you would have to solve a matrix equation p = (LTL)−1LTt of the size 1001 × 1001. A practical scheme for solving this equation is based on the Gauss-Seidel method (Section C.4).

We 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)

 ${\displaystyle d^{\prime }(h,\omega ^{\prime })=\sum _{v}u(v,\omega ^{\prime })\ \exp(-i\omega ^{\prime }4h^{2}/v^{2}),}$ (J-35)

where ω′ is the Fourier dual of t′, with t′ = t2, 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:

 ${\displaystyle \mathbf {L} ={\begin{pmatrix}e^{-i\omega ^{\prime }4h_{1}^{2}/v_{1}^{2}}&e^{-i\omega ^{\prime }4h_{1}^{2}/v_{2}^{2}}&\ldots &e^{-i\omega ^{\prime }4h_{1}^{2}/v_{p}^{2}}\\e^{-i\omega ^{\prime }4h_{2}^{2}/v_{1}^{2}}&e^{-i\omega ^{\prime }4h_{2}^{2}/v_{2}^{2}}&\ldots &e^{-i\omega ^{\prime }4h_{2}^{2}/v_{p}^{2}}\\{\boldsymbol {\vdots }}&{\boldsymbol {\vdots }}&{\boldsymbol {\ddots }}&{\boldsymbol {\vdots }}\\e^{-i\omega ^{\prime }4h_{m}^{2}/v_{1}^{2}}&e^{-i\omega ^{\prime }4h_{m}^{2}/v_{2}^{2}}&\ldots &e^{-i\omega ^{\prime }4h_{m}^{2}/v_{p}^{2}}\end{pmatrix}},}$ (J-36)

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 = (LT*L)−1LT*d, where the asterisk denotes complex conjugate. The size of the matrix LTL is m × p. For m = 120, p = 120, and 500 Fourier components, you would have to solve a matrix equation p = (LTL)−1LTt of the size 120 × 120 for each of the 500 Fourier components. The elements of the matrix 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

${\displaystyle h=\sum _{k=1}^{n}\Delta x_{k}}$

or

 ${\displaystyle h=\sum _{k=1}^{n}\tan \ \theta _{k}\Delta z_{k}.}$ (J-37a)

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

 ${\displaystyle p={\frac {\sin \theta _{k}}{v_{k}}},}$ (J-37b)

where vk is the interval velocity for the kth layer. It follows that

 ${\displaystyle \tan \ \theta _{k}={\frac {pv_{k}}{\sqrt {1-p^{2}v_{k}^{2}}}}.}$ (J-37c)

Also note that

 ${\displaystyle \Delta z_{k}={\frac {v_{k}\Delta \tau _{k}}{2}},}$ (J-37d)

where Δτk is the two-way zero-offset time along the vertical raypath within the kth layer. Substitute equations (J-37c) and (J-37d) into equation (J-37a), and for convenience, define the full offset x = 2h

 ${\displaystyle x=\sum _{k=1}^{n}{\frac {pv_{k}^{2}\Delta \tau _{k}}{\sqrt {1-p^{2}v_{k}^{2}}}}.}$ (J-38)
Figure J-1  Geometry of a horizontally layered earth model with constant velocities to develop the theory for the normal moveout and Dix equations (Section J.4).

From the geometry of Figure J-1, note that the two-way traveltime tn(x) along the raypath from the source S to the reflection point R on the nth interface back to the reflector G is given by the sum of the traveltimes Δtk along the raypath segment within each layer

${\displaystyle t_{n}(x)=2\sum _{k=1}^{n}\Delta t_{k},}$

or

 ${\displaystyle t_{n}(x)=2\sum _{k=1}^{n}{\frac {\Delta z_{k}}{v_{k}\cos \theta _{k}}}.}$ (J-39a)

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

 ${\displaystyle \cos \ \theta _{k}={\sqrt {1-p^{2}v_{k}^{2}}}.}$ (J-39b)

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

 ${\displaystyle t_{n}(x)=\sum _{k=1}^{n}{\frac {\Delta \tau _{k}}{\sqrt {1-p^{2}v_{k}^{2}}}}.}$ (J-40)

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 ${\displaystyle t_{n}^{2}(x)}$ can be represented by a curve that is symmetric with respect to offset x (Taner and Koehler, 1969)

 ${\displaystyle t_{n}^{2}(x)=C_{0}+C_{1}x^{2}+C_{2}x^{4}+\ldots ,}$ (J-41)

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 t2(x) and the even powers of x displayed in equation (J-41) [6]. The series expansions are

 ${\displaystyle x=\sum _{k=1}^{n}pv_{k}^{2}\Delta \tau _{k}\left(1+{\frac {1}{2}}p^{2}v_{k}^{2}+{\frac {3}{8}}p^{4}v_{k}^{4}+\ldots \right)}$ (J-42a)

and

 ${\displaystyle t_{n}(x)=\sum _{k=1}^{n}\Delta \tau _{k}\left(1+{\frac {1}{2}}p^{2}v_{k}^{2}+{\frac {3}{8}}p^{4}v_{k}^{4}+\ldots \right).}$ (J-42b)

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

 ${\displaystyle x=\sum _{k=1}^{n}\left(\sum _{j=1}^{\infty }p^{2j-1}a_{j}v_{k}^{2j}\Delta \tau _{k}\right)}$ (J-43a)

and

 ${\displaystyle t_{n}(x)=\sum _{k=1}^{n}\left(\sum _{j=0}^{\infty }p^{2j}b_{j}v_{k}^{2j}\Delta \tau _{k}\right).}$ (J-43b)

In equation (J-43a), the terms (a1, a2, a3, …) are the fractional coefficients (1, 1/2, 3/8, …) displayed in equation (J-42a). Similarly, in equation (J-43b), the terms (b0, b1, b2, …) 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

 ${\displaystyle x=\sum _{j=1}^{\infty }p^{2j-1}\left(\sum _{k=1}^{n}a_{j}v_{k}^{2j}\Delta \tau _{k}\right)}$ (J-44a)

and

 ${\displaystyle t_{n}(x)=\sum _{j=0}^{\infty }p^{2j}\left(\sum _{k=1}^{n}b_{j}v_{k}^{2j}\Delta \tau _{k}\right).}$ (J-44b)

Define two terms, Aj and Bj

 ${\displaystyle A_{j}=\sum _{k=1}^{n}a_{j}v_{k}^{2j}\Delta \tau _{k},}$ (J-45a)

 ${\displaystyle B_{j}=\sum _{k=1}^{n}b_{j}v_{k}^{2j}\Delta \tau _{k},}$ (J-45b)

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

 ${\displaystyle x=\sum _{j=1}^{\infty }p^{2j-1}A_{j}}$ (J-46a)

and

 ${\displaystyle t_{n}(x)=\sum _{j=0}^{\infty }p^{2j}B_{j}.}$ (J-46b)

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)

 ${\displaystyle x^{2m}=\sum _{j=1}^{\infty }p^{2(j+m-1)}A_{j,2m},}$ (J-47a)

where m = 1, 2, 3, …, and Aj, 2m is a recursive combination of Aj of equation (J-45a) [6].

Similarly, ${\displaystyle t_{n}^{2}(x)}$ that we need to substitute into equation (J-41) can be expressed as a double summation of the form given by equation (J-46b)

 ${\displaystyle t_{n}^{2}(x)=\sum _{j=0}^{\infty }p^{2j}B_{j,2},}$ (J-47b)

where Bj, 2 is a recursive combination of Bj 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, C0, C1, C2, …, can be determined in terms of the model parameters — interval velocities vk and layer thicknesses defined in terms of two-way vertical times Δτk. While details can be found in Hubral and Krey [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)

 ${\displaystyle t_{n}^{2}(x)\approx C_{0}+C_{1}x^{2}.}$ (J-48)

Expand the series in equations (J-44a) and (J-44b) up to the second power of the ray parameter p and set a1 = 1, b0 = 1, and b1 = 1/2 to obtain

 ${\displaystyle (x)\approx p\sum _{k=1}^{n}v_{k}^{2}\Delta \tau _{k}}$ (J-49a)

and

 ${\displaystyle t_{n}(x)\approx \sum _{k=1}^{n}\Delta \tau _{k}+{\frac {p^{2}}{2}}\sum _{k=1}^{n}v_{k}^{2}\Delta \tau _{k}.}$ (J-49b)

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

 ${\displaystyle \tau _{n}=\sum _{k=1}^{n}\Delta \tau _{k}}$ (J-50a)

and

 ${\displaystyle V_{n}^{2}={\frac {1}{\tau _{n}}}\sum _{k=1}^{n}v_{k}^{2}\Delta \tau _{k},}$ (J-50b)

where τn is the total two-way zero-offset elapsed time along the vertical raypath from the surface down to the nth interface, and Vn is the rms velocity at the nth 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

 ${\displaystyle x\approx pV_{n}^{2}\tau _{n}}$ (J-51a)

and

 ${\displaystyle t_{n}(x)\approx \tau _{n}+{\frac {p^{2}V_{n}^{2}\tau _{n}}{2}}.}$ (J-51b)

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

 ${\displaystyle x^{2}\approx p^{2}V_{n}^{4}\tau _{n}^{2}}$ (J-52a)

and

 ${\displaystyle t_{n}^{2}(x)\approx \tau _{n}^{2}+p^{2}V_{n}^{2}\tau _{n}^{2}.}$ (J-52b)

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

 ${\displaystyle \tau _{n}^{2}+p^{2}V_{n}^{2}\tau _{n}^{2}=C_{0}+C_{1}p^{2}V_{n}^{4}\tau _{n}^{2},}$ (J-53)

and set the terms with like powers of p on both sides of the equation equal to one another to determine the coefficients C0 and C1

 ${\displaystyle C_{0}=\tau _{n}^{2}}$ (J-54a)

and

 ${\displaystyle C_{1}={\frac {1}{V_{n}^{2}}}.}$ (J-54b)

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

 ${\displaystyle t_{n}^{2}(x)\approx \tau _{n}^{2}+{\frac {x^{2}}{V_{n}^{2}}}.}$ (J-55)

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

 ${\displaystyle V_{n}^{2}={\frac {1}{\tau _{n}}}\left(\sum _{k=1}^{n-1}v_{k}^{2}\Delta \tau _{k}+v_{n}^{2}\Delta \tau _{n}\right).}$ (J-56a)

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

 ${\displaystyle V_{n-1}^{2}={\frac {1}{\tau _{n-1}}}\sum _{k=1}^{n-1}v_{k}^{2}\Delta \tau _{k}.}$ (J-56b)

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

 ${\displaystyle V_{n}^{2}={\frac {V_{n-1}^{2}\tau _{n-1}+v_{n}^{2}\Delta \tau _{n}}{\tau _{n}}}.}$ (J-56c)

Finally, rearrange the terms and note that Δτn = τn − τn−1 to obtain the Dix equation

 ${\displaystyle v_{n}={\sqrt {\frac {V_{n}^{2}\tau _{n}-V_{n-1}^{2}\tau _{n-1}}{\tau _{n}-\tau _{n-1}}}},}$ (J-57)

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 nth layer boundaries along the vertical raypaths and the associated two-way zero-offset times tn−1 and tn, 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 Vn and Vn−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.

Figure J-2  A time map of a seismic horizon. The faults have been smoothed out.

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.)

Figure J-3  2-D amplitude spectrum of the map in Figure J-2.

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 mi points over the ith ring, the average value gi of the quantity gij that is being mapped is

 ${\displaystyle {\bar {g}}_{i}={\frac {1}{m_{i}}}\sum _{j=1}^{m_{i}}g_{ij}.}$ (J-58a)

Thus, the cumulative average g over n rings is

 ${\displaystyle {\bar {g}}={\frac {1}{n}}\sum _{i=1}^{n}{\bar {g}}_{i}.}$ (J-58b)

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

 ${\displaystyle {\bar {g}}=\sum _{i=1}^{n}w_{i}{\bar {g}}_{i},}$ (J-58c)

where wi are the weights. The residual anomaly is defined by

 ${\displaystyle g=g_{0}-{\bar {g}},}$ (J-59)

where g0 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.

Figure J-4  A smoothed version of the contour map in Figure J-2.

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 (kx, ky) as the Fourier duals of (x, y). Anomaly 3, which has infinite wavelength in the north-south direction (y-axis) and a wavelength of 2AA′ in the east-west direction (x-axis), lies on the kx-axis in the transform domain. Anomaly 7 also lies on the kx-axis, but is farther from the origin since it has a shorter wavelength component (2BB′) in the x direction. Anomaly 8 tilted counterclockwise by an angle θ degrees from the y direction has an apparent wavelength component 2CC′ in the y direction and 2DD′ in the x direction.

Note the following relations:

 ${\displaystyle \tan \theta ={\frac {\lambda _{x}}{\lambda _{y}}},}$ (J-60a)

where λx = DD′ and λy = CC′ are the wavelength components. By definition

${\displaystyle \lambda _{x}={\frac {2\pi }{k_{x}}}}$

and

${\displaystyle \lambda _{y}={\frac {2\pi }{k_{y}}},}$

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

 ${\displaystyle \tan \theta ={\frac {k_{y}}{k_{x}}}.}$ (J-60b)
Figure J-5  A regional anomaly map based on equation (J-58c). The input is the contour map in Figure J-2.

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 kx − ky 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:

 ${\displaystyle {\tilde {g}}(x,y)=a_{0}+a_{1}x+a_{2}y.}$ (J-61)

The least-squares error is

 ${\displaystyle L=\sum _{i=1}^{M}(g_{i}-{\tilde {g}}_{i})^{2},}$ (J-62)

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 (a0, a1, a2) for which L is minimum;

 ${\displaystyle {\frac {\partial L}{\partial a_{0}}}={\frac {\partial L}{\partial a_{1}}}={\frac {\partial L}{\partial a_{2}}}=0.}$ (J-63)
Figure J-6  Anomalies of various shapes in (a) the space and (b) the wavenumber domains.

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

 ${\displaystyle L=\sum _{i=1}^{M}(g_{i}-a_{0}-a_{1}x_{i}-a_{2}y_{i})^{2}.}$ (J-64)

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

${\displaystyle {\begin{array}{l}\sum a_{0}+\sum a_{1}x+\sum a_{2}y=\sum g,\\\sum a_{0}x+\sum a_{1}x^{2}+\sum a_{2}xy=\sum xg,\\\sum a_{0}y+\sum a_{1}xy+\sum a_{2}y^{2}=\sum yg.\end{array}}}$

When put into matrix form, we obtain

 ${\displaystyle {\begin{pmatrix}M&\sum x&\sum y\\\sum x&\sum x^{2}&\sum xy\\\sum y&\sum xy&\sum y^{2}\\\end{pmatrix}}{\begin{pmatrix}a_{0}\\a_{1}\\a_{2}\end{pmatrix}}={\begin{pmatrix}\sum g\\\sum xg\\\sum yg\end{pmatrix}}.}$ (J-65)

Equation (J-65) is solved for the set of coefficients (a0, a1, a2).

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:

1. Define the earth model parameters in terms of slowness functions and depths at the boundaries of the layers included in the initial model.
2. Assume that the initial model is made up of horizontal layers with laterally invariant model parameters.
3. 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.
4. Estimate the changes in the model parameters, rather than the parameters themselves.
5. Perturb the initial model in two parts — slowness perturbation and depth perturbation.
6. 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 Ri at the interface zn in the subsurface to the receiver location Gj at the surface 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 Ri.

Figure J-10  Geometry of a nonzero-offset ray to develop the theory for the reflection traveltime tomography (Section J.6).

The modeled traveltime segment ${\displaystyle {t_{AC}^{\prime }}}$ from A to C along the raypath RiGj within the kth layer is

 ${\displaystyle t_{AC}^{\prime }={\frac {AC}{v_{k}}}}$ (J-66a)

or

 ${\displaystyle t_{AC}^{\prime }=(z_{k}-z_{k-1})s_{k}\sec \theta _{k},}$ (J-66b)

since AB = zk − zk−1 and AC = AB sec θk. In equation (J-66a), vk is the interval velocity for the kth layer and x is the lateral distance from the midpoint location y. Also, in equation (J-66b), sk = 1/vk is the slowness of the kth layer and θk is the angle of incidence at the kth layer boundary.

The modeled traveltime ${\displaystyle {t_{n}^{\prime }}}$ from the reflection point Ri on the nth interface to the receiver Gj is then the sum of the traveltime segments given by equation (J-66b) within n layers between the points Ri and Gj:

 ${\displaystyle t_{n}^{\prime }=\sum _{k=1}^{n}(z_{k}-z_{k-1})s_{k}\sec \theta _{k}.}$ (J-67)

We now derive the expression for the half-offset hn between the midpoint location y and the receiver location Gj. Refer to Figure J-10 once again and note that the offset segment BC is given by

 ${\displaystyle BC=AB\tan \theta _{k}}$ (J-68a)

or

 ${\displaystyle BC=(z_{k}-z_{k-1})\tan \theta _{k}.}$ (J-68b)

The half-offset hn is the sum of the lateral segments within the n layers above the reflection point Ri given by equation (J-68b):

 ${\displaystyle h_{n}=\sum _{k=1}^{n}(z_{k}-z_{k-1})\tan \theta _{k}.}$ (J-69)

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

 ${\displaystyle s=s_{k}\sin \theta _{k},}$ (J-70)

which is the horizontal component of the slowness.

Consider an initial estimate of the parameter vector p : (…, sm, …, zm, …) along the raypath from Ri to Gj in Figure J-10, where 1 ≤ mn. We want to minimize the difference between the observed times tn and the modeled times ${\displaystyle {t_{n}^{\prime }}}$ by iteratively perturbing the initial estimate of the parameter vector. A change Δp in the parameter vector will change the modeled times as follows:

 ${\displaystyle \left[t_{n}^{\prime }\right]_{modeled}=\left[t_{n}^{\prime }\right]_{initial}+\left[{\frac {\partial t_{n}^{\prime }}{\partial p}}\right]_{modeled}\Delta p.}$ (J-71)

The error in modeling the traveltimes is given by

 ${\displaystyle e_{n}=\left[t_{n}\right]_{observed}-\left[t_{n}^{\prime }\right]_{modeled}.}$ (J-72)

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

 ${\displaystyle e_{n}=\left[t_{n}\right]_{observed}-\left[t_{n}^{\prime }\right]_{initial}-\left[{\frac {\partial t_{n}^{\prime }}{\partial p}}\right]_{modeled}\Delta p.}$ (J-73)

Now define the difference Δtn between the observed traveltimes tn and the initial estimate of the modeled traveltimes ${\displaystyle {t_{n}^{\prime }},}$ and rewrite equation (J-73) as

 ${\displaystyle e_{n}=\Delta t_{n}-{\frac {\partial t_{n}^{\prime }}{\partial p}}\Delta p.}$ (J-74a)

The second term on the right is the amount of change in Δtn as a result of the change in the parameter Δp. Define this term as ${\displaystyle \Delta t_{n}^{\prime },}$ and rewrite equation (J-74a) once more as

 ${\displaystyle e_{n}=\Delta t_{n}-\Delta t_{n}^{\prime },}$ (J-74b)

where

 ${\displaystyle \Delta t_{n}^{\prime }={\frac {\partial t_{n}^{\prime }}{\partial p}}\Delta p.}$ (J-75a)

The parameter vector p : (…, sm, …, zm, …), with 1 ≤ mn, is composed of the slowness variable sm and the depth variable ${\displaystyle z_{m},\ \Delta t_{n}^{\prime }}$ in equation (J-75a) actually has two parts — one part associated with the slowness perturbation and the other associated with the depth perturbation, given by

 ${\displaystyle \Delta t_{n}^{\prime }=\sum _{m=1}^{n}\Delta t_{n}^{\prime }(s_{m})+\sum _{m=1}^{n}\Delta t_{n}^{\prime }(z_{m}),}$ (J-75b)

where ${\displaystyle \Delta t_{n}^{\prime }(s_{m})\ {\text{and}}\ \Delta t_{n}^{\prime }(z_{m})}$ 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 Δsm in layer m, where 1 ≤ mn, yields a traveltime perturbation ${\displaystyle \Delta t_{n}^{\prime }(s_{m})}$ given by

 ${\displaystyle \Delta t_{n}^{\prime }(s_{m})={\frac {\partial t_{n}^{\prime }}{\partial s_{m}}}\Delta s_{m}+\sum _{k=1}^{n}{\frac {\partial t_{n}^{\prime }}{\partial \theta _{k}}}\Delta \theta _{k}.}$ (J-76a)

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

 ${\displaystyle \Delta t_{n}^{\prime }(s_{m})=(z_{m}-z_{m-1})\ \sec \theta _{m}\ \Delta s_{m}+\sum _{k=1}^{n}(z_{k}-z_{k-1})s_{k}\ \sec _{k}^{2}\theta _{k}\sin \theta _{k}\ \Delta \theta _{k}.}$ (J-76b)

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

 ${\displaystyle \Delta t_{n}^{\prime }(s_{m})=(z_{m}-z_{m-1})\ \sec \theta _{m}\ \Delta s_{m}+s\sum _{k=1}^{n}(z_{k}-z_{k-1})\ \sec _{k}^{2}\theta _{k}\ \Delta \theta _{k}.}$ (J-76c)

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

 ${\displaystyle \Delta h_{n}(s_{m})={\frac {\partial h_{n}}{\partial s_{m}}}\Delta s_{m}+\sum _{k=1}^{n}{\frac {\partial h_{n}}{\partial \theta _{k}}}\Delta \theta _{k}.}$ (J-77a)

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

 ${\displaystyle \Delta h_{n}(s_{m})=\sum _{k=1}^{n}(z_{k}-z_{k-1})\ \sec _{k}^{2}\theta _{k}\Delta \theta _{k}.}$ (J-77b)

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 Δhn(sm) = 0 in equation (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:

 ${\displaystyle \Delta t_{n}^{\prime }(s_{m})=(z_{m}-z_{m-1})\ \sec \theta _{m}\ \Delta s_{m}.}$ (J-78)

We now evaluate the traveltime perturbation ${\displaystyle \Delta t_{n}^{\prime }(z_{m})}$ caused by a depth perturbation Δzm in layer m, where 1 ≤ mn. This perturbation is given by

 ${\displaystyle \Delta t_{n}^{\prime }(z_{m})={\frac {\partial t_{n}^{\prime }}{\partial z_{m}}}\Delta z_{m}+\sum _{k=1}^{n}{\frac {\partial t_{n}^{\prime }}{\partial \theta _{k}}}\Delta \theta _{k}.}$ (J-79a)

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

 ${\displaystyle \Delta t_{n}^{\prime }(z_{m})=(s_{m}\ \sec \theta _{m}-s_{m+1}\ \sec \theta _{m+1})\Delta z_{m}+\sum _{k=1}^{n}(z_{k}-z_{k-1})s_{k}\ \sec _{k}^{2}\theta _{k}\ \sin \theta _{k}\ \Delta \theta _{k}.}$ (J-79b)

A perturbation in the thickness defined by (zm − zm−1) also causes a change in the thickness (zm+1zm), but in the opposite direction. That is why we have the two parts in the first term of the right-hand side of equation (J-79b) — (sm sec θm) and (sm+1 sec θm+1). Substitute equation (J-70) into the second term on the right-hand side of this equation to obtain

 ${\displaystyle \Delta t_{n}^{\prime }(z_{m})=(s_{m}\ \sec \theta _{m}-s_{m+1}\ \sec \theta _{m+1})\Delta z_{m}+s\sum _{k=1}^{n}(z_{k}-z_{k-1})\ \sec _{k}^{2}\theta _{k}\ \Delta \theta _{k}.}$ (J-79c)

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

 ${\displaystyle \Delta h_{n}(z_{m})={\frac {\partial h_{n}}{\partial z_{m}}}\Delta z_{m}+\sum _{k=1}^{n}{\frac {\partial h_{n}}{\partial \theta _{k}}}\Delta \theta _{k}.}$ (J-80a)

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

 ${\displaystyle \Delta h_{n}(z_{m})=(\tan \theta _{m}-\tan \theta _{m+1})\Delta z_{m}+\sum _{k=1}^{n}(z_{k}-z_{k-1})\ \sec _{k}^{2}\theta _{k}\ \Delta \theta _{k}.}$ (J-80b)

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 Δhn(zm) = 0 in equation (J-80b). This then leads to

 ${\displaystyle \sum _{k=1}^{n}(z_{k}-z_{k-1})\ \sec _{k}^{2}\theta _{k}\ \Delta \theta _{k}=-(\tan \theta _{m}-\tan \theta _{m+1})\Delta z_{m},}$ (J-81a)

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

 ${\displaystyle \Delta t_{n}^{\prime }(z_{m})=(s_{m}\ \sec \theta _{m}-s_{m+1}\ \sec \theta _{m+1})\Delta z_{m}-s\ (\tan \theta _{m}-\tan \theta _{m+1})\Delta z_{m}.}$ (J-81b)

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

 ${\displaystyle {\begin{array}{l}\Delta t_{n}^{\prime }(z_{m})=(s_{m}\ \sec \theta _{m}-s_{m+1}\ \sec \theta _{m+1})\Delta z_{m}\\\quad \quad \quad \quad -(s_{m}\ \sin ^{2}\theta _{m}\ \sec \theta _{m}-s_{m+1}\ \sin ^{2}\theta _{m+1}\ \sec \theta _{m+1})\Delta z_{m}.\end{array}}}$ (J-81c)

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

 ${\displaystyle \Delta t_{n}^{\prime }(z_{m})=(s_{m}\ \cos \theta _{m}-s_{m+1}\ \cos \theta _{m+1})\Delta z_{m}.}$ (J-82)

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

 ${\displaystyle \Delta t_{n}^{\prime }=\sum _{m=1}^{n}(z_{m}-z_{m-1})\ \sec \theta _{m}\Delta s_{m}+\sum _{m=1}^{n}(s_{m}\ \cos \theta _{m}-s_{m+1}\ \cos \theta _{m+1})\Delta z_{m}.}$ (J-83)

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 : (Δsm, …, Δzm, …), where 1 ≤ mn.

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

 ${\displaystyle {\begin{pmatrix}\vdots \\\vdots \\\vdots \\\Delta t_{n}^{\prime }\\\vdots \\\vdots \\\vdots \\\end{pmatrix}}={\begin{pmatrix}\\\ldots \quad Z_{m}\quad \ldots \quad S_{m}\quad \ldots \\\\\end{pmatrix}}{\begin{pmatrix}\vdots \\\Delta s_{m}\\\vdots \\\Delta z_{m}\\\vdots \\\end{pmatrix}},}$ (J-84)

where

 ${\displaystyle Z_{m}=(z_{m}-z_{m-1})\ \sec \theta _{m}}$ (J-85a)

and

 ${\displaystyle S_{m}=s_{m}\ \cos \theta _{m}-s_{m+1}\cos \theta _{m+1}.}$ (J-85b)

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

 ${\displaystyle {\boldsymbol {\Delta }}\mathbf {t^{\prime }} =\mathbf {L} {\boldsymbol {\Delta }}\mathbf {p} ,}$ (J-86)

where Δt′ is the column vector on the left-hand side of equation (J-84), L is the sparse matrix with nonzero elements Zm and Sm given by equations (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

 ${\displaystyle \mathbf {e} ={\boldsymbol {\Delta }}\mathbf {t} -{\boldsymbol {\Delta }}\mathbf {t^{\prime }} }$ (J-87)

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

 ${\displaystyle {\boldsymbol {\Delta }}\mathbf {p} =(\mathbf {L} ^{\mathbf {T} }\mathbf {L} )^{-1}\mathbf {L} ^{\mathbf {T} }\ {\boldsymbol {\Delta }}\mathbf {t} ,}$ (J-88)

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:

1. Perform prestack depth migration using the initial model and generate a set of image gathers.
2. 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.
3. 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 Zm and Sm given by equations (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).
4. Estimate the change in parameters vector Δp, by way of the GLI solution given by equation (J-88).
5. Update the parameter vector p + Δp.
6. 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 2h is given by the hyperbolic moveout equation

 ${\displaystyle t(v,z)={\frac {2{\sqrt {z^{2}+h^{2}}}}{v}}.}$ (J-89)

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

 ${\displaystyle \Delta t={\frac {\partial t}{\partial v}}\Delta v+{\frac {\partial t}{\partial z}}\Delta z.}$ (J-90)

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

 ${\displaystyle {\frac {\partial t}{\partial v}}=-2{\sqrt {z^{2}+h^{2}}}{\frac {1}{v^{2}}},}$ (J-91a)

and simplify by way of equation (J-89)

 ${\displaystyle {\frac {\partial t}{\partial v}}=-{\frac {t}{v}}.}$ (J-91b)

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

 ${\displaystyle {\frac {\partial t}{\partial z}}={\frac {2z}{v{\sqrt {z^{2}+h^{2}}}}},}$ (J-92a)

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

 ${\displaystyle {\frac {\partial t}{\partial z}}={\frac {4z}{v^{2}t}}.}$ (J-92b)

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

 ${\displaystyle {\frac {\Delta t}{t}}=-{\frac {\Delta v}{v}}+{\frac {4z\Delta z}{v^{2}t^{2}}}.}$ (J-93)

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

 ${\displaystyle {\frac {\Delta t}{t}}=-{\frac {\Delta v}{v}}+{\frac {\Delta z}{z}}.}$ (J-94)

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 t0 associated with the two models are identical. Rewrite equation (J-89) in terms of the constant zero-offset time t0 = 2z/v

 ${\displaystyle t(v)=t_{0}{\sqrt {1+{\frac {4h^{2}}{t_{0}^{2}v^{2}}}.}}}$ (J-95)

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

 ${\displaystyle \Delta t={\frac {\partial t}{\partial v}}\Delta v.}$ (J-96)

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

 ${\displaystyle \Delta t=-{\frac {4h^{2}\Delta v}{v^{3}t_{0}}}{\frac {1}{\sqrt {1+{\dfrac {4h^{2}}{t_{0}^{2}v^{2}}}}}}.}$ (J-96)

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

${\displaystyle \Delta t=-{\frac {4h^{2}\Delta v}{v^{3}t_{0}}}\left(1-{\frac {2h^{2}}{t_{0}^{2}v^{2}}}\right).}$

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

 ${\displaystyle {\frac {\Delta v}{v}}={\frac {z^{2}}{h^{2}}}{\frac {\Delta t}{t_{0}}},}$ (J-97)

where z = vt0/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/t0 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

1. 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. 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.
3. 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.
4. Claerbout, 1992, Claerbout, J. F., 1992, Earth soundings analysis — processing versus inversion: Blackwell Scientific Publications.
5. 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. Hubral and Krey (1980), Hubral, P. and Krey, T., 1980, Interval velocities from seismic reflection time measurements: Soc. Expl. Geophys.
7. Aldridge, 1994, Aldridge, D. F., 1994, Linearization of the eikonal equation: Geophysics, 59, 1631–1632.
8. 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.
9. Lines, 1993, Lines, L., 1993, Ambiguity in analysis of velocity and depth: Geophysics, 58. 596–597.
10. Bickel, 1990, Bickel, S. H., 1990, Velocity-depth ambiguity of reflection traveltimes: Geophysics, 55, 266–276.
11. Landa, E., Thore, P., Sorin, V., and Koren, Z., 1991, Interpretation of velocity estimates from coherency inversion: Geophysics, 56, 1377–1383.
12. 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.