# Mathematical foundation of deconvolution

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 |

## Contents

- 1 B.1 Synthetic seismogram
- 2 B.2 The inverse of the source wavelet
- 3 B.3 The inverse filter
- 4 B.4 Frequency-domain deconvolution
- 5 B.5 Optimum Wiener filters
- 6 B.6 Spiking deconvolution
- 7 B.7 Predictive deconvolution
- 8 B.8 Surface-consistent deconvolution
- 9 B.9 Inverse
*Q*filtering - 10 References
- 11 See also
- 12 External links

## B.1 Synthetic seismogram

Consider an earth model that consists of homogeneous horizontal layers with thicknesses corresponding to the sampling interval. The seismic impedance associated with a layer is defined as *I* = *ρv*, where *ρ* is density and *v* is the compressional-wave velocity within the layer. The instantaneous value of seismic impedance for the *k*th layer is given by

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle I_k={\rho_k}{v_k}.}****(**)

For a vertically incident plane wave, the pressure amplitude reflection coefficient associated with an interface is given by

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle c_k=\frac{I_{k+1}-I_k}{I_{k+1}+I_k}.}****(**)

Assume that the change of density with depth is negligible compared to the change of velocity with depth. Equation (**B-2**) then takes the form

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle c_k=\frac{v_{k+1}-v_k}{v_{k+1}+v_k}.}****(**)

If the amplitude of the incident wave is unity, then the magnitude of the reflection coefficient corresponds to the fraction of amplitude reflected from the interface.

With knowledge of the reflection coefficients, we can compute the impulse response of a horizontally layered earth model using the Kunetz method ^{[1]}. The impulse response contains not only the primary reflections, but also all the possible multiples. Finally, convolution of the impulse response with a source wavelet yields the synthetic seismogram.

Although random noise can be added to the seismogram, the convolutional model used here to establish the deconvolution filters does not include random noise. We also assume that the source waveform does not change as it propagates down into the earth; hence, the convolutional model does not include intrinsic attenuation. The convolution of a seismic wavelet *w*(*t*) with the impulse response *e*(*t*) yields the seismogram *x*(*t*)

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle x(t)=w(t)\ast e(t).}****(**)

By Fourier transforming both sides, we get

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle X(\omega)=W(\omega)E(\omega),}****(**)

where *X*(*ω*), *W*(*ω*), and *E*(*ω*) represent the complex Fourier transforms of the seismogram, the source waveform, and the impulse response, respectively. In terms of amplitude and phase spectra, the Fourier transforms are expressed as

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle X(\omega)=A_x(\omega)\ \exp[i\phi_x(\omega)],}****(**)

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle W(\omega)=A_w(\omega)\ \exp[i\phi_w(\omega)],}****(**)

and

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle E(\omega)=A_e(\omega)\ \exp[i\phi_e(\omega)].}****(**)

By substituting equations (**B-6**) into equation (**B-5**), we have

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle A_x(\omega) = A_w(\omega)A_e (\omega)}****(**)

and

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \phi_x(\omega) = \phi_w(\omega) + \phi_e (\omega).}****(**)

Hence, in convolving the seismic wavelet with the impulse response, their phase spectra are added, while the amplitude spectra are multiplied.

We assume that the earth’s impulse response can be represented by a white reflectivity series; hence, its amplitude spectrum is flat

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle A_e(\omega) = A_0 = constant.}****(**)

By substituting into equation (**B-7a**), we obtain

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle A_x(\omega) = A_0A_w(\omega).}****(**)

This implies that the amplitude spectrum of the seismogram is a scaled version of the amplitude spectrum of the source wavelet.

We now examine the autocorrelation functions *r*(*τ*) of *x*(*t*), *w*(*t*), and *e*(*t*). Autocorrelation is a measure of similarity between the events on a time series at different time positions. It is a running sum given by the expression

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle r_e(\tau) = \frac{1}{N}\sum_{t=0}^{N-1} e_te_{t+\tau}, \quad \tau = 0,\ \mp1,\ \mp2,\ \ldots,\ \mp(N-1),}****(**)

where *τ* is time lag. A random time series is an uncorrelated series. (Strictly, it is uncorrelated when it is continuous and infinitely long.) Therefore,

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle r_e(\tau)=0,\quad\tau=0}****(**)

and

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle {{r}_{e}}\left( 0 \right)={{r}_{0}}=constant.}****(**)

Equation (**B-11**) states that the autocorrelation of a perfect random series is zero at all lags except at zero lag. The zero-lag value actually is the cumulative energy contained in the time series:

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle {{r}_{0}}={{e}_{0}^{2}}+{{e}_{1}^{2}}+\cdots +{{e}_{N-1}^{2}}.}****(**)

Consider the *z*-transform of the convolutional model in equation (**B-4**):

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle X\left(z\right)=W\left(z\right)E\left(z\right).}****(**)

By putting 1/*z* in place of *z* and taking the complex conjugate, we get

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \overline{X}\left( {1}/{z}\; \right)=\overline{W}\left( {1}/{z}\; \right)\overline{E}\left( {1}/{z}\; \right),}****(**)

where the bar denotes the complex conjugate. By multiplying both sides of equations (**B-13**) and (**B-14**), we get

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle X(z)\overline{X}(1/z)=[W(z)E(z)][\overline{W}(1/z)\overline{E}(1/z)].}****(**)

By rearranging the right side,

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle X(z)\overline{X}(1/z)=[W(z)\overline{W}(1/z)][E(z)\overline{E}(1/z)].}****(**)

Finally, by definition, equation (**B-16**) yields

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle r_x=r_w \ast r_e,}****(**)

where *r _{x}*,

*r*, and

_{w}*r*are the autocorrelations of the seismogram, seismic wavelet, and impulse response, respectively. Based on the white reflectivity series assumption (equation

_{e}**B-11**), we have

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle {{r}_{x}}={{r}_{0}}{{r}_{w}}.}****(**)

Equation (**B-18**) states that the autocorrelation of the seismogram is a scaled version of that of the seismic wavelet. We will see that conversion of the seismic wavelet into a zero-lag spike requires knowledge of the wavelet’s autocorrelation. Equation (**B-18**) suggests that the autocorrelation of the seismogram can be used in lieu of that of the seismic wavelet, since the latter often is not known.

## B.2 The inverse of the source wavelet

A basic purpose of deconvolution is to compress the source waveform into a zero-lag spike so that closely spaced reflections can be resolved. Assume that a filter operator *f*(*t*) exists such that

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle w(t)\ast f(t)=\delta(t),}****(**)

where *δ*(*t*) is the Kronecker delta function. The filter *f*(*t*) is called the inverse filter for *w*(*t*). Symbolically, *f*(*t*) is expressed in terms of the seismic wavelet *w*(*t*) as

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle f(t)=\delta (t)\ast \frac{1}{w(t)}.}****(**)

The *z*-transform of the seismic wavelet with a finite length *m* + 1 is (Appendix A)

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle W(z)=w_0+w_1z+w_2z^2+\cdots+w_mz^m.}****(**)

The *z*-transform of the inverse filter can be obtained by polynomial division:

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle F(z)=\frac{1}{W(z)}.}****(**)

The result is another polynomial whose coefficients are the terms of the inverse filter

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle F(z)=f_0+f_1z+f_2z^2+\cdots+f_nz^n+\cdots.}****(**)

Note that the polynomial *F*(*z*) in equation (**B-23**) has only positive powers of *z*; this means *f*(*t*) is causal. If the coefficients of *F*(*z*) asymptotically approach zero as time goes to infinity, so that the filter has finite energy, we say that the filter *f*(*t*) is realizable. If the coefficients increase without bound, we say that the filter is not realizable. In practice, we prefer to work with a causal and realizable filter. Such a filter, by definition, also is minimum-phase. If *f*(*t*) is minimum-phase, then the seismic wavelet *w*(*t*) also must be minimum-phase. Finally, to apply the filter with a finite length *n* + 1, the polynomial *F*(*z*) must be truncated. Truncation of the filter operator induces some error in spiking the seismic wavelet.

The inverse of the seismic wavelet also can be computed in the frequency domain. By Fourier transforming equation (**B-19**), we get

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle W(\omega)F(\omega)=1.}****(**)

By substituting equation (**B-6b**), we obtain

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle F(\omega)=\frac{1}{A_w(\omega)\ \ \exp[i\phi_w(\omega)]}.}****(**)

We express the Fourier transform of the inverse filter *F*(*ω*) as

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle F(\omega)={A_f(\omega)\ \ \exp[i\phi_f(\omega)]},}****(**)

and compare it with equation (**B-25**) to get

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle A_f(\omega)=\frac{1}{A_w(\omega)}}****(**)

and

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \phi_f(\omega)=-\phi_w(\omega).}****(**)

Equations (**B-27**) show that the amplitude spectrum of the inverse filter is the inverse of that of the seismic wavelet, and the phase spectrum of the inverse filter is negative of that of the seismic wavelet.

## B.3 The inverse filter

Instead of the polynomial division procedure described by equation (**B-22**), consider a different approach to derive the inverse filter. Start with the *z*-transform of the autocorrelation of the seismic wavelet:

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle R_w(z)=W(z)\overline{W}(1/z),}****(**)

and the *z*-transform of equation (**B-19**):

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle W\left(z\right)F\left(z\right)=1,}****(**)

from which we get

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle W(z)=\frac{1}{F(z)}.}****(**)

By substituting into equation (**B-28**), we obtain

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle R_w(z)F(z)=\overline{W}(1/z).}****(**)

Since the wavelet is a real-time function,

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle r_w(\tau)=r_w(-\tau).}****(**)

Consider the special case of a three-point inverse filter (*f*_{0}, *f*_{1}, *f*_{2}). We will assume that the seismic wavelet *w*(*t*) is minimum-phase (causal and realizable); hence, its inverse *f*(*t*) also is minimum-phase. The *z*-transform of *f*(*t*) is

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle F(z)=f_0+f_1z+f_2z^2.}****(**)

The *z*-transform of its autocorrelogram *r _{w}*(

*τ*) is, by way of equation (

**B-32**),

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle R_w(z)=\cdots+r_2z^{-2}+r_1z^{-1}+r_0+r_1z+r_2z^2+\cdots}****(**)

The *z*-transform of *w*(*t*) is

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle W\left( z \right)={{w}_{0}}+{{w}_{1}}z+{{w}_{2}}{{z}^{2}}+\cdots +{{w}_{m}}{{z}^{m}},}**

therefore,

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \overline{W}(1/z)=\bar{w}_0+\bar{w}_1z+\bar{w}_2z^2+\cdots+\bar{w}_mz^m.}****(**)

By substituting equations (**B-33a**), (**B-33b**), and (**B-33c**) into equation (**B-31**), we obtain

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle (r_2z^{-2}+r_1z^{-1}r_0+r_1z+r_2z^2)(f_0+f_1z+f_2z^2)=\bar{w}_0+\bar{w}_1z^{-1}+\bar{w}_2z^{-2}.}****(**)

To solve for (*f*_{0}, *f*_{1}, *f*_{2}), we identify the coefficients of powers of *z*. The coefficient of *z*^{0} yields

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle r_0f_0+r_1f_1+r_2f_2=\bar{w}_0,}**

the coefficient of *z*^{1} yields

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle {{r}_{1}}{{f}_{0}}+{{r}_{0}}{{f}_{1}}+{{r}_{1}}{{f}_{2}}=\text{0},}**

while the coefficient of *z*^{2} yields

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle {{r}_{2}}{{f}_{0}}+{{r}_{1}}{{f}_{1}}+{{r}_{0}}{{f}_{2}}=\text{0}.}**

When put into matrix form, these equations for the coefficients of powers of *z* yield

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \begin{pmatrix} r_0&r_1&r_2\\ r_1&r_0&r_1\\ r_2&r_1&r_0\\ \end{pmatrix} \begin{pmatrix} f_0\\ f_1\\ f_2 \end{pmatrix}= \begin{pmatrix} \bar{w}_0\\ 0\\ 0 \end{pmatrix}}****(**)

Note that **Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle {{\bar{w}}_{0}},}**
which equals *w*_{0} for the usual case of a real source wavelet, is the amplitude of the wavelet at *t* = 0. There are four unknowns and three equations. By normalizing with respect to *f*_{0}, we get

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \begin{pmatrix} r_0&r_1&r_2\\ r_1&r_0&r_1\\ r_2&r_1&r_0\\ \end{pmatrix} \begin{pmatrix} 1\\ a_1\\ a_2 \end{pmatrix}= \begin{pmatrix} L\\ 0\\ 0 \end{pmatrix}}****(**)

where *a*_{1} = *f*_{1}/*f*_{0}, *a*_{2} = *f*_{2}/*f*_{0}, and *L* = *w*_{0}/*f*_{0}. We now have three unknowns, *a*_{1}, *a*_{2}, and *L*, and three equations. The square matrix elements on the left side of the equation represent the autocorrelation lags of the seismic wavelet, which we do not know. However, the autocorrelation lags from equation (**B-18**) of the seismogram that we do know can be substituted.

For the general case of an *n*-point inverse filter (*f*_{0}, *f*_{1}, *f*_{2}, …, *f*_{n}), equation (**B-36a**) takes the form

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \begin{pmatrix} r_0&r_1&r_2&\cdots&r_n\\ r_1&r_0&r_1&\cdots&r_{n-1}\\ r_2&r_1&r_0&\cdots&r_{n-2}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ r_n&r_{n-1}&r_{n-2}&\cdots&r_{0}\\ \end{pmatrix} \begin{pmatrix} 1\\ a_1\\ a_2\\ \vdots\\ a_{n-1} \end{pmatrix}= \begin{pmatrix} L\\ 0\\ 0\\ \vdots\\ 0 \end{pmatrix}}****(**)

The autocorrelation matrix in equation (**B-36b**) is of a special type. First, it is a symmetric matrix; second, its diagonal elements are identical. This type of matrix is called a Toeplitz matrix. For *n* number of normal equations, the standard algorithms require a memory space that is proportional to *n*^{2} and CPU time that is proportional to *n*^{3}. Because of the special properties of the Toeplitz matrix, Levinson devised a recursive scheme ^{[1]} that requires a memory space and CPU time proportional to *n* and *n*^{2}, respectively (Section B.6).

## B.4 Frequency-domain deconvolution

We want to estimate a minimum-phase wavelet from the recorded seismogram. The inverse of the wavelet is the spiking deconvolution operator. We start with the autocorrelation of the seismic wavelet *w*(*t*) in the frequency domain:

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle R_w(\omega)=W(\omega)\overline{W}(\omega),}****(**)

where **Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \bar{W}\left( \omega \right)}**
is the complex conjugate of the Fourier transform of *w*(*t*). Since *w*(*t*) normally is not known, *R _{w}*(

*ω*) is not known. However, based on an assumption of white reflectivity series (equation

**B-18**), we can substitute the autocorrelation function of the seismogram in equation (

**B-37**).

Define a new function *U*(*ω*):

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle U(\omega)=\ln[R_w(\omega)].}****(**)

When both sides of equation (**B-38**) are exponentiated:

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle R_w(\omega)=\exp[U(\omega)].}****(**)

Suppose that another function, *ϕ*(*ω*), is defined and that equation (**B-39**) is rewritten as ^{[1]}

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle R_w(\omega)=\exp\left\{\frac{1}{2}[U(\omega)+i\phi(\omega)]\right\}\ \ \exp\left\{\frac{1}{2}[U(\omega)-i\phi(\omega)]\right\}}****(**)

When equation (**B-40**) is compared with equation (**B-37**), we see that

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle W(\omega)=\exp\left\{\frac{1}{2}[U(\omega)+i\phi(\omega)]\right\}.}****(**)

**Figure B-2**Frequency domain spiking deconvolution. (a) Starting with the mixed-phase wavelet (a), the spiking deconvolution operator (d) is derived via the procedure in Figure B-1. This operator is minimum phase; its inverse (g) also is minimum phase. The inverse has the same amplitude spectrum as that of the original wavelet [compare (h) with (b)]. Its phase spectrum is a sign-reversed phase spectrum of the deconvolution operator [compare (i) with (f)]. Convolve operator (d) with original wavelet (a) to yield output (j). Although the spectrum of the output is nearly flat (k), the output (j) is far from being a spike at*t*= 0. The desired output (m) is obtained only with a wavelet that is minimum phase (g), rather than mixed phase (a). Compare these results with those obtained by the time-domain approach for deconvolution as shown in Figure 2.3-2.**Figure 2.3-2**Starting with wavelet (a), autocorrelogram (d) is computed to derive spiking deconvolution operator (e). This operator and its inverse (h) are minimum-phase. The inverse of the deconvolution operator has the same amplitude spectrum as that of the original wavelet. [Compare (i) to (b).] Its phase spectrum is simply the sign-reversed phase spectrum of the spiking deconvolution operator. [Compare (j) to (g).] If operator*e*were convolved with original wavelet*a*, then output*k*results. Although the output spectrum nearly is flat (l), it is far from being a spike at*t*= 0 (n). This desired output (n) is obtained if the input is a minimum-phase wavelet (h), rather than a mixed-phase wavelet (a).

We know *R _{w}*(

*ω*) from equation (

**B-18**) and therefore we know

*U*(

*ω*) from equation (

**B-38**). To estimate

*W*(

*ω*) from equation (

**B-41**), we also need to know

*ϕ*(

*ω*). If we make the minimum-phase assumption,

*ϕ*(

*ω*) turns out to be the Hilbert transform of

*U*(

*ω*)

^{[1]}. To perform the Hilbert transform, first inverse Fourier transform

*U*(

*ω*) back to the time domain. Then, double the positive time values, leave the zero-lag alone, and set the negative time values to zero. This operation yields a time function

*u*

^{+}(

*t*), which vanishes before

*t*= 0. Then, return to the transform domain to get

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle U^{+}(\omega)=\frac{1}{2}[U(\omega)+i\phi(\omega)],}****(**)

where *U*^{+}(*ω*) is the Fourier transform of *u*^{+}(*t*). Finally, exponentiating *U*^{+}(*ω*) yields the Fourier transform *W*(*ω*) of the minimum-phase wavelet *w*(*t*) (equation **B-41**).

Once the minimum-phase wavelet *w*(*t*) is estimated, its Fourier transform *W*(*ω*) is rewritten in terms of its amplitude and phase spectra,

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle W(\omega)=A(\omega)\ \ \exp[i\phi(\omega)].}****(**)

The deconvolution filter in the Fourier transform domain is

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle F(\omega)=\frac{1}{W(\omega)}.}****(**)

By substituting equation (**B-43**), we obtain the amplitude and phase spectra of this filter:

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle A_f(\omega)=\frac{1}{A(\omega)}}****(**)

and

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \phi_f(\omega)=-\phi(\omega).}****(**)

Since the estimated wavelet *w*(*t*) is minimum phase, it follows that the deconvolution filter [whose Fourier transform is given by equation (**B-44**)] also is minimum phase. By inverse Fourier transforming, equation (**B-44**) yields the deconvolution operator. The frequency-domain method of deconvolution is outlined in Figure B-1.

Figure B-2 illustrates the frequency-domain deconvolution that is described in Figure B-1. The panels in Figure B-2 should be compared with the corresponding results of the Wiener-Levinson (time-domain) method in Figure 2-20. As expected, there is virtually no difference between the two results.

From the discussion so far, we see that the spiking deconvolution operator is the inverse of the minimum-phase equivalent of the seismic wavelet. The process of estimating the minimum-phase equivalent of a seismic wavelet is called *spectral decomposition*. The minimum-phase wavelet, once computed, can be inverted easily to get the deconvolution operator. To avoid dividing by zero in equation (**B-45a**) and to ensure that the filter is stable, a small number often is added to the amplitude spectrum before division. This is called *prewhitening*. Equation (**B-45a**) then takes the form

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle A_f(\omega)=\frac{1}{A(\omega)+\varepsilon}.}****(**)

The amplitude spectrum also can be smoothed to get a more stable operator. Smoothing the spectrum is analogous to shortening the equivalent time-domain deconvolution operator.

Zero-phase deconvolution can be implemented conveniently in the frequency domain. To do this, the phase spectrum given by equation (**B-45b**) is set to zero and thus yields a deconvolution operator that flattens the amplitude spectrum of the input seismogram, but does not alter the phase.

## B.5 Optimum Wiener filters

The following concise discussion of optimum Wiener filters is based on ^{[2]}. Consider the general filter model in Figure B-3. Wiener filtering involves designing the filter *f*(*t*) so that the least-squares error between the actual and desired outputs is minimum. Error *L* is defined as

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle L=\sum\limits_{t}(d_t-y_t)^2.}****(**)

The actual output is the convolution of the filter with the input

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle y_t=f_t\ast x_t.}****(**)

When equation (**B-48**) is substituted into equation (**B-47**), we get

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle L=\sum\limits_t\left(d_t-\sum\limits_{\tau}f_\tau x_{t-\tau}\right)^2.}****(**)

The goal is to compute the filter coefficients (*f*_{0}, *f*_{1}, …, *f*_{n-1}) so that the error is minimum. Filter length *n* must be predefined. The minimum error is attained by setting the variation of *L* with respect to *f _{i}* to zero:

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \frac{\partial L}{\partial f_i}=0,\quad i=0,1,2,\ldots,(n-1).}****(**)

By expanding the square term in equation (**B-49**), we have

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle L=\sum\limits_td^2_t-2\sum\limits_t d_t \sum\limits_\tau f_\tau x_{t-\tau}+\sum\limits_t\left(\sum\limits_\tau f_\tau x_{t-\tau}\right)^2.}****(**)

By taking the partial derivatives and setting them to zero, we get

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \frac{\partial L}{\partial f_i}=-2\sum\limits_t d_t x_{t-i}+2\sum\limits_t \left(\sum\limits_{\tau}f_\tau x_{t-\tau}\right)x_{t-i}=0,}****(**)

or

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \sum\limits_\tau f_\tau \sum\limits_t x_{t-\tau}x_{t-i}=\sum\limits_t d_t x_{t-i},\qquad\quad i=0,1,2,\ldots,(n-1).}****(**)

By using

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \sum\limits_t x_{t-\tau}x_{t-i}=r_{i-\tau}}****(**)

and

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \sum\limits_t d_t x_{t-i}=g_i}****(**)

for each *i*th term, we get

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \sum\limits_\tau f_\tau r_{i-\tau}=g_i, \quad i=0,1,2,\ldots,(n-1).}****(**)

When put into matrix form, equation (**B-55**) becomes

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \begin{pmatrix} r_0&r_1&r_2&\cdots&r_{n-1}\\ r_1&r_0&r_1&\cdots&r_{n-2}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ r_{n-1}&r_{n-2}&r_{n-3}&\cdots&r_{0}\\ \end{pmatrix} \begin{pmatrix} f_0\\ f_1\\ \vdots\\ f_{n-1} \end{pmatrix} = \begin{pmatrix} g_0\\ g_1\\ \vdots\\ g_{n-1} \end{pmatrix}}****(**)

Here, *r _{i}* are the autocorrelation lags of the input and

*g*are the lags of the crosscorrelation between the desired output and input. Since the autocorrelation matrix is Toeplitz, the optimum Wiener filter coefficients

_{i}*f*can be computed by using Levinson recursion

_{i}^{[1]}.

The least-squares error involved in this process now is computed. Expressed again, equation (**B-51**) becomes

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle L_{\mathit{min}}=\sum\limits_t d^2_t -2 \left(\sum\limits_t\sum\limits_\tau d_t f_{\tau}x_{t-\tau}\right)+\sum\limits_t\left(\sum\limits_\tau f_{\tau}x_{t-\tau}\right)^2.}****(**)

By substituting the relationships

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \sum\limits_t d_t x_{t-\tau}=g_{\tau}}****(**)

and

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \sum\limits_t x_t x_{t-\tau}=r_{\tau}}****(**)

into equation (**B-57**), we get

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle L_{\it min}=\sum\limits_t d_t^2 - 2 \sum\limits_\tau f_\tau g_\tau+\sum\limits_\tau f_\tau \sum\limits_i f_i \sum\limits_t x_{t-\tau}x_{t-i}}****(**)

or

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle L_{\it min}=\sum\limits_t d_t^2 - 2 \sum\limits_\tau f_\tau g_\tau+\sum\limits_\tau f_\tau \sum\limits_i f_i r_{\tau-i}.}****(**)

Finally, by using equation (**B-55**), we get

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle L_{\it min}=\sum\limits_t d_t^2- \sum\limits_\tau f_\tau g_\tau.}****(**)

## B.6 Spiking deconvolution

Consider the desired output to be the zero-delay spike *d _{t}* : (1, 0, 0, …). Given the input series

*x*: (

_{t}*x*

_{0},

*x*

_{1},

*x*

_{2}, …), equation (

**B-56**) takes the form

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \begin{pmatrix} r_0&r_1&r_2&\cdots&r_{n-1}\\ r_1&r_0&r_1&\cdots&r_{n-2}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ r_{n-1}&r_{n-2}&r_{n-3}&\cdots&r_{0}\\ \end{pmatrix} \begin{pmatrix} f_0\\ f_1\\ \vdots\\ f_{n-1} \end{pmatrix}= \begin{pmatrix} x_0\\ 0\\ \vdots\\ 0 \end{pmatrix}}****(**)

Divide both sides by *f*_{0} to obtain

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \begin{pmatrix} r_0&r_1&r_2&\cdots&r_{n-1}\\ r_1&r_0&r_1&\cdots&r_{n-2}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ r_{n-1}&r_{n-2}&r_{n-3}&\cdots&r_{0}\\ \end{pmatrix} \begin{pmatrix} 1\\ a_1\\ \vdots\\ a_{n-1} \end{pmatrix}= \begin{pmatrix} v\\ 0\\ \vdots\\ 0 \end{pmatrix}}****(**)

where *a _{i}* =

*f*/

_{i}*f*

_{0},

*i*= 1, 2, …,

*n*− 1, and

*v*=

*x*

_{0}/

*f*

_{0}. Equation (

**B-63**) is solved for the unknown quantity

*v*and the filter coefficients (

*a*

_{1},

*a*

_{2}, …,

*a*

_{n−1}). Since the desired output is a zero-delay spike, the filter (1,

*a*

_{1},

*a*

_{2}, …,

*a*

_{n−1}) describes a spiking deconvolution process.

The solution to equation (**B-63**) can be obtained efficiently using the Levinson recursion ^{[1]}. Start with the solution of equation (**B-63**) for the two-term filter (1, *a*_{1}), then solve for the filter (1, *a*_{1}, *a*_{2}) and so on. Equation (**B-63**) for *n* = 2 takes the form

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \begin{pmatrix} r_0&r_1\\ r_1&r_0 \end{pmatrix} \begin{pmatrix} 1\\ a_1 \end{pmatrix}= \begin{pmatrix} v\\ 0\end{pmatrix}}****(**)

Write out the simultaneous equations:

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle {{r}_{0}}+{{r}_{1}}{{a}_{1}}=v}**

and

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle {{r}_{1}}+{{r}_{0}}{{a}_{1}}=\text{0},}**

which yield the filter coefficient *a*_{1} and the unknown variable *v* from the first iteration:

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle {a_1} = - {r_1}/{r_0}}****(**)

and

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle v ={{r}_{0}}+{{r}_{1}}{{a}_{1}}.}****(**)

Now write equation (**B-63**) for the three-term filter **Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \left( 1,\ {a'_{1}},\ {a'_{2}} \right):}**

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \left(\begin{matrix} {{r}_{0}} & {{r}_{1}} & {{r}_{2}}\\ {{r}_{1}} & {{r}_{0}} & {{r}_{1}}\\ {{r}_{2}} & {{r}_{1}} & {{r}_{0}}\\ \end{matrix}\right)\left(\begin{matrix} 1\\ {{{{a}'}}_{1}}\\ {{{{a}'}}_{2}}\\ \end{matrix}\right)=\left(\begin{matrix} {{v}'}\\ 0\\ 0\\ \end{matrix}\right).}****(**)

We want to solve for the filter coefficients **Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \left({a'_{1}},\ {a'_{2}} \right)}**
by using the results of the previous step (equations **B-65a**,**B-64**) by adding a row in the following manner:

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \begin{pmatrix} r_0&r_1&r_2\\ r_1&r_0&r_1\\ r_2&r_1&r_0\\ \end{pmatrix} \begin{pmatrix} 1\\ a_1\\ 0 \end{pmatrix}= \begin{pmatrix} v\\ 0\\ e \end{pmatrix}}****()**

where

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle e = {r_2} + {r_1}{a_1}.}****(**)

Rewrite equation (**B-67**) by changing the order of the rows:

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \begin{pmatrix} r_0&r_1&r_2\\ r_1&r_0&r_1\\ r_2&r_1&r_0\\ \end{pmatrix} \begin{pmatrix} 0\\ a_1\\ 1 \end{pmatrix}= \begin{pmatrix} e\\ 0\\ v \end{pmatrix}}****(**)

Multiply both sides of equation (**B-69**) with variable *c*, yet to be determined, and subtract the result from equation (**B-67**):

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \begin{pmatrix} r_0&r_1&r_2\\ r_1&r_0&r_1\\ r_2&r_1&r_0\\ \end{pmatrix} \begin{pmatrix} 1\\ a_1-ca_1\\ -c \end{pmatrix}= \begin{pmatrix} v-ce\\ 0\\ e-cv \end{pmatrix}}****(**)

Finally, compare equation (**B-70**) with equation (**B-66**), and note that

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \begin{pmatrix} 1\\ a'_1\\ a'_2 \end{pmatrix}= \begin{pmatrix} 1\\ a_1-ca_1\\ -c \end{pmatrix}}****(**)

and

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \begin{pmatrix} v'\\ 0\\ 0 \end{pmatrix}= \begin{pmatrix} v-ce\\ 0\\ e-cv \end{pmatrix}}****(**)

Solve for *c* and *v'*:

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle c=\frac{e}{v},}****(**)

and

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle v'=v\left[1-\left(\frac{e}{v}\right)^2\right].}****(**)

Hence, the new filter coefficients **Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \left({a'_{1}},\ {a'_{2}} \right)}**
are (equations **B-71a**)

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle a'_1=a_1-\frac{e}{v}a_1}****(**)

and

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle a'_2=-\frac{e}{v}.}****(**)

This recursive scheme is repeated to determine the next set of filter coefficients:

- Compute
*v*and*e*using the autocorrelation lags of the input series and the present filter coefficients (equations**B-67**). - Compute the next set of filter coefficients (equations
**B-71a**). - Compute a new value for
*v*(equation**B-72b**).

This recursive process yields the Wiener filter (1, *a*_{1}, *a*_{2}, …, *a*_{n−1}) of the desired length *n*.

## B.7 Predictive deconvolution

Suppose that the desired output in the filter model of Figure B-3 is a time-advanced version of the input, *d*(*t*) = *x*(*t* + *α*). We want to design a Wiener filter *f*(*t*) that predicts *x*(*t* + *α*) from the past values of the input *x*(*t*). In this special case, the crosscorrelation function *g* becomes

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle g_\tau=\sum\limits_t d_tx_{t-\tau}= \sum\limits_t x_{t+\alpha}x_{t-\tau}= \sum\limits_t x_t x_{t-(a+\tau)}.}****(**)

By definition, we have

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle r_\tau=\sum\limits_t x_tx_{t-\tau}.}****(**)

For the *α* + *τ* lag, equation (**B-75**) becomes

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle r_{\alpha+\tau}=\sum\limits_t x_t x_{t-(\alpha+\tau)}.}****(**)

Combine equations (**B-74**) and (**B-76**) and note that *r _{α + τ}* =

*g*. By substituting this result into equation (

_{τ}**B-56**), we get the set of normal equations that must be solved to find the

*prediction filter*(

*f*

_{0},

*f*

_{1}, …,

*f*

_{n−1}):

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \begin{pmatrix} r_0&r_1&r_2&\cdots&r_{n-1}\\ r_1&r_0&r_1&\cdots&r_{n-2}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ r_{n-1}&r_{n-2}&r_{n-3}&\cdots&r_{0}\\ \end{pmatrix} \begin{pmatrix} f_0\\ f_1\\ \vdots\\ f_{n-1} \end{pmatrix}= \begin{pmatrix} r_\alpha\\ r_{\alpha+1}\\ \vdots\\ r_{\alpha+n-1} \end{pmatrix}}****(**)

For a unit prediction lag, *α* = 1, equation (**B-77**) takes the form:

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \begin{pmatrix} r_0&r_1&r_2&\cdots&r_{n-1}\\ r_1&r_0&r_1&\cdots&r_{n-2}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ r_{n-1}&r_{n-2}&r_{n-3}&\cdots&r_{0}\\ \end{pmatrix} \begin{pmatrix} f_0\\ f_1\\ \vdots\\ f_{n-1} \end{pmatrix}= \begin{pmatrix} r_1\\ r_2\\ \vdots\\ r_n \end{pmatrix}}****(**)

By augmenting the right side to the square matrix on the left side, we have

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \begin{pmatrix} -r_1&r_0&r_1&r_2&\cdots&r_{n-1}\\ -r_2&r_1&r_0&r_1&\cdots&r_{n-2}\\ \vdots&\vdots&\vdots&\vdots&\ddots&\vdots\\ -r_{n}&r_{n-1}&r_{n-2}&r_{n-3}&\cdots&r_{0}\\ \end{pmatrix} \begin{pmatrix} 1\\ f_0\\ f_1\\ \vdots\\ f_{n-1} \end{pmatrix}= \begin{pmatrix} 0\\ 0\\ 0\\ \vdots\\ 0 \end{pmatrix}}****(**)

By adding one row, then putting the negative sign to the filter column, we obtain

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \begin{pmatrix} r_0&r_1&r_2&\cdots&r_{n}\\ r_1&r_0&r_1&\cdots&r_{n-1}\\ r_2&r_1&r_0&\cdots&r_{n-2}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ r_{n}&r_{n-1}&r_{n-2}&\cdots&r_{0}\\ \end{pmatrix} \begin{pmatrix} 1\\ -f_0\\ -f_1\\ \vdots\\ -f_{n-1} \end{pmatrix}= \begin{pmatrix} L\\ 0\\ 0\\ \vdots\\ 0 \end{pmatrix}}****(**)

We now have *n* + 1 equations and *n* + 1 unknowns — (*f*_{0}, *f*_{1}, …, *f*_{n−1}, *L*). From equation (**B-80**)

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle L=r_0-r_1f_0-r_2f_1-\cdots -r_nf_{n-1}.}****(**)

By using equation (**B-61**), the minimum error associated with the unit prediction-lag filter can be computed as follows. Start with

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle {d_t} = {x_{t + 1}}}****(**)

so that

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \sum_t d^2_t=r_0}****(**)

and

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle {g_t} = {r_{t + 1}}.}****(**)

Substituting into equation (**B-61**), we get

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle L_{min}=r_0-(r_1f_0+r_2f_1+\cdots +r_nf_{n-1}),}****(**)

which is identical to quantity *L* given by equation (**B-81**). Therefore, when we solve equation (**B-80**), we compute both the minimum error and the filter coefficients.

Equation (**B-77**) gives us the prediction filter with prediction lag *α*. The desired output is *x _{t+α}*. The actual output is

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle {{\hat{x}}_{t+\alpha }}}**which is an estimate of the desired output. This desired output is the predictable component of the input series; i.e., periodic events such as multiples. The error series contains the unpredictable component of the series and is defined as

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle e_{t+\alpha}=x_{t+\alpha}-\hat{x}_{t+\alpha}=x_{t+\alpha}-\sum\limits_{\tau}f_\tau x_{t-\tau}.}****(**)

We assume that the unpredictable component *e _{t}* is the uncorrelated reflectivity we want to extract from the seismogram. By taking the

*z*-transform of equation (

**B-85**), we have

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle z^{-\alpha}E(z)=z^{-\alpha}X(z)-F(z)X(z)}****(**)

or

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle E(z)=[1-z^{\alpha}F(z)]X(z).}****(**)

Now define a new filter *a*(*t*) whose *z*-transform is

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle A(z)=1-z^{\alpha}F(z)}****(**)

so that, when substituted into equation (**B-87**), we have

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle E\left(z\right)=A\left(z\right)X\left(z\right).}****(**)

The corresponding time-domain relationship is

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle e(t)=a(t)\ast x(t).}****(**)

After defining *e*(*t*) as reflectivity, equation (**B-90**) states that by applying filter *a*(*t*) to the input seismogram *x*(*t*), we obtain the reflectivity series. Since computing the reflectivity series is a goal of deconvolution, prediction filtering can be used for deconvolution. The time-domain form *a*(*t*) of the filter defined by equation (**B-88**) is

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle {{a}_{t}}=\left(1,\ \overbrace{0,\ 0,\ \cdots ,\ 0,}^{\alpha -1}-{{f}_{0}},\ -{{f}_{1}},\ \cdots ,\ -{{f}_{n-1}}\right).}****(**)

The filter *a*(*t*) is obtained from the prediction filter *f*(*t*), which is the solution to equation (**B-77**). We call *a*(*t*) the *prediction error filter*. For unit-prediction lag, the filter coefficients given by equation (**B-91**) are of the form

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle a_t=(1,-f_0,-f_1,\cdots,-f_{n-1}).}****(**)

This is the same as the solution of equation (**B-80**). Moreover, equation (**B-80**) is equivalent to equation (**B-36**) for *n* = 2. Thus, we can conclude that a prediction error filter with unit-prediction lag and with *n* + 1 length is equivalent to an inverse filter of the same length except for a scale factor.

## B.8 Surface-consistent deconvolution

Deconvolution can be formulated as a surface-consistent spectral decomposition ^{[3]}. In such a formulation, the seismic trace is decomposed into the convolutional effects of source, receiver, offset, and the earth’s impulse response, thus ex plicitly accounting for variations in wavelet shape affected by near-source and near-receiver conditions and source-receiver separation. Decomposition is followed by inverse filtering to recover the earth’s impulse response. The assumption of surface-consistency implies that the basic wavelet shape depends only on the source and receiver locations, and not on the details of the raypath from source to reflector to receiver.

The convolutional model discussed in The convolutional model is described by

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle x(t)=w(t)\ast e(t)+n(t),}****(**)

where *x*(*t*) is the recorded seismogram, *w*(*t*) is the source waveform, *e*(*t*) is the earth’s impulse response that we want to estimate, and *n*(*t*) is the noise component.

A postulated surface-consistent convolutional model is given by

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle x'_{ij}(t)=s_j(t)\ast h_l(t)\ast e_k(t)\ast g_i(t)+n(t),}****(**)

where **Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle {{{x}'}_{ij}}\left(t \right)}**
is a model of the recorded seismogram, *s _{j}*(

*t*) is the waveform component associated with source location

*j*,

*g*(

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

*i*, and

*h*(

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

*l*= |

*i − j*|. As in equation (

**B-93**),

*e*(

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

*k*= (

*i*+

*j*)/2. By comparing equations (

**B-93**) and (

**B-94**), we infer that

*w*(

*t*) represents the combined effects of

*s*(

*t*),

*h*(

*t*), and

*g*(

*t*).

^{[4]}offer an alternative to equation (

**B-94**) in which the offset term is ignored.

To illustrate a method of computing *s _{j}*(

*t*),

*h*(

_{l}*t*),

*e*(

_{k}*t*), and

*g*(

_{i}*t*), assume

*n*(

*t*) = 0 and Fourier transform equation (

**B-94**):

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle X'_{ij}(\omega)=S_j(\omega)H_l(\omega)E_k(\omega)G_i(\omega).}****(**)

This equation can be separated into the following amplitude spectral components:

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \overline{X}'_{ij}(\omega)=\overline{S}_j(\omega)\overline{H}_l(\omega)\overline{E}_k(\omega)\overline{G}_i(\omega),}****(**)

and phase spectral components (Section A.1):

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \phi'_{ij}(\omega)=\phi_{sj}(\omega)+\phi_{hl}(\omega)+\phi_{ek}(\omega)+\phi_{ri}(\omega).}****(**)

If the minimum-phase assumption is made, only the amplitude spectra (equation **B-96a**) need to be considered.

Equation (**B-96a**) now can be linearized by taking the logarithm of both sides:

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \tilde{X}'_{ij}(\omega)=\tilde{S}_j(\omega)+\tilde{H}_l(\omega)+\tilde{E}_k(\omega)+\tilde{G}_i(\omega).}****(**)

The left side is the logarithm of the amplitude spectrum **Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle {{\overline{X}}^{\prime }}_{ij}\left( \omega \right)}**
of the modeled input trace as in the left-hand side of equation (**B-96a**), and the right-hand terms are the logarithms of the amplitude spectra of the individual components as in the right-hand side of equation (**B-96a**).

An alternative model equation is given by ^{[5]} based on the work by ^{[6]}:

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \tilde{X}''_{ij}(\omega)=\tilde{S}_j(\omega)+\tilde{H}_l(\omega)+\tilde{E}_k(\omega)+\tilde{G}_i(\omega),}****(**)

where

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \tilde{X}''_{ij}(\omega)=\tilde{X}'_{ij}(\omega)-\tilde{X}'_{avg}(\omega).}****(**)

The spectral component **Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle {{{\tilde{X}}'}_{avg}}\left( \omega \right)}**
is associated with the average amplitude spectrum for the entire data set. The terms on the right-hand side of equation (**B-98**) now correspond to *residual* spectral components for each source, offset, midpoint, and receiver location.

The spectral components on the right-hand side of equation (**B-98a**) can be computed by least-squares error minimization. For each frequency *ω*, equation (**B-98a**) is written for each trace of each CMP gather in the data set. 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

_{c}**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle {{\tilde{X}}_{ij}}}**at frequency

*ω*, 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}

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \begin{pmatrix} \\ \vdots\\ \tilde{X}''_{ij}\\ \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}}****(**)

Write equation (**B-99a**) in matrix notation:

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \tilde{\mathbf X}''={\mathbf Lp},}****(**)

where **Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle {\tilde{\mathbf X}}''}**
is the column vector of *m*-length on the left-hand side in equation (**B-99**), **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.

We want to estimate for each frequency *ω* the model parameters **p** such that the difference between the actual spectral component **Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle {\tilde{\mathbf X}}}**
and the modeled spectral component **Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle {\tilde{\mathbf X}}''}**
is minimum in the least-squares sense.

The error vector **v** is defined as the difference between the modeled and the actual spectral component for each frequency *ω*

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle {\mathbf v}=\tilde{\mathbf X}-\tilde{\mathbf X}''.}****(**)

Substitute equation (**B-99b**) into equation (**B-100a**) to obtain

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle {\mathbf v}=\tilde{\mathbf X}-{\mathbf Lp}.}****(**)

Following ^{[7]}, the least-squares solution for equation (**B-100b**) can be determined. First, the cumulative squared error *C* is expressed as

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle C={\mathbf v}^{\mathbf T*}{\mathbf v}.}****(**)

where *T* is for transpose and * is for complex conjugate. By substituting for **v** from equation (**B-100b**), we get

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle C=(\tilde{\mathbf X}-{\mathbf Lp})^{\mathbf T*}(\tilde{\mathbf X}-{\mathbf Lp}).}****(**)

Minimization of *C* with respect to **p** requires that

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \frac{\partial C}{\partial \tilde{S}_j}=\frac{\partial C}{\partial \tilde{H}_l}=\frac{\partial C}{\partial \tilde{E}_k}=\frac{\partial C}{\partial \tilde{G}_i}=0.}**

This requirement yields the desired least-squares solution:

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle {\mathbf p}=({\mathbf L}^{\mathbf T\ast}{\mathbf L})^{-1}{\mathbf L}^{\mathbf T*}\tilde{\mathbf X}.}****(**)

Application of the least-squares minimization to surface-consistent prediction-error filtering is given by ^{[8]}. A practical scheme for solving equation (**B-98a**) is based on the Gauss-Seidel method. In this scheme, each term on the right-hand side of equation (**B-98a**) is computed by the following set of recursive equations:

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \tilde{S}^m_j=\frac{1}{n_r}\sum\limits^{n_r}_i\left\{\tilde{X}_{ij}-\tilde{H}^{m-1}_{l}-\tilde{E}^{m-1}_{k}-\tilde{G}^{m-1}_{i}\right\},}****(**)

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \tilde{G}^m_i=\frac{1}{n_s}\sum\limits^{n_s}_j\left\{\tilde{X}_{ij}-\tilde{H}^{m-1}_{l}-\tilde{E}^{m-1}_{k}-\tilde{S}^{m-1}_{j}\right\},}****(**)

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \tilde{H}^m_l=\frac{1}{n_e}\sum\limits^{n_e}_k\left\{\tilde{X}_{ij}-\tilde{S}^{m-1}_{j}-\tilde{E}^{m-1}_{k}-\tilde{G}^{m-1}_{i}\right\},}****(**)

and

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \tilde{E}^m_k=\frac{1}{n_h}\sum\limits^{n_h}_l\left\{\tilde{X}_{ij}-\tilde{S}^{m-1}_{j}-\tilde{H}^{m-1}_{l}-\tilde{G}^{m-1}_{i}\right\},}****(**)

where *m* is the iteration index. The solutions in equations (**B-103**) are based on the orthogonality of the shot and receiver axes, and the orthogonality of the midpoint and offset axes. Equations (**B-103**) can be modified as follows:

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \tilde{S}_j^m=\frac{1}{n_r}\sum^{n_r}_{i}\{\tilde{X}_{ij}\}-\frac{1}{n_r}\sum^{n_r}_i\{\tilde{H}_l^{m-1} -\tilde{E}_k^{m-1}-\tilde{G}^{m-1}_i\},}****(**)

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \tilde{G}_i^m=\frac{1}{n_s}\sum^{n_s}_{j}\{\tilde{X}_{ij}\}-\frac{1}{n_s}\sum^{n_s}_{j}\{\tilde{H}_l^{m-1} -\tilde{E}_k^{m-1} -\tilde{S}^{m-1}_j\},}****(**)

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \tilde{H}_l^m=\frac{1}{n_e}\sum^{n_e}_{k}\{\tilde{X}_{ij}\}-\frac{1}{n_e}\sum^{n_e}_k\{\tilde{S}_j^{m-1} -\tilde{E}_k^{m-1} -\tilde{G}^{m-1}_i\},}****(**)

and

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \tilde{E}_k^m=\frac{1}{n_h}\sum^{n_h}_{l}\{\tilde{X}_{ij}\}-\frac{1}{n_h}\sum^{n_h}_l\{\tilde{S}_j^{m-1} -\tilde{H}_l^{m-1} -\tilde{G}^{m-1}_i\}.}****(**)

This modification enables us to compute and store the sum of the spectral components of input data **Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \sum{{{X}_{ij}}},}**
thus circumventing the need for storing the individual spectral components *X _{ij}*

^{[5]}. The process is iterated until an index

*m*that attains the least-squares minimization.

The parameter vector **p** that contains the spectral components **Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\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 *ω* using equations (**B-104**). Results from all frequency components are then combined to obtain the terms in equation (**B-98a**). The surface-consistent spiking deconvolution operator to be applied to each trace in the data set is then the minimum-phase inverse of **Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle {{s}_{j}}\left(t \right)*{{g}_{i}}\left(t \right)*{{h}_{l}}\left(t \right).}**
In the case of predictive deconvolution with a desired prediction lag, for each source, receiver and, midpoint location, a deconvolution operator is computed by using the autocorrelograms of the terms *s _{j}*(

*t*),

*g*(

_{i}*t*), and

*h*(

_{l}*t*). To each trace in the data set, these operators are then applied in a cascaded manner.

As a by-product of the derivation of the surface-consistent spectral decomposition equation (**B-97**), trace amplitudes themselves can be corrected for in a surface-consistent manner ^{[9]}. Sum the individual terms in equation (**B-97**) over the frequencies:

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \sum_\omega\tilde{X}'_{ij}(\omega) = \sum_\omega\tilde{S}_j(\omega) + \sum_\omega\tilde{H}_l(\omega) + \sum_\omega\tilde{E}_k(\omega) + \sum_\omega\tilde{G}_i(\omega).}****(**)

Each term yields a scalar that is related to the source, receiver, offset and midpoint locations:

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \tilde{x}'_{ij}=\tilde{s}_j+\tilde{h}_l+\tilde{e}_k+\tilde{g}_i.}****(**)

Once computed in the same manner as the terms in equation (**B-98a**), these scalars are then applied to individual traces in the data set for surface-consistent amplitude corrections.

In practice, application of surface-consistent deconvolution to field data usually involves two terms, only — the source term *s _{j}*(

*t*) and the receiver term

*g*(

_{i}*t*). In a transition zone, surface conditions at source and receiver locations may vary significantly from dry to wet surface conditions. Hence, the most likely situation where surface-consistent amplitude corrections and deconvolution may be required is transition-zone data. Figure B-4 shows a field data example of surface-consistent deconvolution. Note the variations in the autocorrelograms from one source location to the next and from one receiver location to the next. Differences in reflection continuity are observed within the first 1 s on the stacked sections created by application of conventional trace-by-trace deconvolution and surface-consistent deconvolution.

## B.9 Inverse *Q* filtering

Consider a 1-D seismogram that represents a compressional plane wave that propagates vertically downward in a homogeneous medium with intrinsic attenuation. This plane wave is expressed as the solution to the scalar wave equation:

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \frac{1}{v^2}\frac{\partial ^2P}{\partial t^2}=\frac{\partial^2P}{\partial{z}^{2}},}****(**)

where *P*(*t, z*) is the plane wave represented by the 1-D seismogram — a CMP-stacked trace, *t* is the traveltime, *z* is the depth variable and *v* is the wave velocity. We shall assume that the wave velocity is constant.

To solve equation (**B-106**), first, Fourier transform in the time direction:

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \frac{\omega^2}{v^2}P=\frac{\partial^2P}{\partial z^2},}****(**)

where *P*(*ω, z*) is the Fourier transform of the wavefield *P*(*t, z*), and *ω* is the angular frequency. The upcoming wave solution is then given by

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle P(\omega,\ z)=P(\omega,\ z=0)\ \ \exp(-i\frac{\omega} {v}z).}****(**)

To include amplitude decay in wave propagation in a medium with intrinsic attenuation, the wave velocity is defined as a complex variable:

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle v=\alpha+i\beta.}****(**)

Substitute equation (**B-109**) into equation (**B-108**) to get

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle P(\omega,\ z)=P(\omega,\ z=0)\ \ \exp(-i\frac {w}{\alpha+i\beta}z).}****(**)

By simple algebra, rewrite equation (**B-110**) as follows:

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle P(\omega,\ z)=P(\omega,\ z=0)\ \ \exp(-i\frac {w\alpha}{\alpha^2+\beta^2}z)\ \ \exp(-\frac{\omega\beta}{\alpha^2+\beta^2}z).}****(**)

For most rocks, the assumption that *β* is much smaller than *α* can be made. As a result, equation (**B-111**) can be simplified as follows:

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle P(\omega,\ z)=P(\omega,\ z=0)\ \ \exp(-i\frac{\omega}{\alpha}z)\ \ \exp(-\frac{\omega\beta}{\alpha^2}z).}****(**)

Now, define a vertical time variable *τ* equivalent to the depth variable *z* via *z* = *ατ*, and rewrite equation (**B-112**):

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle P(\omega,\ \tau)=P(\omega,\ \tau=0)\ \ \exp(-i\omega\tau)\ \exp(-\omega\frac{\beta}{\alpha}\tau).}****(**)

Assume an attenuation constant *Q* that is independent of frequency *ω* ^{[10]}:

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \frac{1}{2Q}=\frac{\beta}{\alpha},}****(**)

and substitute into equation (**B-113**) to obtain

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle P(\omega,\ \tau)=P(\omega,\ \tau=0)\ \ \exp(-i\omega \tau)\ \ \exp(-\frac{\omega \tau}{2Q}).}****(**)

Note from equation (**B-115**) that the higher the frequency, the greater the attenuation.

For a nondissipative medium, *β* = 0; hence, equation (**B-114**) states that *Q* is infinite. As a result, equation (**B-115**) takes the special form

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle P(\omega,\ \tau)=P(\omega,\ \tau=0)\ \ \exp(-i\omega \tau).}****(**)

The amplitude spectrum of the *inverse Q* filter is thus given by the exponential scaling function

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle A(\omega,\ \tau)=\exp(\frac{\omega \tau}{2Q}).}****(**)

The phase spectrum can either be set to zero, or more appropriately, assumed to be minimum-phase. In the latter case, it can be computed by taking the Hilbert transform of the amplitude spectrum given by equation (**B-117a**) (Section B.4):

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \phi(\omega,\ \tau)=\mathbf{H}\{A(\omega,\tau)\},}****(**)

where **H** represents the Hilbert transform.

By combining the amplitude and phase spectra given by equations (**B-117a**,**B-117b**), we define the minimum-phase inverse *Q* filter as

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle W(\omega,\ \tau)=A(\omega,\ \tau)\ \ \exp\{-i\phi(\omega,\tau)\}.}****(**)

The inverse *Q* filtering equation (**B-115**) now takes the form

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle P(\omega,\ \tau)=P(\omega,\ \tau=0)\ \ \exp(-i\omega \tau)W(\omega,\ \tau).}****(**)

Note that the time variable *t* is associated with the input trace *P*(*t, τ* = 0) and the time variable *τ* is associated with the output trace *P*(*t* = 0, *τ*) after the application of the inverse *Q* filter.

To apply the filter *W*(*ω, τ*) to a trace *P*(*ω, τ* = 0), define a time step Δ*τ* and write equation (**B-119**) in its recursive form:

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle P(\omega,\ k\Delta \tau)=P[\omega,\ (k-1)\Delta \tau]\ \ \exp(-i\omega\Delta \tau)W(\omega,\ \Delta \tau),}****(**)

where *k* = 1, 2, …, *n*, with *n* number of time steps (number of samples in the input trace).

Equation (**B-120**) can now be used to describe a procedure for inverse *Q* filtering:

- Fourier transform the input trace
*P*(*t, τ*= 0) to obtain the complex transform function*P*(*ω, τ*= 0). - Define a time step Δ
*τ*and apply the linear phase shift to*P*(*ω, τ*= 0) by multiplying with the exponential exp(−*iω*Δ*τ*). - Specify a constant
*Q*and apply the inverse*Q*filter given by equation (**B-117a**,**B-117b**) that represents the inverse*Q*filter. - Repeat steps (a), (b), and (c) for all frequencies.
- Sum over all frequencies to obtain the inverse
*Q*-filtered wavefield at time step Δ*τ*given by*P*(*t*= 0, Δ*τ*). - Repeat step (d) for all time steps
*k*Δ*τ*,*k*= 1, 2, …,*n*, to obtain the inverse*Q*-filtered wavefield*P*(*t*= 0,*τ*) at all times*τ*.

Inverse *Q* filtering often is applied to data using a constant *Q* factor. An efficient scheme for a vertically varying *Q* factor is described by ^{[12]}.

Figure B-5 shows an example of inverse *Q* filtering applied to field data. Compare the stacked sections and the average amplitude spectra with no deconvolution, inverse *Q* filtering, inverse *Q* filtering followed by deconvolution and deconvolution, only. An inverse *Q* filter restores the high-frequency components of signal subjected to intrinsic attenuation by the propagation medium. The deconvolution that follows the inverse *Q* filtering then easily flattens the spectrum within the passband.

## References

- ↑
^{1.0}^{1.1}^{1.2}^{1.3}^{1.4}^{1.5}Claerbout, 1976, Claerbout, J. F., 1976, Fundamentals of geophysical data processing: McGraw-Hill Book Co. - ↑ Robinson and Treitel, 1980, Robinson, E. A. and Treitel, S., 1980, Geophysical signal analysis: Prentice-Hall Book Co.
- ↑ Taner and Coburn, 1981, Taner, M. T. and Coburn, K., 1981, Surface-consistent deconvolution: Presented at the 51st Ann. Internat. Mtg., Soc. Expl. Geophys.
- ↑ Cambois and Stoffa (1992), Cambois, G. and Stoffa, P. L., 1992, Surface-consistent deconvolution in the log-Fourier domain: Geophysics, 57, 823–840.
- ↑
^{5.0}^{5.1}Cary and Lorentz (1993), Cary, P. W. and Lorentz, G. A., 1993, Four-component surface-consistent deconvolution: Geophysics, 58, 383–392. - ↑ Morley and Claerbout (1983), Morley, L. and Claerbout, J. F., 1983, Predictive deconvolution in shot-receiver space: Geophysics, 48, 515–531.
- ↑ 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.
- ↑ Levin (1989), Levin, S. A., 1989, Surface-consistent deconvolution: Geophysics, 54, 1123–1133.
- ↑ Taner and Koehler, 1981, Taner, M. T. and Koehler, F., Surface-consistent corrections: Geophysics, 46, 17–22.
- ↑ Kjartansson, 1979, Kjartansson, E., 1979, Constant
*Q*-wave propagation and attenuation: J. Geophys. Res., 84, 4737–4748. - ↑ Saatcilar, 1996, Saatcilar, R., 1996, An algorithm for
*Q*-filtering: J. Seis. Expl., 5, 157–168. - ↑ Hargreaves and Calvert, 1991, Hargreaves, N. D. and Calvert, A. J., 1991, Inverse
*Q*filtering by Fourier transform: Geophysics, 56, 519–527.

## See also

- Introduction to deconvolution
- The convolutional model
- Inverse filtering
- Optimum Wiener filters
- Predictive deconvolution in practice
- Field data examples
- The problem of nonstationarity