Digital Imaging and Deconvolution:

The ABCs of Seismic Exploration and Processing

By Enders A. Robinson and Sven Treitel


Digital Imaging and Deconvolution: The ABCs of Seismic Exploration and Processing (SEG Geophysical References Series No. 15), covers the basic ideas and methods used in seismic processing, concentrating on the fundamentals of seismic imaging and deconvolution. Most chapters are followed by problem sets. Some exercises supplement textual material; others are meant to stimulate classroom discussions. Text and exercises deal mostly with simple examples that can be solved with nothing more than pencil and paper. The book covers wave motion; digital imaging; digital filtering; various visualization aspects of the seismic reflection method; sampling theory; the frequency spectrum; synthetic seismograms; wavelets and wavelet processing; deconvolution; the need for continuing interaction between the seismic interpreter and the computer; seismic attributes; phase rotation; and seismic attenuation. The last of the 15 chapters gives a detailed mathematical overview. Digital Imaging and Deconvolution, nominated for the Association of Earth Science Editors award for best geoscience publication of 2008–2009, will interest professional geophysicists, graduate students, and upper-level undergraduates in geophysics. The book also will be helpful to scientists and engineers in other disciplines who use digital signal processing to analyze and image wave-motion data in remote-detection applications. The methods described are important in optical imaging, video imaging, medical and biological imaging, acoustical analysis, radar, and sonar.

  1. Page 1

    The first offshore drilling for oil in Texas occurred along Goose Creek, 21 miles southeast of Houston on Galveston Bay. In 1903, John I. Gaillard noticed bubbles coming to the surface of the water. With a match, he confirmed that the bubbles were natural gas, a strong indication of oil deposits. The discovery well was drilled and hit oil on 2 June 1908, at 1600 ft. In 1916, a well at Goose Creek hit a 10,000-barrels-per-day (bbl/day) gusher at a depth of 2017 ft (Figure 1). Initially, that well produced 8000 bbl/day. The community changed overnight as men rushed to obtain leases, to build derricks, and to drill wells. Within two months, the well leveled off to 300 bbl/day. The largest well of the field was Sweet 16, which came in on 4 August 1917, gushing 35,000 bbl/day from a depth of 3050 ft. This well stayed out of control for three days before the crew could close it.

    The Goose Creek field is on a deep-seated salt dome with slightly arched overlying beds. When a hurricane hit in 1919, the Goose Creek oil field suffered tremendous property damage. The hurricane's relatively mild 39-mph winds destroyed more than 1450 oil derricks.

    At the time of the Goose Creek discoveries, the proper equipment for finding new oil fields included a Brunton compass, a K&E stadia handbook with Jacob's staff, a 7-ft stadia rod, a small bricklayer's hammer, and of course a couple of matches.

    In contrast, this book gives

  2. Page 39

    In the past, most seismic surveys were along surface lines, which yield 2D subsurface images. Because of great strides in computer technology and seismic instrumentation, exploration geophysics has made the transition from 2D to 3D processing.

    The wave equation behaves nicely in one dimension and in three dimensions but not in two dimensions. In one dimension, waves on a uniform string propagate without distortion. In three dimensions, waves in a homogeneous isotropic medium propagate in an undistorted way except for a spherical correction factor. However, in two dimensions, wave propagation is complicated and distorted. By its very nature, 2D processing never can account for events originating outside of the plane. As a result, 2D processing is broken up into a large number of approximate partial steps in a sequence of operations. These steps are ingenious, but they never can give a true image.

    On the other hand, 3D processing accounts for all of the events. It is now cost-effective to lay out seismic surveys over a surface area and to do 3D processing. No longer is the third dimension missing, so consequently, the need for a large number of piecemeal 2D approximations is gone. Prestack depth migration is a 3D imaging process that is computationally extensive but mathematically simple. The resulting 3D images of the interior of the earth surpass all expectations in utility and beauty.

    Reflection seismology is a remote-imaging method used in petroleum exploration. The seismic reflection method was developed in the 1920s. From then until about 1965,

  3. Page 79

    In petroleum exploration, the geophysicist's task is to look beneath the earth's surface in the search for new deposits of oil and natural gas. Subsurface geologic structures of interest can be several miles deep. Geophysicists use the seismic reflection method in their search for oil and natural gas. The exploration geophysicist illuminates the earth's subsurface by means of an energy source that generates seismic waves. The subsurface rock layers transmit and reflect those seismic waves. A seismic survey consists of collecting seismic reflection data over a selected geographic area, called the prospect. The seismic reflection method, when considered as an instrument for remote detection, has much in common with other disciplines that are based on noninvasive techniques for finding the structure of an inaccessible body, such as medical imaging and nondestructive testing.

    What is seismic acquisition? The essential features of seismic data acquisition are as follows: (1) At a fixed point on the surface of the earth, a source of energy — such as an array of dynamite charges or air guns or swept-frequency vibrators (as in vibroseis) — is activated. Such an activated source is called a shot. (2) Seismic waves from the shot propagate downward from the source point and go deep into the earth (Bois, 1968). (3) Eventually, the waves are reflected from geologic interfaces and propagate upward from those interfaces. A primary reflection is a reflection that travels directly down to the interface and then directly back up to the surface. A multiple reflection

  4. Page 99

    Discrete-time signals are sequences of numbers, with each number being identified by a fixed time instant. Such a series of data in a time sequence is called a time series. In other words, a time series xn is a series of data, with each data value xn being associated with a discrete, equally spaced time index n. The time index n is taken to be a whole number, or integer.

    Time series occur in all branches of science (Wold, 1938; Kolmogorov, 1941; Wiener, 1942). Economic data always appear in the form of numerical time series. Some meteorologic data, such as daily temperatures, are numerical time series; other meteorologic data, such as continuous barographic records, are continuous-time signals. Continuous-time signals appear in the engineering, biological, and physical sciences. Such continuous-time signals can be read (or measured, observed, or sampled) at equal intervals of time, thereby generating time series (Robinson and Silvia, 1979, 1980).

    Because a time series represents only the sampled values of a continuous-time signal, it provides only a limited description of the signal. By taking the sampling instants close enough together, the amount of information that is lost by replacing a well-behaved continuous-time function by a time series can be made small. A time spacing that is too gross would mean substantial information loss in the sampling process. At the other extreme, a time spacing that is too fine would mean substantial redundancy in information produced by the sampling

  5. Page 117

    What is digital filtering? The behavior of analog filters ordinarily is studied in the frequency domain. Digital filtering, on the other hand, is treated more fruitfully in the time domain. A digital filter is represented by its impulse response. The impulse response is made up of a sequence of numbers that act as weighting coefficients. The output of a digital filter is obtained by convolving the digitized input signal with the filter's impulse response.

    The mechanics of digital filtering in the time domain can be described with the aid of Z-transform theory. The amplitude spectrum and the phase spectrum represent an important characterization of the filter. A digital filter is said to be causal if its output at time n depends only on its input at time n and on inputs at times before n. In Chapter 6, these ideas are related to the more familiar interpretation of filter behavior in the frequency domain.

    What is a causal digital filter? As we just mentioned, a digital filter is represented by a sequence of numbers called its impulse response or its weighting coefficients. A digital filter is causal if its present output (at time n) depends only on present and past inputs (that is, depends only on inputs at times n, n − 1, n − 2, …, and so on). Another term for a causal filter is a realizable filter.

    What is a constant digital filter? A constant filter is one that has a single constant weighting coefficient

  6. Page 143

    What are the time and frequency domains? A filter's action can be described by its impulse response as well as by its frequency spectrum. The filter's impulse response is in the time domain, as is the input signal itself. The filter's frequency spectrum is in the frequency domain. Both modes of expression are functions of each other – that is, if one is known, the other can be derived from it.

    In digital filtering, either domain can be employed, but generally, both the seismic signal and the characteristics of the filter must be converted into the same form. For example, we need for the operation to be in the time domain, but only the frequency spectrum of the filter is specified. In such a case, the frequency spectrum must be transformed into an impulse response (in the time domain) so that operation can be carried out in the time domain. The Fourier transform and the inverse Fourier transform provide the physical basis for such conversions from one domain to the other.

    What is a Fourier transform? The Fourier transform converts a function of time (the signal) into the corresponding function of frequency (the temporal frequency spectrum). The inverse Fourier transform works in the reverse direction. Specifically, the inverse Fourier transform converts a function of frequency (the temporal frequency spectrum) into the corresponding function of time (the signal).

    The Fourier transform need not apply only to frequency in cycles per second and time in seconds. The Fourier transform also can be

  7. Page 159

    Faraday (1831) investigated the water patterns produced under vibration. Transitory ripples disclosed distinct patterns. From the ripples, Faraday discerned an oscillatory condition that proved useful in his subsequent investigations on light. In the same way, geophysicists must investigate each ripple on a seismic wavelet to unravel the deep secrets of the earth.

    What is a wavelet? A wavelet is a signal that has finite energy (Robinson, 1962, 1964a, 1964b). In other words, a wavelet is a waveform with the bulk of its energy confined to a finite interval on the time scale. A wavelet can be written aswhere bn is the coefficient at discrete time n. The present time instant n = 0 represents the critical point in the classification of a wavelet. Past times would be all the instants n 0 before the present time, and future times would be all the instants n > 0 after the present time.

    For the record, all of the digital filters that we consider (unless otherwise stated) fall under the category of linear time-invariant filters. Note that filter is called linear if it satisfies the additive property and the multiplicative property. Let a given input yield a given output. A filter is called time-invariant if the same input delayed (or advanced) by a given amount yields the same output delayed (or advanced) by the same amount (Robinson and Silvia, 1978).

    Recall our discussion in Chapter 4,

  8. Page 193

    We assume that a seismic trace has been corrected for amplitude decay resulting from spherical spreading over the seismic time scale of interest (say, for example, from 0 to 6 s). However, in reality, other effects also must be considered. One such effect is inelastic absorption — the loss of seismic energy to frictionally generated heat. (We will treat inelastic absorption in Chapter 14.) We must consider the effect of the seismic energy's source. In addition, effects result from the instrumentation; source and instrument effects are man-made at or near the surface of the ground. We lump these surface effects together in the form of a source wavelet, which we denote by s.

    However, the most important effect is that of the deep earth itself. The deep earth is represented by the sequence ε of reflection coefficients. This sequence is called the reflectivity (Robinson, 1957; Ulrych, 1999). In an ideal seismic experiment, an impulsive source (i.e., a sharp spike of energy imposed at time 0) produces the impulse response of the earth, which we denote by h. In the elastic range, the earth is a linear system, so we can write the synthetic trace (also called the synthetic seismogram) x as the convolution

    Equation 1 is called the dynamic convolutional model (Robinson, 1999). The model is called dynamic because the reflection impulse response h is a highly nonlinear function of the reflection coefficients. To deconvolve

  9. Page 217

    What is a wavelet? From a seismic processor's point of view, a wavelet is one of the basic building blocks used to construct the seismic models on which seismic-processing methods are based. In the various processing steps, the wavelets are removed from the seismic data to yield the final sections. Correct estimation and/or measurement of the wavelets allows the good results that are expected in any exploration program. The wavelet is a basic concept, and good wavelet estimation is fundamentally important in exploration geophysics.

    The seismic method is an instrument for remote detection that uses seismic traveling waves to delineate the subsurface structure of the earth. An exploration geophysicist illuminates the earth's subsurface by an energy source that generates those seismic waves. In a 3D earth, the waves travel in all directions, but to keep the present discussion simple, we consider the case of only vertically upgoing waves and vertically downgoing waves.

    Subsurface rock layers transmit and reflect the seismic waves, and seismic theory and practice deal with those traveling seismic waves. The simplest example of a traveling wave is a primary reflection. A primary reflection consists of the downgoing path from the source to the reflection horizon and the returning upgoing path from the reflector to the receiver. A multiple reflection is an event that bounces back and forth among various interfaces as it proceeds on its trip. Directional sensors can be used to record seismic waves, but usually the receiver is either a hydrophone (for exploration at sea

  10. Page 253

    Pythagoras (c. 580-c. 500 B.C.) taught that “all is number.” Pythagoras realized that numbers were hidden in everything, from the harmonies of music to the orbits of the planets. In other words, number and the nature of number make a thing clear either in itself or in its relation to other things. Today's world, with its digital computers, digital pictures, digital animation, digital television, digital telephones, digital regulators, and digital processing, attests to the foresight of Pythagoras. Pythagoras was instrumental in the development of the language of mathematics, which enabled him and others to describe the nature of the universe.

    In additional ways unforeseen by Pythagoras, everything is number. The great mathematician Carl Friedrich Gauss (1777–1855) once said, “Mathematics is the queen of science and number theory the queen of mathematics.” While he was still a teenager, Gauss was intrigued with numbers. At the age of 18, he thought up and justified the numerical method of least squares.

    Gauss's love for numerical calculations stirred his interest in astronomy. On New Year's Eve 1800–1801, Giuseppe Piazzi had discovered what he thought was a new planet (it was the asteroid Ceres). Because observers soon would lose sight of such a small object, it was important to calculate its elliptical orbit as soon as possible. Using only the few observations that had been made of the asteroid, Gauss calculated its orbit (reputedly by least squares) so accurately that astronomers could locate it again late in 1801 and early in 1802.

    In the finest

  11. Page 289

    Let k = (k0, k1,…, kN−1) denote the least-squares prediction filter for prediction distance α. The prediction-error filter results directly from the prediction filter. As we saw in equation 12 of Chapter 10, the prediction-error filter is f = (1, 0, 0,…, 0, − k0, −k1,…, −kN−1). The prediction-error operator has α−1 zeros that lie between the leading coefficient, which is 1, and the negative prediction-filter coefficients. These α−1 zeros constitute the gap.

    Let us begin by examining a prediction-error filter, which is in fact a deconvolution filter. The associated prediction-error series is the deconvolved signal. (See Appendix K, exercise 34, at the end of this chapter for a further description of the prediction filter and the prediction-error filter.) A prediction-error filter must be causal. A successfully deconvolved signal shows improved seismic resolution and provides an estimate of the reflectivity series. Depending on a specified prediction distance α, we distinguish between two types of predictive deconvolution: (1) spiking deconvolution, for which the prediction distance equals one time unit, and (2) gap deconvolution, for which the prediction distance is greater than one time unit.

    Let B(Z) be the Z-transform of a minimum-phase wavelet b. Then A(Z) = 1/B(Z) is the Z-transform of the inverse wavelet a = b−1. For prediction distance α, the head of b is h = (b0, b1, …, bα−1) and the tail is t = (bα, bα+1, …). For both the head and tail, the first coefficient is at time 0. Thus, the

  12. Page 329

    What is petrophysics? Petrophysics (or rock physics) is the branch of geology that is concerned with the physical properties and behavior of rocks. Oil exploration requires direct and indirect measurements of the underground rock strata, the totality of which commonly is referred to as “the geology.” Well cores provide actual rock samples extracted from a well. Downhole instruments placed in oil wells provide direct measurements, or well logs, of the rocks in place. Cores and well logs are used to determine the petrophysics of the strata to help reveal potential oil-bearing reservoirs. Of interest in petrophysics are things such as rock type (e.g., sandstone, shale, limestone, etc.) and rock characteristics (e.g., porosity, permeability, fracturing, etc.).

    The seismic method provides indirect measurements in the form of seismic records. The word indirect is appropriate because the seismic method uses a noninvasive technique (seismic waves and seismic recording instruments) to penetrate and record data from a remote body (the underground rocks). Similar techniques are radar, sonar, and ultrasonic medical imaging. Seismic records are made up of signals called traces, which record reflections of the seismic waves from the boundaries of different rock layers. Such rock-layer boundaries are called interfaces (or horizons). In some cases, horizons are approximately horizontal planar surfaces with low to moderate dips. In other cases, horizons are curved and contorted surfaces with steep dips.

    Most familiar to us are water waves. Ancient Greek philosophers, many of whom were interested in music, hypothesized that a parallel exists between water waves and

  13. Page 349

    Chapter 12 introduced the concept of instantaneous attributes. Instantaneous means that we consider properties that are associated with the trace at one instant of discrete time n. A trace (which is a real-valued signal) can be considered to be the real part of a complex trace

    In this representation, the trace xn can be called the in-phase signal, and the imaginary part yn is called the quadrature signal. The quadrature signal is the Hilbert transform of the in-phase signal. The complex trace can be written as a vector, which in polar form iswhere

    The length An is called the instantaneous amplitude, and the angle θn is called the instantaneous phase lead. For the rest of this chapter, we shall simply say phase instead of phase lead.

    What is phase rotation? Phase rotation refers to rotation of the complex-trace vector. Phase rotation for a specified phase angle is done as follows: A rotation of 90° (π/2 radians) transforms the complex trace, as given by equation 2 of Chapter 12, into the rotated complex trace

    We recognize the real part, namely un = −An sin θn, as the negative of the quadrature component (Figure 1).

    A 180° rotation gives the factor e = −1. Thus, the rotated complex trace is

  14. Page 363

    Attenuation is a reduction in the energy of a traveling wave as it propagates through a medium. Attenuation — the falloff of a wave's energy with distance — has three main causes: (1) transmission loss at interfaces because of reflection, diffraction, mode conversion, and scattering (Bowman, 1955); (2) geometric divergence effects as waves spread out from a source; and (3) absorption, which is the conversion of kinetic energy into heat by friction (note that kinetic energy is the energy of motion). Writers do not always distinguish between the terms attenuation and absorption. Here, we refer to attenuation in the general sense of energy loss to any cause, and we use the term absorption in the special sense of energy loss to heat.

    Transmission loss is a wave's energy loss as the wave travels through an interface. In transmission loss, the energy that is lost is diverted from the traveling wave of interest. There is no loss of total kinetic energy because the lost energy merely travels somewhere else. For example, when a wave meets an interface, some energy is reflected back from the interface, and only part of the wave's energy is transmitted though the interface.

    Mode conversion is the conversion of P-wave energy into S-wave energy or vice versa. Mode conversion occurs when a wave arrives at an interface at an obliquely incident angle to that interface. Converted waves divert energy away from the given wave.

    The energy of a wave in a homogeneous material is proportional

  15. Page 375

    In this chapter, we present the basic mathematical framework for the study of input-output models, along with some simple examples that show various approaches to the description and analysis of those models. The mathematical tools are derived from the general field of operational calculus, so for the most part, we are limited to the consideration of linear systems. Linear digital systems are described by linear difference equations, whereas linear analog systems are described by linear differential equations. Fortunately, many cases occur in engineering and science in which the systems either are linear or can be approximated sufficiently closely by a linear representation.

    Linear methods have been applied very successfully to the analysis of geophysical systems and the processing of geophysical data. The input and output of a system are related by a difference equation (digital case) or a differential equation (analog case), the solution of which gives the output for a given input. This equation provides a complete description of the system, but often it must be converted to other forms to be useful.

    Other modes of description of the system are related to the outputs produced by special types of inputs. Thus, we have

    1. the impulse response of the system, which is the output produced by an impulse input

    2. the frequency response of the system, which relates the outputs produced by sinusoidal inputs

    3. the transfer function (or system function), which is a generalization of the frequency-response function

    These modes of description are related to

Purchase Chapters

Recommended Reading