User:Ifarley/Yilmaz

Seismic Data Analysis

Seismic-data-analysis.jpg
Series	Investigations in Geophysics
Title	Seismic Data Analysis
Author	Öz Yilmaz
DOI	http://dx.doi.org/10.1190/1.9781560801580
ISBN	ISBN 978-1-56080-094-1
Store	SEG Online Store

2.0 Introduction

Deconvolution compresses the basic wavelet in the recorded seismogram, attenuates reverberations and short-period multiples, thus increases temporal resolution and yields a representation of subsurface reflectivity. The process normally is applied before stack; however, it also is common to apply deconvolution to stacked data. Figure 2.0-1 shows a stacked section with and without deconvolution. Deconvolution has produced a section with a much higher temporal resolution. The ringy character of the stack without deconvolution limits resolution, considerably.

Figure 2.0-2 shows selected common-midpoint (CMP) gathers from a marine line before and after deconvolution. Note that the prominent reflections stand out more distinctly on the deconvolved gathers. Deconvolution has removed a considerable amount of ringyness, while it has compressed the waveform at each of the prominent reflections. The stacked sections associated with these CMP gathers are shown in Figure 2.0-3. The improvement observed on the deconvolved CMP gathers also are noted on the corresponding stacked section.

Figure 2.0-4 shows some NMO-corrected CMP gathers from a land line with and without deconvolution. Corresponding stacked sections are shown in Figure 2.0-5. Again, note that deconvolution has compressed the wavelet and removed much of the reverberating energy.

Deconvolution sometimes does more than just wavelet compression; it can remove a significant part of the multiple energy from the section. Note that the stacked section in Figure 2.0-6 shows a marked improvement between 2 and 4 s after deconvolution.

To understand deconvolution, first we need to examine the constituent elements of a recorded seismic trace (The convolutional model). The earth is composed of layers of rocks with different lithology and physical properties. Seismically, rock layers are defined by the densities and velocities with which seismic waves propagate through them. The product of density and velocity is called seismic impedance. The impedance contrast between adjacent rock layers causes the reflections that are recorded along a surface profile. The recorded seismogram can be modeled as a convolution of the earth’s impulse response with the seismic wavelet. This wavelet has many components, including source signature, recording filter, surface reflections, and receiver-array response. The earth’s impulse response is what would be recorded if the wavelet were just a spike. The impulse response comprises primary reflections (reflectivity series) and all possible multiples.

Ideally, deconvolution should compress the wavelet components and eliminate multiples, leaving only the earth’s reflectivity in the seismic trace. Wavelet compression can be done using an inverse filter as a deconvolution operator. An inverse filter, when convolved with the seismic wavelet, converts it to a spike (Inverse filtering). When applied to a seismogram, the inverse filter should yield the earth’s impulse response. An accurate inverse filter design is achieved using the least-squares method (Inverse filtering).

The fundamental assumption underlying the deconvolution process (with the usual case of unknown source wavelet) is that of minimum phase. This issue is dealt with also in Inverse filtering.

The optimum Wiener filter, which has a wide range of applications, is discussed in Optimum wiener filters. The Wiener filter converts the seismic wavelet into any desired shape. For example, much like the inverse filter, a Wiener filter can be designed to convert the seismic wavelet into a spike. However, the Wiener filter differs from the inverse filter in that it is optimal in the least-squares sense. Also, the resolution (spikiness) of the output can be controlled by designing a Wiener prediction error filter — the basis for predictive deconvolution (Optimum wiener filters). Converting the seismic wavelet into a spike is like asking for a perfect resolution. In practice, because of noise in the seismogram and assumptions made about the seismic wavelet and the recorded seismogram, spiking deconvolution is not always desirable. Finally, the prediction error filter can be used to remove periodic components — multiples, from the seismogram. Practical aspects of predictive deconvolution are presented in Predictive deconvolution in practice, and field data examples are provided in Field data examples. Finally, time-varying aspects of the source waveform — nonstationarity, are discussed in The problem of nonstationarity.

• 

  Figure 2.0-1  Interpreters prefer the crisp, finely detailed appearance of the deconvolved section (right) as opposed to the blurred, ringy appearance of the section without deconvolution (left). (Data courtesy Enterprise Oil.)

• 

  Figure 2.0-2  Note the prominent reflections on the deconvolved gathers (b). The reverberations would make it difficult to distinguish prominent reflections on the undeconvolved gathers (a).

• 

  Figure 2.0-3  (a) The section obtained from the undeconvolved gathers of Figure 2.0-2a, and (b) the section obtained from the deconvolved gathers of Figure 2.0-2b.

The mathematical treatment of deconvolution is found in Appendix B. However, several numerical examples, which provide the theoretical groundwork from a heuristic viewpoint, are given in the text. Much of the early theoretical work on deconvolution came from the MIT Geophysical Analysis Group, which was formed in the mid-1950s.

2.1 The convolutional model

A sonic log segment is shown in Figure 2.1-1a. The sonic log is a plot of interval velocity as a function of depth based on downhole measurement using logging tools. Here, velocities were measured between the 1000- to 5400-ft depth interval at 2-ft intervals. The velocity function was extrapolated to the surface by a linear ramp. The sonic log exhibits a strong low-frequency component with a distinct blocky character representing gross velocity variations. Actually, it is this low-frequency component that normally is estimated by velocity analysis of CMP gathers (Velocity analysis).

In many sonic logs, the low-frequency component is an expression of the general increase of velocity with depth due to compaction. In some sonic logs, however, the low-frequency component exhibits a blocky character (Figure 2.1-1a), which is due to large-scale lithologic variations. Based on this blocky character, we may define layers of constant interval velocity (Table 2-1), each of which can be associated with a geologic formation (Table 2-2).

Table 2-1. The interval velocity trend obtained from the sonic log in Figure 2.1-1a.

Layer Number	Interval Velocity,* ft/s	Depth Range, ft
1	21 000	1000 – 2000
2	19 000	2000 – 2250
3	18 750	2250 – 2500
4	12 650	2500 – 3775
5	19 650	3775 – 5400

* The velocity in Layer 2 gradually decreases from the top of the layer to the bottom.

Table 2-2. Stratigraphic identification associated with layering described in Table 2-1.

Layer Number	Lithologic Unit
1	Limestone
2	Shaly limestone with gradual increase in shale content
3	Shaly limestone
4	Sandstone
5	Dolomite

The sonic log also has a high-frequency component superimposed on the low-frequency component. These rapid fluctuations can be attributed to changes in rock properties that are local in nature. For example, the limestone layer can have interbeddings of shale and sand. Porosity changes also can affect interval velocities within a rock layer. Note that well-log measurements have a limited accuracy; therefore, some of the high-frequency variations, particularly those associated with a first arrival that is strong enough to trigger one receiver but not the other in the log tool (cycle skips), are not due to changes in lithology.

Well-log measurements of velocity and density provide a link between seismic data and the geology of the substrata. We now explain the link between log measurements and the recorded seismic trace provided by seismic impedance — the product of density and velocity. The first set of assumptions that is used to build the forward model for the seismic trace follows:

Assumption 1. The earth is made up of horizontal layers of constant velocity.

Assumption 2. The source generates a compressional plane wave that impinges on layer boundaries at normal incidence. Under such circumstances, no shear waves are generated.

Assumption 1 is violated in both structurally complex areas and in areas with gross lateral facies changes. Assumption 2 implies that our forward model for the seismic trace is based on zero-offset recording — an unrealizable experiment. Nevertheless, if the layer boundaries were deep in relation to cable length, we assume that the angle of incidence at a given boundary is small and ignore the angle dependence of reflection coefficients. Combination of the two assumptions thus imply a normal-incidence one-dimensional (1-D) seismogram.

• 

  Figure 2.0-4  Some NMO-corrected gathers associated with the stacked sections in Figure 2.0-5, (a) before, (b) after deconvolution. Deconvolution has removed the ringy character from the data.

• 

  Figure 2.0-5  Deconvolution helps distinguish prominent reflections with ease (b). However, on a section without deconvolution (a), reflections are buried in reverberating energy. Selected CMP gathers for both sections are shown in Figure 2.0-4.

• 

  Figure 2.0-6  CMP stacks with (a) no deconvolution before stack, (b) spiking deconvolution before stack. Deconvolution can remove a significant amount of multiple energy from seismic data. (Data courtesy Elf Aquitane and partners.)

• 

  Figure 2.1-1  (a) A segment of a measured sonic log, (b) the reflection coefficient series derived from (a), (c) the series in (b) after converting the depth axis to two-way time axis, (d) the impulse response that includes the primaries (c) and multiples, (e) the synthetic seismogram derived from (d) convolved with the source wavelet in Figure 2.1-4. One-dimensional seismic modeling means getting (e) from (a). Deconvolution yields (d) from (e), while 1-D inversion means getting (a) from (d). Identify the event on (a) and (b) that corresponds to the big spike at 0.5 s in (c). Impulse response (d) is a composite of the primaries (c) and all types of multiples.

• 

  Figure 2.1-2  A seismic source wavelet after onset takes the form shown at top left. As the wavelet travels into the earth, the amplitude level drops (geometric spreading) and a loss of high frequencies occurs (frequency absorption).

Based on assumptions 1 and 2, the reflection coefficient c (for pressure or stress), which is associated with the boundary between, say, layers 1 and 2, is defined as

$c={\frac {I_{2}-I_{1}}{I_{2}+I_{1}}},$

(1a)

where I is the seismic impedance associated with each layer given by the product of density ρ and compressional velocity v.

From well-log measurements, we find that the vertical density gradient often is much smaller than the vertical velocity gradient. Therefore, we often assume that the impedance contrast between rock layers is essentially due to the velocity contrast, only. Equation (1a) then takes the form:

$c={\frac {v_{2}-v_{1}}{v_{2}+v_{1}}}.$

(1b)

If v₂ is greater than v₁, the reflection coefficient would be positive. If v₂ is less than v₁, then the reflection coefficient would be negative.

The assumption that density is invariant with depth or that it does not vary as much as velocity is not always valid. The reason we can get away with it is that the density gradient usually has the same sign as the velocity gradient. Hence, the impedance function derived from the velocity function only should be correct within a scale factor.

For vertical incidence, the reflection coefficient is the ratio of the reflected wave amplitude to the incident wave amplitude. Moreover, from its definition (equation 1a), the reflection coefficient is seen as the ratio of the change in acoustic impedance to twice the average acoustic impedance. Therefore, seismic amplitudes associated with earth models with horizontal layers and vertical incidence (assumptions 1 and 2) are related to acoustic impedance variations.

The reflection coefficient series c(z), where z is the depth variable, is derived from sonic log v(z) and is shown in Figure 2.1-1b. We note the following:

1. The position of each spike gives the depth of the layer boundary, and
2. the magnitude of each spike corresponds to the fraction of a unit-amplitude downward-traveling incident plane wave that would be reflected from the layer boundary.

To convert the reflection coefficient series c(z) (Figure 2.1-1b) derived from the sonic log into a time series c(t), select a sampling interval, say 2 ms. Then use the velocity information in the log (Figure 2.1-1a) to convert the depth axis to a two-way vertical time axis. The result of this conversion is shown in Figure 2.1-1c, both as a conventional wiggle trace and as a variable area and wiggle trace (the same trace repeated six times to highlight strong reflections). The reflection coefficient series c(t) (Figure 2.1-1c) represents the reflectivity of a series of fictitious layer boundaries that are separated by an equal time interval — the sampling rate . The major events in this reflectivity series are from the boundary between layers 2 and 3 located at about 0.3 s, and the boundary between layers 4 and 5 located at about 0.5 s.

The reflection coefficient series (Figure 2.1-1c) that was constructed is composed only of primary reflections (energy that was reflected only once). To get a complete 1-D response of the horizontally-layered earth model (assumption 1), multiple reflections of all types (surface, intrabed and interbed multiples) must be included. If the source were unit-amplitude spike, then the recorded zero-offset seismogram would be the impulse response of the earth, which includes primary and multiple reflections. Here, the Kunetz method is used to obtain such an impulse response. The impulse response derived from the reflection coefficient series in Figure 2.1-1c is shown in Figure 2.1-1d with the variable area and wiggle display.

The characteristic pressure wave created by an impulsive source, such as dynamite or air gun, is called the signature of the source. All signatures can be described as band-limited wavelets of finite duration — for example, the measured signature of an Aquapulse source in Figure 2.1-2. As this waveform travels into the earth, its overall amplitude decays because of wavefront divergence. Additionally, frequencies are attenuated because of the absorption effects of rocks (see Gain applications). The progressive change of the source wavelet in time and depth also is shown in Figure 2.1-2. At any given time, the wavelet is not the same as it was at the onset of source excitation. This time-dependent change in waveform is called nonstationarity.

Wavefront divergence is removed by applying a spherical spreading function (The 2-D Fourier transform). Frequency attenuation is compensated for by the processing techniques discussed in The problem of nonstationarity. Nevertheless, the simple convolutional model discussed here does not incorporate nonstationarity. This leads to the following assumption:

Assumption 3. The source waveform does not change as it travels in the subsurface — it is stationary.

The convolutional model in the time domain

A convolutional model for the recorded seismogram now can be proposed. Suppose a vertically propagating downgoing plane wave with source signature (Figure 2.1-3a) travels in depth and encounters a layer boundary at 0.2-s two-way time. The reflection coefficient associated with the boundary is represented by the spike in Figure 2.1-3b. As a result of reflection, the source wavelet replicates itself such that it is scaled by the reflection coefficient. If we have a number of layer boundaries represented by the individual spikes in Figures 2.1-3b through 2.1-3f, then the wavelet replicates itself at those boundaries in the same manner. If the reflection coefficient is negative, then the wavelet replicates itself with its polarity reversed, as in Figure 2.1-3c.

Now consider the ensemble of the reflection coefficients in Figure 2.1-3g. The response of this sparse spike series to the basic wavelet is a superposition of the individual impulse responses. This linear process is called the principle of superposition. It is achieved computationally by convolving the basic wavelet with the reflectivity series (Figure 2.1-3g). The convolutional process already was demonstrated by the numerical example in The 1-D Fourier transform.

The response of the sparse spike series to the basic wavelet in Figure 2.1-3g has some important characteristics. Note that for events at 0.2 and 0.35 s, we identify two layer boundaries. However, to identify the three closely spaced reflecting boundaries from the composite response (at around 0.6 s), the source waveform must be removed to obtain the sparse spike series. This removal process is just the opposite of the convolutional process used to obtain the response of the reflectivity series to the basic wavelet. The reverse process appropriately is called deconvolution.

• 

  Figure 2.1-3  A wavelet (a) traveling in the earth repeats itself when it encounters a reflector along its path (b, c, d, e, f). The left column represents the reflection coefficients, while the right column represents the response to the wavelet. Amplitudes of the response are scaled by the reflection coefficient. The resulting seismogram (bottom right) represents the composite response of the earth’s reflectivity (bottom left) to the wavelet (top right).

• 

  Figure 2.1-4  The top frame is the same as in Figure 2.1-1d. The asterisk denotes convolution. The recorded seismogram (bottom frame) is the sum of the noise-free seismogram and the noise trace. This figure is equivalent to equation (2a).

• 

  Figure 2.1-5  A random signal with infinite length has a flat amplitude spectrum and an autocorrelogram that is zero at all lags except the zero lag. The discrete random series with finite length shown here seems to satisfy these requirements. What distinguishes a random signal from a spike (1, 0, 0, …)?

The principle of superposition now is applied to the impulse response derived from the sonic log in Figure 2.1-1d. Convolution of a source signature with the impulse response yields the synthetic seismogram shown in Figure 2.1-4. The synthetic seismogram also is shown in Figure 2.1-1e. This 1-D zero-offset seismogram is free of random ambient noise. For a more realistic representation of a recorded seismogram, noise is added (Figure 2.1-4).

The convolutional model of the recorded seismogram now is complete. Mathematically, the convolutional model illustrated in Figure 2.1-4 is given by

${x(t)=w(t)\ast e(t)+n(t)},$

(2a)

where x(t) is the recorded seismogram, w(t) is the basic seismic wavelet, e(t) is the earth’s impulse response, n(t) is the random ambient noise, and * denotes convolution. Deconvolution tries to recover the reflectivity series (strictly speaking, the impulse response) from the recorded seismogram.

An alternative to the convolutional model given by equation (2a) is based on a surface-consistent spectral decomposition . In such a formulation, the seismic trace is decomposed into the convolutional effects of source, receiver, offset, and the earth’s impulse response, thus explicitly accounting for variations in wavelet shape caused by near-source and near-receiver conditions and source-receiver separation. The following equation describes the surface-consistent convolutional model (Section B.8):

${x\prime _{ij}(t)=s_{j}(t)\ast h_{l}(t)\ast e_{k}(t)\ast g_{i}(t)+n(t)},$

(2b)

where ${{{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 (2a), e_k(t) represents the earth’s impulse response at the source-receiver midpoint location, k = (i + j)/2. By comparing equations (2a) and (2b), we infer that w(t) represents the combined effects of s(t), h(t), and g(t).

The assumption of surface-consistency implies that the basic wavelet shape depends only on the source and receiver locations, not on the details of the raypath from source to reflector to receiver. In a transition zone, surface conditions at the source and receiver locations may vary significantly from dry to wet surface conditions. Hence, the most likely situation where the surface-consistent convolutional model may be applicable is with transition-zone data. Nevertheless, the formulation described in this section is the most accepted model for the 1-D seismogram.

The random noise present in the recorded seismogram has several sources. External sources are wind motion, environmental noise, or a geophone loosely coupled to the ground. Internal noise can arise from the recording instruments. A pure-noise seismogram and its characteristics are shown in Figure 2.1-5. A pure random-noise series has a white spectrum — it contains all the frequencies. This means that the autocorrelation function is a spike at zero lag and zero at all other lags. From Figure 2.1-5, note that these characteristic requirements are reasonably satisfied.

Now examine the equation for the convolutional model. All that normally is known in equation (2a) is x(t) — the recorded seismogram. The earth’s impulse response e(t) must be estimated everywhere except at the location of wells with good sonic logs. Also, the source waveform w(t) normally is unknown. In certain cases, however, the source waveform is partly known; for example, the signature of an air-gun array can be measured. However, what is measured is only the waveform at the very onset of excitation of the source array, and not the wavelet that is recorded at the receiver. Finally, there is no a priori knowledge of the ambient noise n(t).

We now have three unknowns — w(t), e(t), and n(t), one known — x(t), and one single equation (2a). Can this problem be solved? Pessimists would say no. However, in practice, deconvolution is applied to seismic data as an integral part of conventional processing and is an effective method to increase temporal resolution.

To solve for the unknown e(t) in equation (2a), further assumptions must be made.

Assumption 4. The noise component n(t) is zero.

Assumption 5. The source waveform is known.

Under these assumptions, we have one equation,

${x(t)=w(t)\ast e(t)}.$

(3a)

and one unknown, the reflectivity series e(t). In reality, however, neither of the above two assumptions normally is valid. Therefore, the convolutional model is examined further in the next section, this time in the frequency domain, to relax assumption 5.

If the source waveform were known (such as the recorded source signature), then the solution to the deconvolution problem is deterministic. In Inverse filtering, one such method of solving for e(t) is considered. If the source waveform were unknown (the usual case), then the solution to the deconvolution problem is statistical. The Wiener prediction theory (Optimum wiener filters) provides one method of statistical deconvolution.

The convolutional model in the frequency domain

The convolutional model for the noise-free seismogram (assumption 4) is represented by equation (3a). Convolution in the time domain is equivalent to multiplication in the frequency domain (Section A.1). This means that the the amplitude spectrum of the seismogram equals the product of the amplitude spectra of the seismic wavelet and the earth’s impulse response (Section B.1):

${{A}_{x}}\left(\omega \right)={{A}_{w}}\left(\omega \right){{A}_{e}}\left(\omega \right),$

(3b)

where A_x(ω), A_w(ω), and A_e(ω) are the amplitude spectra of x(t), w(t), and e(t), respectively.

Figure 2.1-6 shows the amplitude spectra (top row) of the impulse response e(t), the seismic wavelet w(t), and the seismogram x(t). The impulse response is the same as that shown in Figure 2.1-1d. The similarity in the overall shape between the amplitude spectrum of the wavelet and that of the seismogram is apparent. In fact, a smoothed version of the amplitude spectrum of the seismogram is nearly indistinguishable from the amplitude spectrum of the wavelet. It generally is thought that the rapid fluctuations observed in the amplitude spectrum of a seismogram are a manifestation of the earth’s impulse response, while the basic shape is associated primarily with the source wavelet.

Mathematically, the similarity between the amplitude spectra of the seismogram and the wavelet suggests that the amplitude spectrum of the earth’s impulse response must be nearly flat (Section B.1). By examining the amplitude spectrum of the impulse response in Figure 2.1-6, we see that it spans virtually the entire spectral bandwidth. As seen in Figure 2.1-5, a time series that represents a random process has a flat (white) spectrum over the entire spectral bandwidth. From close examination of the amplitude spectrum of the impulse response in Figure 2.1-6, we see that it is not entirely flat — the high-frequency components have a tendency to strengthen gradually. Thus, reflectivity is not entirely a random process. In fact, this has been observed in the spectral properties of reflectivity functions derived from a worldwide selection of sonic logs .

We now study the autocorrelation functions (middle row, Figure 2.1-6) of the impulse response, seismic wavelet, and synthetic seismogram. Note that the autocorrelation functions of the basic wavelet and seismogram also are similar. This similarity is confined to lags for which the autocorrelation of the wavelet is nonzero. Mathematically, the similarity between the autocorrelogram of the wavelet and that of the seismogram suggests that the impulse response has an autocorrelation function that is small at all lags except the zero lag (Section B.1). The autocorrelation function of the random series in Figure 2.1-5 also has similar characteristics. However, there is one subtle difference. When compared, Figures 2.1-5 and 2.1-6 show that autocorrelation of the impulse response has a significantly large negative lag value following the zero lag. This is not the case for the autocorrelation of random noise. The positive peak (zero lag) followed by the smaller negative peak in the autocorrelogram of the impulse response arises from the spectral behavior discussed above. In particular, the positive peak and the adjacent, smaller negative peak of the autocorrelogram together nearly act as a fractional derivative operator (Section A.1), which has a ramp effect on the amplitude spectrum of the impulse response as seen in Figure 2.1-6.

Figure 2.1-6  Convolution of the earth’s impulse response (a) with the wavelet (b) (equation 2a) yields the seismogram (c) (bottom row). This process also is convolutional in terms of their autocorrelograms (middle row) and multiplicative in terms of their amplitude spectra (top row). Assumption 6 (white reflectivity) is based on the similarity between autocorrelograms and amplitude spectra of the impulse response and wavelet.

The above observations made on the amplitude spectra and autocorrelation functions (Figure 2.1-6) imply that reflectivity is not entirely a random process. Nonetheless, the following assumption almost always is made about reflectivity to replace the statement made in assumption 5.

Assumption 6. Reflectivity is a random process. This implies that the seismogram has the characteristics of the seismic wavelet in that their autocorrelations and amplitude spectra are similar.

This assumption is the key to implementing the predictive deconvolution. It allows the autocorrelation of the seismogram, which is known, to be substituted for the autocorrelation of the seismic wavelet, which is unknown. In Optimum wiener filters, we shall see that as a result of assumption 6, an inverse filter can be estimated directly from the autocorrelation of the seismogram. For this type of deconvolution, Assumption 5, which is almost never met in reality, is not required. But first, we need to review the fundamentals of inverse filtering.

2.2 Inverse filtering

If a filter operator f(t) were defined such that convolution of f(t) with the known seismogram x(t) yields an estimate of the earth’s impulse response e(t), then

${e(t)=f(t)\ast x(t)}.$

(4)

By substituting equation (4) into equation (3a), we get

${x(t)=w(t)\ast f(t)\ast x(t)}.$

(5)

When x(t) is eliminated from both sides of the equation, the following expression results:

${\delta (t)=w(t)\ast f(t)},$

(6)

where δ(t) represents the Kronecker delta function:

$\delta (t)=\left\{{\begin{array}{ll}1,&t=0,\\0,&{\mbox{otherwise}}.\end{array}}\right.$

(7)

By solving equation (6) for the filter operator f(t), we obtain

$f(t)=\delta (t)\ast {\frac {1}{w(t)}}.$

(8)

Therefore, the filter operator f(t) needed to compute the earth’s impulse response from the recorded seismogram turns out to be the mathematical inverse of the seismic wavelet w(t). Equation (8) implies that the inverse filter converts the basic wavelet to a spike at t = 0. Likewise, the inverse filter converts the seismogram to a series of spikes that defines the earth’s impulse response. Therefore, inverse filtering is a method of deconvolution, provided the source waveform is known (deterministic deconvolution). The procedure for inverse filtering is described in Figure 2.2-1.

Figure 2.2-1  A flowchart for inverse filtering.

The inverse of the source wavelet

Computation of the inverse of the source wavelet is accomplished mathematically by using the z-transform (Section A.2). For example, let the basic wavelet be a two-point time series given by w(t) : (1, -  1/2). The z-transform of this wavelet is defined by the following polynomial:

${W(z)=1-{\frac {1}{2}}z}.$

(9)

The power of variable z is the number of unit time delays associated with each sample in the series. The first term has zero delay, so z is raised to zero power. The second term has unit delay, so z is raised to first power. Hence, the z-transform of a time series is a polynomial in z, whose coefficients are the values of the time samples.

A relationship exists between the z-transform and the Fourier transform (Section A.2). The z-variable is defined as

${z=\exp(-iw\Delta t)},$

(10)

where ω is angular frequency and Δt is sampling interval.

The convolutional relation in the time domain given by equation (8) means that the z-transform of the inverse filter, F(z), is obtained by polynomial division of the z-transform of the input wavelet, W(z), given by equation (9) (Section A.2):

${F(z)={\frac {1}{1-{\frac {1}{2}}z}}=1+{\frac {1}{2}}z+{\frac {1}{4}}z^{2}+\ldots }.$

(11)

The coefficients of $F\left(z\right):\left(1,\ {\frac {1}{2}},\ {\frac {1}{4}},\ \ldots \right)$ represent the time series associated with the filter operator f(t). Note that the series has an infinite number of coefficients, although they decay rapidly. As in any filtering process, in practice the operator is truncated.

First consider the first two terms in equation (11) which yield a two-point filter operator (1, 1/2). The design and application of this operator is summarized in Table 2-3. The actual output is (1, 0, -  1/4), whereas the ideal result is a zero-delay spike (1, 0, 0). Although not ideal, the actual result is spikier than the input wavelet, (1, -  1/2).

Can the result be improved by including one more coefficient in the inverse filter? As shown in Table 2-4, the actual output from the three-point filter is (1, 0, 0, -  1/8). This is a more accurate representation of the desired output (1, 0, 0, 0) than that achieved with the output from the two-point filter (Table 2-3). Note that there is less energy leaking into the nonzero lags of the output from the three-point filter. Therefore, it is spikier. As more terms are included in the inverse filter, the output is closer to being a spike at zero lag. Since the number of points allowed in the operator length is limited, in practice the result never is a perfect spike.

Table 2-3. Design and application of the truncated inverse filter (1, 1/2) with the input wavelet (1, -  1/2).

Filter Design
Input Wavelet					$w\left(t\right):\left(1,\ -{\frac {1}{2}}\right)$
The z-Transform					$W\left(z\right)=1-\ {\frac {1}{2}}z$
The Inverse					$F\left(z\right)=1\ +\ {\frac {1}{2}}z\ +\ {\frac {1}{4}}z^{2}\ +\ \cdots$
The Inverse Filter					$f\left(t\right):\left(1,\ {\frac {1}{2}},\ {\frac {1}{4}},\ \cdots \right)$
Filter Application
Truncated Inverse Filter						$\left(1,\ {\frac {1}{2}}\right)$
Input Wavelet						$\left(1,\ -{\frac {1}{2}}\right)$
Actual Output						$\left(1,\ 0,\ -{\frac {1}{4}}\right)$
Desired Output						(1, 0, 0)
Convolution Table:
		1	-  1/2			Output
	1/2	1				1
		1/2	1			0
			1/2	1		-  1/4

Table 2-4. Design and application of the truncated inverse filter (1, 1/2, 1/4) with the input wavelet (1, -  1/2).

Filter Design
Input Wavelet						$w\left(t\right):\left(1,\ -{\frac {1}{2}}\right)$
The z-Transform						$W\left(z\right)=1-\ {\frac {1}{2}}z$
The Inverse						$F\left(z\right)=1\ +\ {\frac {1}{2}}z\ +\ {\frac {1}{4}}z^{2}\ +\ \cdots$
The Inverse Filter						$f\left(t\right):\left(1,\ {\frac {1}{2}},\ {\frac {1}{4}},\ \cdots \right)$
Filter Application
Truncated Inverse Filter							$\left(1,\ {\frac {1}{2}}\ {\frac {1}{4}}\right)$
Input Wavelet							$\left(1,\ -{\frac {1}{2}}\right)$
Actual Output							$\left(1,\ 0,\ 0,\ -{\frac {1}{8}}\right)$
Desired Output							(1, 0, 0, 0)
Convolution Table:
			1	-  1/2			Output
	1/4	1/2	1				1
		1/4	1/2	1			0
			1/4	1/2	1		0
				1/4	1/2	1	-  1/8

The inverse of the input wavelet $w\left(t\right):\left(1,\ -{\frac {1}{2}}\right)$ has coefficients that rapidly decay to zero (equation 11). What about the inverse of the input wavelet $w\left(t\right):\left(-{\frac {1}{2}},\ 1\right)?$ Again, define the z–transform:

${W(z)=-{\frac {1}{2}}+z}.$

(12)

The z–transform of its inverse is given by the polynomial division:

${F(z)={\frac {1}{-{\frac {1}{2}}+z}}=-2-4z-8z^{2}-\cdots }$

(13)

As a result, the inverse filter coefficients are given by the divergent series f(t) : (−2, −4, −8, …). Truncate this series and convolve the two-point operator with the input wavelet (-  1/2, 1) as shown in Table 2-5. The actual output is (1, 0, −4), while the desired output is (1, 0, 0). Not only is the result far from the desired output, but also it is less spiky than the input wavelet (-  1/2, 1). The reason for this poor result is that the inverse filter coefficients increase in time rather than decay (equation 13). When truncated, the larger coefficients actually are excluded from the computation.

If we kept the third coefficient of the inverse filter in the above example (equation 13), then the actual output (Table 2-6) would be (1, 0, 0, −8), which also is a bad approximation to the desired output (1, 0, 0, 0).

Table 2-5. Design and application of the truncated inverse filter (−2, −4) with input wavelet (-  1/2, 1).

Filter Design
Input Wavelet					$w\left(t\right):\left(-{\frac {1}{2}},\ 1\right)$
The z-Transform					$W\left(z\right)=-{\frac {1}{2}}+z$
The Inverse					$F\left(z\right)=-2\ -\ 4z\ -\ 8z^{2}\ -\ \cdots$
The Inverse Filter					f(t): (−2, −4, −8, …)
Filter Application
Truncated Inverse Filter						(−2, −4)
Input Wavelet						$\left(-{\frac {1}{2}},\ 1\right)$
Actual Output						(1, 0, −4)
Desired Output						(1, 0, 0)
Convolution Table:
		-  1/2	1			Output
	−4	−2				1
		−4	−2			0
			−4	−2		−4

Table 2-6. Design and application of the truncated inverse filter (−2, −4, −8) with input wavelet (-  1/2, 1).

Filter Design
Input Wavelet					$w\left(t\right):\left(-{\frac {1}{2}},\ 1\right)$
The z-Transform					$W\left(z\right)=-{\frac {1}{2}}+z$
The Inverse					$F\left(z\right)=-2\ -\ 4z\ -\ 8z^{2}\ -\ \cdots$
The Inverse Filter					f(t): (−2, −4, −8, …)
Filter Application
Truncated Inverse Filter							(−2, −4, −8)
Input Wavelet							$\left(-{\frac {1}{2}},\ 1\right)$
Actual Output							(1, 0, 0, −8)
Desired Output							(1, 0, 0, 0)
Convolution Table:
			-  1/2	1			Output
	−8	−4	−2				1
		−8	−4	−2			0
			−8	−4	−2		0
				−8	−4	−2	−8

Least-squares inverse filtering

A well-behaved input wavelet, such as (1, -  1/2) as opposed to (-  1/2, 1), has a z-transform whose inverse can be represented by a convergent series. Then the inverse filtering described above yields a good approximation to a zero-lag spike output (1, 0, 0). Can we do even better than that?

Formulate the following problem: Given the input wavelet (1, -  1/2), find a two-term filter (a, b) such that the error between the actual output and the desired output (1, 0, 0) is minimum in the least-squares sense.

Compute the actual output by convolving the filter (a, b) with the input wavelet (1, -  1/2) (Table 2-7). The cumulative energy of the error L is defined as the sum of the squares of the differences between the coefficients of the actual and desired outputs:

${L=(a-1)^{2}+(b-{\frac {a}{2}})^{2}+(-{\frac {b}{2}})^{2}}.$

(14)

The task is to find coefficients (a, b) so that L takes its minimum value. This requires variation of L with respect to the coefficients (a, b) to vanish (Section B.5). By simplifying equation (14), taking the partial derivatives of quantity L with respect to a and b, and setting the results to zero, we get

${{\frac {5}{2}}a-b=2},$

(15a)

and

$-a+{{\frac {5}{2}}b=0}.$

(15b)

We have two equations and two unknowns; namely, the filter coefficients (a, b). The so-called normal set of equations (15a) and (15b) can be put into the following convenient matrix form

${\begin{pmatrix}{5}/{2}&-1\\-1&{5}/{2}\end{pmatrix}}{\begin{pmatrix}a\\b\end{pmatrix}}={\begin{pmatrix}2\\0\end{pmatrix}}.$

(16)

By solving for the filter coefficients, we obtain (a, b): (0.95, 0.38). Design and application of this least-squares inverse filter are summarized in Table 2-7.

To quantify the spikiness of this result and compare it with the result from the inverse filter in Table 2-3, compute the energy of the errors made in both (Table 2-8). Note that the least-squares filter yields less error when trying to convert the input wavelet (1, -  1/2) to a spike at zero lag (1, 0, 0).

We now examine the performance of the least-squares filter with the input wavelet (-  1/2, 1). Note that the inverse filter produced unstable results for this wavelet (Table 2-5). We want to find a two-term filter (a, b) that, when convolved with the input wavelet (-  1/2, 1), yields an estimate of the desired spike output (1, 0, 0) (Table 2-9). As before, the least-squares error between the actual output and the desired output should be minimal.

The cumulative energy of the error is given by

${L=(-{\frac {a}{2}}-1)^{2}+(-{\frac {b}{2}}+a)^{2}+b^{2}}.$

(17)

Table 2-7. Design and application of a two-term least-squares inverse filter (a, b).

Filter Design
Convolution of the filter (a, b) with input wavelet (1, -  1/2):
	1	-  1/2		Actual Output	Desired Output
b	a			a	1
	b	a		b − a/2	0
		b	a	−b/2	0
Filter Application
Least-Squares Filter				(0.95, 0.38)
Input Wavelet				(1, −0.5)
Actual Output				(0.95, −0.09, −0.19)
Desired Output				(1, 0, 0)

Table 2-8. Error in two-term inverse and least-squares filtering.

Input: $(1,\ -{\frac {1}{2}})$
Desired Output: (1, 0, 0)
	Actual Output	Error Energy
Inverse Filter	(1, 0, −0.25)	0.063
Least-Squares Filter	(0.95, −0.09, −0.19)	0.048

Table 2-9. Design and application of a two-term least-squares inverse filter (a, b).

Filter Design
Convolution of the filter (a, b) with input wavelet $(-{\frac {1}{2}},\ 1)$:
	-  1/2	1		Actual Output	Desired Output
b	a			−a/2	1
	b	a		−b/2 + a	0
		b	a	b	0
Filter Application
Least-Squares Filter				(−0.95, −0.19)
Input Wavelet				(−0.5, 1)
Actual Output				(0.24, −0.38, −0.19)
Desired Output				(1, 0, 0)

Table 2-10. Error in two-term inverse and least-squares filtering.

Input: $(-{\frac {1}{2}},\ 1)$
Desired Output: (1, 0, 0)
	Actual Output	Error Energy
Inverse Filter	(1, 0, −4)	16
Least-Squares Filter	(0.24, −0.38, −0.19)	0.792

By simplifying equation (17), taking the partial derivatives of quantity L with respect to a and b, and setting the results to zero, we obtain

${\frac {5}{2}}a-b=-1,$

(18a)

and

$-a+{\frac {5}{2}}b=0.$

(18b)

Combine equations (18a,18b) into a matrix form

${\begin{pmatrix}5/2&-1\\-1&5/2\end{pmatrix}}{\begin{pmatrix}a\\b\end{pmatrix}}={\begin{pmatrix}-1\\0\end{pmatrix}}.$

(19)

By solving for the filter coefficients, we obtain (a, b): (−0.95, −0.19). The design and application of this filter are summarized in Table 2-9.

Table 2-10 shows the results from the inverse filter and least-squares filter quantified. The error made by the least-squares filter is, again, much less than the error made by the truncated inverse filter. However, both filters yield larger errors for input wavelet (-  1/2, 1) (Table 2-10) as compared to errors for wavelet (1, -  1/2) (Table 2-8). The reason for this is discussed next.

Minimum phase

Two input wavelets, wavelet 1: (1, -  1/2) and wavelet 2: (-  1/2, 1), were used for numerical analyses of the inverse filter and least-squares inverse filter in this section. The results indicate that the error in converting wavelet 1 to a zero-lag spike is less than the error in converting wavelet 2 (Tables 2-8 and 2-10).

Is this also true when the desired output is a delayed spike (0, 1, 0)? The cumulative energy of the error L associated with the application of a two-term least-squares filter (a, b) (Table 2-11) to convert the input wavelet (1, -  1/2) to a delayed spike (0, 1, 0) is

${L=a^{2}+\left[\left(b-{\frac {a}{2}}\right)-1\right]^{2}+\left(-{\frac {b}{2}}\right)^{2}}.$

(20)

Table 2-11. Design and application of a two-term least-squares inverse filter (a, b).

Filter Design
Convolution of the filter (a, b) with input wavelet $(1,\ -{\frac {1}{2}})$:
	1	-  1/2		Actual Output	Desired Output
b	a			a	0
	b	a		b − a/2	1
		b	a	−b/2	0
Filter Application
Least-Squares Filter				(−0.09, 0.76)
Input Wavelet				$(1,\ -{\frac {1}{2}})$
Actual Output				(−0.09, 0.81, −0.38)
Desired Output				(0, 1, 0)

Table 2-12. Error in least-squares filtering.

Input Wavelet: $(1,\ -{\frac {1}{2}})$
Desired Output	Actual Output	Error Energy
(1, 0, 0)	(0.95, −0.09, −0.19)	0.048
(0, 1, 0)	(−0.09, 0.81, −0.38)	0.190

By simplifying equation (20), taking the partial derivatives of quantity L with respect to a and b, and setting the results to zero, we obtain

${{\frac {5}{2}}a-b=-1},$

(21a)

and

${-a+{\frac {5}{2}}b=2}.$

(21b)

Combine equations (21a,21b) into a matrix form

${\begin{pmatrix}5/2&-1\\-1&5/2\end{pmatrix}}{\begin{pmatrix}a\\b\end{pmatrix}}={\begin{pmatrix}-1\\2\end{pmatrix}}.$

(22)

By solving for the filter coefficients, we obtain (a, b): (−0.09, 0.76). The design and application of this filter are summarized in Table 2-11.

Table 2-12 shows the results of the least-squares filtering to convert the input wavelet (1, -  1/2) to zero-lag (Table 2-7) and delayed spikes (Table 2-11). Note that the input wavelet is converted to a zero-lag spike with less error, and the corresponding actual output more closely resembles a zero-lag spike desired output.

Table 2-13. Design and application of a two-term least-squares inverse filter (a, b).

Filter Design
Convolution of the filter (a, b) with input wavelet $(-{\frac {1}{2}},\ 1)$:
	-  1/2	1		Actual Output	Desired Output
b	a			−a/2	0
	b	a		−b/2 + a	1
		b	a	b	0
Filter Application
Least-Squares Filter				(0.76, −0.09)
Input Wavelet				(−0.5, 1)
Actual Output				(−0.38, 0.81, −0.09)
Desired Output				(0, 1, 0)

We now examine the performance of the least-squares filter with the input wavelet (-  1/2, 1). The cumulative energy of the error L associated with the application of a two-term least-squares filter (a, b) (Table 2-13) to convert the input wavelet (-  1/2, 1) to a delayed spike (0, 1, 0) is

$L=\left(-{\frac {a}{2}}\right)^{2}=\left[\left(-{\frac {b}{2}}+a\right)-1\right]^{2}+b^{2}.$

(23)

By simplifying equation (23), taking the partial derivatives of quantity L with respect to a and b, and setting the results to zero, we obtain

${{\frac {5}{2}}a-b=2},$

(24a)

and

${-a+{\frac {5}{2}}b=-1}.$

(24b)

Combine equations (24a,24b) into a matrix form

${\begin{pmatrix}5/2&-1\\-1&5/2\end{pmatrix}}{\begin{pmatrix}a\\b\end{pmatrix}}={\begin{pmatrix}2\\-1\end{pmatrix}}.$

(25)

By solving for the filter coefficients, we obtain (a, b): (0.76, −0.09). The design and application of this filter are summarized in Table 2-13.

Table 2-14 shows the results of the least-squares filtering to convert the input wavelet (-  1/2, 1) to zero-lag (Table 2-9) and delayed spikes (Table 2-13). Note that the input wavelet is converted to a delayed spike with less error, and the corresponding actual output more closely resembles a delayed spike desired output.

Table 2-14. Error in least-squares filtering.

Input Wavelet: $(-{\frac {1}{2}},\ 1)$
Desired Output	Actual Output	Error Energy
(1, 0, 0)	(0.24, −0.38, 0.19)	0.762
(0, 1, 0)	(−0.38, 0.81, −0.09)	0.190

Now, evaluate the results of the least-squares inverse filtering summarized in Tables 2-12 and 2-14. Wavelet 1: (1, -  1/2) is closer to being a zero-delay spike (1, 0, 0) than wavelet 2: (-  1/2, 1). On the other hand, wavelet 2 is closer to being a delayed spike (0, 1, 0) than wavelet 1. We conclude that the error is reduced if the desired output closely resembles the energy distribution in the input series. Wavelet 1 has more energy at the onset, while wavelet 2 has more energy concentrated at the end.

Figure 2.2-2 shows three wavelets with the same amplitude spectrum, but with different phase-lag spectra. As a result, their shapes differ. (From the 1-D Fourier transform, we know that the shape of a wavelet can be altered by changing the phase spectrum without modifying the amplitude spectrum.) The wavelet on top has more energy concentrated at the onset, the wavelet in the middle has its energy concentrated at the center, and the wavelet at the bottom has most of its energy concentrated at the end.

We say that a wavelet is minimum phase if its energy is maximally concentrated at its onset. Similarly, a wavelet is maximum phase if its energy is maximally concentrated at its end. Finally, in all in-between situations, the wavelet is mixed phase. Note that a wavelet is defined as a transient waveform with a finite duration — it is realizable. A minimum-phase wavelet is one-sided — it is zero before t = 0. A wavelet that is zero for t < 0 is called causal. These definitions are consistent with intuition — physical systems respond to an excitation only after that excitation. Their response also is of finite duration. In summary, a minimum-phase wavelet is realizable and causal.

These observations are quantified by considering the following four, three-point wavelets :

Wavelet A : (4, 0, -1)

Wavelet B : (2, 3, -2)

Wavelet C : (-2, 3, 2)

Wavelet D : (-1, 0, 4)

Compute the cumulative energy of each wavelet at any one time. Cumulative energy is computed by adding squared amplitudes as shown in Table 2-15. These values are plotted in Figure 2.2-3. Note that all four wavelets have the same amount of total energy — 17 units. However, the rate at which the energy builds up is significantly different for each wavelet. For example, with wavelet A, the energy builds up rapidly close to its total value at the very first time lag. The energy for wavelets B and C builds up relatively slowly. Finally, the energy accumulates at the slowest rate for wavelet D. From Figure 2.2-3, note that the energy curves for wavelets A and D form the upper and lower boundaries. Wavelet A has the least energy delay, while wavelet D has the largest energy delay.

Table 2-15. Cumulative energy of wavelets A, B, C, and D at time samples 0, 1 and 2.

Wavelet	0	1	2
A	16	16	17
B	4	13	17
C	4	13	17
D	1	1	17

• 

  Figure 2.2-2  A wavelet has a finite duration. If its energy is maximally front-loaded, then it is minimum-phase (top). If its energy is concentrated mostly in the middle, then it is mixed-phase (middle). Finally, if its energy is maximally end-loaded, then the wavelet is maximum-phase. A quantitative analysis of this phase concept is provided in Figure 2.2-3.

• 

  Figure 2.2-3  A quantitative analysis of the minimum- and maximum-phase concept. The fastest rate of energy buildup in time occurs when the wavelet is minimum-phase (A). The slowest rate occurs when the wavelet is maximum-phase (D).

• 

  Figure 2.2-4  All wavelets referred to in Figure 2.2-3 (A, B, C, and D) have the same amplitude spectrum as shown above Adapted from .

• 

  Figure 2.2-5  Phase-lag spectra of the wavelets referred to in Figure 2.2-3. They have the common amplitude spectrum of Figure 2.2-4 Adapted from .

Given a fixed amplitude spectrum as in Figure 2.2-4, the wavelet with the least energy delay is called minimum delay, while the wavelet with the most energy delay is called maximum delay. This is the basis for Robinson’s energy delay theorem: A minimum-phase wavelet has the least energy delay.

Time delay is equivalent to a phase-lag. Figure 2.2-5 shows the phase spectra of the four wavelets. Note that wavelet A has the least phase change across the frequency axis; we say it is minimum phase. Wavelet D has the largest phase change; we say it is maximum phase. Finally, wavelets B and C have phase changes between the two extremes; hence, they are mixed phase.

Since all four wavelets have the same amplitude spectrum (Figure 2.2-4) and the same power spectrum, they should have the same autocorrelation. This is verified as shown in Table 2-16, where only one side of the autocorrelation is tabulated, since a real time series has a symmetric autocorrelation (The 1-D Fourier transform).

Note that zero lag of the autocorrelation (Table 2-16) is equal to the total energy (Table 2-15) contained in each wavelet — 17 units. This is true for any wavelet. In fact, Parseval’s theorem states that the area under the power spectrum is equal to the zero-lag value of the autocorrelation function (Section A.1).

The process by which the seismic wavelet is compressed to a zero-lag spike is called spiking deconvolution. In this section, filters that achieve this goal were studied — the inverse and the least-squares inverse filters. Their performance depends not only on filter length, but also on whether the input wavelet is minimum phase.

Table 2-16. Autocorrelation lags of wavelets A, B, C, and D.

Wavelet A
	4	0	−1			Output
	4	0	−1			17
		4	0	−1		0
			4	0	−1	−4
Wavelet B
	2	3	−2			Output
	2	3	−2			17
		2	3	−2		0
			2	3	−2	−4
Wavelet C
	−2	3	2			Output
	−2	3	2			17
		−2	3	2		0
			−2	3	2	−4
Wavelet D
	−1	0	4			Output
	−1	0	4			17
		−1	0	4		0
			−1	0	4	−4

The spiking deconvolution operator is strictly the inverse of the wavelet. If the wavelet were minimum phase, then we would get a stable inverse, which also is minimum phase. The term stable means that the filter coefficients form a convergent series. Specifically, the coefficients decrease in time (and vanish at t = ∞); therefore, the filter has finite energy. This is the case for the wavelet (1, -  1/2) with an inverse $\left(1,\ {\frac {1}{2}},\ {\frac {1}{4}},\ \cdots \right).$ The inverse is a stable spiking deconvolution filter. On the other hand, if the wavelet were maximum phase, then it does not have a stable inverse. This is the case for the wavelet (-  1/2, 1), whose inverse is given by the divergent series (−2, −4, −8, …). Finally, a mixed-phase wavelet does not have a stable inverse. This discussion leads us to assumption 7.

Assumption 7. The seismic wavelet is minimum phase. Therefore, it has a minimum-phase inverse.

Now, a summary of the implications of the underlying assumptions for deconvolution stated in Sections 2.1 and 2.2 is appropriate.

1. Assumptions 1, 2, and 3 allow formulating the convolutional model of the 1-D seismogram by equation (2).
2. Assumption 4 eliminates the unknown noise term in equation (2a) and reduces it to equation (3a).
3. Assumption 5 is the basis for deterministic deconvolution — it allows estimation of the earth’s reflectivity series directly from the 1-D seismogram described by equation (3a).
4. Assumption 6 is the basis for statistical deconvolution — it allows estimates for the autocorrelogram and amplitude spectrum of the normally unknown wavelet in equation (3a) from the known recorded 1-D seismogram.
5. Finally, assumption 7 provides a minimum-phase estimate of the phase spectrum of the seismic wavelet from its amplitude spectrum, which is estimated from the recorded seismogram by way of assumption 6.

Once the amplitude and phase spectra of the seismic wavelet are statistically estimated from the recorded seismogram, its least-squares inverse — spiking deconvolution operator, is computed using optimum Wiener filters. When applied to the wavelet, the filter converts it to a zero-delay spike. When applied to the seismogram, the filter yields the earth’s impulse response (equation 4). In optimum wiener filters, we show that a known wavelet can be converted into a delayed spike even if it is not minimum phase.

2.3 Optimum wiener filters

Return to the desired output — the zero-delay spike (1, 0, 0), that was considered when studying inverse and least-squares filters (Inverse filtering). Rewrite equation (16), which we solved to obtain the least-squares inverse filter, as follows:

$2{\begin{pmatrix}5/4&-1/2\\-1/2&5/4\end{pmatrix}}{\begin{pmatrix}a\\b\end{pmatrix}}={\begin{pmatrix}2\\0\end{pmatrix}}.$

(26)

Divide both sides by 2 to obtain

${\begin{pmatrix}5/4&-1/2\\-1/2&5/4\end{pmatrix}}{\begin{pmatrix}a\\b\end{pmatrix}}={\begin{pmatrix}1\\0\end{pmatrix}}.$

(27)

The autocorrelation of the input wavelet (1, -  1/2) is shown in Table 2-17. Note that the autocorrelation lags are the same as the first column of the 2 × 2 matrix on the left side of equation (27).

Now compute the crosscorrelation of the desired output (1, 0, 0) with the input wavelet (1, -  1/2) (Table 2-18). The crosscorrelation lags are the same as the column matrix on the right side of equation (27).

Table 2-17. Autocorrelation lags of input wavelet (1, -  1/2).

	1	-  1/2		Output
	1	-  1/2		5/4
		1	-  1/2	-  1/2

Table 2-18. Crosscorrelation lags of desired output (1, 0, 0) with input wavelet (1, -  1/2).

	1	0	0	Output
	1	-  1/2		1
		1	-  1/2	0

In general, the elements of the matrix on the left side of equation (27) are the lags of the autocorrelation of the input wavelet, while the elements of the column matrix on the right side are the lags of the crosscorrelation of the desired output with the input wavelet.

Now perform similar operations for wavelet (-  1/2, 1). By rewriting the matrix equation (19), we obtain

$2{\begin{pmatrix}5/4&-1/2\\-1/2&5/4\end{pmatrix}}{\begin{pmatrix}a\\b\end{pmatrix}}={\begin{pmatrix}-1\\0\end{pmatrix}}.$

(28)

Divide both sides by 2 to obtain

${\begin{pmatrix}5/4&-1/2\\-1/2&5/4\end{pmatrix}}{\begin{pmatrix}a\\b\end{pmatrix}}={\begin{pmatrix}-1/2\\0\end{pmatrix}}.$

(29)

The autocorrelation of wavelet (-  1/2, 1) is given in Table 2-19. The elements of the matrix on the left side of equation (29) are the autocorrelation lags of the input wavelet. Note that autocorrelation of wavelet (-  1/2, 1) is identical to that of wavelet (1, -  1/2) (Table 2-17). As discussed in Inverse filtering, an important property of a group of wavelets with the same amplitude spectrum is that they also have the same autocorrelation.

The crosscorrelation of the desired output (1, 0, 0) with input wavelet (-  1/2, 1) is given in Table 2-20. Note that the right side of equation (29) is the same as the crosscorrelation lags.

Table 2-19. Autocorrelation lags of input wavelet (-  1/2, 1).

	-  1/2	1		Output
	-  1/2	1		5/4
		-  1/2	1	-  1/2

Table 2-20. Crosscorrelation lags of desired output (1, 0, 0) with input wavelet (-  1/2, 1).

	1	0	0	Output
	-  1/2	1		-  1/2
		-  1/2	1	0

Matrix equations (27) and (29) were used to derive the least-squares inverse filters (Inverse filtering). These filters then were applied to the input wavelets to compress them to zero-lag spike. The matrices on the left in equations (27) and (29) are made up of the autocorrelation lags of the input wavelets. Additionally, the column matrices on the right are made up of lags of the crosscorrelation of the desired output — a zero-lag spike, with the input wavelets. These observations were generalized by Wiener to derive filters that convert the input to any desired output .

The general form of the matrix equation such as equation (29) for a filter of length n is (Section B.5):

${\begin{pmatrix}r_{0}&r_{1}&r_{2}&\cdots &r_{n-1}\\r_{1}&r_{0}&r_{1}&\cdots &r_{n-2}\\r_{2}&r_{1}&r_{0}&\cdots &r_{n-3}\\\vdots &\vdots &\vdots &\ddots &\vdots \\r_{n-1}&r_{n-2}&r_{n-3}&\cdots &r_{0}\end{pmatrix}}{\begin{pmatrix}a_{0}\\a_{1}\\a_{2}\\\vdots \\a_{n-1}\\\end{pmatrix}}={\begin{pmatrix}g_{0}\\g_{1}\\g_{2}\\\vdots \\g_{n-1}\end{pmatrix}}$

(30)

Here r_i, a_i, and g_i, i = 0, 1, 2, …, n − 1 are the autocorrelation lags of the input wavelet, the Wiener filter coefficients, and the crosscorrelation lags of the desired output with the input wavelet, respectively.

The optimum Wiener filter (a₀, a₁, a₂, …, a_n−1) is optimum in that the least-squares error between the actual and desired outputs is minimum. When the desired output is the zero-lag spike (1, 0, 0, …, 0), then the Wiener filter is identical to the least-squares inverse filter. In other words, the least-squares inverse filter really is a special case of the Wiener filter.

The Wiener filter applies to a large class of problems in which any desired output can be considered, not just the zero-lag spike. Five choices for the desired output are:

Type 1: Zero-lag spike,

Type 2: Spike at arbitrary lag,

Type 3: Time-advanced form of input series,

Type 4: Zero-phase wavelet,

Type 5: Any desired arbitrary shape.

These desired output forms will be discussed in the following sections.

The general form of the normal equations (30) was arrived at through numerical examples for the special case where the desired output was a zero-lag spike. Section B.5 provides a concise mathematical treatment of the optimum Wiener filters. Figure 2.3-1 outlines the design and application of a Wiener filter.

Determination of the Wiener filter coefficients requires solution of the so-called normal equations (30). From equation (30), note that the autocorrelation matrix is symmetric. This special matrix, called the Toeplitz matrix, can be solved by Levinson recursion, a computationally efficient scheme (Section B.6). To do this, compute a two-point filter, derive from it a three-point filter, and so on, until the n-point filter is derived . In practice, filtering algorithms based on the optimum Wiener filter theory are known as Wiener-Levinson algorithms.

Spiking deconvolution

The process with type 1 desired output (zero-lag spike) is called spiking deconvolution. Crosscorrelation of the desired spike (1, 0, 0, …, 0) with input wavelet (x₀, x₁, x₂, …, x_n−1) yields the series (x₀, 0, 0, …, 0).

Figure 2.3-1  A flowchart for Wiener filter design and application.

The generalized form of the normal equations (30) takes the special form:

${\begin{pmatrix}r_{0}&r_{1}&r_{2}&\cdots &r_{n-1}\\r_{1}&r_{0}&r_{1}&\cdots &r_{n-2}\\r_{2}&r_{1}&r_{0}&\cdots &r_{n-3}\\\vdots &\vdots &\vdots &\ddots &\vdots \\r_{n-1}&r_{n-2}&r_{n-3}&\cdots &r_{0}\end{pmatrix}}{\begin{pmatrix}a_{0}\\a_{1}\\a_{2}\\\vdots \\a_{n-1}\\\end{pmatrix}}={\begin{pmatrix}1\\0\\0\\\vdots \\0\end{pmatrix}}$

(31)

Equation (31) was scaled by (1/x₀). The least-squares inverse filter, which was discussed in Inverse filtering, has the same form as the matrix equation (31). Therefore, spiking deconvolution is mathematically identical to least-squares inverse filtering. A distinction, however, is made in practice between the two types of filtering. The autocorrelation matrix on the left side of equation (31) is computed from the input seismogram (assumption 6) in the case of spiking deconvolution (statistical deconvolution), whereas it is computed directly from the known source wavelet in the case of least-squares inverse filtering (deterministic deconvolution).

Figure 2.3-2 is a summary of spiking deconvolution based on the Wiener-Levinson algorithm. Frame (a) is the input mixed-phase wavelet. Its amplitude spectrum shown in frame (b) indicates that the wavelet has most of its energy confined to a 10- to 50-Hz range. The autocorrelation function shown in frame (d) is used in equation (31) to compute the spiking deconvolution operator shown in frame (e). The amplitude spectrum of the operator shown in frame (f) is approximately the inverse of the amplitude spectrum of the input wavelet shown in frame (b). (The approximation improves as operator length increases.) This should be expected, since the goal of spiking deconvolution is to flatten the output spectrum. Application of this operator to the input wavelet gives the result shown in frame (k).

Ideally, we would like to get a zero-lag spike, as shown in frame (n). What went wrong? Assumption 7 was violated by the mixed-phase input wavelet shown in frame (a). Frame (h) shows the inverse of the deconvolution operator. This is the minimum-phase equivalent of the input mixed-phase wavelet in frame (a). Both wavelets have the same amplitude spectrum shown in frames (b) and (i), but their phase spectra are significantly different as shown in frames (c) and (j). Since spiking deconvolution is equivalent to least-squares inverse filtering, the minimum-phase equivalent is merely the inverse of the deconvolution operator. Therefore, the amplitude spectrum of the operator is the inverse of the amplitude spectrum of the minimum-phase equivalent as shown in frames (f) and (i), and the phase spectrum of the operator is the negative of the phase spectrum of the minimum-phase wavelet as shown in frames (g) and (j). One way to extract the seismic wavelet, provided it is minimum phase, is to compute the spiking deconvolution operator and find its inverse.

In conclusion, if the input wavelet is not minimum phase, then spiking deconvolution cannot convert it to a perfect zero-lag spike as in frame (k). Although the amplitude spectrum is virtually flat as shown in frame (l), the phase spectrum of the output is not minimum phase as shown in frame (m). Finally, note that the spiking deconvolution operator is the inverse of the minimum-phase equivalent of the input wavelet. This wavelet may or may not be minimum phase.

Prewhitening

From the preceding section, we know that the amplitude spectrum of the spiking deconvolution operator is (approximately) the inverse of the amplitude spectrum of the input wavelet. This is sketched in Figure 2.3-3. What if we had zeroes in the amplitude spectrum of the input wavelet? To study this, apply a minimum-phase band-pass filter (Exercise 2-10) with a wide passband (3-108 Hz) to the minimum-phase wavelet of Figure 2.3-2, as shown in frame (h). Deconvolution of the filtered wavelet does not produce a perfect spike; instead, a spike accompanied by a high-frequency pre-and post-cursor results (Figure 2.3-4). This poor result occurs because the deconvolution operator tries to boost the absent frequencies, as seen from the amplitude spectrum of the output. Can this problem occur in a recorded seismogram? Situations in which the input amplitude spectrum has zeroes rarely occur. There is always noise in the seismogram and it is additive in both the time and frequency domains. Moreover, numerical noise, which also is additive in the frequency domain, is generated during processing. However, to ensure numerical stability, an artificial level of white noise is added to the amplitude spectrum of the input seismogram before deconvolution. This is called prewhitening and is referred to in Figure 2.3-3.

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

If the percent prewhitening is given by a scalar, 0 ≤ ε < 1, then the normal equations (31) are modified as follows:

${\begin{pmatrix}\beta r_{0}&r_{1}&r_{2}&\cdots &r_{n-1}\\r_{1}&\beta r_{0}&r_{1}&\cdots &r_{n-2}\\r_{2}&r_{1}&\beta r_{0}&\cdots &r_{n-3}\\\vdots &\vdots &\vdots &\ddots &\vdots \\r_{n-1}&r_{n-2}&r_{n-3}&\cdots &\beta r_{0}\end{pmatrix}}{\begin{pmatrix}a_{0}\\a_{1}\\a_{2}\\\vdots \\a_{n-1}\\\end{pmatrix}}={\begin{pmatrix}1\\0\\0\\\vdots \\0\end{pmatrix}},$

(32)

where β = 1 + ε. Adding a constant εr₀ to the zero lag of the autocorrelation function is the same as adding white noise to the spectrum, with its total energy equal to that constant. The effect of the prewhitening level on performance of deconvolution is discussed in Predictive deconvolution in practice.

Figure 2.3-3  Prewhitening amounts to adding a bias to the amplitude spectrum of the seismogram to be deconvolved. This prevents dividing by zero since the amplitude spectrum of the inverse filter (middle) is the inverse of that of the seismogram (left). Convolution of the filter with the seismogram is equivalent to multiplying their respective amplitude spectra — this yields nearly a white spectrum (right).

Wavelet processing by shaping filters

Spiking deconvolution had trouble compressing wavelet (-  1/2, 1) to a zero-lag spike (1, 0, 0) (Table 2-14). In terms of energy distribution, this input wavelet is more similar to a delayed spike, such as (0, 1, 0), than it is to a zero-lag spike, (1, 0, 0). Therefore, a filter that converts wavelet (-  1/2, 1) to a delayed spike would yield less error than the filter that shapes it to a zero-lag spike (Table 2-14).

Recast the filter design and application outlined in Table 2-13 in terms of optimum Wiener filters by following the flowchart in Figure 2.3-1. First, compute the crosscorrelation (Table 2-21). From Table 2-19, we know the autocorrelation of the input wavelet. By substituting the results from Tables 2-19 and 2-21 into the matrix equation (30), we get

${\begin{pmatrix}5/4&-1/2\\-1/2&5/4\\\end{pmatrix}}{\begin{pmatrix}a\\b\\\end{pmatrix}}={\begin{pmatrix}1\\-1/2\end{pmatrix}}.$

(33)

By solving for the filter coefficients, we obtain $\left(a,\ b\right):\left({\frac {16}{21}},\ -{\frac {2}{21}}\right).$ This filter is applied to the input wavelet as shown in Table 2-22. As we would expect, the output is the same as that of the least-squares filter (Table 2-13). Note that, from Table 2-14, the energy of the least-squares error between the actual and desired outputs was 0.190 and 0.762 for a delayed spike and a zero-lag spike desired output, respectively. This shows that there is less error when converting wavelet (-  1/2, 1) to the delayed spike (0, 1, 0) than to zero-lag spike (1, 0, 0).

In general, for any given input wavelet, a series of desired outputs can be defined as delayed spikes. The least-squares errors then can be plotted as a function of delay. The delay (lag) that corresponds to the least error is chosen to define the desired delayed spike output. The actual output from the Wiener filter using this optimum delayed spike should be the most compact possible result.

Table 2-21. Crosscorrelation lags of desired output (0, 1, 0) with input wavelet (-  1/2, 1).

	0	1	0	Output
	-  1/2	1		1
		-  1/2	1	-  1/2

Table 2-22. Convolution of input wavelet $(-{\frac {1}{2}},\ 1)$ with filter coefficients $\left({\frac {16}{21}},\ -{\frac {2}{21}}\right).$

		-  1/2	1		Output
	-  2/21	-  16/21			−0.38
		-  2/21	-  16/21		0.81
			-  2/21	-  16/21	−0.09

The process that has a type 5 desired output (any desired arbitrary shape) is called wavelet shaping. The filter that does this is called a Wiener shaping filter. In fact, type 2 (delayed spike) and type 4 (zero-phase wavelet) desired outputs are special cases of the more general wavelet shaping.

Figure 2.3-5 shows a series of wavelet shapings that use delayed spikes as desired outputs. The input is a mixed-phase wavelet. Filter length was held constant in all eight cases. Note that the zero-delay spike case (spiking deconvolution) does not always yield the best result (Figure 2.3-5a). A delay in the neighborhood of 60 ms (Figure 2.3-5e) seems to yield an output that is closest to being a perfect spike. Typically, the process is not very sensitive to the amount of delay once it is close to the optimum delay. If the input wavelet were minimum-phase, then the optimum delay of the desired output spike generally is zero. On the other hand, if the input wavelet were mixed-phase, as illustrated in Figure 2.3-5, then the optimum delay is nonzero. Finally, if the input wavelet were maximum-phase, then the optimum delay is the length of that wavelet .

Can we not delay the desired spike output (Figure 2.3-5) and obtain a better result than we obtained from spiking deconvolution? This goal is achieved by applying a constant-time shift (60 ms in Figure 2.3-5) to a delayed spike result. Better yet, the same result can be obtained by shifting the shaping filter operator as much as the delay in the spike and applying it to the input wavelet. Such a filter operator is two-sided (non-causal), since it has coefficients for negative and positive time values. The one-sided filter defined along the positive time axis has an anticipation component, while the filter defined along the negative time axis has a memory component . The two-sided filter has an anticipation component and a memory component. Figure 2.3-6 shows a series of shaping filterings with two-sided Wiener filters for various spike delay values.

• 

  Figure 2.3-4  (a) Minimum-phase wavelet, (b) after band-pass filtering, (c) followed by deconvolution. The amplitude spectrum of the band-pass filtered wavelet is zero above 108 Hz (middle row); therefore, the inverse filter derived from it yields unstable results (bottom row). The time delays on the wavelets in the left frames of the middle and bottom rows are for display purposes only.

• 

  Figure 2.3-5  Shaping filtering. (0) Input wavelet, (1) desired output, (2) shaping filter operator, (3) actual output. Here, the purpose is to convert the mixed-phased wavelet (0) to a series of delayed spikes as shown in (a) through (h) by using a one-sided operator (anticipation component only). The best result is with a 60-ms delay (e).

• 

  Figure 2.3-6  Shaping filtering. (0) Input wavelet, (1) desired output, (2) shaping filter operator, (3) actual output. Here, the purpose is to convert the mixed-phase wavelet (0) to a series of delayed spikes as shown in (a) through (h) using a two-sided operator (with memory and anticipation components). The best result is obtained with a zero-delay spike using a two-sided filter (a).

Figure 2.3-7 shows examples of wavelet shaping. The input wavelet represented by trace (b) is the same mixed-phase wavelet as in Figure 2.3-6 (top left frame). This wavelet is shaped into zero-phase wavelets with three different bandwidths represented by traces (c), (d) and (e). This process commonly is referred to as dephasing. Figure 2.3-7 shows another wavelet shaping in which the input wavelet is converted to its minimum-phase equivalent represented by trace (f). This conversion is often applied to recorded air-gun signatures.

Figure 2.3-8 shows examples of a recorded air-gun signature that was shaped into its minimum-phase equivalent and into a spike. When the input is the recorded signature, then the wavelet shapings in Figure 2.3-8 are called signature processing.

Wavelet shaping requires knowledge of the input wavelet to compute the crosscorrelation column on the right side of equation (30). If it is unknown, which is the case in reality, then the minimum-phase equivalent of the input wavelet can be estimated statistically from the data. This minimum-phase estimate then is shaped to a zero-phase wavelet.

Wavelet processing is a term that is used with flexibility. The most common meaning refers to estimating (somehow) the basic wavelet embedded in the seismogram, designing a shaping filter to convert the estimated wavelet to a desired form, usually a broad-band zero-phase wavelet (Figure 2.3-8), and finally, applying the shaping filter to the seismogram. Another type of wavelet processing involves wavelet shaping in which the desired output is the zero-phase wavelet with the same amplitude spectrum as that of the input wavelet (Figure 2.3-9). Note that this type of wavelet processing does not try to flatten the spectrum, but only tries to correct for the phase of the input wavelet, which sometimes is assumed to be minimum-phase.

Predictive deconvolution

The type 3 desired output, a time-advanced form of the input series, suggests a prediction process. Given the input x(t), we want to predict its value at some future time (t + α), where α is prediction lag. Wiener showed that the filter used to estimate x(t + α) can be computed by using a special form of the matrix equation (30) . Since the desired output x(t + α) is the time-advanced version of the input x(t), we need to specialize the right side of equation (30) for the prediction problem.

Figure 2.3-7  Shaping filtering with various desired outputs. (a) Impulse response, (b) input seismogram. Here, (c), (d) and (e) show three possible desired outputs that are band-limited zero-phase wavelets, while (f) shows a desired output that is the minimum-phase equivalent of the input wavelet (b). Finally, (g) and (h) are desired outputs that are band-pass filtered versions of (f).

Consider a five-point input time series x(t): (x₀, x₁, x₂, x₃, x₄), and set α = 2. The autocorrelation of the input series is computed in Table 2-23, and the crosscorrelation between the desired output x(t + 2) and the input x(t) is computed in Table 2-24. Compare the results in Tables 2-23 and 2-24, and note that g_i = r_i+α for α = 2 and i = 0, 1, 2, 3, 4.

Equation (30), for this special case, is rewritten as follows:

${\begin{pmatrix}r_{0}&r_{1}&r_{2}&r_{3}&r_{4}\\r_{1}&r_{0}&r_{1}&r_{2}&r_{3}\\r_{2}&r_{1}&r_{0}&r_{1}&r_{2}\\r_{3}&r_{2}&r_{1}&r_{0}&r_{1}\\r_{4}&r_{3}&r_{2}&r_{1}&r_{0}\\\end{pmatrix}}{\begin{pmatrix}a_{0}\\a_{1}\\a_{2}\\a_{3}\\a_{4}\\\end{pmatrix}}={\begin{pmatrix}r_{2}\\r_{3}\\r_{4}\\r_{5}\\r_{6}\\\end{pmatrix}}.$

(34)

The prediction filter coefficients a(t) : (a₀, a₁, a₂, a₃, a₄) can be computed from equation (34) and applied to the input series x(t) : (x₀, x₁, x₂, x₃, x₄) to compute the actual output y(t) : (y₀, y₁, y₂, y₃, y₄) (Table 2-25). We want to predict the time-advanced form of the input; hence, the actual output is an estimate of the series x(t + α) : (x₂, x₃, x₄), where α = 2. The prediction error series e(t) = x(t + α) − y(t) : (e₂, e₃, e₄, e₅, e₆) is given in Table 2-26.

The results in Table 2-26 suggest that the error series can be obtained more directly by convolving the input series x(t) : (x₀, x₁, x₂, x₃, x₄) with a filter with coefficients (1, 0, −a₀, −a₁, −a₂, −a₃, −a₄) (Table 2-27). The results for (e₂, e₃, e₄, e₅, e₆) are identical (Tables 2-26 and 2-27). Since the series (a₀, a₁, a₂, a₃, a₄) is called the prediction filter, it is natural to call the series (1, 0, −a₀, −a₁, −a₂, −a₃, −a₄) the prediction error filter. When applied to the input series, this filter yields the error series in the prediction process (Table 2-27).

Figure 2.3-8  Signature processing: (a) Recorded signature, (b) desired output, (c) shaping operator, (d) shaped signature. The desired output is a zero-delay spike (top) and the minimum-phase equivalent of the recorded signature (bottom).

Table 2-23. Autocorrelation lags of input series x(t) : (x₀, x₁, x₂, x₃, x₄).

$r_{0}=x_{0}^{2}+x_{1}^{2}+x_{2}^{2}+x_{3}^{2}+x_{4}^{2}$

r1 = x0x1 + x1x2 + x2x3 + x3x4

r2 = x0x2 + x1x3 + x2x4

r3 = x0x3 + x1x4

r4 = x0x4

r5 = 0

r6 = 0

Table 2-24. Crosscorrelation of desired output x(t + α): (x₂, x₃, x₄), α = 2, with input x(t) : (x₀, x₁, x₂, x₃, x₄)

g0 = x0x2 + x1x3 + x2x4

g1 = x0x3 + x1x4

g2 = x0x4

g3 = 0

g4 = 0

Table 2-25. Convolution of prediction filter a(t): (a₀, a₁, a₂, a₃, a₄) with input series x(t): (x₀, x₁, x₂, x₃, x₄) to compute actual output y(t) : (y₀, y₁, y₂, y₃, y₄).

y0 = a0x0

y1 = a1x0 + a0x1

y2 = a2x0 + a1x1 + a0x2

y3 = a3x0 + a2x1 + a1x2 + a0x3

y4 = a4x0 + a3x1 + a2x2 + a1x3 + a0x4

Table 2-26. The error series e(t) = x(t + α) − y(t) : (e₂, e₃, e₄, e₅, e₆), α = 2. For y(t), see Table 2-25.

e2 = x2 − a0x0

e3 = x3 − a1x0 − a0x1

e4 = x4 − a2x0 − a1x1 − a0x2

e5 = 0 − a3x0 − a2x1 − a1x2 − a0x3

e6 = 0 − a4x0 − a3x1 − a2x2 − a1x3 − a0x4

Table 2-27. Convolution of prediction error filter coefficients (1, 0, −a₀, −a₁, −a₂, −a₃, −a₄) with input series x(t) : (x₀, x₁, x₂, x₃, x₄).

e0 = x0

e1 = x1

e2 = x2 − a0x0

e3 = x3 − a1x0 − a0x1

e4 = x4 − a2x0 − a1x1 − a0x2

e5 = 0 − a3x0 − a2x1 − a1x2 − a0x3

e6 = 0 − a4x0 − a3x1 − a2x2 − a1x3 − a0x4

Figure 2.3-9  Wavelet processing. An autocorrelogram (a), estimated from the seismic trace, is used after smoothing (b) to compute the spiking deconvolution operator (d). Here (c) is just a one-sided version of (b). The inverse of the operator (d) is the minimum-phase wavelet (e), which is sometimes assumed to be the basic wavelet contained in the original seismic trace. It is easy to compute its zero-phase equivalent (f) and design a shaping filter (g) that converts the minimum-phase wavelet (e) to the zero-phase wavelet (f). The actual output is (h), which should be compared with (f). The zero-phase equivalent (f) has the same amplitude spectrum as the minimum-phase wavelet (e).

Why place so much emphasis on the error series? Consider the prediction process as it relates to a seismic trace. From the past values of a time series up to time t, a future value can be predicted at time t + α, where α is the prediction lag. A seismic trace often has a predictable component (multiples) with a periodic rate of occurrence. According to assumption 6, anything else, such as primary reflections, is unpredictable.

Figure 2.3-10  A flowchart for predictive deconvolution using a prediction filter.

Some may claim that reflections are predictable as well; this may be the case if deposition is cyclic. However, this type of deposition is not often encountered. While the prediction filter yields the predictable component (the multiples) of a seismic trace, the remaining unpredictable part, the error series, is essentially the reflection series.

Equation (34) can be generalized for the case of an n-long prediction filter and an α-long prediction lag.

${\begin{pmatrix}r_{0}&r_{1}&r_{2}&\cdots &r_{n-1}\\r_{1}&r_{0}&r_{1}&\cdots &r_{n-2}\\r_{2}&r_{1}&r_{0}&\cdots &r_{n-3}\\\vdots &\vdots &\vdots &\ddots &\vdots \\r_{n-1}&r_{n-2}&r_{n-3}&\cdots &r_{0}\\\end{pmatrix}}{\begin{pmatrix}a_{0}\\a_{1}\\a_{2}\\\vdots \\a_{n-1}\\\end{pmatrix}}={\begin{pmatrix}r_{\alpha }\\r_{\alpha +1}\\r_{\alpha +2}\\\vdots \\r_{\alpha +n-1}\\\end{pmatrix}}$

(35)

Note that design of the prediction filters requires only autocorrelation of the input series.

There are two approaches to predictive deconvolution:

1. The prediction filter (a₀, a₁, a₂, …, a_n−1) may be designed using equation (35) and applied on input series as described in Figure 2.3-10.
2. Alternatively, the prediction error filter (1, 0, 0, …, 0, −a₀, −a₁, −a₂, …, −a_n−1) can be designed and convolved with the input series as described in Figure 2.3-11.

Now consider the special case of unit prediction lag, α = 1. For n = 5, equation (35) takes the following form:

${\begin{pmatrix}r_{0}&r_{1}&r_{2}&r_{3}&r_{4}\\r_{1}&r_{0}&r_{1}&r_{2}&r_{3}\\r_{2}&r_{1}&r_{0}&r_{1}&r_{2}\\r_{3}&r_{2}&r_{1}&r_{0}&r_{1}\\r_{4}&r_{3}&r_{2}&r_{1}&r_{0}\\\end{pmatrix}}{\begin{pmatrix}a_{0}\\a_{1}\\a_{2}\\a_{3}\\a_{4}\\\end{pmatrix}}={\begin{pmatrix}r_{1}\\r_{2}\\r_{3}\\r_{4}\\r_{5}\\\end{pmatrix}}.$

(36)

Figure 2.3-11  A flowchart for predictive deconvolution using a prediction error filter.

By augmenting the right side to the left side, we obtain:

${\begin{pmatrix}-r_{1}&r_{0}&r_{1}&r_{2}&r_{3}&r_{4}\\-r_{2}&r_{1}&r_{0}&r_{1}&r_{2}&r_{3}\\-r_{3}&r_{2}&r_{1}&r_{0}&r_{1}&r_{2}\\-r_{4}&r_{3}&r_{2}&r_{1}&r_{0}&r_{1}\\-r_{5}&r_{4}&r_{3}&r_{2}&r_{1}&r_{0}\\\end{pmatrix}}{\begin{pmatrix}1\\a_{0}\\a_{1}\\a_{2}\\a_{3}\\a_{4}\\\end{pmatrix}}={\begin{pmatrix}0\\0\\0\\0\\0\\0\\\end{pmatrix}}.$

(37)

Add one row and move the negative sign to the column matrix that represents the filter coefficients to get:

${\begin{pmatrix}r_{0}&r_{1}&r_{2}&r_{3}&r_{4}&r_{5}\\r_{1}&r_{0}&r_{1}&r_{2}&r_{3}&r_{4}\\r_{2}&r_{1}&r_{0}&r_{1}&r_{2}&r_{3}\\r_{3}&r_{2}&r_{1}&r_{0}&r_{1}&r_{2}\\r_{4}&r_{3}&r_{2}&r_{1}&r_{0}&r_{1}\\r_{5}&r_{4}&r_{3}&r_{2}&r_{1}&r_{0}\\\end{pmatrix}}{\begin{pmatrix}1\\-a_{0}\\-a_{1}\\-a_{2}\\-a_{3}\\-a_{4}\\\end{pmatrix}}={\begin{pmatrix}L\\0\\0\\0\\0\\0\\\end{pmatrix}}.$

(38)

where L = r₀ − r₁a₀ − r₂a₁ − r₃a₂ − r₄a₃ − r₅a₄. Note that there are six unknowns, (a₀, a₁, a₂, a₄, a₅, L), and six equations. Solution of these equations yields the unit-delay prediction error filter (1, −a₀, −a₁, −a₂, −a₄, −a₅), and the quantity L — the error in the filtering process (Section B.5). We can rewrite equation (38) as follows:

${\begin{pmatrix}r_{0}&r_{1}&r_{2}&r_{3}&r_{4}&r_{5}\\r_{1}&r_{0}&r_{1}&r_{2}&r_{3}&r_{4}\\r_{2}&r_{1}&r_{0}&r_{1}&r_{2}&r_{3}\\r_{3}&r_{2}&r_{1}&r_{0}&r_{1}&r_{2}\\r_{4}&r_{3}&r_{2}&r_{1}&r_{0}&r_{1}\\r_{5}&r_{4}&r_{3}&r_{2}&r_{1}&r_{0}\\\end{pmatrix}}{\begin{pmatrix}b_{0}\\b_{1}\\b_{2}\\b_{3}\\b_{4}\\b_{5}\\\end{pmatrix}}={\begin{pmatrix}L\\0\\0\\0\\0\\0\\\end{pmatrix}}.$

(39)

where b₀ = 1, b_i = −a_i, and i = 1, 2, 3, 4, 5. This equation has a familiar structure. In fact, except for the scale factor L, it has the same form as equation (31), which yields the coefficients for the least-squares zero-delay inverse filter. This inverse filter is therefore the same as the prediction error filter with unit prediction lag, except for a scale factor. Hence, spiking deconvolution actually is a special case of predictive deconvolution with unit prediction lag.

We now know that predictive deconvolution is a general process that encompasses spiking deconvolution. In general, the following statement can be made: Given an input wavelet of length (n + α), the prediction error filter contracts it to an α-long wavelet, where α is the prediction lag . When α = 1, the procedure is called spiking deconvolution.

Figure 2.3-12 interrelates the various filters discussed in this chapter and indicates the kind of process they imply. From Figure 2.3-12, note that Wiener filters can be used to solve a wide range of problems. In particular, predictive deconvolution is an integral part of seismic data processing that is aimed at compressing the seismic wavelet, thereby increasing temporal resolution. In the limit, it can be used to spike the seismic wavelet and obtain an estimate for reflectivity.

Figure 2.3-12  A flowchart for interrelations between various deconvolution filters.

2.4 Predictive deconvolution in practice

It now is appropriate to review the implications of the assumptions stated in Sections 2.1 and 2.2 that underlie the process of deconvolution within the context of predictive deconvolution.

1. Assumptions 1, 2, and 3 are the basis for the convolutional model of the recorded seismogram (The convolutional model). In practice, deconvolution often yields good results in areas where these three assumptions are not strictly valid.
2. Assumption 3 can be relaxed in practice by considering a time-variant deconvolution (The problem of nonstationarity). In this technique, a seismogram is divided into a number of time gates, typically three or more. Deconvolution operators then are designed from each gate and convolved with data within that gate. Alternatively, time-variant spectral whitening can be used to account for nonstationarity (The problem of nonstationarity).
3. Not much can be done about assumption 4. However, noise can be minimized in the recording process. Deconvolution operators can be designed using time gates and frequency bands with low noise levels. Poststack deconvolution can be used in an effort to take advantage of the noise reduction inherent in the stacking process.
4. If the source wavelet were minimum-phase and known (assumption 5), then a perfect result could be obtained from deconvolution in the noise-free case as in trace (c) of Figures 2.4-1 and 2.4-2.
5. If assumption 6 were violated and if the source waveform were not known, then you would have problems as in trace (d) of Figures 2.4-1 and 2.4-2.
6. The quality of the output from spiking deconvolution is degraded further when the source wavelet is not minimum-phase as in Figures 2.4-3 and 2.4-4; that is, when assumption 7 is violated.
7. Finally, in addition to violating assumptions 5 and 7, if there were noise in the data, that is, when assumption 4 is violated, then the result of the deconvolution would be unacceptable as in Figure 2.4-5.

Figures 2.4-1 through 2.4-5 test our confidence in the usefulness of predictive deconvolution. In reality, deconvolution has been applied to billions of seismic traces; most of the time it has yielded satisfactory results. Figures 2.4-1 through 2.4-5 emphasize the critical assumptions that underlie predictive deconvolution. When deconvolution does not work on some data, the most probable reason is that one or more of the above assumptions has been violated. In the remaining part of this section, a series of numerical experiments will be performed to examine the validity of these assumptions. The purpose of these experiments is to gain a basic understanding of deconvolution from a practical point of view.

Operator length

We start with a single, isolated minimum-phase wavelet as in trace (b) of Figure 2.4-6. Assumptions 1 through 5 are satisfied for this wavelet. The ideal result of spiking deconvolution is a zero-lag spike, as indicated by trace (a). In this and the following numerical analyses, we refer to the autocorrelogram and amplitude spectrum (plotted with linear scale) of the output from each deconvolution test to better evaluate the results. In Figure 2.4-6 and the following figures, n, α, and ε refer to operator length of the prediction filter, prediction lag, and percent prewhitening, respectively. The length of the prediction error filter then is n + α.

In Figure 2.4-6, prediction lag is unity and equal to the 2-ms sampling rate, prewhitening is 0%, and operator length varies as indicated in the figure. Short operators yield spikes with small-amplitude and relatively high-frequency tails. The 128-ms-long operator gives an almost perfect spike output. Longer operators whiten the spectrum further, bringing it closer to the spectrum of the impulse response.

The action of spiking deconvolution on the seismogram derived by convolving the minimum-phase wavelet with a sparse-spike series is similar (Figure 2.4-7) to the case of the single isolated wavelet (Figure 2.4-6). Recall that spiking deconvolution basically is inverse filtering where the operator is the least-squares inverse of the seismic wavelet. Therefore, an increasingly better result should be obtained when more and more coefficients are included in the inverse filter.

Now consider the real situation of an unknown source wavelet. Based on assumption 6, autocorrelation of the input seismogram rather than that of the seismic wavelet is used to design the deconvolution operator. The result of using the trace rather than the wavelet autocorrelation is shown in Figure 2.4-8. Deconvolution recovers the gross aspects of the spike series, trace (a). However, note that the deconvolved traces have spurious small-amplitude spikes trailing each of the real spikes. We see that increasing operator length does not indefinitely improve the results; on the contrary, more and more spurious spikes are introduced.

Very short operators produce the same type of noise spikes as in Figures 2.4-7 and 2.4-8. Examine the series of deconvolution tests in Figure 2.4-8 and note that the 94-ms operator does the best job. Compare the autocorrelogram of trace (b) in Figure 2.4-8 with that of trace (b) in Figure 2.4-6. Note that only the first 100-ms portion represents the autocorrelation of the source wavelet. This explains why the 94-ms operator worked best; that is, the autocorrelation lags of trace (b) in Figure 2.4-8 beyond 94 ms do not represent the seismic wavelet.

• 

  Figure 2.4-1  (a) Impulse response, (b) seismogram, (c) spiking deconvolution using known, minimum-phase wavelet, (d) deconvolution assuming an unknown, minimum-phase source wavelet. Impulse response (a) is a sparse-spike series. For an unknown source wavelet (in violation of assumption 4), spiking deconvolution yields a less than perfect result (compare (c) and (d)).

• 

  Figure 2.4-2  (a) Impulse response, (b) seismogram, (c) spiking deconvolution using known, minimum-phase source wavelet, (d) deconvolution assuming an unknown, minimum-phase source wavelet. Impulse response (a) is based on a sonic log (Figure 2.1-1a). For the unknown source wavelet (in violation of assumption 4), spiking deconvolution yields a less than perfect result. (Compare (c) and (d).

• 

  Figure 2.4-3  (a) Impulse response, (b) seismogram, (c) deconvolution using a known, mixed-phase source wavelet, (d) deconvolution assuming an unknown, mixed-phase source wavelet. Impulse response (a) is a sparse-spike series. For a mixed-phase source wavelet (in violation of assumption 5), spiking deconvolution yields a degraded output (d), even when the wavelet is known (c).

• 

  Figure 2.4-4  (a) Impulse response, (b) seismogram, (c) deconvolution using a known, mixed-phase source wavelet, (d) deconvolution assuming an unknown, mixed-phase source. Impulse response (a) is based on the sonic log of Figure 2.1-1a. For the mixed-phase source wavelet (in violation of assumption 5), spiking deconvolution yields a degraded output (d) even when the wavelet is known (c).

• 

  Figure 2.4-5  (a) Impulse response, (b) seismogram with noise, (c) deconvolution assuming an unknown, mixed-phase source wavelet. Impulse response (a) is based on the sonic log of Figure 2.1-1a. In the presence of random noise (in violation of assumption 3), spiking deconvolution can produce a result with nelation to the earth’s reflectivity (compare (a) to (c)).

Consider the seismogram in Figure 2.4-9, where the wavelet is assumed to be unknown. Deconvolution has restored the spikes that correspond to major reflections in the impulse response as in trace (b) with some success. The 64-ms operator is a good choice.

The mixed-phase wavelet in Figure 2.4-10 shows what can happen when assumption 7 is violated. The wavelet in Figure 2.4-6 is the minimum-phase equivalent of the mixed-phase wavelet in Figure 2.4-10. Both wavelets have the same autocorrelograms and amplitude spectra. Hence, the deconvolution operators for both wavelets are identical. Because the minimum-phase assumption was violated, deconvolution does not convert the mixed-phase wavelet to a perfect spike. Instead, the deconvolved output is a complicated high-frequency wavelet. Also note that the dominant peak in the output is negative, while the impulse response has a positive spike. This difference in the sign can happen when a mixed-phase wavelet is deconvolved. Increasing the operator length further whitens the spectrum; however, the 128-ms operator yields a result that cannot be improved further by longer operators.

The seismogram obtained from the mixed-phase wavelet and the sparse-spike series (used in the preceding figures) is shown in Figure 2.4-11. The 94-ms operator gives the best result. This also is the case in Figure 2.4-12, where both assumptions 6 and 7 are violated. The situation with the seismogram in Figure 2.4-13 is not very good. The spikes that correspond to major reflections in the impulse response were restored; however, there are some timing errors and polarity reversals. (Compare these results with those in Figure 2.4-9 for the events between 0.2 and 0.3 s and 0.6 and 0.7 s.) The 64-ms operator gives an output that cannot be improved by longer operators.

What kind of operator length should be used for spiking deconvolution? To select an operator length, ideally we want to use the autocorrelation of the unknown seismic wavelet. Fortunately, the autocorrelation of the input seismogram has the characteristics of the wavelet autocorrelation (assumption 6). Therefore, it seems appropriate that we should use part of the autocorrelation obtained from the input seismogram that most resembles the autocorrelation of the unknown seismic wavelet. That part is the first transient zone in the autocorrelation, as seen by comparing the autocorrelations of trace (b) in Figure 2.4-6 and trace (c) in Figure 2.4-9. The autocorrelations of trace (b) in Figure 2.4-10 and trace (c) in Figure 2.4-13 suggest the same principle.

Prediction lag

So far, we have learned that predictive deconvolution has two uses: (a) spiking deconvolution — the case of unit prediction lag, and (b) predicting the input seismogram at a future time defined by the prediction lag. Case (b) is used to predict and attenuate multiples.

• 

  Figure 2.4-6  Test of operator length for a single, isolated input wavelet, where n = operator length, α = prediction lag, and ε = percent prewhitening. (a) Impulse response, (b) seismogram with minimum-phase source wavelet.

• 

  Figure 2.4-7  Test of operator length where n = operator length, α = prediction lag, and ε = percent prewhitening. (a) Impulse response, (b) seismogram with known, minimum-phase source wavelet.

• 

  Figure 2.4-8  Test of operator length where n = operator length, α = prediction lag, and ε = percent prewhitening. (a) Impulse response, (b) seismogram with known, minimum-phase source wavelet.

• 

  Figure 2.4-9  Test of operator length where n = operator length, α = prediction lag, and ε = percent prewhitening. (a) Reflectivity, (b) impulse response, (c) seismogram with unknown, minimum-phase source wavelet.

• 

  Figure 2.4-10  Test of operator length for a single, isolated input wavelet where n = operator length, α = prediction lag, and ε = percent prewhitening. (a) Impulse response, (b) seismogram with mixed-phase source wavelet.

• 

  Figure 2.4-11  Test of operator length where n = operator length, α = prediction lag, and ε = percent prewhitening. (a) Impulse response, (b) seismogram with known, mixed-phase source wavelet.

• 

  Figure 2.4-12  Test of operator length where n = operator length, α = prediction lag, and ε = percent prewhitening. (a) Impulse response, (b) seismogram with unknown, mixed-phase source wavelet.

• 

  Figure 2.4-13  Test of operator length where n = operator length, α = prediction lag, and ε = percent prewhitening. (a) Reflectivity, (b) impulse response, (c) seismogram with unknown, mixed-phase source wavelet.

Now, the effect of the prediction lag parameter is examined from an interpretive point of view. Consider the single, isolated minimum-phase wavelet in Figure 2.4-14. Here, operator length and percent prewhitening are kept constant, while prediction lag is varied. When prediction lag is equal to the sampling rate, then the result is equivalent to spiking deconvolution. Predictive deconvolution using a prediction lag greater than unity yields a wavelet of finite duration instead of a spike. Given an input wavelet of α + n samples, predictive deconvolution using prediction filter with length n and prediction lag α converts this wavelet into another wavelet that is α samples long. The first α lags of the autocorrelation are preserved, while the next n lags are zeroed out. Additionally, the amplitude spectrum of the output increasingly resembles that of the input wavelet as prediction lag is increased (Figure 2.4-14). At a 94-ms prediction lag, predictive deconvolution does nothing to the input wavelet because almost all the lags of its autocorrelation have been left untouched. This experiment has an important practical implication: Under the ideal, noise-free conditions, resolution on the output from predictive deconvolution can be controlled by adjusting the prediction lag. Unit prediction lag implies the highest resolution, while a larger prediction lag implies less than full resolution. However, in reality, these assessments are dictated by the signal-to-noise ratio. The deconvolved output using a unit prediction lag contains high frequencies; nevertheless, resolution may be degraded if the high-frequency energy is mostly noise, not signal.

In Figure 2.4-14, prediction lags of 8 and 22 ms correspond to the first and second zero crossings on autocorrelation of the input wavelet, respectively. The first zero crossing produces a spike with some width, while the second zero crossing lag produces a wavelet with a positive and negative lobe.

The relationship between prediction lag and whitening also holds for the sparse-spike series in Figure 2.4-15 and when the input wavelet is unknown (Figure 2.4-16).

The effect of prediction lag on the output from predictive deconvolution of a synthetic seismogram, which was obtained from the sonic log (Figure 2.1-1a), is demonstrated in Figures 2.4-17 and 2.4-18. As the prediction lag is increased, the output spectrum becomes increasingly less broadband. Predictive deconvolution of seismograms constructed from the mixed-phase wavelet again demonstrates that output resolution can be controlled by adjusting prediction lag (Figures 2.4-19 through 2.4-23).

If prediction lag is increased, then the output amplitude spectrum becomes increasingly band-limited. The output also can be band-limited by applying a band-pass filter on the spiking deconvolution output. Are these two ways of band-limiting equivalent? Refer to the results from both the minimum- and mixed-phase wavelets in Figures 2.4-14 and 2.4-19, respectively. Note that the output of the 22-ms prediction lag has an amplitude spectrum that is band-limited to approximately 0 to 100 Hz. However, the spectral shape within this bandwidth is not a boxcar, but rather similar to that of the input wavelet. The boxcar shape would be the case if a band-pass filter (0 to 100 Hz) were applied to the output of the spiking deconvolution (2-ms prediction lag). Hence, spiking deconvolution followed by band-pass filtering is not equivalent to predictive deconvolution with a prediction lag greater than unity.

• 

  Figure 2.4-14  Test of prediction lag for a single, isolated input wavelet where n = operator length, α = prediction lag, and ε = percent prewhitening. (a) Impulse response, (b) seismogram with minimum-phase source wavelet.

• 

  Figure 2.4-15  Test of prediction lag where n = operator length, α = prediction lag, and ε = percent prewhitening. (a) Impulse response, (b) seismogram with known, minimum-phase source wavelet.

• 

  Figure 2.4-16  Test of prediction lag where n = operator length, α = prediction lag, and ε = percent prewhitening. (a) Impulse response, (b) seismogram with unknown, minimum-phase source wavelet.

• 

  Figure 2.4-17  Test of prediction lag where n = operator length, α = prediction lag, and ε = percent prewhitening. (a) Reflectivity, (b) impulse response, (c) seismogram with known minimum-phase source wavelet.

• 

  Figure 2.4-18  Test of prediction lag where n = operator length, α = prediction lag, and ε = percent prewhitening. (a) Reflectivity, (b) impulse response, (c) seismogram with unknown, minimum-phase source wavelet.

• 

  Figure 2.4-19  Test of prediction lag for a single, isolated input wavelet where n = operator length, α = prediction lag, and ε = percent prewhitening. (a) Impulse response, (b) seismogram with mixed-phase source wavelet.

• 

  Figure 2.4-20  Test of prediction lag where n = operator length, α = prediction lag, and ε = percent prewhitening. (a) Impulse response, (b) seismogram with known, mixed-phase source wavelet.

• 

  Figure 2.4-21  Test of prediction lag where n = operator length, α = prediction lag, and ε = percent prewhitening. (a) Impulse response, (b) seismogram with unknown, mixed-phase source wavelet.

• 

  Figure 2.4-22  Test of prediction lag where n = operator length, α = prediction lag, and ε = percent prewhitening. (a) Reflectivity, (b) impulse response, (c) seismogram with known, mixed-phase source wavelet.

• 

  Figure 2.4-23  Test of prediction lag where n = operator length, α = prediction lag, and ε = percent prewhitening. (a) Reflectivity, (b) impulse response, (c) seismogram with unknown, mixed-phase source wavelet.

In conclusion, if prediction lag is increased, the output from predictive deconvolution becomes less spiky. This effect can be used to our advantage, since it allows the bandwidth of deconvolved output to be controlled by adjusting the prediction lag. The application of spiking deconvolution to field data is not always desirable, since it boosts high-frequency noise in the data. The most prominent effect of the nonunity prediction lag is suppression of the high-frequency end of the spectrum and preservation of the overall spectral shape of the input data. This effect is seen in Figures 2.4-18 and 2.4-23, which correspond to the minimum-and mixed-phase seismic wavelets. If prediction lag is increased further, then the low-frequency end of the spectrum is affected as well, making the output more band-limited.

Percent prewhitening

The reasons for prewhitening were discussed in Optimum wiener filters. Consider the single, isolated minimum-phase wavelet in Figure 2.4-24. Keep the operator length and prediction lag constant and vary the percent prewhitening. Note that the effect of varying prewhitening is similar to that of varying the prediction lag; that is, the spectrum increasingly becomes less broadband as the percent prewhitening is increased. Compare Figure 2.4-14 with Figure 2.4-24. Note that prewhitening narrows the spectrum without changing much of the flatness character, while larger prediction lag narrows the spectrum and alters its shape, making it look more like the spectrum of the input seismic wavelet. These characteristics also can be inferred from the shapes of the output wavelets. Prewhitening preserves the spiky character of the output, although it adds a low-amplitude, high-frequency tail (Figure 2.4-24). On the other hand, increasing prediction lag produces a wavelet with a duration equal to the prediction lag (Figure 2.4-14).

The effect of prewhitening on the sparse-spike train seismogram with a known and unknown minimum-phase wavelet is shown in Figures 2.4-25 and 2.4-26, respectively. The effect of prewhitening on deconvolution of the synthetic seismogram obtained from the sonic log (Figure 2.1-1a) is shown in Figures 2.4-27 and 2.4-28 for known and unknown minimum-phase wavelets. Prewhitening tests using the mixed-phase wavelet are shown in Figure 2.4-29. Finally, the combined effects of a prediction lag that is greater than unity and prewhitening for the single, isolated wavelet are shown in Figure 2.4-30. These figures demonstrate that prewhitening narrows the output spectrum, making it band-limited. In particular, the tests in Figures 2.4-24 and 2.4-29 using the single, isolated minimum- and mixed-phase wavelets suggest that spiking deconvolution with some prewhitening is somewhat equivalent to spiking deconvolution without prewhitening followed by post-deconvolution broad band-pass filtering. However, this is not exactly true, for prewhitening still leaves some relatively suppressed energy at the high-frequency end of the spectrum. From Figure 2.4-30, we infer that predictive deconvolution with a prediction lag greater than unity and with some prewhitening yields a result somewhat equivalent to a spiking deconvolution followed by band-pass filtering.

• 

  Figure 2.4-24  Test of percent prewhitening for a single, isolated input wavelet where n = operator length, α = prediction lag, and ε = percent prewhitening. (a) Impulse response, (b) seismogram with minimum-phase source wavelet.

• 

  Figure 2.4-25  Test of percent prewhitening where n = operator length, α = prediction lag, and ε = percent prewhitening. (a) Impulse response, (b) seismogram with known, minimum-phase source wavelet.

• 

  Figure 2.4-26  Test of percent prewhitening where n = operator length, α = prediction lag, and ε = percent prewhitening. (a) Impulse response, (b) seismogram with unknown, minimum-phase source wavelet.

• 

  Figure 2.4-27  Test of percent prewhitening where n = operator length, α = prediction lag, and ε = percent prewhitening. (a) Reflectivity, (b) impulse response, (c) seismogram with known, minimum-phase source wavelet.

• 

  Figure 2.4-28  Test of percent prewhitening where n = operator length, α = prediction lag, and ε = percent prewhitening. (a) Reflectivity, (b) impulse response, (c) seismogram with unknown, minimum-phase source wavelet.

• 

  Figure 2.4-29  Test of percent prewhitening for a single, isolated input wavelet where n = operator length, α = prediction lag, and ε = percent prewhitening. (a) Impulse response, (b) seismogram with mixed-phase source wavelet.

• 

  Figure 2.4-30  Test of percent prewhitening for a single, isolated input wavelet where n = operator length, α = prediction lag, and ε = percent prewhitening. (a) Impulse response, (b) seismogram with minimum-phase source wavelet.

In conclusion, we can say that prewhitening yields a band-limited output. However, the effect is less controllable when compared to varying the prediction lag. By varying prediction lag, we have some idea of the output bandwidth, since it is related to prediction lag. The smaller the prediction lag, the broader the output bandwidth. Prewhitening is used only to ensure that numerical instability in solving for the deconvolution operator (equation 32) is avoided. In practice, typically 0.1 to 1% prewhitening is standard.

Effect of random noise on deconvolution

We assume that the noise component in the recorded seismogram is zero (assumption 4). The autocorrelation of ideal random noise is zero at all lags except the zero lag (Figure 2.1-5). Therefore, the effect of random noise on deconvolution operators should be somewhat similar to the effect of prewhitening. Both effects modify the diagonal of the autocorrelation matrix, making it more dominant [equation (32)]. However, the noise component also slightly modifies the nonzero lags of the autocorrelation. Compare the autocorrelograms of traces (b) in Figures 2.4-24 and 2.4-31. In Figure 2.4-24, an isolated minimum-phase wavelet was considered, while in Figure 2.4-31, random noise was added to the same wavelet. The output wavelet shape from spiking deconvolution of the noisy wavelet using a 128-ms operator is similar to the output from spiking deconvolution of the wavelet without noise, using the same operator length but with, say, 20 percent prewhitening. This result has practical importance — prewhitening is equivalent to adding perfect random noise to the system. Since a recorded seismogram always contains some amount of random noise, only a minute amount, say 0.1 percent, of the white noise needs to be added to the seismogram for numerical stability.

The effect of random noise on the performance of deconvolution is examined further in Figures 2.4-32 and 2.4-33. These results should be compared with their noiseless counterparts in Figures 2.4-9 and 2.4-13, respectively. Observe that the noise component has a harmful effect on deconvolution. For example, when comparing Figures 2.4-9 and 2.4-32, note that the deconvolution result from the noisy seismogram has spurious spikes (for instance, between 0.5 and 0.6 s), which could be interpreted as genuine reflections. Noisy field data, which yield better stack when not treated by deconvolution, have been noted. Only by testing can we determine whether deconvolution performs satisfactorily on data with a severe noise problem.

• 

  Figure 2.4-31  Effect of random noise on deconvolution performance. Input seismogram (b) associated with reflectivity (a) contains a single, isolated wavelet (at around 0.2 s) buried in random noise. Here, n = operator length, α = prediction lag, and ε = percent prewhitening.

• 

  Figure 2.4-32  Effect of random noise on deconvolution performance. (a) Reflectivity, (b) impulse response. Input seismogram (c) with unknown, minimum-phase source wavelet is noise contaminated. Here, n = operator length, α = prediction lag, and ε = percent prewhitening.

• 

  Figure 2.4-33  Effect of random noise on deconvolution performance. (a) Reflectivity, (b) impulse response. Input seismogram (c) with unknown, mixed-phase source wavelet is noise contaminated. Here, n = operator length, α = prediction lag, and ε = percent prewhitening.

Multiple attenuation

We have learned that a prediction filter predicts periodic events, like multiples, in the seismogram. The prediction error filter yields the unpredictable component of the seismogram — the reflectivity series. For example, consider the simple case of water-bottom multiples. If the reflection coefficient of the water bottom is c_w and if water depth is equivalent to a two-way time t_w, then the time series is

$\left(1,\ 0,\ \ldots ,\ 0,\ -{{c}_{w}},\ 0,\ \ldots ,\ 0,\ c_{w}^{2},\ 0,\ \ldots ,\ 0,\ -c_{w}^{3},\ 0,\ \ldots \right)$

as represented by trace (b) in Figure 2.4-34. The separation between the spikes is t_w in trace (b). Note that the periodicity in the time series (trace (b) or (c)) manifests itself in the amplitude spectrum as periodic peaks (or notches). The greater the spike separation in time, the closer the peaks (or notches) in the amplitude spectrum.

The noise-free convolutional model for the seismogram that contains the water-bottom multiples can be written as

$x(t)=w(t)\ast m(t)\ast e(t),$

(40)

where m(t) represents the water-layer reverberation spike series as in trace (b) of Figure 2.4-34, and e(t) now represents the earth’s impulse response excluding multiples associated with the water bottom. Predictive deconvolution can suppress the periodic component m(t) in the seismogram as demonstrated by trace (d) in Figure 2.4-3.

Note the two distinct goals for predictive deconvolution: (a) spiking the seismic wavelet w(t), and (b) predicting and attenuating multiples m(t). The first goal is achieved using an operator with unit prediction lag, while the second is achieved using an operator with a prediction lag greater than unity.

The autocorrelation of the input trace can be used to determine the appropriate prediction lag for multiple suppression. Periodicity associated with multiples is evident in the autocorrelogram of trace (c) in Figure 2.4-34, as an isolated series of energy lobes in the neighborhood of 0.2 and 0.4 s. Prediction lag should be chosen to bypass the first part of the autocorrelogram that represents the seismic wavelet. Operator length should be chosen to include the first isolated energy packet in the autocorrelogram. After applying predictive deconvolution, we are left with only the water-bottom primary reflection. Isolated bursts in the autocorrelogram have been suppressed, while periodic peaks in the amplitude spectrum have been eliminated as shown in trace (d). If desired, the basic wavelet can be compressed into a spike as shown in trace (e) by applying spiking deconvolution to the output of predictive deconvolution as shown in trace (d). The sequence can be interchanged by first applying spiking deconvolution as shown in trace (f) followed by predictive deconvolution as shown in trace (g).

• 

  Figure 2.4-34  (a) Reflectivity, (b) impulse response, (c) seismogram. Two-step deconvolution aimed at attenuating the multiples, then spiking the remaining primary wavelet (d) to (e). The process can be performed in the reverse order (f) to (g). Here, n = operator length and α = prediction lag.

• 

  Figure 2.4-35  Predictive deconvolution for multiple attenuation. (a) Reflectivity, (b) impulse response, (c) seismogram. Two-step deconvolution: Predictive deconvolution (d) followed by spiking deconvolution (e). Traces (f), (g), and (h) result when single-step deconvolution is applied to the input trace (c), using the operator lengths n and prediction lags α, as indicated.

By using a sufficiently long spiking deconvolution operator, two goals are achieved in one step as seen in trace (h). However, this approach can be dangerous if primary reflections are unintentionally suppressed. This is the case in Figure 2.4-35. Here, the water-bottom reflection is followed by a deeper event at about 0.28 s as seen in trace (a). The impulse response contains water-bottom multiples and the peg-leg multiples that are associated with the deeper reflector as seen in trace (b). The amplitude spectrum has peaks that come in pairs, indicating the presence of two different periodic components in the seismogram. Careful choice of predictive deconvolution parameters yields an output with only the wavelets associated with the water bottom and the deeper reflector as seen in trace (d). This is followed by a spiking deconvolution that yields two spikes representative of the water bottom and the deep primary as seen in trace (e). Spiking deconvolution alone produces the reflection coefficient series and the spikes that represent the multiples as seen in trace (f). If a longer spiking deconvolution operator is used, then the primary reflection easily can be eliminated as in trace (g). If a predictive deconvolution operator is used with an improper parameter choice, then again, the primary reflection can be eliminated easily as in trace (h).

How can we ensure that no primaries are destroyed by deconvolution? Examine the autocorrelogram of trace (c) in Figure 2.4-35. The first 50-ms portion represents the seismic wavelet. This is followed by a burst between 50 to 170 ms that represents the correlation of the water bottom and primary. The isolated burst between 170 to 340 ms represents the actual multiple series (both the peg-legs and water-bottom multiples). The prediction lag must be chosen to bypass the first part of the autocorrelogram, which represents the seismic wavelet and possible correlation between the primaries. The operator length must be chosen to include the first isolated burst, in this case between 170 to 340 ms.

It is only with vertical incidence and zero-offset recording that periodicity of the multiples is preserved. Therefore, predictive deconvolution aimed at multiple suppression may not be entirely effective when applied to nonzero-offset data, such as common-shot or common-midpoint data. Figure 2.4-36a shows a common-shot gather with its autocorrelogram and average amplitude spectrum. The field record is prepared for deconvolution by first applying t²–scaling (Figure 2.4-36b) and muting the first arrivals associated with largely guided waves (Figure 2.4-36c). The autocorrelogram in Figure 2.4-36c indicates the presence of multiples. Note on the shot record the water-bottom reflection is at 0.4 s at near offset. Additionally, there are two strong primary reflectors at 0.6 and 1.4 s at near offset; these primaries give rise to a first-order and peg-leg multiples. Despite the flatness of the spectrum and the attenuation of nonzero lags of the autocorrelogram after deconvolution, the first-order multiples associated with the water-bottom reflection and the peg-leg multiples associated with the primary reflection at 0.6 s still persist in the record (Figures 2.4-36d). This occurs because these events have large moveout which causes significant departure from periodicity at nonzero offsets. On the other hand, note that the peg-leg multiples associated with the primary reflection at 1.4 s have been attenuated significantly by deconvolution. This event has a very small moveout; thus, its perodicity is much preserved.

Predictive deconvolution sometimes is applied to CMP stacked data in an effort to suppress multiples. The performance of such an approach can be unsatisfactory, because the amplitude relationships between multiples often are grossly altered by the stacking process, primarily because of velocity differences between primaries and multiples. Also, geometric spreading compensation by using primary velocity function adversely affects the amplitudes of multiples on nonzero-offset data.

There is one domain in which the periodicity and amplitudes of multiples are preserved — the slant-stack domain. In The slant-stack transform, the application of predictive deconvolution to data in the slant-stack domain for multiple suppression is discussed.

2.5 Field data examples

The deconvolution parameters now are examined using field data examples. We shall discuss application of statistical deconvolution to pre- and poststack data. Additionally, we shall discuss application of deterministic deconvolution to marine data to convert the recorded source signature to its minimum-phase equivalent, and to land data recorded using a vibroseis source to convert the autocorrelogram of the sweep signal to its minimum-phase equivalent.

• 

  Figure 2.4-36  (a) A common-shot gather, (b) after t²–scaling and (c) after muting first arrivals which are largely guided waves, and (d) after spiking deconvolution. The amplitude spectra averaged over the shot record are shown at the top and the autocorrelograms are shown at the bottom.

• 

  Figure 2.5-1  An autocorrelation window test used to design deconvolution operators. The solid bars indicate the window boundaries. The entire 6-s length was included in (a). The autocorrelograms are displayed beneath the records.

Prestack deconvolution

Figure 2.5-1 shows a CMP gather that contains five prominent reflections at around 1.1, 1.35, 1.85, 2.15, and 3.05 s. The gather also contains strong reverberations associated with these reflections. The examination of the deconvolution parameters will begin with an analysis of the time gate to estimate the autocorrelation function. A first gate selected may be the entire length (6 s) of the record as seen in panel (a). The solid lines on the CMP gathers refer to the gate start and end times. The autocorrelogram of the record is shown at the bottom of each panel. A second choice might be to exclude the deeper part of the record where ambient noise dominates. The start of the gate is chosen as the first arrival path as shown in panel (b). A third choice may be to exclude not only the deeper portion, but also the early part of the record that contains energy corresponding to the guided waves as shown in panel (c). These waves travel within the water layer and are not part of the signal reflected from the substrata.

• 

  Figure 2.5-2  Test of operator length. The corresponding autocorrelogram is beneath each record. The window used in autocorrelation estimation is shown in Figure 2.5-1c. (a) Input gather. Deconvolution using prediction lag = 4 ms (spiking deconvolution), 0.1% prewhitening, and prediction filter operator lengths (b) 40 ms, (c) 80 ms, (d) 160 ms, (e) 240 ms.

• 

  Figure 2.5-3  Test of prediction lag. The corresponding autocorrelogram is beneath each record. The window used in autocorrelation estimation is shown in Figure 2.5-1c. (a) Input gather. Deconvolution using prediction filter operator length = 160 ms, 0.1 percent prewhitening, and prediction lags (b) 12 ms, (c) 32 ms, (d) 64 ms, (e) 128 ms.

By comparing the autocorrelograms from these different windows, note that the third choice best represents the reverberatory character of the data as shown in panel (c) over most of the offsets. All of the traces in the autocorrelogram within approximately the first 150 ms have a common appearance. This early portion of the autocorrelogram characterizes the basic seismic wavelet contained in the data.

In general, the autocorrelation window should include part of the record that contains useful reflection signal, and should exclude coherent or incoherent noise. An autocorrelation function contaminated by noise is undesirable since the deconvolution process is most effective on noise-free data (assumption 4).

Another aspect of the autocorrelation window is length. Panel (d) of Figure 2.5-1 shows the autocorrelogram estimated from a narrow window. The autocorrelogram estimated from the narrower part of the time gate (the right side of the record) in some data cases may lack the characteristics of the reverberations, and even those of the basic seismic wavelet.

Figure 2.5-4  Test of percent prewhitening. The corresponding autocorrelogram is beneath each record. The window used in autocorrelation estimation is shown in Figure 2.5-1c. (a) Input gather. Deconvolution using prediction filter operator length = 160 ms, prediction lag = 4 ms (spiking deconvolution), and percent prewhitening (b) 1 percent, (c) 4 percent, (d) 16 percent, (e) 32 percent.

In general, any autocorrelation function is biased; that is, the first lag value is computed from, say, n nonzero samples, the second lag value is computed from n – 1 nonzero samples, and so on. If n is not large enough, then there can be an undesirable biasing effect. How large should the data window be to avoid such biasing? If the largest autocorrelation lag used in designing the deconvolution operator were m, an accepted rule of thumb is that the number of data samples should be no less than 8m.

Figure 2.5-5  (a) A common-shot gather, (b) after muting guided waves, (c) after t²-scaling, and (d) after spiking deconvolution using an operator length of 320 ms. The amplitude spectra (top) averaged over the shot record, and the autocorrelograms (bottom) are used to to choose deconvolution parameters and evaluate the data after the application of deconvolution.

Now that the autocorrelation window is determined, we examine operator length. In Figure 2.5-2, prediction lag (4 ms, the same as the sampling rate) and percent prewhitening (0.1%) are fixed. The autocorrelograms (at the bottom of each gather) are displayed for diagnostic purposes. From the analyses of the single spike, sparse spike, and reflectivity models (Predictive deconvolution in practice), recall that the short (40-ms) operator leaves some residual energy that corresponds to the basic wavelet and reverberating wavetrain in the record. For a spiking deconvolution with a 160-ms-long operator, no remnant of the energy is associated with the basic wavelet and reverberations. Any operator longer than 160 ms does not change the result, significantly. From the output of the 160-ms operator, note that the prominent reflections (at 1.1, 1.35, 1.85, and 2.15 s at the near-offset) have been uncovered, the seismic wavelet has been compressed, and the reverberations have been significantly suppressed.

The effect of prediction lag now is examined. In Figure 2.5-3, the 160-ms operator length and 0.1% prewhitening are fixed, while prediction lag is varied. If prediction lag were increased, the deconvolution process would be increasingly less effective in broadening the spectrum, and the autocorrelograms would contain increasingly more energy at nonzero lags. In the extreme, the deconvolution process is ineffective with a 128-ms prediction lag. In practice, common values for the prediction lag are unity (spiking deconvolution) or the first or second zero crossing of the autocorrelation function (predictive deconvolution).

Finally, the percent of prewhitening is varied, while the 4-ms prediction lag and 160-ms operator length are fixed. These tests are shown in Figure 2.5-4. By increasing the percent prewhitening, the deconvolution process becomes less effective. The high end of the spectrum is not flattened as much as the rest of the spectrum (Figure 2.4-24). Note that the autocorrelograms contain increasingly more energy at nonzero lags with increasing percent prewhitening. In practice, it is not advisable to assign a large percent of prewhitening. Typically, a value between 0.1 and 1 percent is sufficient to ensure stability in designing the deconvolution operator (equation 32).

We now examine the effect of operator length and prediction lag on amplitude spectrum. Figure 2.5-5a shows a common-shot gather with its autocorrelogram and average amplitude spectrum. The field record is prepared for deconvolution by first muting the guided waves (Figure 2.5-5b) and applying t²-scaling (Figure 2.5-5c). Figure 2.5-5d shows the same record after spiking deconvolution. Note the flattening of the spectrum within the passband of the recorded data and attenuation of the energy at nonzero lags of the autocorrelogram. With prediction lag greater than unity (Figure 2.5-6), for the same operator length, we note insufficient flattening at the high-frequency end of the spectrum. At larger prediction lags, note the insufficient flattening at the low-frequency end of the spectrum. A very large prediction lag causes the amplitude spectrum of the deconvolved data remain similar to that of the input data (compare Figures 2.5-6c with 2.5-5c).

Figure 2.5-6  The shot record in Figure 2.5-5c after predictive deconvolution using an operator length of 320 ms and a prediction lag of: (a) unit-prediction, (b) 8 ms, and (c) 24 ms. The amplitude spectra (top) averaged over the shot record, and the autocorrelograms (bottom) are used to to choose deconvolution parameters and evaluate the data after the application of deconvolution.

The data sometimes must be preconditioned for deconvolution. If the data were too noisy, then a wide band-pass filter could be necessary before deconvolution. If there is significant coherent noise in the data, dip filtering (Sections 6.2 and 6.3) can be applied before deconvolution so that coherent noise is not included in the autocorrelation estimate. Alternatively, time-variant spectral whitening (The problem of nonstationarity) can be applied to balance the spectrum before deconvolution.

Signature deconvolution

In marine seismic exploration, the far-field signature of the source array can be recorded. The idea is to apply a deterministic deconvolution to remove the source signature, then to apply predictive deconvolution. The convolutional model is given by

$x(t)=s(t)\ast w(t)\ast e(t),$

(41)

where s(t) is the source signature recorded in the far-field just before it travels down into the earth, which has an impulse response e(t). Since s(t) is recorded, an inverse filter can be deterministically designed, as discussed in Inverse filtering, then applied to the recorded seismogram to remove it from equation (41). The unknown wavelet w(t) includes the propagating effects in the earth and the response of the recording system. This remaining wavelet then is removed by the statistical method of spiking deconvolution as discussed in Optimum wiener filters. Compare equation (41) with equation (3a) and note that the old w(t) of equation (3a) is split into two parts — the source signature s(t), which is the known component, and the new w(t), which is the unknown component.

There are two ways to handle s(t). One way is to convert it to its minimum-phase equivalent followed by predictive deconvolution (Figure 2.5-7). Another way is to convert s(t) into a spike followed by predictive deconvolution (Figure 2.5-8). The process involves the following steps:

1. Estimate the minimum-phase equivalent of the recorded source signature by computing the spiking deconvolution operator (equation 39) and taking its inverse.
2. (b) Design a shaping filter to convert the source signature to its minimum-phase equivalent or a zero-delay spike (equation 30).
3. Apply the shaping filter to each trace in each recorded shot record.
4. Apply predictive deconvolution to output data from step (c).

The results shown in Figures 2.5-7 and 2.5-8 (panels (c)) should be compared with single-step statistical deconvolution (Figure 2.5-9). Since the source was not minimum-phase in this case, Figure 2.5-9b should be better than Figure 2.5-9d. Is it?

Actually, results of signature processing depend on the accuracy of the recorded signature. One should avoid signature processing of old marine data unless there exists a concurrently recorded source signature. Contemporary marine data almost always include the recorded source signature at each shot location. Figure 2.5-10a shows a recorded water-gun signature with its amplitude and phase spectra. The minimum-phase equivalent is shown in Figure 2.5-10b, and the result of signature deconvolution to convert the recorded waveform to its minimum-phase equivalent is shown in Figure 2.5-10c. While the amplitude spectrum is unaltered, the phase spectrum is minimum phase.

Application of signature processing described in Figure 2.5-10 to a recorded shot gather is shown in Figure 2.5-11. Again, note that signature processing aimed at converting the recorded source signature to its minimum-phase equivalent leaves the amplitude spectrum and autocorrelogram of the shot record unaltered (Figures 2.5-11a,b). We can also observe from the corresponding CMP-stacked sections that the process has not made any impact on the degree of vertical resolution; only a change in phase has taken place (Figure 2.5-12). Following the deterministic step of Figure 2.5-11b, statistical deconvolution is applied to flatten the spectrum (Figure 2.5-11c).

• 

  Figure 2.5-7  Signature processing. A shaping filter is designed to convert the recorded signature [s(t) in equation (41)] to its minimum-phase equivalent and applied to the input record (a). The output (b) has the same bandwidth as the input (a). The output (b) then has been processed by predictive deconvolution using operator length of 160 ms, and prediction lags of (c) 4 ms (spiking deconvolution), (d) 12 ms, (e) 32 ms. Examine the autocorrelograms (bottom) and note that those of (a) and (b) should be identical.

• 

  Figure 2.5-8  Signature processing. A shaping filter is designed to convert the recorded signature (s(t) in equation 41) to a spike and applied to the input record (a). The output (b) then has been processed by predictive deconvolution using an operator length of 160 ms and prediction lags of (c) 4 ms (spiking deconvolution), (d) 12 ms, and (e) 32 ms. Note from the autocorrelograms (bottom) that the output from signature processing (b) still contains a wavelet component — the w(t) component in equation (41), that still needs to be removed.

• 

  Figure 2.5-9  Signature processing compared with statistical deconvolution. (a) Input shot record. (b) A shaping filter is designed to convert the recorded signature (s(t) in equation 41) to its minimum-phase equivalent and applied to the input record followed by spiking deconvolution. (This panel is the same as Figure 2.5-7c.) (c) A shaping filter is designed to convert the recorded signature (s(t) in equation 41) to a spike and applied to the input record followed by spiking deconvolution. (This panel is the same as Figure 2.5-8c.) (d) Spiking deconvolution of the input record (a). The autocorrelograms (bottom) suggest that the wavelet compression is achieved in all three cases, (b), (c) and (d).

• 

  Figure 2.5-10  (a) A water-gun source signature and its amplitude and phase spectra (top and middle graphs), (b) its minimum-phase equivalent, and (c) the output of signature processing to convert (a) to (b).

Vibroseis deconvolution

The vibroseis source is a long-duration sweep signal in the form of a frequency-modulated sinusoid that is tapered on both ends. Just as a convolutional model was proposed for the marine seismogram given by equation (41), a similar convolutional model can be proposed for the vibroseis seismogram

$x(t)=s(t)\ast w(t)\ast e(t),$

(42)

where x(t) is the recorded seismogram, s(t) is the sweep signal, w(t) is the seismic wavelet with the same meaning as in equation (41), and e(t) is the earth’s impulse response. Convolutions in equation (42) become multiplications in the frequency domain:

$X(\omega )=S(\omega )W(\omega )E(\omega ).$

(43)

In terms of amplitude A(ω) and phase ϕ(ω) spectra, equation (43) yields

$A_{x}(\omega )=A_{s}(\omega )A_{w}(\omega )A_{e}(\omega )$

(44a)

and

$\phi _{x}(\omega )=\phi _{s}(\omega )+\phi _{w}(\omega )+\phi _{e}(\omega ).$

(44b)

Crosscorrelation of the recorded seismogram x(t) with the sweep signal s(t) is equivalent to multiplying equation ('44a) by A_s(ω) and subtracting ϕ_s(ω) from equation (44b). The correlated vibroseis seismogram x(t) therefore would have the following amplitude and phase spectra:

$A'(\omega )=A_{s}^{2}(\omega )A_{w}(\omega )A_{e}(\omega )$

(45a)

and

$\phi \prime (\omega )=\phi _{w}(\omega )+\phi _{e}(\omega ).$

(45b)

The inverse Fourier transform of $A_{s}^{2}\left(\omega \right)$ yields the autocorrelation of the sweep signal, which is called the Klauder wavelet k(t). Returning to the time domain, equations (45a,45b) yield

$x\prime (t)=k(t)\ast w(t)\ast e(t).$

(46)

Figure 2.5-13 outlines the process of vibroseis correlation where the seismic wavelet w(t) has been omitted for convenience. Note that, following vibroseis correlation, the sweep s(t) contained in the recorded seismogram x(t) is replaced with its autocorrelogram — the Klauder wavelet k(t).

Since it is an autocorrelation, the Klauder wavelet is zero-phase. Convolution of k(t) with the assumingly minimum-phase wavelet w(t) yields a mixed-phase wavelet. Because spiking deconvolution is based on the minimum-phase assumption, it cannot recover e(t) properly from vibroseis data.

One approach to deconvolution of vibroseis data is to apply a zero-phase inverse filter to remove k(t), followed by a minimum-phase deconvolution to remove w(t). The amplitude spectrum of the inverse filter is defined as ${1}/{A_{s}^{2}\left(\omega \right)}$. In practice, problems arise because of zeroes in the spectrum that are caused by the band-limited nature of the Klauder wavelet. Inversion of an amplitude spectrum, which has zeroes, yields an unstable operator (Optimum wiener filters). To circumvent this problem, a small percent of white noise, say 0.1%, usually is added before inverting the Klauder wavelet spectrum.

Another approach is to design a filter that converts the Klauder wavelet to its minimum-phase equivalent . A technique to compute the minimum-phase spectrum from a given amplitude spectrum is described in Section B.4 and is included in the discussion on frequency-domain deconvolution. If the Klauder wavelet were converted to its minimum-phase equivalent, then equation (45b) would take the form:

$\phi \prime (\omega )=\phi _{k}(\omega )+\phi _{w}(\omega )+\phi _{e}(\omega ).$

(47)

• 

  Figure 2.5-11  (a) A shot record and its amplitude spectrum and autocorrelogram (top and bottom, respectively), (b) after signature processing to convert the recorded water-gun signature to its minimum-phase equivalent as illustrated in Figure 2.5-10, and (c) followed by spiking deconvolution.

• 

  Figure 2.5-12  (a) A portion of a CMP-stacked section that corresponds to the data in Figure 2.5-11a, and (b) after signature processing of shot records as in Figure 2.5-11b to convert the recorded water-gun signature to its minimum-phase equivalent as described in Figure 2.5-10. Note that signature processing was not aimed at wavelet compression; instead, it was done to convert the source signature to its minimum-phase equivalent.

• 

  Figure 2.5-13  An outline of vibroseis correlation.

If we assume that w(t) is minimum-phase and if we make k(t) minimum-phase, then the result of their convolution also is minimum-phase. Spiking deconvolution now is applicable since the minimum-phase assumption is satisfied.

There is a 90-degree phase difference in some vibrator systems between the control sweep signal and the baseplate response. As an option, we may want to subtract out this phase difference. Figure 2.5-14 shows the recommended sequence of operations for vibroseis processing.

Figure 2.5-14  A flowchart for vibroseis deconvolution.

Figure 2.5-15 shows how the flowchart in Figure 2.5-14 is used with a synthetic reflectivity series. By including the step to convert the Klauder wavelet into its minimum-phase equivalent before spiking deconvolution, a closer representation of the impulse response is produced as seen by comparing steps (k) and (l) with (m).

Although sound in theory, the above scheme may have problems in practice. Fundamental issues, such as whether the convolutional model given in equation (42) really represents what goes on in the earth, are not resolved.

Vibroseis data often are deconvolved as dynamite data, without converting the Klauder wavelet to its minimum-phase equivalent. An example of deconvolution of a correlated vibroseis record is shown in Figure 2.5-16.

Despite the fact that the basic minimum-phase assumption is violated for vibroseis data as compared to explosive data, spiking deconvolution without conversion of the Klauder wavelet to its minimum-phase equivalent seems to work for most field data. Figure 2.5-17 shows a set of correlated vibroseis records before and after spiking deconvolution. Prominent reflections after deconvolution are enhanced and reverberations are attenuated significantly. Nonetheless, the problem of tying vibrator lines to lines recorded with other sources, say dynamite, is more difficult if the vibrator data have not been phase-corrected. Field systems now exist to do minimum-phase vibroseis correlation in the field.

Poststack deconvolution

Poststack deconvolution often is considered for several reasons. First, a residual wavelet almost always is present on the stacked section. This is because none of the underlying assumptions for deconvolution is completely met in real data; therefore, deconvolution never can completely compress the basic wavelet contained in prestack data to a spike. Second, since a CMP stack is an approximation to the zero-offset section, predictive deconvolution aimed at removing multiples may be a viable process after stack. Figure 2.5-18 is an example of poststack deconvolution applied to marine data. After deconvolution, the spectrum is flattened, albeit incompletely, the wavelet is compressed further and the marker horizons are better characterized. Again, as the prediction lag is increased, the flatness character of the spectrum and thus vertical resolution is increasingly compromised (Figure 2.5-19).

Figure 2.5-20 shows poststack deconvolution applied to land data. Note the significant improvement in vertical resolution as it can be verified by the autocorrelogram and average amplitude spectrum of the data.

2.6 The problem of nonstationarity

Figure 2.6-1 shows a CMP gather and its filtered versions before deconvolution. The filter scans show that there is signal between 10 to 40 Hz. Note that higher frequencies are confined to the shallower part of the gather. The same field record after spiking deconvolution is shown in Figure 2.6-2. Filter scans of the deconvolved record also are shown in this figure. A comparison of the records before and after deconvolution (Figures 2.6-1 and 2.6-2) demonstrates the effects of the process; in particular, compression of the wavelet and broadening of the spectrum. The input signal level above 40 Hz is relatively weaker than that below 40 Hz. Deconvolution has attempted to reduce the differences between the signal levels within different frequency bands by flattening the spectrum. The flattening, however, was more effective in the shallow part of the record than in the deeper part.

As discussed in Sections 1.4 and 2.1, the source wavelet is not stationary — its shape and bandwidth change with traveltime (Figure 2.1-2). Specifically, attenuation of high frequencies in the wavelet increases as waves travel deeper in the subsurface. Although multiwindow deconvolution was used in this case (Figure 2.6-2), spectral flattening was not achieved over the entire length of the record because of the large degree of non-stationarity of the data. Nonstationarity is primarily the result of the effects of wavefront divergence and frequency attenuation.

• 

  Figure 2.5-15  An illustration of vibroseis deconvolution as described in Figure 2.5-14. (a) Impulse response spike series; (b) vibroseis sweep signal (10 to 120-Hz bandwidth); (c) convolution of (a) with (b); (d) minimum-phase source wavelet; (e) convolution of (d) with (c); (f) crosscorrelation of the sweep signal (b) with (e); (g) Klauder wavelet; (h) minimum-phase equivalent of (g); (i) the operator that converts (g) to (h); (j) the operator (i) applied to (f); (k) spiking deconvolution of (f); (l) spiking deconvolution of (j); (m) 10 to 120-Hz zero-phase band-pass filtered version of (a).

• 

  Figure 2.5-16  Deconvolution of a vibroseis record. Three windows were used. (a) Correlated vibroseis record and its autocorrelograms (b); (d) spiking deconvolution output and its autocorrelograms (c).

• 

  Figure 2.5-17  (a) Common-shot gathers recorded with vibroseis source. Geometric spreading correction and trace balancing have been applied. (b) After spiking deconvolution.

• 

  Figure 2.5-18  (a) A portion of a CMP-stacked section, and after spiking deconvolution using an operator length of: (b) 120 ms, (b) 160 ms, (c) 220 ms, and (d) 320 ms. Note from the autocorrelograms (bottom) that much of the reverberating energy is attenuated using a 320-ms operator length. The amplitude spectra averaged over the CMP stack (top) indicate that, irrespective of how long the operator length is, spiking deconvolution in this case has failed to flatten the spectrum completely within the passband. This is because of nonstationarity of the signal.

• 

  Figure 2.5-19  A portion of a CMP-stacked section as in Figure 2.5-7a after predictive deconvolution using an operator length of 320 ms and a prediction lag of: (a) 8 ms, (b) 12 ms, (c) 24 ms, (d) 32 ms, and (e) 48 ms. The amplitude spectra (top) averaged over the CMP stack, and the autocorrelograms (bottom) are used to to choose deconvolution parameters and evaluate the data after the application of deconvolution.

• 

  Figure 2.5-20  (a) A portion of a CMP-stacked section, and (b) after spiking deconvolution using an operator length of 240 ms. The amplitude spectra (top) averaged over the CMP stack, and the autocorrelograms (bottom) are used to choose deconvolution parameters and evaluate the data after the application of deconvolution.

The phenomenon of nonstationarity on stacked data is exemplified in Figure 2.6-3. Following a single-gate poststack spiking deconvolution, note that the average amplitude spectrum of the data does not indicate a flat spectrum (Figure 2.6-3b). Instead, the spectral behavior is similar to data with predictive deconvolution with a large prediction lag (Figure 2.6-3c). The strong reflector in the neighborhood of 2.5 s separates the zone with two different bandwidths — a broad-band signal zone above, and a relatively narrow-band zone below.

The usual approach to reduce nonstationarity is to apply processes designed to compensate for the effects wavefront divergence and frequency attenuation before deconvolution. Wavefront divergence is corrected for by applying a geometric spreading function (Gain applications). As yet, a method to compensate for attenuation has not been discussed. Attenuation is measured by a quantity called the quality factor Q. An infinite Q means that there is no attenuation. This factor can change in depth and in the lateral direction. If we had an analytic form for an attenuation function, then it would be easy to compensate for its effect. Several models for Q have been proposed. The constant-Q model is quite plausible and the easiest to deal with . However, the big problem of estimating Q from seismic data still remains. If a reliable Q value is available, say from borehole measurements, then, as will be discussed later in this section, inverse Q filtering can be applied to data to compensate for the frequency attenuation.

Time-variant deconvolution

Nonstationarity was discussed in Sections 1.4 and 2.1. The time-variant character of the seismic wavelet (Figure 2.1-2) often requires a multiwindow deconvolution. Figure 2.6-4 is a field record that was deconvolved by using three time gates. The autocorrelograms from gates 1, 2, and 3 are shown in Figure 2.6-5. Note the difference in character of the reverberatory energy from one gate to another. The shallow gate (1) has more high-frequency signal than the middle gate (2); while the middle gate has more high-frequency signal than the deeper gate (3). For best results, we must design different deconvolution operators from different parts of the record and apply them to the corresponding time gates. Up to three windows usually are sufficient to handle the nonstationary character of the seismic signal.

• 

  Figure 2.6-1  A field record (far left panel) with its band-pass filtered versions.

• 

  Figure 2.6-2  Spiking deconvolution applied to the field record in Figure 2.6-1 (far left panel), followed by application of a series of band-pass filters.

• 

  Figure 2.6-3  (a) A portion of a CMP-stacked section with its amplitude spectrum averaged over the CMP range (top) and autocorrelogram (bottom); (b) after time-invariant spiking deconvolution; and (c) after time-invariant predictive deconvolution with a prediction lag of 24 ms.

• 

  Figure 2.6-4  Three-window deconvolution. The solid bars indicate window boundaries. With data from each window, a deconvolution operator is designed and applied to the data in that window. The operators are blended across the window boundaries. (a) Input gather. Deconvolution using operator length of 160 ms and prediction lags of (b) 4 ms, (c) 12 ms, (d) 32 ms.

• 

  Figure 2.6-5  Autocorrelograms associated with the records shown in Figure 2.6-4. Note the difference in character between windows 1, 2, and 3 (from top to bottom in the far left panel). Windows are indicated in Figure 2.6-4a.

• 

  Figure 2.6-6  Spiking deconvolution (c) on a shot record (a) followed by band-pass filtering (d). (b) Autocorrelograms before and after deconvolution.

Another example of single- and multiwindow deconvolution is shown in Figures 2.6-6 and 2.6-7. Here, autocorrelograms from different gates do not show significant variations. Therefore, it probably does not make any difference whether a single or multigate deconvolution is used. In Figures 2.6-6 and 2.6-7, the record is shown after deconvolution followed by a wide bandpass filter application. Since the amplitude spectrum of the input data is flattened as a result of spiking deconvolution, both the high-frequency ambient noise as much as the high-frequency components of the signal are boosted. Therefore, the output of spiking deconvolution often is filtered with a wide band-pass operator.

A practical problem with time-variant deconvolution is limiting design gates to small time windows. Consider, for instance, a three-window deconvolution of a 5s data. This means that at best an average gate kength of 1.5 s at near offset and less than 1 s at far offset can be used to design a deconvolution operator. To attain good statistics in an autocorrelation estimate, an operator length no more than one-eighth to one-tenth of the design gate, say 150 ms, should be considered. Hence, if you need to use a longer operator, time-variant deconvolution may have limited effectiveness in attenuating reverberations and short-period multiples. A way to account for nonstationarity while avoiding the short-operator effect of multiwindow deconvolution follows.

Time-variant spectral whitening

Frequency attenuation and a way to compensate for it are illustrated in Figure 2.6-8. Let us assume that we have an input seismogram with amplitudes decaying in time, as depicted. Now apply a series of narrow bandpass filters to this trace. Examine the field record in Figure 2.6-1 and associate the filter panels with the traces sketched in Figure 2.6-8. Note that the low-frequency component of trace F_L has a lower decay rate than the moderate-frequency component F_M. Likewise, the moderate-frequency component F_M has a lower decay rate than the high-frequency component of the signal F_H. A series of gain functions, such as G₁, G₂, G₃, can be computed to describe the decay rates for each frequency band. This is done by computing the envelope of the band-pass filtered traces (Figure 2.6-8). The inverses of these gain functions then are applied to each frequency band and the results are summed. The amplitude spectrum of the resulting trace has thus been whitened in a time-variant manner. This time-variant spectral whitening process is outlined in Figure 2.6-9. The number of the filter bands, the width of each band, and the overall bandwidth of application of time-variant spectral whitening are parameters that can be prescribed for a particular application.

• 

  Figure 2.6-7  Three-window deconvolution on the same shot record as in Figure 2.6-6. In this case, there is no significant difference between the characters of the autocorrelograms estimated from three windows. (a) Input gather; (b) autocorrelograms before and after spiking deconvolution; (c) three-window spiking deconvolution on (a); (d) band-pass filtering on (c).

• 

  Figure 2.6-8  A schematic illustration of the rate of decay of the frequencies in a seismic trace .

• 

  Figure 2.6-9  A flowchart for time-variant spectral whitening.

Figures 2.6-10 and 2.6-11 show some field records before and after spiking deconvolution, respectively. Note that this process not only has compressed the wavelet, but also has tried to suppress any reverberations in the data. Refer to Figure 2.6-12 and note that time-variant spectral whitening mainly has compressed the wavelet without changing much of the ringy character of the data. Also note that little was done explicitly to the phase. Therefore, the action of time-variant spectral whitening may be close to a zero-phase deconvolution, although there is no rigorous theoretical proof of this.

In practice, one of the main differences between time-variant spectral whitening and conventional deconvolution is that the former seems to be able to do a better job of flattening the amplitude spectrum. This difference can be significant for broad-band data with large dynamic range.

Time-variant spectral whitening sometimes helps attenuate ground roll on land records by way of its spectral balancing effect. Note that, in Figure 2.6-13, spiking deconvolution with different operator lengths has failed to flatten the spectrum, sufficiently. On the other hand, following spiking deconvolution, application of time-variant spectral whitening has balanced the spectrum and thus attenuated the ground-roll energy.

The ability of time-variant spectral whitening in flattening the spectrum within the passband of stacked data is observed in Figure 2.6-14. Note that spiking deconvolution is fairly effective, but not sufficient, for attaining a flat spectrum (Figure 2.6-14b). Following time-variant spectral whitening, the spectrum is flattened within the passband of the data as seen in Figure 2.6-14c.

It is a requirement to prepare stacked data input to amplitude inversion with the broadest possible bandwidth and flattest possible spectrum. Hence, a processing sequence tailored for amplitude inversion almost always includes poststack deconvolution and time-variant spectral flattening steps.

Frequency-domain deconvolution

Spectral flattening can be achieved by an alternate approach in the frequency domain. As discussed in Section B.4, minimum-phase spiking deconvolution can be formulated in the frequency domain. Alternatively, we can flatten the amplitude spectrum without modifying the phase. This is called zero-phase frequency-domain deconvolution. When performed over multiple time gates down the trace, it is essentially equivalent to time-variant spectral flattening. If we only want to flatten the spectrum, then the approach shown in Figure 2.6-15 can be taken.

Although the domains of operations may differ, both minimum-phase frequency-domain and Wiener-Levinson deconvolution techniques should yield equivalent results. Differences between the results shown in Figures 2.6-11 and 2.6-16 mainly are due to their computational aspects.

The zero-phase frequency-domain deconvolution aimed at achieving time-variant spectral whitening requires partitioning the input seismogram into small time gates, as well as designing and applying the process described in Figure 2.6-15 to each gate, individually. Figure 2.6-17 shows the field records after zero-phase frequency-domain deconvolution. The output is comparable to the time-variant spectral flattening output shown in Figure 2.6-12.

Figure 2.6-10  CMP gathers with no deconvolution.

Inverse Q filtering

Frequency attenuation caused by the intrinsic attenuation in rocks was discussed in Sections 1.4 and 2.1. Attenuation causes loss of high frequencies in the propgating waveform with increasing traveltime. This gives rise to a nonstationary behavior in the shape of the wavelets associated with reflection events at different times.

Wave attenuation usually is described by a dimensionless factor Q, which is defined by the ratio of the mean stored energy to the energy loss over a period of time that is equivalent to one cycle of a frequency component of the waveform . Time-variant deconvolution and time-variant spectral whitening discussed in this section are processes that can correct for the time-varying effects of attenuation by spectral flattening. A deterministic alternative to compensate for frequency-dependent attenuation is provided by inverse Q filtering.

The amplitude spectrum of the inverse Q filter is given by (Section B.9)

$A(\omega ,\tau )=\exp({\frac {\omega \tau }{2Q}}),$

(48)

where ω is the angular frequency component associated with the input trace and τ is the time variable associated with the output trace from inverse Q filtering. The phase spectrum of the inverse Q filter usually is assumed to be minimum-phase, which can be computed by taking the Hilbert transform of the amplitude spectrum (Section B.4).

• 

  Figure 2.6-11  Spiking deconvolution applied to the CMP gathers in Figure 2.6-10.

• 

  Figure 2.6-12  Time-variant spectral whitening (TVSW) applied to the CMP gathers in Figure 2.6-10. Compare this with Figures 2.6-11, 2.6-16, and 2.6-17.

• 

  Figure 2.6-13  (a) A land record with its average amplitude spectrum (top) and autocorrelogram (bottom); (b) t²-scaling, (c) AGC scaling of amplitudes. Spiking deconvolution with operator length of (d) 160 ms, (e) 240 ms, and (f) 320 ms. (g) Spiking deconvolution using an operator length of 240 ms followed by time-variant spectral whitening.

• 

  Figure 2.6-14  (a) A portion of a CMP stack and its amplitude spectrum averaged over the CMP range (top) and autocorrelogram (bottom); (b) after spiking deconvolution, (c) followed by time-variant spectral whitening, and (d) time-variant filtering.

• 

  Figure 2.6-15  A flowchart for a frequency-domain deconvolution.

Application of the inverse Q filter requires knowledge of the attenuation factor Q, which usually is assumed to be constant. A compilation of laboratory measurements of Q for some rock samples is given by Table 2-28.

Note from Table 2-28 that most measurements have been made at extremely high frequencies compared to the typical frequency band of seismic waves. Nevertheless, by assuming frequency-independent Q factor , these measurements can still be considered useful. Also note that the Q factor can vary significantly for limestone, sandstone, and shale rock samples of different origin.

The inverse of the amplitude spectrum defined by equation (48) can be used to obtain a quantitative measure of attenuation. In terms of frequency f, wave velocity v and depth z = vτ, the inverse is

$A^{-1}(f,v,z,Q)=\exp(-{\frac {\pi fz}{Qv}}).$

(49)

To determine how far in depth the wave has to travel before its amplitude reduces to, say, one-tenth of its amplitude at the surface z = 0, rewrite equation (49) as follows:

$z={\frac {2.3Qv}{\pi f}}.$

(50)

• 

  Figure 2.6-16  Minimum-phase frequency-domain deconvolution applied to the CMP gathers in Figure 2.6-10. Compare this with Figures 2.6-11, 2.6-12, and 2.6-17.

• 

  Figure 2.6-17  Zero-phase frequency-domain deconvolution applied to the CMP gathers in Figure 2.6-10. Compare this with Figures 2.6-11, 2.6-12, and 2.6-16.

Table 2-28. Intrinsic attenuation measurements in rocks adapted from .

Rock Type	Attenuation Constant, Q	Frequency Range (Hz)
Basalt	550	3,000-4,000
Granite	300	20,000-200,000
Quartzite	400	3,000-4,000
Limestone I	200	10,000-15,000
Limestone II	50	2-120
Limestone III	650	4-18,000
Chalk	2	150
Sandstone I	25	550-4,000
Sandstone II	125	20,000
Sandstone III	75	2,500-5,000
Sandstone IV	100	2-40
Shale I	15	75-550
Shale II	75	3,300-12,800

Note that the smaller the Q factor, the lower the velocity and the higher the frequency, the shallower the depth at which the wave amplitude decays to a fraction of the wave amplitude at z = 0. Table 2-29 lists the z values for a frequency of 30 Hz and a velocity of 3000 m/s for a range of Q values.

Table 2-29. Depth at which wave amplitude drops to one-tenth of its original at the surface for a range of Q values (equation 50).
v = 3000 m/s
f = 30 Hz

Q Factor	Depth in m
25	1,830
50	3,660
100	7,325
250	18,312
500	36,625

Note that the smaller the Q factor the shallower the depth at which the amplitude drops to the specified value of one-tenth of the original value at the surface z = 0. For very large Q values, that is, for small attenuation, the amplitude reduction to the specified value does not take place until the wave reaches very large depths beyond the exploration objectives.

Unfortunately, there is no reliable way to estimate the attenuation factor Q directly from seismic data. At best, inverse Q filtering can be applied to post- or prestack data (Section B.9) using a range of constant Q factors to create a Q panel, much like a filter panel (The 1-D Fourier transform). The factor that yields the flattest spectrum in combination with other signal processing applications — deconvolution and time-variant spectral whitening, is chosen as the optimum Q value.

Deconvolution strategies

Throughout the development of deconvolution theory, several alternatives have been proposed to better solve the deconvolution problem. Still, predictive deconvolution is used more than the other methods, although the minimum-phase and white reflectivity sequence assumptions have been key issues of concern.

Follow the common sequence for deconvolution of marine data in Figures 2.6-18 through 2.6-22. Note the presence of reverberations and short-period multiples in the CMP-stacked data with no deconvolution applied (Figure 2.6-18). Signature processing, in this case, was done to convert the recorded source signature to its minimum-phase equivalent (Figure 2.6-19). Therefore, aside from phase, there is no difference between the sections in Figures 2.6-18 and 2.6-19. Deconvolution before stack has helped attenuate reverberations and short-period multiples and, to some extent, compressed the basic wavelet (2.6-20). The additional step of poststack deconvolution has restored much of the high frequencies attenuated during stacking (Figure 2.6-21). Finally, time-variant spectral whitening has flattened the spectrum within the passband of the data (Figure 2.6-22) and yielded a crisp section with high resolution.

The same sequence can be followed in Figures 2.6-23 through 2.6-27. The CMP-stacked section includes reflections associated with a shallow, low-relief sedimentary strata (Figure 2.6-23). Following the signature processing (Figure 2.6-24), observe the gradual increase in the vertical resolution by prestack deconvolution (Figure 2.6-25), poststack deconvolution (Figure 2.6-26) and time-variant spectral whitening (Figure 2.6-27).

The following formal processing sequence for deconvolution theoretically should yield optimum results:

1. Apply a geometric spreading compensation function to remove the amplitude loss due to wavefront divergence.
2. Apply an exponential gain or minimum-phase inverse Q filter ; . This compensates for frequency attenuation.
3. Optionally apply signature processing to marine data. For vibroseis data, apply the filter that converts the Klauder wavelet to its minimum-phase equivalent.
4. Apply predictive deconvolution to compress the basic wavelet and attenuate reverberations and short-period multiples. If required, apply surface-consistent deconvolution (Section B.8). This accounts for the near-surface variations effect on the wavelet because of inhomogeneities in the vicinity of sources and receivers. In step (b), the vertical variations effect on the wavelet is handled.
5. Apply predictive deconvolution after stack to broaden the spectrum and attenuate short-period multiples.
6. Apply time-variant spectral whitening. This provides further flattening of the spectrum within the signal bandwidth without affecting phase.

• 

  Figure 2.6-18  A portion of a CMP-stacked section with no deconvolution.

• 

  Figure 2.6-19  The section in Figure 2.6-18 with signature processing to convert the recorded source signature to its minimum-phase equivalent.

• 

  Figure 2.6-20  The section in Figure 2.6-19 with prestack spiking deconvolution.

• 

  Figure 2.6-21  The section in Figure 2.6-20 with poststack spiking deconvolution.

• 

  Figure 2.6-22  The section in Figure 2.6-21 with time-variant spectral whitening.

• 

  Figure 2.6-23  A portion of a CMP-stacked section with no deconvolution.

• 

  Figure 2.6-24  The section in Figure 2.6-23 with signature processing to convert the recorded source signature to its minimum-phase equivalent.

• 

  Figure 2.6-25  The section in Figure 2.6-24 with prestack spiking deconvolution.

• 

  Figure 2.6-26  The section in Figure 2.6-25 with poststack spiking deconvolution.

• 

  Figure 2.6-27  The section in Figure 2.6-26 with time-variant spectral whitening.

The idea is to do as much deterministic deconvolution as possible. Inverse Q filtering, signature deconvolution, and the filter that converts the Klauder wavelet to its minimum-phase equivalent are deterministic operators. Any remaining issues then are handled by statistical means. However, for most data cases, just doing the geometric spreading correction followed by predictive prestack and poststack deconvolution is adequate.

Exercises

Exercise 2-1. Write the z-transform of wavelet $\left(1,\ 0,\ -{\frac {1}{4}}\right).$ Design a three-term inverse filter and apply it to the original. Hint: The z-transform of the wavelet can be written as a product of two doublets, (1, -  1/2) and (1, 1/2).

Exercise 2-2. Consider the following set of wavelets:

Wavelet A : (3, -2, 1)

Wavelet B : (1, -2, 3)

1. Plot the percent of cumulative energy as a function of time for each wavelet. Use Robinson’s energy delay theorem to determine the minimum- and maximum-phase wavelet.
2. Set up matrix equation (31) for each wavelet, compute the spiking deconvolution operators, then apply them.
3. Let the desired output be (0, 0, 1, 0). Set up matrix equation (30) for each wavelet, compute the shaping filters, and apply them. Find that the error for wavelet B with the delayed spike is smaller.

Exercise 2-3. Consider wavelet A in Exercise 2-2. Set up matrix equation (32), where ε = 0.01, 0.1. Note that ε = 0 already is assigned in Exercise 2-2. As the percent prewhitening increases, the spikiness of the deconvolution output decreases.

Exercise 2-4. Consider a multiple series associated with a water bottom with a reflection coefficient c_w and two-way time t_w. Design an inverse filter to suppress the multiples. [This is called the Backus filter ].

Exercise 2-5. Consider an earth model that comprises a water-bottom reflector and a deep reflector at two-way times of 500 and 750 ms, respectively. What prediction lag and operator length should you choose to suppress (a) water-bottom multiples, and (b) peg-leg multiples?

Exercise 2-6. Refer to Figure 2.6-9. Consider the following three bandwidths — low (F_L), medium (F_M) and high (F_H), for TVSW application:

F_L : 10 to 30 Hz

F_M : 30 to 50 Hz

F_H : 50 to 70 Hz

What kind of slopes should you assign to each bandwidth so that the output trace has an amplitude spectrum that is unity over the 10-to-70-Hz bandwidth?

Exercise 2-7. If the signal character down the trace changes rapidly (strong nonstationarity), should you consider narrow or broad bandwidths for the filters used in TVSW?

Exercise 2-8. Consider a minimum-phase wavelet and the following two processes applied to it:

1. Spiking deconvolution followed by 10-to-50-Hz zero-phase band-pass filtering.
2. Shaping filter to convert the minimum-phase wavelet to a 10-to-50-Hz zero-phase wavelet.

What is the difference between the two outputs?

Exercise 2-9. How would you design a minimum-phase band-pass filter operator?

Exercise 2-10. Consider (a) convolving a minimum-phase wavelet with a zero-phase wavelet, (b) convolving a minimum-phase wavelet with a minimum-phase wavelet, and (c) adding two minimum-phase wavelets. Are the resulting wavelets minimum-phase?

Exercise 2-11. Consider the sinusoid shown in Figure 1-1 (frame 1) as input to spiking deconvolution. What is the output?

Exercise 2-12. Order the panels in Figure 2.E-1 with increasing prediction lag.

Figure 2.E-1  The shot record shown in Figure 2.4-36c after predictive deconvolution using an operator length of 480 ms, and four different prediction lags (Exercise 2-12). The amplitude spectra averaged over the shot record are shown at the top and the autocorrelograms are shown at the bottom.

Appendix B: Mathematical foundation of deconvolution

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 kth layer is given by

$I_{k}={\rho _{k}}{v_{k}}.$

(B-1)

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

$c_{k}={\frac {I_{k+1}-I_{k}}{I_{k+1}+I_{k}}}.$

(B-2)

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

$c_{k}={\frac {v_{k+1}-v_{k}}{v_{k+1}+v_{k}}}.$

(B-3)

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

$x(t)=w(t)\ast e(t).$

(B-4)

By Fourier transforming both sides, we get

$X(\omega )=W(\omega )E(\omega ),$

(B-5)

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

$X(\omega )=A_{x}(\omega )\ \exp[i\phi _{x}(\omega )],$

(B–6a)

$W(\omega )=A_{w}(\omega )\ \exp[i\phi _{w}(\omega )],$

(B–6b)

and

$E(\omega )=A_{e}(\omega )\ \exp[i\phi _{e}(\omega )].$

(B–6c)

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

$A_{x}(\omega )=A_{w}(\omega )A_{e}(\omega )$

(B–7a)

and

$\phi _{x}(\omega )=\phi _{w}(\omega )+\phi _{e}(\omega ).$

(B–7b)

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

$A_{e}(\omega )=A_{0}=constant.$

(B-8)

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

$A_{x}(\omega )=A_{0}A_{w}(\omega ).$

(B-9)

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

$r_{e}(\tau )={\frac {1}{N}}\sum _{t=0}^{N-1}e_{t}e_{t+\tau },\quad \tau =0,\ \mp 1,\ \mp 2,\ \ldots ,\ \mp (N-1),$

(B-10)

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

$r_{e}(\tau )=0,\quad \tau =0$

(B–11a)

and

${{r}_{e}}\left(0\right)={{r}_{0}}=constant.$

(B–11b)

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:

${{r}_{0}}={{e}_{0}^{2}}+{{e}_{1}^{2}}+\cdots +{{e}_{N-1}^{2}}.$

(B-12)

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

$X\left(z\right)=W\left(z\right)E\left(z\right).$

(B-13)

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

${\overline {X}}\left({1}/{z}\;\right)={\overline {W}}\left({1}/{z}\;\right){\overline {E}}\left({1}/{z}\;\right),$

(B-14)

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

$X(z){\overline {X}}(1/z)=[W(z)E(z)][{\overline {W}}(1/z){\overline {E}}(1/z)].$

(B-15)

By rearranging the right side,

$X(z){\overline {X}}(1/z)=[W(z){\overline {W}}(1/z)][E(z){\overline {E}}(1/z)].$

(B-16)

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

$r_{x}=r_{w}\ast r_{e},$

(B-17)

where r_x, r_w, and r_e are the autocorrelations of the seismogram, seismic wavelet, and impulse response, respectively. Based on the white reflectivity series assumption (equation B-11), we have

${{r}_{x}}={{r}_{0}}{{r}_{w}}.$

(B-18)

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

$w(t)\ast f(t)=\delta (t),$

(B-19)

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

$f(t)=\delta (t)\ast {\frac {1}{w(t)}}.$

(B-20)

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

$W(z)=w_{0}+w_{1}z+w_{2}z^{2}+\cdots +w_{m}z^{m}.$

(B-21)

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

$F(z)={\frac {1}{W(z)}}.$

(B-22)

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

$F(z)=f_{0}+f_{1}z+f_{2}z^{2}+\cdots +f_{n}z^{n}+\cdots .$

(B-23)

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

$W(\omega )F(\omega )=1.$

(B-24)

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

$F(\omega )={\frac {1}{A_{w}(\omega )\ \ \exp[i\phi _{w}(\omega )]}}.$

(B-25)

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

$F(\omega )={A_{f}(\omega )\ \ \exp[i\phi _{f}(\omega )]},$

(B-26)

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

$A_{f}(\omega )={\frac {1}{A_{w}(\omega )}}$

(B–27a)

and

$\phi _{f}(\omega )=-\phi _{w}(\omega ).$

(B–27b)

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:

$R_{w}(z)=W(z){\overline {W}}(1/z),$

(B-28)

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

$W\left(z\right)F\left(z\right)=1,$

(B-29)

from which we get

$W(z)={\frac {1}{F(z)}}.$

(B-30)

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

$R_{w}(z)F(z)={\overline {W}}(1/z).$

(B-31)

Since the wavelet is a real-time function,

$r_{w}(\tau )=r_{w}(-\tau ).$

(B-32)

Consider the special case of a three-point inverse filter (f₀, f₁, f₂). 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

$F(z)=f_{0}+f_{1}z+f_{2}z^{2}.$

(B–33a)

The z-transform of its autocorrelogram r_w(τ) is, by way of equation (B-32),

$R_{w}(z)=\cdots +r_{2}z^{-2}+r_{1}z^{-1}+r_{0}+r_{1}z+r_{2}z^{2}+\cdots$

(B–33b)

The z-transform of w(t) is

$W\left(z\right)={{w}_{0}}+{{w}_{1}}z+{{w}_{2}}{{z}^{2}}+\cdots +{{w}_{m}}{{z}^{m}},$

therefore,

${\overline {W}}(1/z)={\bar {w}}_{0}+{\bar {w}}_{1}z+{\bar {w}}_{2}z^{2}+\cdots +{\bar {w}}_{m}z^{m}.$

(B–33c)

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

$(r_{2}z^{-2}+r_{1}z^{-1}r_{0}+r_{1}z+r_{2}z^{2})(f_{0}+f_{1}z+f_{2}z^{2})={\bar {w}}_{0}+{\bar {w}}_{1}z^{-1}+{\bar {w}}_{2}z^{-2}.$

(B-34)

To solve for (f₀, f₁, f₂), we identify the coefficients of powers of z. The coefficient of z⁰ yields

$r_{0}f_{0}+r_{1}f_{1}+r_{2}f_{2}={\bar {w}}_{0},$

the coefficient of z¹ yields

${{r}_{1}}{{f}_{0}}+{{r}_{0}}{{f}_{1}}+{{r}_{1}}{{f}_{2}}={\text{0}},$

while the coefficient of z² yields

${{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

${\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}}$

(B-35)

Note that ${{\bar {w}}_{0}},$ which equals w₀ 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₀, we get

${\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}}$

(B–36a)

where a₁ = f₁/f₀, a₂ = f₂/f₀, and L = w₀/f₀. We now have three unknowns, a₁, a₂, 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₀, f₁, f₂, …, f_n), equation (B-36a) takes the form

${\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}}$

(B–36b)

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² and CPU time that is proportional to n³. Because of the special properties of the Toeplitz matrix, Levinson devised a recursive scheme that requires a memory space and CPU time proportional to n and n², 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:

$R_{w}(\omega )=W(\omega ){\overline {W}}(\omega ),$

(B-37)

where ${\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(ω):

$U(\omega )=\ln[R_{w}(\omega )].$

(B-38)

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

$R_{w}(\omega )=\exp[U(\omega )].$

(B-39)

Suppose that another function, ϕ(ω), is defined and that equation (B-39) is rewritten as

$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\}$

(B-40)

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

$W(\omega )=\exp \left\{{\frac {1}{2}}[U(\omega )+i\phi (\omega )]\right\}.$

(B-41)

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

$U^{+}(\omega )={\frac {1}{2}}[U(\omega )+i\phi (\omega )],$

(B-42)

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,

$W(\omega )=A(\omega )\ \ \exp[i\phi (\omega )].$

(B-43)

The deconvolution filter in the Fourier transform domain is

$F(\omega )={\frac {1}{W(\omega )}}.$

(B-44)

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

$A_{f}(\omega )={\frac {1}{A(\omega )}}$

(B–45a)

and

$\phi _{f}(\omega )=-\phi (\omega ).$

(B–45b)

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-1  Flowchart for frequency-domain deconvolution.

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

$A_{f}(\omega )={\frac {1}{A(\omega )+\varepsilon }}.$

(B-46)

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

$L=\sum \limits _{t}(d_{t}-y_{t})^{2}.$

(B-47)

• 

  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 B-3  Wiener filter model.

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

$y_{t}=f_{t}\ast x_{t}.$

(B-48)

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

$L=\sum \limits _{t}\left(d_{t}-\sum \limits _{\tau }f_{\tau }x_{t-\tau }\right)^{2}.$

(B-49)

The goal is to compute the filter coefficients (f₀, f₁, …, 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:

${\frac {\partial L}{\partial f_{i}}}=0,\quad i=0,1,2,\ldots ,(n-1).$

(B-50)

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

$L=\sum \limits _{t}d_{t}^{2}-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}.$

(B-51)

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

${\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,$

(B-52)

or

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

(B-53)

By using

$\sum \limits _{t}x_{t-\tau }x_{t-i}=r_{i-\tau }$

(B–54a)

and

$\sum \limits _{t}d_{t}x_{t-i}=g_{i}$

(B–54b)

for each ith term, we get

$\sum \limits _{\tau }f_{\tau }r_{i-\tau }=g_{i},\quad i=0,1,2,\ldots ,(n-1).$

(B-55)

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

${\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}}$

(30)

Here, r_i are the autocorrelation lags of the input and g_i are the lags of the crosscorrelation between the desired output and input. Since the autocorrelation matrix is Toeplitz, the optimum Wiener filter coefficients f_i can be computed by using Levinson recursion .

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

$L_{\mathit {min}}=\sum \limits _{t}d_{t}^{2}-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}.$

(B-57)

By substituting the relationships

$\sum \limits _{t}d_{t}x_{t-\tau }=g_{\tau }$

(B–58a)

and

$\sum \limits _{t}x_{t}x_{t-\tau }=r_{\tau }$

(B–58b)

into equation (B-57), we get

$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}$

(B-59)

or

$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}.$

(B-60)

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

$L_{\it {min}}=\sum \limits _{t}d_{t}^{2}-\sum \limits _{\tau }f_{\tau }g_{\tau }.$

(B-61)

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₀, x₁, x₂, …), equation (B-56) takes the form

${\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}}$

(B-62)

Divide both sides by f₀ to obtain

${\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}}$

(B-63)

where a_i = f_i/f₀, i = 1, 2, …, n − 1, and v = x₀/f₀. Equation (B-63) is solved for the unknown quantity v and the filter coefficients (a₁, a₂, …, a_n−1). Since the desired output is a zero-delay spike, the filter (1, a₁, a₂, …, a_n−1) describes a spiking deconvolution process.

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

${\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}}$

(B-64)

Write out the simultaneous equations:

${{r}_{0}}+{{r}_{1}}{{a}_{1}}=v$

and

${{r}_{1}}+{{r}_{0}}{{a}_{1}}={\text{0}},$

which yield the filter coefficient a₁ and the unknown variable v from the first iteration:

${a_{1}}=-{r_{1}}/{r_{0}}$

(B–65a)

and

$v={{r}_{0}}+{{r}_{1}}{{a}_{1}}.$

(B–65b)

Now write equation (B-63) for the three-term filter $\left(1,\ {a'_{1}},\ {a'_{2}}\right):$

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

(B-66)

We want to solve for the filter coefficients $\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:

${\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

$e={r_{2}}+{r_{1}}{a_{1}}.$

(B-68)

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

${\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}}$

(B-69)

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

${\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}}$

(B-70)

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

${\begin{pmatrix}1\\a'_{1}\\a'_{2}\end{pmatrix}}={\begin{pmatrix}1\\a_{1}-ca_{1}\\-c\end{pmatrix}}$

(B–71a)

and

${\begin{pmatrix}v'\\0\\0\end{pmatrix}}={\begin{pmatrix}v-ce\\0\\e-cv\end{pmatrix}}$

(B–71b)

Solve for c and v':

$c={\frac {e}{v}},$

(B–72a)

and

$v'=v\left[1-\left({\frac {e}{v}}\right)^{2}\right].$

(B–72b)

Hence, the new filter coefficients $\left({a'_{1}},\ {a'_{2}}\right)$ are (equations B-71a)

$a'_{1}=a_{1}-{\frac {e}{v}}a_{1}$

(B–73a)

and

$a'_{2}=-{\frac {e}{v}}.$

(B–73b)

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

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

This recursive process yields the Wiener filter (1, a₁, a₂, …, 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

$g_{\tau }=\sum \limits _{t}d_{t}x_{t-\tau }=\sum \limits _{t}x_{t+\alpha }x_{t-\tau }=\sum \limits _{t}x_{t}x_{t-(a+\tau )}.$

(B-74)

By definition, we have

$r_{\tau }=\sum \limits _{t}x_{t}x_{t-\tau }.$

(B-75)

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

$r_{\alpha +\tau }=\sum \limits _{t}x_{t}x_{t-(\alpha +\tau )}.$

(B-76)

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₀, f₁, …, f_n−1):

${\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}}$

(B-77)

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

${\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}}$

(B-78)

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

${\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}}$

(B-79)

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

${\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}}$

(B-80)

We now have n + 1 equations and n + 1 unknowns — (f₀, f₁, …, f_n−1, L). From equation (B-80)

$L=r_{0}-r_{1}f_{0}-r_{2}f_{1}-\cdots -r_{n}f_{n-1}.$

(B-81)

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

${d_{t}}={x_{t+1}}$

(B-82)

so that

$\sum _{t}d_{t}^{2}=r_{0}$

(B–83a)

and

${g_{t}}={r_{t+1}}.$

(B–83b)

Substituting into equation (B-61), we get

$L_{min}=r_{0}-(r_{1}f_{0}+r_{2}f_{1}+\cdots +r_{n}f_{n-1}),$

(B-84)

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 ${{\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

$e_{t+\alpha }=x_{t+\alpha }-{\hat {x}}_{t+\alpha }=x_{t+\alpha }-\sum \limits _{\tau }f_{\tau }x_{t-\tau }.$

(B-85)

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

$z^{-\alpha }E(z)=z^{-\alpha }X(z)-F(z)X(z)$

(B-86)

or

$E(z)=[1-z^{\alpha }F(z)]X(z).$

(B-87)

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

$A(z)=1-z^{\alpha }F(z)$

(B-88)

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

$E\left(z\right)=A\left(z\right)X\left(z\right).$

(B-89)

The corresponding time-domain relationship is

$e(t)=a(t)\ast x(t).$

(B-90)

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

${{a}_{t}}=\left(1,\ \overbrace {0,\ 0,\ \cdots ,\ 0,} ^{\alpha -1}-{{f}_{0}},\ -{{f}_{1}},\ \cdots ,\ -{{f}_{n-1}}\right).$

(B-91)

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

$a_{t}=(1,-f_{0},-f_{1},\cdots ,-f_{n-1}).$

(B-92)

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

$x(t)=w(t)\ast e(t)+n(t),$

(B-93)

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

$x'_{ij}(t)=s_{j}(t)\ast h_{l}(t)\ast e_{k}(t)\ast g_{i}(t)+n(t),$

(B-94)

where ${{{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). 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):

$X'_{ij}(\omega )=S_{j}(\omega )H_{l}(\omega )E_{k}(\omega )G_{i}(\omega ).$

(B-95)

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

${\overline {X}}'_{ij}(\omega )={\overline {S}}_{j}(\omega ){\overline {H}}_{l}(\omega ){\overline {E}}_{k}(\omega ){\overline {G}}_{i}(\omega ),$

(B–96a)

and phase spectral components (Section A.1):

$\phi '_{ij}(\omega )=\phi _{sj}(\omega )+\phi _{hl}(\omega )+\phi _{ek}(\omega )+\phi _{ri}(\omega ).$

(B–96b)

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:

${\tilde {X}}'_{ij}(\omega )={\tilde {S}}_{j}(\omega )+{\tilde {H}}_{l}(\omega )+{\tilde {E}}_{k}(\omega )+{\tilde {G}}_{i}(\omega ).$

(B-97)

The left side is the logarithm of the amplitude spectrum ${{\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 based on the work by :

${\tilde {X}}''_{ij}(\omega )={\tilde {S}}_{j}(\omega )+{\tilde {H}}_{l}(\omega )+{\tilde {E}}_{k}(\omega )+{\tilde {G}}_{i}(\omega ),$

(B–98a)

where

${\tilde {X}}''_{ij}(\omega )={\tilde {X}}'_{ij}(\omega )-{\tilde {X}}'_{avg}(\omega ).$

(B–98b)

The spectral component ${{{\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_c channels, so that the total number of traces is n_s × n_c. For n_s × n_c values of the actual spectral components ${{\tilde {X}}_{ij}}$ at frequency ω, and n_s shot locations, n_r receiver locations, n_e midpoint locations, and n_h offsets, we have the following set of model equations: