User:Ageary/Yilmaz 11

From SEG Wiki
Jump to: navigation, search
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


Contents

11.1 Introduction

In earth imaging in depth, we reviewed the two-dimensional (2-D) and three-dimensional (3-D), post- and prestack migration strategies for imaging the earth’s interior in depth. In earth modeling in depth, we learned traveltime inversion techniques for estimating a structural model of the earth that is needed to obtain an accurate image in depth. In structural inversion, structural inversion case studies for earth modeling and imaging in depth were presented. By structural inversion, we define the geometry of the reservoir unit, and the overlying and underlying depositional units. Traveltimes, however, are only one of the two components of recorded seismic wavefields; amplitudes are the other component.

In this chapter, we shall turn our attention to inversion of reflection amplitudes to infer petrophysical properties within the depositional unit associated with the reservoir rocks. The petrophysical properties include porosity, permeability, pore pressure, and fluid saturation. Specifically, we shall discuss prestack amplitude inversion to derive the amplitude variation with offset (AVO) attributes (analysis of amplitude variation with offset) and poststack amplitude inversion to estimate an acoustic impedance model of the earth (acoustic impedance estimation). The processes of estimating the acoustic impedance and AVO attributes by way of inversion of amplitudes may be appropriately referred to as stratigraphic inversion. Our goal ultimately is reservoir characterization based on structural and stratigraphic inversion of seismic data with calibration to well data.

We appropriately begin this chapter by investigating the resolution we can achieve from seismic data in defining vertical and lateral variations in the geometry of the reservoir unit. Resolution is the ability to separate two events that are very close together (seismic resolution). There are two aspects of seismic resolution: vertical (or temporal) and lateral (or spatial). Seismic resolution becomes especially important in mapping small structural features, such as subtle sealing faults, and in delineating thin stratigraphic features that may have limited areal extent.

Reservoir characterization involves calibration of the results of analysis of surface seismic data — both from structural and stratigraphic inversion, to well data. One category of well data includes various types of logs recorded in the borehole. Logs that are most relevant to seismic data are sonic, shear, and density. Another category of well data is a vertical seismic profile (VSP) (vertical seismic profiling).

Just as we can seismically characterize a reservoir, we also can seismically monitor its depletion. This is achieved by recording 3-D seismic data over the field that is being developed and produced at appropriate time intervals and detecting changes in the reservoir conditions; specifically, changes in petrophysical properties of the reservoir rocks, such as fluid saturation and pore pressure. Specifically, such changes may be related to changes in the seismic amplitudes from one 3-D survey to the next. Time-lapse 3-D seismic monitoring of reservoirs is referred to as the 4-D seismic method (4-D seismic method). The fourth dimension represents the calendar time over which the reservoir is being monitored.

Some reservoirs can be better identified and monitored by using shear-wave data. For instance, acoustic impedance contrast at the top-reservoir boundary may be too small to detect, whereas shear-wave impedance contrast may be sufficiently large to detect. By recording multicomponent data at the ocean bottom, P-wave and S-wave images can be derived. Commonly, four data components are recorded — the pressure wavefield, and inline, crossline, and vertical components of particle velocity. Thus, the multicomponent seismic data recording and analysis is often referred to as the 4-C seismic method (4-C seismic method).

This chapter ends with a brief discussion on anisotropy. While exploration seismology at large is based on the assumption of an isotropic medium, the earth in reality is anisotropic. This means that elastic properties of the earth vary from one recording direction to another. Seismic anisotropy often is associated with directional variations in velocities. For instance, in a vertically fractured limestone reservoir, velocity in the fracture direction is lower than velocity in the direction perpendicular to the plane of fracturing (azimuthal anisotropy). Another directional variation of velocities involves horizontal layering and fracturing of rocks parallel to the layering. In this case, velocity in the horizontal direction is higher than the vertical direction (transverse isotropy). In seismic anisotropy, we shall review seismic anisotropy in relation to velocity analysis, migration, DMO correction, and AVO analysis.

Elastic waves and rock properties

Seismic waves induce elastic deformation along the propagation path in the subsurface. The term elastic refers to the type of deformation that vanishes upon removal of the stress which has caused it. To study seismic amplitudes and thus investigate their use in exploration seismology, it is imperative that we review wave propagation in elastic solids. This gives us the opportunity to appreciate the underlying assumptions in estimating acoustic impedance and AVO attributes.

A summary of the elastic wave propagation theory is provided in Appendix L. To facilitate the forthcoming discussion on the link between elastic waves and rock properties, first, we summarize the definitions of elastic wave theory that should always be remembered.

  1. Stress is force per unit area. Imagine a particle represented by an infinitesimally small volume around a point within a solid body with dimensions (dx, dy, dz) as depicted in Figure L-1. The stress acting upon one of the surfaces, say dy − dz, can in general be at some arbitrary direction. It can, however, be decomposed into three components — one which is normal to the surface, and two which are tangential to the surface. The normal component of the stress is called the normal stress and the tangential components are called the shear stress. A normal stress component is tensional if it is positive and compressional if it is negative. Fluids cannot support shear stress. In a fluid medium, only one independent stress component exists — the hydrostatic pressure.
  2. Strain is deformation measured as the fractional change in dimension or volume induced by stress. Strain is a dimensionless quantity. The stress field away from the typical seismic source is so small that it does not cause any permanent deformation on rock particles along the propagation path. Hence, the strain induced by seismic waves is very small, usually around 10−6. Consider two points, P and Q, within a solid body as indicated in Figure L-2. Subject to a stress field, the solid is deformed in some manner and the particles at points P and Q are displaced to new locations P′ and Q′. Consider deformations of specific types illustrated in Figure L-3. The simplest deformation is the extension in one direction as a result of a tensional stress (Figure L-3a). The fractional change in length in a given direction is defined as the principal strain component. A positive strain refers to an extension and a negative strain refers to a contraction. Other types of deformation are caused by shearing (Figure L-3b), rotation (Figure L-3c), and a combination of the two (Figure L-3d). These angular deformations are called shear strains since they result in a shearing of the volume around a point within a solid body (Figure L-3b).
  3. Elastic deformation is a deformation in solid bodies that vanishes once the stress is released.
  4. Hooke’s law for elastic deformations states that the strain at any point is directly proportional to the stresses applied at that point.
  5. Elastic moduli are material constants that describe stress-strain relations:
    • Bulk modulus is the ratio of hydrostatic stress to volumetric strain; hence, it is a measure of incompressibility.
    • Modulus of rigidity is the ratio of shear stress to shear strain; hence, it is a measure of resistance to shear stress.
    • Young’s modulus is the ratio of the longitudinal stress to the longitudinal strain associated with a cylindrical rod that is subjected to a longitudinal extension in the axial direction. Since strain is a dimensionless quantity, Young’s modulus has the dimensions of stress.
    • Poisson’s ratio is the ratio of the lateral contraction to longitudinal extension associated with a cylindrical rod that is subjected to a longitudinal extension in the axial direction. Since strain is a dimensionless quantity, Poisson’s ratio is a pure number.
  6. Seismic waves are elastic waves that propagate in the earth.
  7. P-waves (or equivalently, compressional waves, longitudinal waves, or dilatational waves) are waves with particle motion in the direction of wave propagation.
  8. S-waves (or equivalently, shear waves, transverse waves, or rotational waves) are waves with particle motion in the direction perpendicular to the direction of wave propagation.
  9. Reflection is the wavefield phenomenon associated with the fraction of incident wave energy that is returned from an interface that separates two layers with different elastic moduli.
  10. Refraction is the the wavefield phenomenon associated with the fraction of incident wave energy that is transmitted into the next layer.
  11. Diffraction is the wavefield phenomenon associated with energy that propagates outward from a sharp discontinuity in the subsurface.

In exploration seismology, we are primarily interested in compressional and shear waves that travel through the interior of solid layers, and thus are characterized as body waves. Whereas in earthquake seismology, we also make use of Love and Rayleigh waves which travel along layer boundaries, and thus are characterized as surface waves.

Both body waves and surface waves are different forms of elastic waves, each associated with a specific type of particle motion. In the case of compressional waves, the particle motion induced by a compressional stress is in the direction of wave propagation. The compressional stress causes a change in the particle dimension or volume. The more the rock resists to the compressional stress, the higher the compressional wave velocity. In the case of shear waves, the particle motion induced by a shear stress is in the direction perpendicular to the direction of wave propagation. The shear stress does not cause a change in the particle dimension or volume; instead, it changes particle shape. The more the rock resists shear stress, the higher the shear wave velocity. Under the assumption that both wave types are elastic, whatever the change induced by the wave motion — the elastic deformation in particle shape, dimension or volume, vanishes once the wave motion on the particle vanishes and is propagated onto the neighboring particle.

Figure 11.0-1 outlines the interrelationships between the various elastic parameters. Starting with Young’s modulus — the ratio of principal stress to principal strain, and Poisson’s ratio — the ratio of shear strain to principal strain, Lamé’s constants — λ and μ are defined. The Lamé constant μ is indeed the modulus of rigidity and the Lamé constant λ = κ − (2/3)μ, where κ is the bulk modulus (Section L.3). Then, the two wave velocities — compressional (P-waves) and shear (S-waves), are derived in terms of the Lamé constants, or the bulk modulus and modulus of rigidity, and density.

From the definitions of the P- and S-wave velocities in Figure 11.0-1, note that both are inversely proportional to density ρ. At first thought, this means that the lower the rock density the higher the wave velocity. A good example is halite which has low density (1.8 gr/cm3) and high P-wave velocity (4500 m/s). In most cases, however, the higher the density the higher the velocity (Figure 11.0-2). This is because an increase in density usually is accompanied by an increase in the ability of the rock to resist compressional and shear stresses. So an increase in density usually implies an increase in bulk modulus and modulus of rigidity. Returning to the expressions for the P- and S-wave velocities in Figure 11.0-1, note that the greater the bulk modulus or the modulus of rigidity, the higher the velocity. Based on field and laboratory measurements, Gardner [1] established an empirical relationship between density ρ and P-wave velocity α. Known as Gardner’s formula for density, this relationship given by ρ = 0.25, where c is a constant that depends on the rock type, is useful to estimate density from velocity when the former is unknown. With the exception of anhydrites, most rock types — sandstones, shales, and carbonates, tend to obey Gardner’s equation for density.

In introduction to velocity analysis and statics corrections, a brief review of the results of some of the key laboratory experiments on seismic velocities was made. For a given lithologic composition, seismic velocities in rocks are influenced by porosity, pore shape, pore pressure, pore fluid saturation, confining pressure, and temperature. It is generally accepted that confining pressure, and thus the depth of burial, has the most profound effect on seismic velocities (Figure 3.0-3). For instance, the P-wave velocities for clastics can vary from 2 km/s at the surface up to 5 km/s and for carbonates from 3 km/s at the surface up to 6 km/s at depths greater than 5 km.

Because of the large variations in P-wave velocities caused by all these factors, P-wave velocity alone is not adequate to infer the lithology, unambiguously. The ambiguity in lithologic identification can be resolved to some extent if we have the additional knowledge of S-wave velocities. Here, we examine the ratio of the S-wave velocity to the P-wave velocity, β/α, which only depends on Poisson’s ratio σ (Figure 11.0-1). In some instances, we refer to the inverse ratio α/β. The higher the Poisson’s ratio, the higher the velocity ratio α/β. This relationship is supported by the physical meaning of Poisson’s ratio — the ratio of shear strain to principle strain. A way to describe the physical meaning of Poisson’s ratio is to consider a metal rod that is subject to an extensional strain. As the rod is stretched, its length increases while its thickness decreases. Hence, the less rigid the rock, the higher the Poisson’s ratio. This is exactly what is implied by the expression in Figure 11.0-1 that relates the modulus of rigidity μ to Poisson’s ratio σ. Unconsolidated sediments or fluid-saturated reservoir rocks have low rigidity, hence high Poisson’s ratio and high velocity ratio α/β. Here is the first encounter with a direct hydrocarbon indicator — the P- to S-wave velocity ratio. Ostrander [3] was the first to publish the link between a change in Poisson’s ratio and change in reflection amplitude as a function of offset.

Aside from the direct measurement of S-wave velocities down the borehole, there are three indirect ways to estimate the S-wave velocities. The first approach is to perform prestack amplitude inversion to estimate the P- and S-wave reflectivities and thus compute the corresponding acoustic impedances (analysis of amplitude variation with offset). The second approach is to record multicomponent seismic data and estimate the S-wave velocities from the P-to-S converted-wave component (4-C seismic method). The third approach is to generate and record S-waves themselves.

Figure 11.0-3 shows a plot of the S-wave slowness (inverse of the S-wave velocity) versus the P-wave slowness (inverse of the P-wave velocity) based on laboratory measurements [2]. Figure 11.0-4a shows a plot of the P-wave velocity to S-wave velocity based on full-waveform sonic logs from a producing oil field [4]. The key observation made from these results is that a lithologic composition may be associated with a reasonably distinctive velocity ratio α/β. The shale and limestone samples fall on a linear trend that corresponds to a velocity ratio of 1.9, whereas the dolomite samples have a velocity ratio of 1.8. The sandstone samples have a range of velocity ratio of 1.6 to 1.7. Lithologic distinction sometimes is more successful with a crossplot of P-wave to S-wave velocity ratio versus the P-wave velocity itself [4]. This is illustrated in Figure 11.0-4b which shows the same sample points as in Figure 11.0-4a.

Effect of shale and clay content on the velocity ratio α/β is an important factor in lithologic identification. Field and laboratory data from sandstone cores indicate that the velocity ratio α/β increases with increasing shale and clay content as a result of a decrease in S-wave velocity (Figure 11.0-5).

Finally, effect of porosity on the velocity ratio α/β is generally dictated by the pore shape. For limestones with their pores in the form of microcracks, the velocity ratio increases as the percent porosity increases [5]. For sandstones with their rounded pores, the velocity ratio does not increase as much with increasing porosity [4]. The difference between the rounded pores and microcracks lies in the fact that it is easier to collapse a rock with microcracks, hence lower modulus of rigidity.

11.2 Seismic resolution

Resolution relates to how close two points can be, yet still be distinguished. Two types of resolution are considered — vertical and lateral, both of which are controlled by signal bandwidth. The yardstick for vertical resolution is the dominant wavelength, which is wave velocity divided by dominant frequency. Deconvolution tries to increase the vertical resolution by broadening the spectrum, thereby compressing the seismic wavelet. The yardstick for lateral resolution is the Fresnel zone, a circular area on a reflector whose size depends on the depth to the reflector, the velocity above the reflector and, again, the dominant frequency. Migration improves the lateral resolution by decreasing the width of the Fresnel zone, thus separating features that are blurred in the lateral direction.

Vertical resolution

For two reflections, one from the top and one from the bottom of a thin layer, there is a limit on how close they can be, yet still be separable. This limit depends on the thickness of the layer and is the essence of the problem of vertical resolution.

The dominant wavelength of seismic waves is given by


(1)

where v is velocity and f is the dominant frequency. Seismic wave velocities in the subsurface range between 2000 and 5000 m/s and generally increase in depth. On the other hand, the dominant frequency of the seismic signal typically varies between 50 and 20 Hz and decreases in depth. Therefore, typical seismic wavelengths range from 40 to 250 m and generally increase with depth. Since wavelength determines resolution, deep features must be thicker than the shallow features to be resolvable. A graph of wavelength as a function of velocity for various values of frequency is plotted in Figure 11.1-1. The wavelength is easily determined from this graph, given the velocity and dominant frequency.

Table 11-1. Threshold for vertical resolution.
λ/4 = v/4f
v (m/s) f (Hz) λ/4 (m)
2000 50 10
3000 40 18
4000 30 33
5000 20 62

The acceptable threshold for vertical resolution generally is a quarter of the dominant wavelength. This is subjective and depends on the noise level in the data. Sometimes the quarter-wavelength criterion is too generous, particularly when the reflection coefficient is small and no reflection event is discernable. Sometimes the criterion may be too stringent, particularly when events do exist and their amplitudes can be picked with ease.

Table 11-1 contains the wavelength threshold values for vertical resolution, considering the realistic velocity and frequency ranges. For example, a shallow feature with a 2000-m/s velocity and 50-Hz dominant frequency potentially can be resolved if it is as thin as 10 m. A thinner feature cannot be resolved. Similarly, for a deep feature with a velocity as high as 5000 m/s and dominant frequency as low as 20 Hz, the thickness must be at least 62 m for it to be resolvable.

It is now appropriate to ask whether a thin stratigraphic unit must be resolved to be mapped. The answer is no. Resolution as defined here and in the geophysical literature implies that reflections from the top and bottom of a thin bed are seen as separate events or wavelet lobes. Using this definition, resolution does not consider amplitude effects. The thickness and areal extent of beds below the resolution limit often can be mapped on the basis of amplitude changes. This amplitude-based analysis can be especially precise when used for mapping gas-generated bright spots in tertiary rocks. Thus, in many stratigraphic plays, resolution in the strict sense is not an issue. For these plays, detection, not resolution, is the problem.

Vertical resolution is a concern when discontinuities are inferred along a reflection horizon because of faults. Figure 11.1-2 shows a series of faults with vertical throws that are equal to 1, 1/2, 1/4, 1/8, and 1/16th of the dominant wavelength. It is when the throw is equal or greater then one-fourth of the dominant wavelength that the presence of the fault can be inferred easily. Perhaps a smaller throw can be inferred by using diffractions from faults along the reflection horizon, provided the noise level is low in the data.

Figure 11.1-1  The relationship between velocity, dominant frequency, and wavelength. Here, wavelength = velocity/frequency. (Adapted from Sheriff, 1976; courtesy American Association of Petroleum Geologists.)

Clearly the ability to resolve or detect small targets can be increased by increasing the dominant frequency of the stacked data. The dominant frequency of a stacked section from a given area is governed by the physical properties of the subsurface, processing quality, and recording parameters. Since we cannot control the subsurface properties, the high-frequency signal level can only be influenced by the effort put into recording and processing.

The emphasis in recording should be to preserve high frequencies and suppress noise. The sampling rate and antialiasing filters should be adequate to record the desired frequencies. Receiver arrays should be small enough to prevent the significant loss of high-frequency signal because of intragroup moveout and statics. However, the arrays should not be too small, since small arrays are not as effective at suppressing random, high-frequency noise (wind noise) as large arrays. Finally, the source effort should be high enough to provide adequate signal level relative to noise level within the desired frequency band. Unless the signal-to-noise ratio of the field data is above some minimal level, say 0.25, processing algorithms have difficulty in recovering the signal. The signal has to be detectable before it can be enhanced.

The emphasis in processing should be to preserve and display the high-frequency signal present in the input data. Filters with good high-frequency response should be used for interpolation processes such as NMO removal, datum and statics corrections, and multiplex skew corrections. Extra care should be taken to ensure that small-scale residual moveout or statics, which might cause loss of high-frequency signal during stacking, are removed before stack. Nonsurface-consistent alignment programs (often called trim statics programs) sometimes are used for this purpose. Finally, care must be taken to ensure that all the high-frequency signal is displayed on the final stack. Poststack deconvolution is a useful tool for this purpose.

Figure 11.1-2  Faults with different amounts of vertical throws expressed in fractions of the dominant wavelength.

Lateral resolution

Lateral resolution refers to how close two reflecting points can be situated horizontally, yet be recognized as two separate points rather than one. Consider the spherical wavefront that impinges on the horizontal planar reflector AA′ in Figure 11.1-3. This reflector can be visualized as a continuum of point diffractors. For a coincident source and receiver at the earth’s surface (location S), the energy from the subsurface point (0) arrives at t0 = 2z0/v. Now let the incident wavefront advance in depth by the amount λ/4. Energy from subsurface location A, or A′, will reach the receiver at time t1 = 2(z0 + λ/4)/v. The energy from all the points within the reflecting disk with radius OA′ will arrive sometime between t0 and t1. The total energy arriving within the time interval (t1t0), which equals half the dominant period (T/2), interferes constructively. The reflecting disk AA′ is called a half-wavelength Fresnel zone [6] or the first Fresnel zone [7]. Two reflecting points that fall within this zone generally are considered indistinguishable as observed from the earth’s surface.

Since the Fresnel zone depends on wavelength, it also depends on frequency. For example, if the seismic signal riding along the wavefront is at a relatively high frequency, then the Fresnel zone is relatively narrow. The smaller the Fresnel zones, the easier it is to differentiate between two reflecting points. Hence, the Fresnel-zone width is a measure of lateral resolution. Besides frequency, lateral resolution also depends on velocity and the depth of the reflecting interface — the radius of the wavefront is approximated by (Exercise 11-1)


(2a)

In terms of dominant frequency f (equation 1), the Fresnel-zone width is


(2b)
Table 11-2. Threshold for lateral resolution (first Fresnel zone).
t0 (s) v (m/s) f (Hz) r (m)
1 2000 50 141
2 3000 40 335
3 4000 30 632
4 5000 20 1118

Table 11-2 shows the Fresnel zone radius, where r = OA′ in Figure 11.1-3 for a range of frequency and velocity combinations at various depths t0 = 2z/v. From Table 11-2, note that the shallower the event (and the higher the dominant frequency), the smaller the Fresnel zone. Since the Fresnel zone generally increases with depth, spatial resolution also deteriorates with depth.

Figure 11.1-4 shows reflections from four interfaces, each with four nonreflecting segments. The actual sizes of these segments are indicated by the solid bars on top. On the seismic section, the reflections appear to be continuous across some of these segments. This is because the size of some of the nonreflecting segments is much less than the width of the Fresnel zone; they are beyond the lateral resolution limit.

Spatial resolution is better understood in terms of diffractions. Note that in Figure 11.1-4, the diffraction energy is smeared across the nonreflecting segments on the deeper reflectors. Since migration is the process that collapses diffractions, it is reasonable to think that migration increases spatial resolution. Remember that migration can be achieved by downward continuation of receivers from the surface to the reflecting horizons. As a result of downward continuation, the observation points get closer and closer to the reflection points and, therefore, the Fresnel zone gets smaller and smaller. A smaller Fresnel zone means a higher spatial resolution (equation 2).

Migration tends to collapse the Fresnel zone to approximately the dominant wavelength (equation 1) [8]. Therefore, we anticipate that migration will not resolve the horizontal limits of some of the nonreflecting segments along the deeper reflectors in Figure 11.1-4. Tables 11-1 and 11-2 can be used to estimate the potential resolution improvement that may result from migration. Unless three-dimensional (3-D) migration (3-D poststack migration) is performed, the actual resolution will be less than that indicated. Two-dimensional (2-D) migration only shortens the Fresnel zone in the direction parallel to the line orientation. Resolution in the perpendicular direction is not affected by 2-D migration.

Figure 11.1-5 indicates how vertical and lateral resolution problems are inter-related. We want to determine the edge of the pinchout. The basis of the pinchout model is a wedge of material represented at a given midpoint location by a two-term reflectivity sequence, one term associated with the top and one with the bottom of the wedge. The true thickness of the wedge at various locations is indicated on top of Figure 11.1-5a. The velocity within the wedge is 2500 m/s.

We first consider the reflectivity sequence with two spikes of equal amplitude and identical polarity. The vertical-incidence seismic response (Figure 11.1-5a) is obtained by convolving the sequences with a zero-phase wavelet with a 20-Hz dominant frequency. (The zero-phase response simplifies event tracking from the top and bottom of the wedge.) Based on this response, the edge of the wedge can be inferred as left of location B, where the waveform reduces to a single wavelet (Figure 11.1-5a). From the resolution threshold criterion, the smallest thickness that can be resolved is (2500 m/s)/(4Hz) = 31.25 m. Figures 11.1-5a, b, and c show the same pinchout modeled using three different zero-phase wavelets with increasing dominant frequency (20, 30, and 40 Hz). Separation between the true location of pinchout A and the position of the minimum resolvable wedge thickness B decreases with increasing wavelet bandwidth.

While the resolution threshold criterion allows us to say only that the thickness of the wedge is less than 31.25 m left of location B, an amplitude-based criterion can provide a substantially more accurate location for the edge of the wedge. Again, refer to Figure 11.1-5a and observe the sudden change in amplitude at location A where the true edge of the pinchout is located. Hence, the edge still can be reliably detected, even though it may not be resolved, provided the signal-to-noise (signal-to-noise) ratio is favorable. Assuming that the relative size of the top and bottom reflection coefficients is known, amplitudes also can be used to estimate the wedge thickness between locations B and A.

Figures 11.1-5a, b, and c show an apparent lateral variation in layer thickness. To see the difference between the true thickness and the apparent thickness (peak-to-peak time), refer to Figure 11.1-5d. This figure shows the data of Figure 11.1-5b with the actual geometry of the wedge superimposed on the seismic response. Since the composite wavelet has only one positive peak, the apparent thickness between locations A and B is nearly zero. At location B, the composite wavelet has a flat top. Immediately to the right of location B, the flat top disappears and the composite wavelet splits. The flat-top character can be identified as the limit of vertical resolution [9]. A short distance to the right of the point at which the composite wavelet first splits into two peaks, the apparent thickness becomes equal to the true thickness. This thickness is called the tuning thickness and is equal to peak-to-trough separation (one half the dominant period) of the convolving wavelet [10]. Beyond the point of tuning thickness, note the apparent thickening of the layer between locations B and C. To the right of location C, the apparent and true thicknesses become equal.

Figure 11.1-5  (a) The result of convolving a zero-phase wavelet of 20-Hz dominant frequency with a wedge reflectivity model. The reflection coefficients associated with the top and bottom of the wedge are of equal amplitude and identical polarity. The true edge of the wedge is beneath location A and the true thickness of the wedge is indicated by the numbers on top; (b) same as (a) except the dominant frequency of the wavelet is 30 Hz; (c) same as (a) except the dominant frequency of the wavelet is 40 Hz; (d) same as (b) with the actual geometry of the wedge superimposed on the seismic response; (e) same as (b) except the reflection coefficients from the top and bottom of the wedge have opposite polarity; (f) same as (e) with the actual geometry of the wedge superimposed on the seismic response.

Besides apparent thickness, the maximum absolute amplitude of the composite wavelet along the pinchout also changes [10]. To the left of location A in Figure 11.1-5b, note the single isolated zero-phase wavelet. Immediately to the right of location A, the response of the two closely spaced spikes with identical polarity results in the maximum absolute amplitude. This amplitude gradually decreases to a minimum exactly at the tuning thickness. It then increases and reaches the amplitude value of the original single wavelet to the right of location C.

Maximum amplitude and apparent thickness change in reverse when the reflectivity model consists of reflection coefficients with equal amplitude and opposite polarity (Figure 11.1-5e). The composite waveform resulting from this reflectivity model is discussed by Widess [11]. Two spikes of opposite polarity with a small separation between them act as a derivative operator. When applied to a zero-phase wavelet, this operator causes a 90-degree phase shift. This phase shift can be seen in Figure 11.1-5e on the wavelets between locations A and B. Widess [11] observed that the composite wavelet within this zone basically retains its shape while its amplitude changes.

Figure 11.2-1  A moveout-corrected CMP gather with a reflection event at 1.25 s that exhibits amplitude variations with offset. (Courtesy Western Geophysical.)

Figure 11.1-5f shows data of Figure 11.1-5e with the actual geometry of the wedge superimposed on the seismic response. Note that the wedge appears thicker than it actually is between locations A and B. Also note the apparent thinning of the layer between locations B and C. Beyond location C, the apparent and true thicknesses become equal. Immediately to the right of location A in Figure 11.1-5e, the response of the two closely spaced spikes with opposite polarity results in the cancellation of the amplitudes. The largest absolute amplitude of the composite wavelet gradually increases to a maximum immediately to the right of location B. It gradually decreases and reaches the amplitude value of the original single wavelet to the right of location C.

From the above discussion, we see that peak-to-peak time measurements and amplitude information can aid in detecting pinchouts that may otherwise be unresolvable. If the size of the reflection coefficients were known, then the amplitudes could be used to map the thickness below the resolution limit.

Nevertheless, the reliability of the analysis depends to some extent on the signal-to-noise ratio. Finally, the deceptive changes in amplitude and apparent thickness must be noted during the mapping of the top and bottom of the pinchout.

11.3 Analysis of amplitude variation with offset

In Appendix L, we review the theory of seismic wave propagation in an elastic continuum. The earth’s upper crust that is of interest in seismic prospecting, however, is made up of rock layers with different elastic moduli. When seismic waves travel down in the earth and encounter layer boundaries with velocity and density contrast, the energy of the incident wave is partitioned at each boundary. Specifically, part of the incident energy associated with a compressional source is mode-converted to a shear wave; then, both the compressional- and shear-wave energy are partly reflected from and partly transmitted through each of these layer boundaries.

The fraction of the incident energy that is reflected depends upon the angle of incidence. Analysis of reflection amplitudes as a function of incidence angle can sometimes be used to detect lateral changes in elastic properties of reservoir rocks, such as change in Poisson’s ratio. This may then suggest a change in the ratio of P-wave velocity to S-wave velocity, which in turn may imply a change in fluid saturation within the reservoir rocks.

By way of CMP recording geometry, reflection amplitudes are not measured as a function of angle; instead, they are measured as a function of source-receiver offset. Nevertheless, a range of incidence angles is spanned by a range of offsets. Amplitude-versus-offset analysis, therefore, provides the information on amplitude-versus-angle.

Figure 11.2-1 shows a moveout-corrected CMP gather with a strong reflection at 1.25 s. Note the amplitude variations with offset — in this case, amplitudes increasing with offset. By picking the peak amplitudes and plotting them against offset, an amplitude variation with offset (AVO) curve is derived for a target horizon at each CMP location. Then, an AVO pattern may emerge, which can then be used to infer reservoir parameters.

The pattern with which amplitudes vary with offset depends on the combination of reservoir rock and fluid properties. Detection of a pattern is dictated primarily by the signal-to-noise ratio and the range of incidence angle that is spanned by the offset range of the CMP gather for a target horizon. The shallower the reflector, the wider the range of incidence angle; hence, AVO indicators are best determined for shallow targets.

The following discussion on reflection and refraction of seismic waves is based on flat layer boundaries. Reflection amplitudes also depend upon the dip of the reflecting boundary and its curvature. We can remove the dip and curvature effects by performing prestack time migration. The resulting CMP gathers are associated with reflectors in their migrated positions and reflection amplitudes can then be associated with a locally flat earth model.

Reflection and refraction

For simplicity, consider a monochromatic compressional plane wave that impinges at normal incidence upon a flat layer boundary at z = 0. The incident energy is partitioned between a reflected and transmitted compressional plane wave. For this special case, there is only one stress component, Pzz, and one displacement component, w, which is only a function of z. The equation of wave motion (L-29c) for this special case takes the form


(3a)

where ρ is density of the medium, and λ and μ are Lamé’s constants (equations L-19a,L-19b) associated with an isotropic solid. They are directly related to the compressional-wave velocity α by equation (L-35) which is rewritten below as


(3b)

A solution to equation (3a) can be written as


(4a)

where w0 is the wave function for the incident compressional wave, A0 is its amplitude, ω is the angular frequency, α1 is the compressional-wave velocity in the upper layer (equation 3b), and z-axis is downward positive. Similar wave functions can be written for the reflected plane wave:


(4b)

and for the transmitted plane wave:


(4c)

where α2 is the compressional-wave velocity in the lower layer.

Given the incident wave amplitude A0, we want to compute the reflected and refracted wave amplitudes, A1 and A2, respectively. Equations (4a,4b,4c) must be accompanied with boundary conditions at z = 0 that satisfy the continuity of displacement and stress. The continuity of displacement condition, w0 + w1 = w2 at the interface z = 0 gives the relation


(5)

The stress component Pzz can be specialized from Hooke’s law (equation L-18c) for the present case of a normal-incident compressional plane wave


(6)

The continuity of stress condition at the interface z = 0 gives the relation


(7)

Now, differentiate the wave functions of equations (4a,4b,4c) with respect to z, substitute into equation (7) and set z = 0 to obtain the following expression:


(8)

where ρ1 and ρ2 are the densities of the upper and lower layers, respectively.

We now have two equations, (5) and (8), and two unknowns, A1 and A2. By combining equations (5) and (8), we can derive the ratio of the reflected wave amplitude to the incident wave amplitude, which is called the reflection coefficient c, associated with the layer boundary as


(9)

Define the product of density and velocity as the seismic impedance, I = ρα. If there is a difference between the seismic impedances of the two layers, then a reflection occurs at the interface. If the upper layer has a higher impedance than the lower layer, the reflection coefficient becomes negative causing a phase reversal on the reflected waveform.

Figure 11.2-2  A long, deep-water seismic section that shows internal reflections within the water layer caused by density contrast associated with temperature and salinity changes in the water layer. (Data courtesy IFP.)

An impedance contrast at a layer boundary often is largely caused by velocity contrast. Nevertheless, conditions exist for which density contrast can be significant in giving rise to reflections. Figure 11.2-2 illustrates one such case. The internal reflections within the water layer occur because of changes in water temperature and salinity that causes variations in water density.

A typical reflection coefficient for a strong reflector is about 0.2. The reflection coefficient associated with a hard water bottom is about 0.3. Note that it is the impedance contrast, and not the density or velocity contrasts, that gives rise to reflection energy.

In the foregoing discussion, a normal-incident compressional plane wave was considered. If the same compressional plane wave was incident at an oblique angle to the interface, then derivation of the reflected and transmitted wave functions gets complicated, and we find that the reflection coefficient changes with angle of incidence. Moreover, at non-normal incidence, the incident compressional-wave energy is partitioned at the interface into four components: reflected compressional, reflected shear, transmitted compressional, and transmitted shear waves.

For simplicity, consider a 2-D compressional plane wave that impinges on a layer boundary with an angle of incidence φ0 as depicted in Figure 11.2-3a. The incident wavefront is denoted by AC. Point A at the layer boundary acts as a Huygens’ secondary source and generates its own spherical wavefronts associated with compressional and shear waves propagating in both upper and lower media with the corresponding velocities. In Figure 11.2-3a, only the reflected compressional wavefront and raypath are shown. By the time the incident wavefront at C reaches the reflecting interface at B, the spherical wavefront associated with the reflected compressional wave reaches D, so that AD = CB and the tangential line DB becomes the wavefront for the reflected compressional wave. Since the incident and the reflected compressional waves travel with the same velocity, the angle of incidence φ0 is equal to the angle of reflection φ1.

Figure 11.2-3  Reflection and refraction of an incident P-wave at a layer boundary. Medium parameters: ρ is density, α is P-wave velocity, β is S-wave velocity. (a) Reflected P-wave; (b) reflected S-wave; (c) refracted P-wave; (d) refracted S-wave; (e) raypaths associated with the incident P wave, and reflected and refracted P- and S-waves. The radius of the circular wavefront associated with Huygens’ secondary source at A on the layer boundary is CB for the reflected wave, (β1/α1) for the reflected S-wave, (α2/α1)CB for the refracted P-wave, and (β2/α1) for the refracted S-wave. The relationship between the angles in (e) is given by Snell’s law (equation 10).

Now consider the reflected shear wave as shown in Figure 11.2-3b. By the time the incident wavefront at C reaches the reflecting interface at B, the spherical wavefront associated with the reflected shear wave reaches D, so that AD = (β1/α1)CB, and the tangential line DB becomes the wavefront for the reflected shear wave. The angle for the reflected shear wave ψ1 is no longer equal to the angle of incidence φ0.

Huygens’ principle can also be applied to describe the refracted wave at the interface. Refer to Figures 11.2-3c and 11.2-3d and note that the compressional incident plane wave at A acts as a Huygens’ secondary source and generates its own compressional wavefront that travels into the lower medium with velocity α2 and shear wavefront that travels into the lower medium with velocity β2. By the time the incident wavefront at C in Figure 11.2-3c reaches the layer boundary at B, the spherical wavefront associated with the refracted compressional wave reaches D, so that AD = (α2/α1)CB, and the tangential line DB becomes the wavefront for the refracted compressional wave. Similarly, by the time the incident wavefront at C in Figure 11.2-3d reaches the layer boundary at B, the spherical wavefront associated with the refracted shear wave reaches D, so that AD = (β2/α1)CB, and the tangential line DB becomes the wavefront for the refracted shear wave.

From the geometry of reflected and refracted ray-paths shown in Figure 11.2-3e, Snell’s law of refraction can be deducted as


(10)

Note that for all four cases in Figure 11.2-3, the horizontal distance AB at the interface is the same for the incident and reflected or the refracted wave. If AC is set to the wavenumber along the propagation path of the incident wave, then AB is the horizontal wavenumber which is invariant as a result of reflection or refraction. Actually, Snell’s law given by equation (10) is a direct consequence of this physical observation.

If the compressional velocity α1 of the upper layer is less than the compressional velocity α2 of the lower layer, then there exists an angle of incidence such that no refracted compressional energy is transmitted into the lower layer. Instead, the refracted energy travels along the interface and is refracted back to the upper layer with an angle equal to the angle of incidence. This angle is called the critical angle of incidence for compressional waves and is given by


(11a)

The critically refracted wave is often called the head wave and is the basis for refraction statics (refraction statics corrections).

If the compressional velocity α1 of the upper layer is less than the shear velocity β2 of the lower layer, then there exists an angle of incidence such that no refracted shear energy is transmitted into the lower layer. Again, the refracted energy travels along the interface and is refracted back to the upper layer with an angle which is called the critical angle of P-to-S conversion given by


(11b)

For the general case of non-normal incidence, boundary conditions at the interface involve not only principal stress and strain components but also the shear stress and strain components. Again, by using the requirement that the stress and displacement must be continuous at the interface, a set of equations can be derived to compute the amplitudes of the reflected and refracted P- and S-wave components associated with an incident compressional source (Section L.5):


(12a)


(12b)


(12c)

and


(12d)

These are the Zoeppritz equations which can be solved for the four unknowns, the reflected compressional-wave amplitude A1, the reflected shear-wave amplitude B1, the refracted compressional-wave amplitude A2, and the refracted shear-wave amplitude B2. Equations (12a,12b,12c,12d) have been normalized by the incident-wave amplitude A0 = 1 (Section L.5).

Figure 11.2-4  Framework for derivation of the Zoeppritz equations. For details see Section L.5.

From Snell’s law (equation 10), given the incident angle for the compressional wave and specifying the compressional- and shear-wave velocities, the angles of reflection and refraction can be computed. Substitution into equation (12) yields the required wave amplitudes. These wave amplitudes, of course, depend on the angle of incidence (Figure 11.2-3e).

Figure 11.2-4 outlines the framework for deriving the Zoeppritz equations based on an earth model that comprises two layers separated by a horizontal interface. Details are left to Section L.5. Starting with the equations of motion and Hooke’s law, derive the wave equation for elastic waves in isotropic media in which elastic properties are invariant in any spatial direction at any given location. Then, use the equations of continuity which state that the vertical and tangential stress and stress components coincide at layer boundary, plane-wave solutions to the wave equation and Snell’s law that relates propagation angles to wave velocities to derive the equations for computing the amplitudes of the reflected and transmitted P- and S-waves.

Refer to Figure 11.2-5 for a specific case of the partition of energy of an incident compressional-wave amplitude into four components. Note the significant changes at critical angles of refraction for the compressional and shear wave. It is important to keep in mind that the shape of these curves varies greatly with different situations of medium parameters. Also, note that at normal incidence no P-to-S conversion takes place, and equations (12a,12b,12c,12d) reduce to the special case described by equations (5) and (8).

From a practical standpoint, the angle-dependency of reflection coefficients implies that the reflection amplitude associated with a reflecting boundary varies with source-receiver separation as well as the depth of the reflector. For sufficiently deep reflectors and the typical source-receiver separations used in practice, the amplitude for the reflected compressional wave is nearly constant or slowly varying with offset (left of the critical angle on the curve corresponding to the reflected compressional-wave energy in Figure 11.2-5). It is this slowly varying portion of the P-to-P reflection curve that we have to detect from CMP data with limited offset range and in the presence of noise. Note also from Figure 11.2-5 that the largest P-to-S conversion occurs beyond the critical angle, corresponding to large source-receiver separations.

Figure 11.2-5  Partitioning of a unit-amplitude incident P-wave energy into four components — reflected and refracted P- and S-waves [12].).

Reflection amplitude variations with angle of incidence, and therefore with offset, can be modeled using the Zoeppritz equations. For modeling the reflection amplitudes, you need well-log curves for the P-wave velocity, S-wave velocity, and density. Then, using equation (13a) compute the P-to-P reflection amplitudes as a function offset. Figure 11.2-6 shows a modeled CMP gather using well data. A sonic log was first converted to a blocky form to simplify the modeled amplitudes on the gather. Then, using the Gardner relation (equation 34a), a density profile was derived. Also, using a ratio of P-wave to S-wave velocity, a shear-wave velocity profile was generated. The objective in this modeling exercise was to see the effect of a change in Poisson’s ratio at some depth on the amplitude variation with offset. The Poisson’s ratio profile shows a change at approximately 650 ms at which time we observe a marked variation of amplitudes with offset on the CMP gather.

Reflector curvature

Define the ratio CE of the reflection amplitude at normal incidence from a curved boundary to the reflection amplitude from a flat boundary at the same depth [13] as


(13a)

where z is the depth to the reflector, and A is the reflector curvature — negative for synclines and positive for anticlines. Note that for a flat reflector, whatever its depth, CE = 1. But for a curved reflector, the reflection amplitude at normal incidence to a synclinal interface is greater than the case of a flat reflector, and the reflection amplitude at normal incidence to an anticlinal interface is smaller than the case of a flat reflector. The physical basis of equation (13a) is that a synclinal interface focuses the energy associated with the reflecting wave, whereas an anticlinal interface defocuses it.

Figure 11.2-6  An example of modeling of a CMP gather based on Zoeppritz equations. (Courtesy Hampson-Russell.)

The effect of reflector curvature on amplitude variation with offset is quantified as [14]


(13b)

where θ is the angle of incidence. Note that equation (13b) reduces to equation (13a) for the case of normal incidence (θ = 0). Both equations are for the case of a 2-D reflecting interface. Equation (13b) was extended by Bernitsas [15] to the case of a 3-D reflecting interface with curvature in both the inline and crossline directions.

To understand the effect of reflector curvature on amplitude variation with offset, it is convenient to study the ratio CE(θ)/CE(θ = 0) as a function of angle of incidence [14] [16]. Figure 11.2-7 shows the behavior of this ratio for the cases of a syncline and an anticlne with varying degrees of curvature. Note that, for an anticlinal interface, the curvature effect defined by the ratio CE(θ)/CE(θ = 0) decreases with angle of incidence or offset (Figure 11.2-7c). In practice, this means that the reflection amplitudes from an anticlinal interface at far offsets are lower than those from a flat interface. Also, note that the larger the curvature defined by the ratio z/x denoted in Figure 11.2-7, the more prominent the curvature effect. For a synclinal interface with a mild curvature (0 > z/x > −1), the curvature effect on amplitudes increases with angle of incidence, and it decreases with angle of incidence for a tight syncline (z/x < −1). Again, note that the larger the curvature, the more prominent the curvature effect.

Figure 11.2-7  Effect of reflector curvature on angular dependence of reflected wave amplitudes. The incident wave is of compressional type. (a) A synclinal reflector (x < 0), (b) an anticlinal reflector (x > 0), (c) effect of curvature on the angle-dependent reflection amplitudes associated with the reflector as in (a), and (d) the same as in (c) for the reflector as in (b). Adapted from [14].

It is clear from the model studies summarized above that reflection amplitudes are influenced by the reflector curvature. This is true for both normal incidence (equation 13a) and non-normal incidence (equation 13b). Reflection amplitudes must be corrected for the reflector curvature by prestack time migration before prestack amplitude inversion to derive the AVO attributes which are discussed in this section. Similarly, reflection amplitudes must be corrected for the reflector curvature by poststack time migration before poststack amplitude inversion to derive the acoustic impedance attribute which is discussed in the next section.

AVO equations

Consider the two elastic half-space layers in Figure 11.2-3e. The Zoeppritz equations (14) can be solved for the reflected and refracted P- and S-wave amplitudes, A1, B1, A2, and B2. However, our interest in exploration seismology is largely the angle-dependency of the P-to-P reflections given by the coefficient A1. Specifically, we wish to infer or possibly estimate elastic parameters of reservoir rocks from reflection amplitudes and relate these parameters to reservoir fluids.

The exact expression for A1 derived from the solution of the Zoeppritz equations (13) is complicated and not intuitive in terms of its practical use for inferring petrophysical properties of reservoir rocks. The first approximation to the Zoeppritz equation for P-to-P reflection amplitude is given by Bortfeld [17] as


(14)

The arrangement of the two terms on the right-hand side of equation (14) is based on separating the acoustic (the first term) and the elastic (the second term) effects on reflection amplitudes. As such, equation (14) does not explicitly indicate angle- or offset-dependence of reflection amplitudes; therefore, its practical implementation for AVO analysis has not been considered.

Instead, we shall use the approximation provided by Aki and Richards [18] as the starting point for deriving a series of practical AVO equations. Now that we only need to deal with the P-to-P reflection amplitude A1, we shall switch to the conventional notation by replacing A1 with R(θ) as the angle-dependent reflection amplitude for AVO analysis.

By assuming that changes in elastic properties of rocks across the layer boundary are small and propagation angles are within the subscritical range, the exact expression for R(θ) given by the Zoeppritz equation can be approximated by [18]


(15)

where α = (α1 + α2)/2, average P-wave velocity and Δα = (α2α1), β = (β1 + β2)/2, average S-wave velocity and Δβ = β2β1, ρ = (ρ1 + ρ2)/2, average density and Δρ = ρ2ρ1, and θ = (φ1 + φ2)/2, average of the incidence and transmission angles for the P-wave (Figure 11.0-2e).

Figure 11.2-8 shows the angle-dependent reflection amplitude associated with an interface with contrast in P- and S-wave velocities and densities based on the exact Zoeppritz equation, the Bortfeld approximation described by equation (14), and the Aki-Richards approximation described by equation (15) [19]. Note that these approximations closely follow the exact solution within the range of angles of incidence that are achievable by the recording of seismic data used in exploration. For very shallow reflectors, however, the approximate solutions would deviate from the exact solution significantly at very wide angles of incidence.

Note that the Aki-Richards approximation to the angle-dependent reflection amplitude R(θ) given by equation (15) has three parts in terms of Δα/α which describes the fractional change in P-wave velocity across the layer boundary and hence may be referred to as the P-wave reflectivity, Δβ/β which describes the fractional change in the S-wave velocity across the layer boundary and hence may be referred to as the S-wave reflectivity, and Δρ/ρ which describes the fractional change in density across the layer boundary.

In practice, we do not observe the separate effects of P-wave reflectivity Δα/α, S-wave reflectivity Δβ/β and fractional change in density Δρ/ρ on the reflection amplitudes R(θ). Instead, we observe changes in reflection amplitudes as a function of angle of incidence. In fact, it is the elastic parameters such as the P-wave reflectivity Δα/α, S-wave reflectivity Δβ/β and fractional change in density Δρ/ρ that we wish to estimate from the observed angle-dependent reflection amplitudes. To use the Aki-Richards equation (15) in the inversion of reflection amplitudes for these elastic parameters, we first need to recast it in successive ranges of angle of incidence. This change of philosophy in arranging the terms in the Aki-Richards equation (15) was first introduced by Shuey [20] and led to practical developments in AVO analysis. The new arrangement in terms of successive ranges of angle of incidence is given by


(16)
Figure 11.2-8  P-to-P reflection amplitude as a function of angle of incidence computed by using the exact Zoeppritz equation, and Bortfeld (equation 14) and Aki-Richards (equation 15) approximations [19].

Another practical matter of concern is that the Aki-Richards equation (15) or any of its modifications that we shall derive in this section describe the modeled reflection amplitudes as a function of angle of incidence. However, the observed reflection amplitudes are available from CMP data as a function of offset. A need then arises either to transform the model equation for the reflection amplitudes from angle to offset coordinates [21] or to actually transform the CMP data from offset to angle coordinates. While the first approach is theoretically appealing, the practical schemes are based on the latter approach. We have already discussed such a transformation in the radon transform — the Radon transform using the linear moveout equation or its robust variation in the form of slant stacking. Figure 11.2-9 shows Zoeppritz amplitude curves as a function of angle of incidence and offset [21].

Based on the theoretical conjecture made earlier by Koefoed [22] that the elastic property that is most directly related to angular dependence of reflection coefficient R(θ) is Poisson’s ratio σ, Shuey [20] introduced a variable transformation from S-wave velocity β to σ. The relationship between the two variables is given by equation (L-49) which we rewrite below as


(17a)

to perform the necessary differentiation


(17b)

The compressional-wave velocity α is given by equation (3b) and the shear-wave velocity β is given by equation (L-47) which is rewritten below as


(17c)

where μ is Lamé’s constant.

We also define the P-wave reflection amplitude RP at normal incidence as


(18)

Substitute equations (17a), (17b), and (18) into the Aki-Richards equation (16) and perform some algebraic simplification to obtain


(19)

Define a new term H


(20)

and by way of equation (18) note that


(21a)

Next combine equations (20) and (21a) to obtain


(21b)

Finally, substitute equations (21a) and (21b) into the second term on the right-hand side of equation (19) to obtain


(22)

where


(23)

Equation (22) is known as Shuey’s three-term AVO equation. The first term RP is the reflection amplitude at normal incidence. At intermediate angles (0 < θ < 30 degrees), the third term may be dropped, thus leading to a two-term approximation


(24)

where


(25)
Figure 11.2-9  Reflection and refraction amplitude curves as a function of offset and angle of incidence. P1, S1: reflected P- and S-waves, P2, S2: transmitted P- and S-waves. [21].

Equation (24) is known as Shuey’s two-term AVO equation. In practice, amplitudes picked along a moveout-corrected event on a CMP gather plotted against sin2 θ can be fitted to a straight line. The slope of the line gives the AVO gradient attribute and the ordinate at zero angle gives the AVO intercept attribute. The AVO gradient given by equation (25) is directly related to change in Poisson’s ratio Δσ, which in turn, is directly related to fluid saturation in reservoir rocks. The AVO intercept attribute represents the reflectivity RP at normal incidence. Therefore, the AVO intercept attribute, in lieu of conventional stack, can be used as input to derive the acoustic impedance attribute (acoustic impedance estimation), which is indirectly related to porosity in reservoir rocks.

Shown in Figure 11.2-10a is a portion of a section derived from 2-D prestack time migration. The objective is to identify fluid-saturated reservoir zones at the apex and the flanks of the structural closure. This image section is derived from the stacking of common-reflection-point (CRP) gathers associated with the prestack time-migrated data. The CRP gather in Figure 11.2-10b shows three events with amplitude variations with offset which are plotted in Figure 11.2-10c. By using Shuey’s equation (24), the AVO gradient and intercept sections are computed from the CRP gathers as shown in Figures 11.2-11 and 11.2-12, respectively. Note that the gradient section exhibits a group of AVO anomalies in the vicinity of the structural apex, possibly indicating fluid-saturated reservoir rocks.

At large angles of incidence beyond 30 degrees, the third term in equation (22) gradually becomes dominant. Note that this term is related directly to fractional change in P-wave velocity, Δα/α. So, not only the reflection traveltimes at far offsets (normal moveout) corresponding to large angles of incidence, but also the reflection amplitudes at large angles of incidence make the biggest contribution to the resolution needed to estimate the changes in P-wave velocities.

The two-term equation (24) can be specialized for a specific value of Poisson’s ratio, σ = 1/3 and H0 = −1 so that equation (25) takes the form


(26)

which can be solved for the change in Poisson’s ratio across a layer boundary


(27)

This is the AVO attribute equation for estimating changes in Poisson’s ratio [23]. Actually, as described by equation (27), this attribute is the scaled sum of the AVO intercept RP and AVO gradient G attributes.

By recasting the first-order approximation to the Zoeppritz equation, Wiggins [24] and Spratt [25] derived a practical expression for S-wave reflectivity. First, drop the term with tan2 θ in equation (16) and rearrange the remaining terms to obtain


(28)

such that, much like the definition for the P-wave reflectivity RP given by equation (18), an expression for the S-wave reflectivity RS


(29)

can be explicitly inserted back into equation (28) to get


(30)

Set β/α = 0.5 to make the last term on the right-hand side vanish and obtain [25]


(31)

This equation is of the form given by equation (24) where


(32)

which can be rewritten explicitly in terms of the shear-wave reflectivity RS


(33)

This is the AVO attribute equation for estimating the shear-wave reflectivity. Given the AVO intercept RP and AVO gradient G attributes, simply take half of the difference between the two attributes to derive the shear-wave reflectivity RS as described by equation (33).

Figure 11.2-13 shows the reflection amplitudes as a function of angle as predicted by equation (31) for three combinations of RP and RS. Note that equation (31) is a good approximation to the P-to-P reflection amplitudes as predicted by the exact Zoeppritz equation up to nearly 30 degrees of angles of incidence.

Return to the Aki-Richards equation (15) and consider the case of N-fold CMP data represented in the domain of angle of incidence. Note that the reflection amplitude R(θ) is a linear combination of three elastic parameters — P-wave reflectivity Δα/α, S-wave reflectivity Δβ/β and fractional change in density Δρ/ρ.

Smith and Gidlow [19] argue in favor of solving for only two of the three parameters by making use of the empirical relation between density ρ and P-wave velocity α [1]:


(34a)

where k is a scalar. This relation holds for most water-saturated rocks. Differentiate to get


(34b)

Now substitute equation (34b) into the original Aki-Richards equation (15) and combine the terms with Δα/α


(35)

Simplify the algebra to obtain the desired two-parameter model equation for prestack amplitude inversion (Section L.6)


(36)

Redefine the coefficients ai and bi and write the discrete form of equation (36)


(37)

where


(38a)

and


(38b)

Consider the case of N-fold CMP data represented in the domain of angle of incidence. We have, for each CMP location and for a specific reflection event associated with a layer boundary, N equations of the form as equation (37) and two unknowns — Δα/α and Δβ/β. We have more equations than unknowns; hence, we encounter yet another example of a generalized linear inversion (GLI) problem (Section J.1). The objective is to determine the two parameters such that the difference between the modeled reflection amplitudes Ri represented by equation (37) and the observed reflection amplitudes Xi is minimum in the least-squares sense [19].

The GLI solution in matrix form is given by (Section L.6)


(39)

This matrix equation is solved for the two parameters — P-wave reflectivity Δα/α and S-wave reflectivity Δβ/β. Note from the definitions of the coefficients ai and bi given by equations (38a,38b) that you have to choose a value for the ratio β/α to compute the two reflectivities.

Shown in Figure 11.2-14 are the P- and S-wave reflectivity sections that were derived from prestack amplitude inversion applied to the CRP gathers associated with the data shown in Figure 11.2-10a. Note the AVO anomalies along the flanks of the structural closure. It is apparent that some of the AVO anomalies in the P-wave reflectivity section stand out more distinctively as compared with those in the S-wave reflectivity section. This may be related to the fact that changing from brine to gas causes a small change in S-wave reflectivity, while it causes a significant change in P-wave reflectivity. Such inference is supported by the laboratory measurements of P- and S-wave reflectivities at a gas-brine interface in sandstones [25]. Figure 11.2-15 shows P-wave reflectivities for a group of reservoir sandstones in which the pore fluid changes from brine to gas. For rocks with low induration (weak rocks), there is a large change in P-wave reflectivity, while for rocks with high induration the change becomes less significant (Figure 11.2-15a). The change in S-wave reflectivity is negligible for rocks with low and high induration (Figure 11.2-15b). A demonstrative field data example of P- and S-wave reflectivity contrasts that arise from a change in lithology and fluid saturation is shown in Figure 11.2-16.

Figure 11.2-18  Crossplot of P- to S-wave velocity measured down a borehole [27].

Refer to equation (17b) and note that the difference between the P-wave and S-wave reflectivities is related to change in Poisson’s ratio Δσ — a direct hydrocarbon indicator. In fact, Smith and Gidlow [19] have coined the term pseudo-Poisson reflectivity to describe the difference between the P-wave and S-wave reflectivities


(40)

Note from equation (17b) that the pseudo-Poisson reflectivity is not exactly the same as what may be referred to as the proper Poisson reflectivity Δσ/σ. If we define the ratio


(41a)

by differentiation, we can derive the difference relation given by equation (40) and thus show that


(41b)

Castagna [26] defined a straight line in the plane of S-wave velocity versus P-wave velocity as shown in Figure 11.2-17. This is called the mudrock line and is represented by the equation


(42)

where the scalar coefficients c0 and c1 are empirically determined for various types of rocks. Suggested values for these scalars are c0 = 1360 and c1 = 1.16 for water-saturated clastics [26]. Gas-bearing sandstones lie above the mudrock line, and carbonates lie below the mudrock line as shown on a crossplot of P- and S-wave velocities sketched in Figure 11.2-17. The crossplot shown in Figure 11.2-18 is based on log measurements at a gas-producing well [27]. Note the separation of gas-sandstone cluster from the water-sandstone and shale clusters in the manner as sketched in Figure 11.2-17.

A way to quantify the prospectivity of the reservoir rock of interest is by defining a fluid factor attribute that indicates the position of the rock property with respect to the mudrock line. First, apply differentiation to both sides of equation (42) and note that


(43)

Then, define the fluid factor ΔF


(44)

If ΔF is close to zero, it means that you have water-saturated rock. If you have a gas-saturated sandstone, ΔF will be negative at the top and positive at the base of the reservoir unit [19].

Figure 11.2-19 shows the pseudo-Poisson and fluid-factor sections that were derived from prestack amplitude inversion applied to the CRP gathers associated with the data shown in Figure 11.2-10a. Note the distinctive AVO anomalies along the flanks of the structural closure.

The two parameters Δα/α and Δβ/β, estimated by using the least-squares solution given by equation (39), represent fractional changes in P- and S-wave velocities. As such, they are related to P- and S-wave reflectivities, ΔIP/IP and ΔIS/IS, respectively, where IP and IS are the P- and S-wave impedances given by


(45a)

and


(45b)

From Section L.6, the P- and S-wave reflectivities are given by


(46a)

and


(46b)

Assuming that the density obeys Gardner’s relation given by equation (34a), the P-wave reflectivity ΔIP/IP is simply the fractional change of the P-wave velocity Δα/α scaled by a constant, whereas the S-wave reflectivity ΔIS/IS is a linear combination of the fractional changes of the P- and S-wave velocities, Δα/α and Δβ/β, respectively.

Figure 11.2-22  Framework for derivation of the various AVO equations in analysis of amplitude variation with offset.

A direct estimation of the P- and S-wave reflectivities given by equations (L-46a,L-46b) can be made by using an alternative formulation of prestack amplitude inversion which is based on transforming the Aki-Richards equation (15) to the new variables ΔIP/IP and ΔIS/IS [28]. Solve equation (46a) for Δα/α and equation (46b) for Δβ/β, and substitute into the Aki-Richards equation (15). Following the algebraic simplification, we obtain (Section L.6)


(47)

Goodway [28] have implemented a specific form of equation (47) to derive the AVO attributes ΔIP/IP and ΔIS/IS. For a specific value of α/β = 2 and small angles of incidence for which tan θ ≈ sin θ, the third term in equation (47) vanishes. Compare equations (18) and (46a), and equations (29) and (46b), and note that equation (47) with the remaining terms takes the form


(48)

The resulting equation then is solved for the P- and S-wave reflectivities, ΔIP/IP = 2RP and ΔIS/IS = 2RS, respectively. Again, consider the case of N-fold CMP data represented in the domain of angle of incidence. In discrete form, equation (48) can be rewritten as


(49)

where i is the trace index and the coefficients ai and bi are given by


(50a)

and


(50b)

The reflection amplitude Ri in equation (49) is a linear combination of the two parameters, RP and RS. As for the Smith-Gidlow equation (37), the two parameters, RP and RS, can be estimated for a specific event from CMP data using the least-squares solution given by


(51)

Following the estimation of the P-wave reflectivity RP and the S-wave reflectivity RS using the least-squares solution given by equation (51), the P-wave impedance IP and the S-wave impedance IS can be computed by integration.

By using the impedance attributes IP and IS, Goodway [28] compute two additional AVO attributes in terms of Lamé’s constants scaled by density — λρ and μρ. Substitute equation (3b) into (45a) to get the relation


(52a)

and substitute equation (17c) into equation (45b) to get the relation


(52b)

Note from equations (52a,52b) that


(52c)

Figure 11.2-20 shows the μρ and λρ AVO attribute sections which were derived from prestack amplitude inversion applied to the CRP gathers associated with the data shown in Figure 11.2-10a. As in the case of the pseudo-Poisson and fluid-factor AVO attribute sections shown in Figure 11.2-19, note the distinctive AVO anomalies along the flanks of the structural closure. Goodway [28] convincingly demonstrates that separation of gas-bearing sands from tight sands and carbonates is much better with the crossplot of the Lamé attributes (μρ, λρ) in contrast with the crossplot of the P- and S-wave reflectivities or impedances (Figure 11.2-21).

Figure 11.2-22 outlines a summary of the AVO equations that are described in this section based on the various approximations to the Aki-Richards equation (15). Start with Shuey’s arrangement of the terms in the Aki-Richards equation, which itself is an approximation to the exact expression for the P-to-P reflection amplitudes given by the Zoeppritz equations. Then, apply a transformation from S-wave velocity to Poisson’s ratio, and drop the third term to get a simple expression for the P-to-P reflection amplitude that is a linear function of sin2 θ. This linear approximation yields the AVO intercept and gradient attributes.

Alternatively, you may drop the third term in the original Shuey arrangement given by equation (16) and apply Wiggins’ rearrangement of the terms. Then, assume the specific case of the P-wave velocity to be twice the S-wave velocity to obtain the equation that yields, once again, the AVO intercept and gradient attributes. The Wiggins AVO intercept and gradient attributes, unlike the Shuey counterparts, directly lead to deriving the P-wave and S-wave reflectivity attributes, albeit only under the assumption that the P-wave velocity is twice the S-wave velocity.

Returning to the three-term Aki-Richards equation (15), you may reduce it to the two-term Smith-Gidlow equation (36) by using Gardner’s relation of density and P-wave velocity. Then, for a specified ratio of P- to S-wave velocity, perform inversion of prestack amplitudes associated with target events on CMP or CRP gathers to obtain the P- and S-wave reflectivity attributes. In addition, using these attributes to derive two more AVO attributes — psuedo-Poisson and fluid factor.

Finally, you may recast the Aki-Richards equation (15) in terms of P- and S-wave reflectivities (equations 46a,46b) and assume that α/β = 2 to obtain the two-term Goodway et al. equation (48). Again, perform inversion of prestack amplitudes associated with target events on CMP or CRP gathers to obtain the P- and S-wave impedance attributes. Subsequently, by using the relations between impedances and Lamé’s constants, you can compute two additional AVO attributes — λρ and μρ.

The AVO equations derived in this section are all expressed in terms of angle of incidence (equations 24, 31, and 36). In practice, the mapping of amplitudes associated with a reflection event on a CMP gather from offset to angle of incidence needs to be performed. Assume that the CMP gather is associated with a horizontally layered earth model so that the event moveout on the CMP gather is given by the hyperbolic equation


(53)

where t is the two-way traveltime from the source to the flat reflecting interface back to the receiver, t0 is the two-way zero-offset time, x is the offset, and vrms is the rms velocity down to the reflector. The ray parameter p is given by the stepout dt/dx measured along the hyperbolic moveout trajectory (the slant-stack transform). Apply differentiation to equation (53) to get


(54)

Since the ray parameter p is also expressed as the ratio of sin θ/vint, where vint is the interval velocity above the reflecting interface, it follows that


(55)

By using equation (55), reflection amplitudes on a CMP gather can be transformed from offset to angle domain, and subsequently used in the AVO equations (24), (31), (36), and (48).

Processing sequence for AVO analysis

There are three important aspects of a processing sequence tailored for AVO analysis.

  1. The relative amplitudes of the seismic data must be preserved throughout the analysis so as to recognize amplitude variation with offset. This requirement often leads to an application of a parsimonious sequence of signal processing to avoid distortion of amplitudes by undesirable effects of some processing algorithms.
  2. The processing sequence must retain the broadest possible signal band in the data with a flat spectrum within the passband.
  3. Prestack amplitude inversion to derive the AVO attributes must be applied to common-reflection-point (CRP) gathers (migration velocity analysis), not to common-midpoint (CMP) gathers. This is because all the AVO equations described in this section are based on a locally flat earth model that can be related to CRP raypaths but not to CMP raypaths. Specifically, the CRP gathers are associated with events in their migrated positions, whereas the CMP gathers are associated with events in their unmigrated positions. The AVO equations described in this section are all based on a horizontally layered earth model that can be related to CRP raypaths but not to CMP raypaths. Prestack time migration compensates for the effect of reflector curvature (equation 13b) on amplitude variation with offset so that the reflection amplitudes in each of the resulting CRP gathers can be associated with a locally flat earth model.

We shall consider a marine 2-D seismic data set from the North Sea recorded over a producing gas field. Data specifications for the line are listed in Table 11-3.

The prestack signal processing that complies with the above requirements includes the following steps:

  1. If a reliable source signature is available, perform signature processing (field data examples) to convert the source waveform to its minimum-phase equivalent. In the present study, no source signature was available. Apply a mute to shot records to remove the guided waves associated with the dispersive, normal-mode propagation within the water layer (Section F.1).
  2. Apply t2-scaling to correct for geometric spreading. A gain function for geometric spreading correction that depends on primary velocities will overcorrect the amplitudes of multiples. Figure 11.2-23 shows a recorded shot record and 11.2-24 shows the same record after guided-wave muting geometric spreading correction.
  3. Apply time-invariant spiking deconvolution with a sufficiently long operator length (Figure 11.2-25). For land data, it may be necessary to apply surface-consistent amplitude corrections and surface-consistent deconvolution to account for the near-surface effects on amplitudes (Section B.8).
  4. When it is required, apply time-variant spectral whitening to account for the nonstationarity of the waveform.

Figure 11.2-26 shows the amplitude spectrum of the shot record in Figure 11.2-25 computed after each of the steps in the sequence prescribed above. To begin with, note that much of the energy in the raw shot record is confined to the shallow portion and is mostly associated with the guided waves. The t2-scaling restores the amplitudes in the deeper portion of the record. Spiking deconvolution followed by a wide-bandpass filter gives the desired shape to the spectrum with a flat, broadband character. Figure 11.2-27 shows the autocorrelograms of the shot record in Figure 11.2-25, again, computed after each of the steps in the sequence prescribed above. Note that the autocorrelogram of the raw shot record is dominated by the guided-wave energy at far offsets on the left and the energy associated with short period multiples and reverberations at near offsets on the right. The guided-wave muting followed by t2-scaling exposes the energy associated with multiples that range from a short- to a moderate-period. Note that spiking deconvolution with an operator length of 480 ms (about one-tenth of the total length of the shot record) has removed a significant amount of the multiple energy from the record. Both the amplitude spectrum (Figure 11.2-26c) and the autocorrelogram (Figure 11.2-27c) of the deconvolved shot record indicate that time-variant spectral whitening does not have to be included in the sequence.

Table 11-3. Data specifications for the line used in AVO analysis.
Line Length 23.5 km
Shot Spacing 37.5 m
Receiver Spacing 18.75 m
CMP Spacing 9.375 m
Minimum Offset 150 m
Maximum Offset 3050 m
Number of Receivers 156
Fold of Coverage 78
Sampling interval 2 ms

The prestack data after signal processing are ready now for prestack time migration to generate the CRP gathers (prestack time migration).

  1. Sort the signal-processed shot records to CMP gathers and perform velocity analysis at sparse intervals along the line. To increase the accuracy of prestack amplitude inversion and thus increase the confidence level for the AVO attributes, the fold of coverage and the offset range of the recorded data must not be reduced during the processing for any reason.
  2. Apply NMO correction using the flat-event velocities. In this case, only three velocity functions were used in creating the velocity field for NMO correction.
  3. When required, attenuate those multiples that have not been handled by deconvolution during the signal processing stage described above. Data suitable for AVO analysis are often associated with near-horizontal earth models or low-relief structures. Hence, techniques based on the discrete Radon transform (the radon transform) often are best suited for amplitude-preserving multiple attenuation.
  4. Sort the data to common-offset sections and apply DMO correction.
  5. Perform zero-offset, frequency-wavenumber time-migration on each of the common-offset sections after NMO and DMO corrections using a single, vertically varying velocity function. Again, data suitable for AVO analysis are generally associated with near-horizontal earth models or low-relief structures. Hence, use of a single, vertically varying velocity function and a zero-offset phase-shift algorithm for common-offset migration often is appropriate.
  6. Sort the data to CRP gathers and apply inverse NMO correction using the velocity field from step (b).
  7. Perform velocity analysis at frequent intervals along the line, and create a velocity field associated with the migrated data.
  8. Apply NMO correction using the velocity field from step (g) and, when required, as it was in the present case study, perform residual moveout velocity analysis to remove any residual moveout errors that may be present in the CRP gathers. The residual moveout analysis actually is a good practice whatever the degree of residual moveout errors in the data, since prestack amplitude inversion is extremely sensitive to how flat the events are in the CRP gathers. Figure 11.2-28 shows selected CRP gathers that exhibit mostly flat events. A close-up view of the CRP gather (Figure 11.2-29) at the well location (CRP 1134) shows the reservoir zone highlighted by the rectangle situated between 3.1-3.3 s.
  9. Stack the CRP gathers to obtain the image section from prestack time migration. Multiples that may have remained in the data will exhibit some residual moveout that can be exploited to further attenuate them during stacking.
  10. For completeness of the prestack time migration sequence, unmigrate the CRP stack from step (i) using the same vertically varying velocity function as in step (e) and a phase-shift zero-offset wavefield modeling algorithm (Figure 11.2-30).
  11. Next, apply poststack spiking deconvolution, often using the same parameters as for prestack deconvolution (Figure 11.2-31).
  12. Remigrate the modeled zero-offset section from step (k) using the migration velocity field derived from the CRP gathers in step (g) (Figure 11.2-32). In this instance, any lateral velocity variations implied by the migration velocity field need to be accounted for by the migration algorithm of choice. Poststack amplitude inversion to derive the acoustic impedance attribute (acoustic impedance estimation) may then be applied to the final product from prestack time migration (Figure 11.2-32). Figure 11.2-33 shows the amplitude spectrum of the CRP stack after each step of the poststack processing. The amplitude spectrum of the CRP stack from step (i) shows the attenuation of high frequencies (Figure 11.2-33a), which are restored by deconvolution (Figure 11.2-33b), but are slightly shifted to lower frequencies by migration (Figure 11.2-33c). The autocorrelogram of the CRP stack from step (i) exhibits remnants of the water-bottom multiples (Figure 11.2-34a), which are largely attenuated by poststack deconvolution (Figure 11.2-34b).

Derivation of AVO attributes by prestack amplitude inversion

Following the prestack signal processing and prestack time migration designed for AVO analysis, prestack amplitude inversion is applied to the CRP gathers from prestack time migration. Whatever the basis for the AVO equations (Figure 11.2-22), derivation of the AVO attributes requires knowledge of the angle of incidence. This in turn requires, strictly, an estimated earth model in depth that can be used to trace rays associated with the CRP geometry. Figure 11.2-35 shows a velocity-depth model associated with the CRP stack shown in Figure 11.2-32. Often, the following procedure is adequate for a velocity-depth model estimation for the purpose of AVO analysis (models with horizontal layers):

  1. Interpret a set of time horizons from the CRP stack (Figure 11.2-32).
  2. Intersect the migration velocity field from step (g) of the prestack time migration sequence described above, which is equivalent to the rms velocity field associated with the vertical raypaths, with the time horizons from step (a) and create a set of horizon-consistent rms velocity profiles.
  3. Perform Dix conversion of the rms velocity profiles to derive a set of interval velocity profiles.
  4. Convert the time horizons from step (a) to depth horizons by using the interval velocity profiles from step (c) along image rays.
  5. Combine the depth horizons from step (d) with the interval velocity profiles from step (c) and build a velocity-depth model as shown in Figure 11.2-35.

Before we proceed with prestack amplitude inversion to derive the AVO attributes, we first need to verify the existence of amplitude variation with offset on recorded data by modeling the amplitudes on a CRP gather coincident with a well location. Actually, we model a moveout-corrected CMP gather at the well location using a locally flat earth model, which is consistent with the earth model associated with the CRP gather. Figure 11.2-36 shows a sonic and a density log at a well location that coincides with the CMP location 1134 (Figure 11.2-32). By using the log curves and the S-wave velocity calculated from the mudrock relation (equation 42), compute the synthetic CMP gather associated with a horizontally layered earth model at the well location. Figure 11.2-37 shows a portion of the CRP gather from prestack time migration as in Figure 11.2-28 and the synthetic CMP gather computed by using the Zoeppritz equation for P-to-P reflection amplitudes associated with a horizontally layered earth model (Section L.5). Shown also in this figure is the angle of incidence within the same time window as for the CRP gather derived from the velocity-depth model in Figure 11.2-35. It is evident that given the offset range of 150-3050 m, the maximum angle of incidence at the target zone (3.1-3.3 s) is nearly 30 degrees. This means that the AVO attributes based on the Shuey approximation (equation 24) can be considered accurate.

A close-up view of the real CRP gather and the synthetic gather of Figure 11.2-37 at the reservoir zone is shown in Figure 11.2-38. Note that the amplitude variations in the two gathers are fairly consistent. In fact, examine in Figure 11.2-39 the amplitude variation with offset at three time levels that coincide with top-gas, top-oil, and top-water boundaries. We may draw two conclusions from the comparison of the real and modeled amplitudes as a function of offset shown in Figure 11.2-39. First, note that the Zoeppritz-modeled (red) and the best-fit actual (dotted green) amplitude variation with offset are in very good agreement — it appears that if you scale one AVO curve by a constant, you will get the other AVO curve. Second, the amplitudes do indeed vary with offset — there is a decrease in amplitude with increasing offset at the top-gas and top-water boundaries, and there is an increase in amplitudes with increasing offset at the top-oil boundary.

The top-gas AVO curve shown in Figure 11.2-39 extracted from the real CRP gather is consistent with the observations made by Rutherford and Williams [29] who classified the AVO anomalies associated with a shale formation over a gas-bearing sand formation. A Class 1 AVO anomaly indicates a large positive P-to-P normal-incidence amplitude (large positive AVO intercept attribute value) and a marked decrease in amplitudes with increasing offset, with a possible phase change at large offsets. According to the Rutherford-Williams classification, the top-gas AVO curve in Figure 11.2-39 can be labeled as Class 1. A Class 2 AVO anomaly indicates a small positive P-to-P normal-incidence amplitude (small positive AVO intercept attribute value) and a marked decrease in amplitudes with increasing offset, with a possible phase change at small-to-moderate offsets. Finally, a Class 3 AVO anomaly indicates a large negative P-to-P normal-incidence amplitude (large negative AVO intercept attribute value) and a marked decrease in amplitudes with increasing offset, typically manifested as bright-spot anomalies.

Another AVO indicator is the difference between a near-angle stack (say, up to 15 degrees of angle of incidence as shown in Figure 11.2-37c) and a wide-angle stack (say, between 15-30 degrees of angles of incidence as shown in Figure 11.2-37c). If there were no amplitude variation with offset, there would not be a difference between the two sections shown in Figure 11.2-40, except for noise and residual moveout errors.

To calibrate the sections in Figure 11.2-40 to well data, compute the reflectivity series using the two log curves in Figure 11.2-36. Then, by using a zero-phase band-limited wavelet with a passband that corresponds to the signal band associated with the prestack time-migrated data (Figure 11.2-32), compute the zero-offset synthetic seismogram that contains primaries only. Insert the synthetic seismogram into the CRP stack at the well location as shown in Figure 11.2-40. Apply constant time shift and constant phase shift until a good match at the reservoir level is attained between the real and the synthetic data. In the present case, a constant time shift only, without a constant phase shift, was adequate to achieve a good match at the reservoir level.

Interpretation of AVO attributes

We are now ready to perform the actual prestack amplitude inversion to derive the AVO attributes. Figure 11.2-41 shows close-up portions of the AVO attribute sections at the vicinity of the well location. In each of the AVO attribute sections, the red represents the high and the green represents the low attributes values. The reservoir zone is between the two red horizons in the CRP stack section shown in Figure 11.2-40. For reference, use Figure 11.2-22 to identify the equations associated with the attributes shown in Figure 11.2-41.

Compare the intercept section (Figure 11.2-41) and the CRP stack section (Figure 11.2-40), and note that the former better matches the synthetic attribute derived from the well data. In theory, the intercept attribute is associated with zero-offset, normal-incidence reflectivity. As such, amplitude inversion for acoustic impedance estimation, whenever possible, should be made from the intercept attribute derived from the amplitude inversion of CRP gathers.

Now compare each attribute section in Figure 11.2-41 with the synthetic equivalent at the well location and note that the match between the two within the reservoir zone generally is very good. While the intercept section shows the sandstone (green) with shale interbeddings (pale orange) within the reservoir unit, the gradient section indicates variations in fluid saturation (indicated by the different tones of orange) within the reservoir rocks. The fluid saturation can also be inferred from the Poisson reflectivity (equation 27) or the pseudo-Poisson (equation 40) section shown in Figure 11.2-41. The orange color in the fluid-factor section represents the gas-saturated reservoir sandstone.

AVO attributes are usually interpreted by crossplotting one attribute against another. Figure 11.2-42 shows three different zones — postreservoir (blue rectangle), within reservoir (red rectangle) and prereservoir (green rectangle), labeled on the intercept and gradient sections. The crossplots of the intercept and gradient attribute values within each of these zones are shown in Figure 11.2-43. The slopes of the best-fit lines for the postreservoir (-4.5) and the prereservoir (-5.5) zones are comparable. Whereas, the slope from the crossplot is significantly larger in the reservoir zone (-7.5). The increase in slope within the reservoir zone indicates an increase in the gradient attribute value, which in turn, suggests an increase in the change in Poisson’s ratio (equation 25), hence an indication of the gas-saturated reservoir sands.

3-D AVO analysis

The ultimate goal in seismic interpretation is to derive a reservoir model. The geometry of the reservoir model is provided by structural interpretation (interpretation of 3-D seismic data), the result of which is downscaled in the vertical direction and upscaled in the lateral direction to build the structural framework of the reservoir model. This means that the reservoir layer is sliced into layers as thin as 1 m in the vertical direction, and each thin layer is compartmentalized into cells with lateral dimensions as big as a few hundred meters.

The internal constitution of the reservoir layer itself is specified by using the result of a stratigraphic interpretation (interpretation of 3-D seismic data). The interior of the reservoir layer eventually needs to be populated by the petrophysical properties, such as porosity, permeability, and fluid saturation. When derived from 3-D seismic data, AVO attributes can be used to infer these properties within the reservoir zone in inline, crossline, and vertical directions.

The prestack signal processing sequence tailored for AVO analysis described in this section is applicable directly to 3-D seismic data. To generate the CRP gathers used in prestack amplitude inversion to derive the volume AVO attributes, however, you need to follow the workflow for 3-D prestack time migration (3-D prestack time migration). Figures 11.2-44 and 11.2-45 show the gradient and intercept volumes, respectively, derived from the CRP gathers as in Figure 7.4-21, following the application of opacity removal. While the opacity test applied to the image volume as shown in Figure 7.5-33 enabled us to identify the presence of a channel system within the reservoir zone, the opacity tests applied to the attribute volumes shown in Figures 11.2-44 and 11.2-45 indicate characteristics of fluid distribution within and away from the channel. This information is needed to choose the most desirable drilling locations and thus better manage the development of the reservoir. Note from Figure 11.2-45b that the amplitudes labeled in red that define the channel also are present along one of the faults. This suggests that the fault may represent the migration path for fluids. Also note the presence of a red spot at the center of the channel between the two faults. This may represent a thicker zone in the flood plain.

11.4 Acoustic impedance estimation

The goal in amplitude inversion is to obtain a broadband impedance model of a reservoir zone from a band-limited time-migrated CMP-stacked data. Alternatively, amplitude inversion can be applied to the CRP stack derived from prestack time migration or the intercept AVO attribute section. The trace amplitudes on a migrated CMP-stack, CRP stack, or the intercept section represent the one-dimensional (1-D), primary, P-wave reflectivity series associated with vertical incidence. Care must be given in processing to preserve relative amplitudes and increase vertical and lateral resolution of the seismic data. Reflection amplitude variations on any of these sections may indicate changes in some rock parameters, such as porosity or fluid saturation, within the reservoir unit.

By considering a sparse-spike reflectivity model of the earth, we first obtain a broad-band reflectivity section. Mathematically, the process involves minimizing the quantity that is equal to the sum of the absolute values of the trace amplitudes on the migrated CMP stack. Computation of the impedance series at a CMP location then involves simple integration of the broadband reflectivity series.

Recall from analysis of amplitude variation with offset that we assume a seismic source that generates a compressional plane wave propagating down into the earth at an angle from the vertical. When this incident plane wave encounters a layer boundary, which is assumed to be flat with no curvature, it is partitioned into four components — reflected and refracted P- and S-waves. The Zoeppritz equations (12a,12b,12c,12d) describe the amplitudes of the partitioned wave components. Aside from the layer parameters — density, P- and S-wave velocities, the wave amplitudes depend upon the angle of incidence (analysis of amplitude variation with offset). Implicit to CMP stacking, we must assume that the reflected P-wave amplitude variation with angle of incidence measured along the normal-moveout trajectory associated with the reflection event itself on the CMP gather is negligible. This assumption becomes critical for shallow targets at large source-receiver offsets. A test of the offset range used in stacking for the time window of interest is crucial for the meaningful interpretation of a poststack amplitude inversion result. By making the vertical-incidence assumption, we also assume that there is no P-to-S conversion within the offset range used in CMP stacking.

The special case of the Zoeppritz equation for vertically incident P-to-P reflection amplitude takes a simple form in which the amplitude depends only on the impedance values of the layer above and the layer below the reflecting boundary (equation 9). Based on this special case of the Zoeppritz equation, computing the impedance series from a reflectivity series involves a simple integration of the latter. The assumption of vertical incidence also requires that the stacked section input to inversion must be migrated.

Figure 11.3-2  A synthetic sonic log section derived from the stacked section in Figure 11.3-1.

Results of poststack amplitude inversion must be viewed within the realm of the underlying critical assumptions we make about seismic amplitudes. Amplitude inversion of poststack time-migrated data should be limited to earth models with low-relief structures and stratigraphic targets. Many of the assumptions indicated above are violated in the presence of reflectors with dip and curvature. To account for the dip and curvature effects, at least kinematically, amplitude inversion must be applied to prestack time-migrated data. Moreover, limitations in the removal of multiple reflections, linear and random noise must be taken into consideration when interpreting anomalies observed on an acoustic impedance section.

Synthetic sonic logs

From the discussion on the convolutional model of a seismic trace in the convolutional model, recall that the starting point was the sonic log (Figure 2.1-8a). The reflection coefficient is defined as the ratio of the reflection amplitude to the incident wave amplitude. In terms of acoustic impedance I = ρα, where ρ is rock density and α is the P-wave velocity, the reflection coefficient c is given by equation (9).

From equation (9), the reflection coefficient is interpreted as the ratio of the change in acoustic impedance to twice the average acoustic impedance. If we assume that density is constant, then equation (9) takes the form


(56)

where i is the sampling index. Therefore, reflection coefficients can be computed by differentiating the sonic log. The inverse process of getting the interval velocities for synthetic sonic logs from the reflection coefficients, which are assumed proportional to the stacked trace amplitudes, involves integration. In practice, only the high-frequency component of the interval velocity function can be obtained from this inversion. The low-frequency trend must be obtained from other sources of information such as conventional velocity analysis or nearby sonic logs. In many practical situations, there is insufficient well control and the bands of the low- and high-frequency information derived from the seismic data do not overlap. In these cases, problems arise in merging the various types of information; hence, the fidelity of the synthetic sonic logs is degraded.

Lindseth [30] first introduced the concept of synthetic sonic logs and used it in the interpretation of stratigraphic prospects. Figure 11.3-1 illustrates a stacked section with a bright-spot prospect. The corresponding synthetic sonic log section is shown in Figure 11.3-2 wherein the seismic data are superimposed as wiggle traces.

Measured sonic logs generally contain much higher frequencies than seismic data. Integration of seismic traces to get synthetic sonic logs implies a further lowering of the frequency content. A synthetic sonic log and a measured log can only be compared when a high-cut filter is applied to the measured log to make the two appear to have equivalent bandwidths.

The inversion described here is based on the recursive relation in equation (56). Specifically, given values for ci and αi, αi+1 can be computed. Model-based, iterative inversion schemes also exist [31]. These begin with an initial impedance function (specified at a CMP location) from which a synthetic seismogram is computed. This seismogram is crosscorrelated with the actual stack trace at that CMP location. The initial impedance model then is perturbed until the best match is attained between the estimated and actual seismic traces. Finally, inversion techniques, which are based on a certain type of characterization of the reflection coefficient series estimated from stacked data, also exist. For example, the estimated reflection coefficient series can be characterized by a sparse spike series. The sparseness requirement can be satisfied by the condition that the sum of the absolute values of the estimated reflection coefficients is minimum [32].

Note that regardless of the method, any inversion process suffers from the uncertainty (nonuniqueness) associated with the frequency components outside the dominant seismic signal bandwidth. Thus, the results of the inversion may be questionable for the low-and high-frequency ends of the spectrum at which noise dominates the signal. Constraints (information other than the seismic data being inverted, such as geologic or well control) often are incorporated into inversion schemes used to handle the low-and high-frequency ends of the spectrum.

Processing sequence for acoustic impedance estimation

The parsimonious signal processing sequence described in analysis of amplitude variation with offset tailored for AVO analysis also should be used for acoustic impedance estimation. Again, it is important to preserve the signal bandwidth of there-corded data during processing and attain a broad-band spectrum with a flat passband for the data input to poststack amplitude inversion. Also, the processing sequence must be designed to preserve relative amplitudes in the data. We shall analyze a 2-D marine line for acoustic impedance estimation. Data specifications for the line are listed in Table 11-4. In addition to the seismic data, sonic and density logs from a well located on the line traverse were used in estimating and removing the constant and linear phase associated with the residual wavelet in the poststack time-migrated CMP stacked data prior to amplitude inversion.

Table 11-4. Data specifications for the line used in acoustic impedance estimation.
Line Length 23.2 km
Shot Spacing 25 m
Receiver Spacing 16.67 m
CMP Spacing 8.33 m
Minimum Offset 80 m
Maximum Offset 3064 m
Number of Receivers 180
Fold of Coverage 60
Sampling interval 2 ms

The processing sequence included the following steps:

  1. Design a Wiener shaping filter (field data examples) to convert the recorded far-field source signature to its minimum-phase equivalent. Then, perform signature deconvolution by applying the filter to the recorded data.
  2. Mute guided waves since they are confined to the shallow water layer and contain no information about the subsurface reflectivity.
  3. Apply t2-scaling to compensate for geometric spreading. As mentioned in analysis of amplitude variation with offset, a velocity-dependent scaling function should be avoided to prevent overcorrection of amplitudes of multiples.
  4. Perform predictive deconvolution with unit-prediction lag to increase the vertical resolution, and remove short-period multiples and reverberations. With increasing prediction lag, the amplitude spectrum departs from a flat character (predictive deconvolution in practice). No data-dependent scaling should be applied to ensure preservation of relative amplitudes.
  5. Test CMP stacking for optimum fold (Figure 11.3-3). The principle criterion is to ensure minimal dependency of reflection amplitudes on the angle of incidence within the time window of interest. By using the NMO-corrected CMP gathers at velocity analysis locations, derive a spatially varying mute pattern such that much of the far-offset data are excluded from stacking. The CMP stack with optimum mute (Figure 11.3-3a) appears to best preserve the flatness character of the average amplitude spectrum of the data. The amplitude spectra and the autocorrelograms imply that there needs to be further flattening of the spectrum within the passband and removal of reverberations.
  6. Perform velocity analyses at frequent intervals, and obtain an optimum CMP stack using a time- and spatially varying mute function.
  7. Apply a frequency-space complex Wiener filter with unit-prediction lag (f − x deconvolution) to attenuate random noise uncorrelated from trace to trace (Figure 11.3-4a).
  8. Apply minimum-phase deconvolution to restore flat spectrum within the passband (Figure 11.3-4b).
  9. Apply time-variant spectral whitening to account for nonstationarity (Figure 11.3-4c).
  10. Design and apply a time-varying filter, while maintaining the broadest possible signal band. Again, any data-dependent scaling must not be applied to preserve relative amplitudes.
  11. Migrate the optimum CMP stacked data using a phase-shift algorithm. Migration compensates for the effect of reflector curvature on reflection amplitudes (equation 13a) so that the resulting amplitudes can be related directly to acoustic impedance variations. The reason for using the phase-shift algorithm is that, compared to finite-difference or integral methods, it causes the least phase and amplitude errors. Additionally, amplitude inversion often is applied to data associated with a velocity field that varies very mildly in the lateral direction.

Poststack time migration may precede the poststack signal processing, particularly deconvolution. As a result of migration, we transform a 2-D or a 3-D wavefield to a 1-D wavefield along the seismic traverse. Hence, we map the normal-incident energy to vertical-incident energy and obtain 1-D zero-offset seismograms at each CMP location. This satisfies the underlying assumption for poststack deconvolution.

Figure 11.3-5 shows the test panel for poststack processing in which migration precedes the signal processing. The sequence includes f − x deconvolution, phase-shift migration (Figure 11.3-5a), spiking deconvolution (Figure 11.3-5b), time-variant spectral whitening (Figure 11.3-5c), and time-variant filtering (Figure 11.3-5d). Finally, a mild (3-trace) Karhunen-Loeve filter was applied to the data to attenuate coherent linear noise and any remaining random noise uncorrelated from trace to trace (Figure 11.3-5e). The unmigrated CMP stack with the poststack processing as in Figure 11.3-4 is shown in Figure 11.3-6 and the migrated CMP stack with the poststack processing as in Figure 11.3-5 is shown in Figure 11.3-7.

Below is a summary of key issues regarding the processing sequence appropriate for amplitude inversion applied to poststack time-migrated CMP data.

  1. Avoid f − k filtering. This type of filtering invariably causes amplitude distortions. For multiple attenuation, use predictive deconvolution with a long operator length; the process will be effective at near offsets where periodicity of multiples are nearly preserved.
  2. Do not apply a velocity-dependent scaling function to compensate for geometric spreading. The amplitudes of multiples will be overcorrected with a velocity-dependent scaling function. A t2-scaling is favored for data to be used in amplitude inversion.
  3. Use the phase-shift method for migration. Differencing approximations to differential operators used in finite-difference migration induce phase errors and cause amplitude distortions. The usual implementation of the Kirchhoff migration (equation H-17, Section H.1) does not include all the terms of the integral solution (equation H-16) to the scalar wave equation. As such, the missing terms can influence the amplitude and phase accuracy of the resulting migrated data.
  4. You may be required to apply a harsh mute in stacking. Design a mute function at each velocity analysis location and apply a spatially varying mute pattern to the data for stacking. The design criterion is such that reflection amplitudes should be uniform within the offset range included in the stacking for an event.
  5. You may want to migrate before, and not after, poststack deconvolution. This is to conform with the fundamental assumption that deconvolution acts upon a 1-D zero-offset reflection seismogram recorded over a horizontally layered earth model. Amplitudes associated with any phantom diffractions should not be present on the traces input to deconvolution.

Derivation of acoustic impedance attribute

As input to poststack amplitude inversion, you may use one of the following data types:

  1. A time-migrated CMP-stacked section or volume of data,
  2. A CRP-stacked section or volume of data derived from prestack time migration,
  3. The intercept AVO attribute section or volume of data (equation 24), or
  4. The P-wave reflectivity section or volume of data (equation 37) derived from prestack amplitude inversion.

Poststack deconvolution and time-variant spectral whitening (Figures 11.3-5b,c) flatten the spectrum within the passband under the minimum-phase assumption. Nevertheless, there may be residual phase remaining in the stacked data; this residual phase typically is considered to have a constant and a linear phase component. Use of well data makes it possible to estimate the residual phase. After the completion of the conventional processing phase, the residual wavelet in the data input to amplitude inversion needs to be estimated and removed. This requires computing synthetic seismograms at well locations and comparing them with the stacked traces to determine constant-phase and linear-phase components associated with the residual phase left in the migrated data input to amplitude inversion.

  1. To start with, inspect and edit available sonic and density logs, and use check-shot information for depth-to-time conversion of the sonic logs (Figure 11.3-8a).
  2. Compute the acoustic impedance curve and obtain synthetic seismograms (Figure 11.3-9a) using zero-phase wavelets with a range of bandwidths that cover the passband of the time-migrated stacked data.
  3. Apply a range of constant phase rotations (Figure 11.3-9b) to the traces from the time-migrated CMP stacked section at the vicinity of the well and compare the results with the synthetic seismograms. In the present case study, we observed that the time-migrated stack did not require any phase correction within the upper portion of the zone of interest (2.5-3.5 s). Nevertheless, a better match was observed between the well data and the surface seismic data with a constant phase rotation of −90 degrees within the lower portion of the zone of interest (3.5-4.5 s). Figures 11.3-10 and 11.3-11 show portions of the time-migrated stacked section as in Figure 11.3-7 with the application of a constant phase rotation of −90 degrees.
  4. Compute a sparse-spike reflectivity model of the earth by a constrained L1-norm minimization applied to the time-migrated sections in Figures 11.3-10b and 11.3-11b.
  5. Apply a mild (3-trace) Karhunen-Loeve filter to the reflectivity sections to remove any spurious, geologically implausable variations in the reflectivity model (Figures 11.3-12a and 11.3-13a).
  6. Construct the acoustic impedance model from the broad-band sparse-spike reflectivity model by integrating the reflectivity series at each CMP location (Figures 11.3-12b and 11.3-13b) within the zone of interest (2.5 — 4.5 s).
  7. Repeat steps (d), (e), and (f) for the phase-rotated sections in Figures 11.3-10c and 11.3-11c. Results are shown in Figures 11.3-14 and 11.3-15.

We applied a range of constant phase rotations to a group of traces at two different locations on the migrated CMP stack and computed the acoustic impedance panels shown in Figures 11.3-16 and 11.3-17. Comparing these panels with the actual impedance curve at the well location (Figure 11.3-8b), we observed that the match between the two impedances in the upper portion of the zone of interest (2.5-3.5 s) is satisfactory without phase rotation applied to the data, while a better match is attained with a 90-degree phase rotation in the lower portion of the zone of interest (3.5-4.5 s). Any discrepancies may be attributed primarily to inaccuracies in the check-shot survey information.

It is important to evaluate the results of amplitude inversion with a cautious perspective. Specifically, one cannot and should not infer individual reservoir properties — pore pressure, confining pressure, porosity, permeability, and fluid saturation, from amplitude inversion. One can only hope to infer a composite effect associated with these properties. Given additional borehole information, however, one may be able to infer that a high impedance value may correspond to, for example, low porosity. In brief, results of amplitude inversion should be used as auxiliary information in support of other tools used in exploration and development of hydrocarbons.

3-D acoustic impedance estimation

As for the AVO attributes, when derived from 3-D seismic data, an acoustic impedance attribute volume can be used to infer lithology or porosity within a reservoir zone. We shall analyze a land 3-D data set for acoustic impedance estimation. Data specifications for the 3-D seismic data are listed in Table 11-5. The fold of coverage over the survey area is fairly uniform (Figure 11.3-18). In addition to the seismic data, sonic and density logs from a well from within the survey area were used in estimating and removing the constant and linear phase associated with the residual wavelet in the time-migrated CMP stacked data prior to amplitude inversion.

Note the high level of ground-roll energy on the selected shot records with the display gain shown in Figure 11.3-19. The noise characteristics vary from shot to shot. The processing sequence included the following steps:

  1. Apply t2-scaling to compensate for geometric spreading. Figure 11.3-20a shows a raw shot record and Figure 11.3-20b shows the same record after the application of t2-scaling.
  2. Perform spiking deconvolution (Figure 11.3-21a).
  3. Note from the spectral estimate shown in Figure 11.3-22 that deconvolution alone does not flatten the spectrum. Apply time-variant spectral whitening (Figure 11.3-21b) to attain a flat spectrum within the signal passband (Figure 11.3-22d) and attenuate the ground-roll energy.
  4. Sort the data into common-cell gathers and perform velocity analysis at coarse grid spacing.
  5. Create a 3-D stacking velocity field and apply NMO correction to the cell gathers.
  6. Apply residual statics corrections (processing of 3-D seismic data) and inverse NMO correction using the velocity field from step (e).
  7. Perform velocity analysis at tight grid spacing and create a new 3-D velocity field.
  8. Apply NMO correction to the cell gathers from step (f), mute and stack.
  9. Perform poststack deconvolution and a wide band-pass filter.
  10. Apply f − x deconvolution (linear uncorrelated noise attenuation) to attenuate random noise uncorrelated from trace to trace.
  11. Finally, perform 3-D poststack phase-shift migration. Note from Figure 11.3-23 that the processing sequence described above yields the spectral shape of the data that we wish to input to poststack amplitude inversion — a broadband spectrum with nearly a flat passband.
Table 11-5. Data specifications for the 3-D seismic data used in acoustic impedance estimation.
Survey size 27.3 km2
Inline dimension 7.35 km
Crossline dimension 3.75 km
Receiver line spacing 60 m
Bin size 30 × 30 m
Number of inlines 124
Number of crosslines 245
Number of bins 30,380
Fold of coverage 16
Sampling interval 2 ms

Figure 11.3-24 shows selected inlines from the image volume derived from 3-D poststack time migration. Following the principle of parsimony in processing for AVO analysis and acoustic impedance estimation, no DMO correction was deemed necessary in the present case study since the subsurface geology can be described by a horizontally layered earth model with no significant faulting or structural distortions.

The next phase in the analysis involves the estimation and removal of the residual phase contained in the time-migrated volume of data. Figure 11.3-25 shows a sonic and a density log measured at a well located at the intersection of inline 76 and crossline 79. Compute the acoustic impedance and the reflectivity series; then, use a zero-phase wavelet with a passband comparable to the signal bandwidth of the seismic data to compute the zero-offset, vertical-incidence synthetic seismogram also shown in Figure 11.3-25. Note from Figure 11.3-26 that the sonic log inserted into a portion of inline 76 at the well location exhibits a very good match with the time-migrated data within the reservoir level indicated by the arrow. To achieve this match, a constant time shift, which is equivalent to a linear phase shift, was applied to the sonic log. A similarly good match is observed in Figure 11.3-27 between the synthetic seismogram and the time-migrated data. To achieve this match, in addition to the constant time shift, a -90-degree phase rotation was applied to the synthetic seismogram. To remove the residual phase in the migrated data set and match it to the zero-phase synthetic seismogram, you therefore need to apply a 90-degree phase rotation to the former.

Figure 11.3-28a shows inline 76 of the image volume from 3-D poststack time migration, and Figure 11.3-28b shows the same section after the 90-degree phase rotation. The phase-rotated image volume is then used in amplitude inversion to create a sparse-spike reflectivity volume. Figure 11.3-28c shows inline 76 from the sparse-spike reflectivity volume. Finally, each sparse-spike reflectivity trace is integrated to estimate the acoustic impedance attribute (Figure 11.3-28d). Selected inlines from the acoustic impedance volume are shown in Figure 11.3-29. These inline traverses coincide with those shown in Figure 11.3-24 from the image volume derived from 3-D poststack time migration before the application of phase rotation.

The acoustic impedance volume is then interpreted to identify the spatial distribution of high-impedance and low-impedance areas at the reservoir level. First, pick the time horizon at the reservoir level from the image volume and isolate a horizon-consistent thin slab of subvolume. Then, apply transparency to the subvolume to obtain the interpretation result shown in Figure 11.3-30. The red represents the areas of high-amplitude reflections at the reservoir level. Apply the same interpretation strategy to the acoustic impedance volume and obtain a consistent result as shown in Figure 11.3-31. The orange tones represent the zones of high acoustic impedance contrast.

Instantaneous attributes

Aside from the AVO attributes, the instantaneous attributes can be helpful in inferring reservoir parameters. When considered as an analytic signal (in the mathematical sense), a seismic trace can be expressed as a complex function [33]. An analytic signal is expressed by a time-dependent complex variable u(t)


(57)

where x(t) is the recorded trace itself and y(t) is its quadrature. The quadrature is a 90-degree phase-shifted version of the recorded trace. It is obtained by taking the Hilbert transform of x(t) [34]


(58)

When substituted into equation (57), we have


(59)

Thus, to get the analytic signal u(t) for a seismic trace x(t), the complex operator [δ(t) + i/πt] is applied to the trace. When analyzed in the Fourier transform domain, this operator is zero for negative frequencies. Therefore, the complex trace u(t) does not contain negative frequency components.

Once u(t) is computed, it can be expressed in polar form as


(60)

where


(61a)

and


(61b)

R(t) represents instantaneous amplitude and ϕ(t) represents instantaneous phase at time t. In practice, instantaneous phase is computed by first taking the logarithm of both sides of equation (60)


(62)

and extracting the imaginary part


(63)

where Im denotes imaginary.

The instantaneous frequency ω(t) is the temporal rate of change of the instantaneous phase function


(64a)

When the derivative of equation (63) is taken with respect to time,


(64b)

For practical implementation, equation (64b) is written as the difference equation


(65)

The instantaneous amplitude measures the reflectivity strength, which is proportional to the square root of the total energy of the seismic signal at an instant of time. The instantaneous phase is used to emphasize the continuity of events on a seismic section. The temporal rate of change of the instantaneous phase is the instantaneous frequency. The instantaneous frequency may have a high degree of variation, which may be related to stratigraphy. However, it also may be difficult to interpret all this variation. Therefore, the instantaneous frequency values often are smoothed in time.

The instantaneous measurements that are related to an analytic signal are associated with an instant of time, rather than an average over a time interval. These measurements are reliable when the seismic signal is recorded and processed so that the CMP stack closely represents the subsurface. In other words, to deduce any stratigraphic meaning from the seismic data before estimating the instantaneous parameters, the amplitude and frequency content of the seismic signal must be preserved in each processing step. Any variation in the shape of the basic waveform that is not attributable to the subsurface geology must be eliminated. Multiples and all types of random noise limit the reliability of the results.

Reflectivity strength is an effective tool to identify bright and dim spots. Phase information is useful in delineating such interesting features as pinchouts, faults, onlaps, and prograding reflections. Instantaneous frequency information can help to identify condensate reservoirs and gas reservoirs, which tend to attenuate high frequencies.

Instantaneous attributes often are displayed in color for interpretation. Figure 11.3-32 shows the instantaneous attributes that correspond to the seismic section in Figure 11.3-1. The seismic data are displayed on top of the attributes as wiggle traces. Note the distinctive anomaly on the instantaneous amplitude section and enhancement of the continuity of reflections on the instantaneous phase section.

11.5 Vertical seismic profiling

Vertical seismic profiling data often provide more reliable correlation of well control to seismic data than synthetic seismograms derived from sonic logs. There are two reasons for this. First, the VSP data have a signal bandwidth closer to the seismic data than the sonic logs. More importantly, VSP data usually are not as sensitive to borehole conditions such as washouts.

VSP acquisition geometry

Acquisition of VSP data involves a surface source that is located either close to the well head (zero-offset case) or away from the well head (offset VSP) and a geophone in the well bore. Several traces are recorded at the same geophone depth, then edited and summed. (Repeatable sources such as vibrators or air guns are used.) The geophone then is moved to a new depth location and recording is repeated. The resulting section is a profile displayed in depth and time. A comprehensive reference on VSP is the book by Hardage [35].

Figure 11.4-1a shows VSP acquisition geometry. Consider a few of the raypaths that are intercepted by the receivers down the borehole: the direct arrivals from shot to receivers AC and AE, reflection ABC, and refraction ABF. Each receiver location yields a trace on the VSP that is in the depth-time domain (Figure 11.4-1b). Note that trace C has both the direct arrival (1) and the reflection from the first interface (3). The direct arrival (1 on trace C) and the refraction (4 on trace F) have only downgoing paths; therefore, they are called downgoing waves. On the other hand, reflection path ABC has a final upcoming segment BC (Figure 11.4-1a) and therefore constitutes an upcoming wave (Figure 11.4-1b). Note that direct arrival 2 on trace E coincides with a reflection arrival, provided the receiver is situated on the interface that causes the reflection. At this location, the upcoming and downgoing waves coincide (arrival 2 on trace E in Figure 11.4-1b).

Processing of VSP data

We now discuss the basic steps in VSP processing. After some trace editing, VSP processing starts with the separation of the downgoing waves from the upcoming waves (reflections). One separation technique is based on f − k filtering. Examine the zero-offset VSP in Figure 11.4-2a and note that the upcoming and downgoing waves have opposite dips. Because of this, each wave type should map in a different half plane in the f − k domain. Hence, downgoing waves can be suppressed with an f − k dip filter, thereby leaving only the reflection and associated multiples that constitute the upcoming waves (Figure 11.4-2b).

A VSP data set may not always have uniform receiver sampling in depth — a prerequisite for f − k filtering. The problems of edge effects and amplitude smearing often are observed on data after f − k filtering (frequency-wavenumber filtering).

An alternate approach to extracting upcoming waves is to use median filtering [35]. To start, apply first-arrival time shifts to traces in a VSP data set VSP(z, t) to flatten the downgoing waves. (First arrivals now are at t = 0.) Then apply a median filter to each horizontal array of samples, VSP(z, t = constant). Median filtering is best explained by an example. Consider the following array of numbers: (-1,2,1,2.5,1.5). Reorder its elements from small to large values: (-1,1,1.5,2,2.5). The median of the series then is the middle sample, 1.5. Median filtering rejects noise bursts and any event that is not flat. Apply the median filter to the VSP data set with flattened downgoing waves to yield the downgoing waves. This result can be subtracted from the input to obtain the upcoming waves. The last step involves unflattening the data.

The next VSP processing step involves datuming all receivers to the well head (D in Figure 11.4-1a). From Figure 11.4-1, this static correction is the same as correcting each trace by an amount that is equal to the one-way traveltime down to the corresponding receiver location. For example, trace C is corrected by an amount equal to the traveltime associated with raypath DC.

The static corrections are followed by deconvolution and filtering (Figure 11.4-2c). In principle, the deconvolution operators can be designed from downgoing or upcoming waves. These deconvolution operators then are applied to traces of the upcoming wave profile. The common practice is to use the downgoing waves to design the deconvolution operators. This is because downgoing waves on a VSP record are much stronger than upcoming waves. Hence, designing deconvolution operators from downgoing waves has an advantage in that the operators are based on stronger signal, with more emphatically represented multiples [35].

The last step involves stacking the traces in Figure 11.4-2c. Stacking normally includes a narrow corridor along the region in which upcoming and downgoing waves coincide. The resulting trace, repeated a few times, is shown in Figure 11.4-2d. To a large degree, corridor stacking prevents the multiples that do not merge with the downgoing wave path from being stacked in. The trace in Figure 11.4-2d can be considered an alternative to a zero-offset synthetic seismogram derived from the sonic log; thus, it can be compared to the CMP stack (not shown here) at the well location.

Figure 11.4-3a shows raw data from an offset VSP data set. By using the median filtering scheme described earlier, downgoing waves are extracted from the raw data (Figure 11.4-3b). The subtraction of downgoing waves from original data (followed by another median filtering for suppression of noise bursts) yields upcoming waves (Figure 11.4-3c). The upcoming waves then are deconvolved (Figure 11.4-3d) using operators designed from the downgoing waves of Figure 11.4-3b. This step normally is followed by static corrections. For nonzero-offset data, we also must correct for moveout resulting from the offset separation between the well head and the shot location. From Figure 11.4-1, this dynamic correction involves mapping the traveltime associated with raypath ABCD to 2DE. An NMO-corrected upcoming-wave profile can be compared with the surface seismic at the well location, provided the subsurface consists of horizontal layers with no lateral velocity variations.

VSP-CDP transform

When there are dipping interfaces, the upcoming-wave profile needs to be migrated; that is, the energy must be mapped to the actual subsurface reflection points. This is true even for zero-offset VSP data. A ray-tracing procedure for reflector mapping is illustrated in Figure 11.4-4. Note that reflection points D, E, and F have different lateral displacements OA, OB, and OC, respectively, from borehole Oz (Figure 11.4-4a). However, upcoming wave energy from all three reflection points is recorded on the same VSP trace at the receiver location R. The reflection times RG, RH, and RK (Figure 11.4-4b) are associated with raypaths SDR, SER, and SFR, respectively. Mapping this energy to reflection points involves a coordinate transformation [36] [37]. In this transformation, the amplitudes on a single VSP trace are mapped onto several traces on the x − t plane, where x is the lateral distance of reflection points from the borehole (Figure 11.4-4b). The event times RG, RH, and RK are mapped onto two-way vertical times AL, BM, and CN, respectively. These vertical times are associated with raypaths AD, BE, and CF in Figure 11.4-4a. The resulting (x, t) section consists of traces similar to the traces of a migrated zero-offset section. Hence, the described ray-tracing procedure often is called VSP-CDP transformation.

The VSP-CDP transform of the data in Figure 11.4-3d is shown in Figure 11.4-3e [38]. Comparison with the migrated surface seismic section at the well location indicates a good correlation of events. The difference in frequency content is partly attributed to differences in processing these two sections and is partly attributed to less high-frequency attenuation effects because of the shorter traveltimes associated with VSP recording.

The VSP-CDP transformation requires knowledge of the velocity-depth model around the borehole, since we must determine the location of the reflection points in the subsurface to perform mapping. The velocity-depth model can be derived using an iterative approach [37]. Starting with an initial velocity-depth model and the recording geometry for the VSP data, traveltimes for upcoming waves are computed. When these estimated traveltimes are compared with the observed traveltimes, discrepancies are noted and the velocity-depth model is modified accordingly. The process is repeated until a good match is attained between the estimated and observed traveltimes.

Figure 11.4-4  (a) Source-receiver geometry for offset VSP; (b) VSP-CDP transformation: Traveltimes RG, RH, and RK, associated with raypaths SDR, SER, and SFR are mapped to traveltimes AL, BM, and CN, associated with the two-way vertical raypaths 2AD, 2BE, and 2CF. Also note that the amplitudes on the VSP trace are positioned on traces after transformation with x-coordinate values the same as those of the reflection points OA, OB, and OC.

Finally, note that the VSP-CDP transform is not exactly a migration process. It handles neither diffractions nor curved interfaces. To handle these features, VSP data must be migrated [39]. The VSP geometry is like the geometry of a common-shot gather, except the shot axis is perpendicular to the receiver axis. Migration of VSP data can be viewed as mapping amplitudes along semielliptical trajectories with their focal points being the source and receiver locations. Superposition of all these trajectories yields the migrated section. The aperture width for VSP data often is inadequate to obtain a migrated section without much smearing.

11.6 4-D seismic method

There is a strong resemblence between the techniques used in clinical medicine and geophysical prospecting. Here, I will refer to this resemblence within the context of the 4-D seismic method. A heart patient is monitored from year to year by the patient’s cardiologist to observe and detect changes in the parameters that describe the heart itself and its condition to sustain the patient’s life, and the composition of the patient’s blood. By using an echocardiogram derived from ultrasound waves, the cardiologist measures the size of the heart and observes whether the valves have any leakage. By recording an electrocardiogram, the cardiologist observes the systolic pressure associated with the rythmical contraction of the heart during which the blood is pumped out and the diastolic pressure associated with the rythmical dilatation of the heart during which the blood is pumped in. It is not just the heart’s parameters themselves that are important to the cardiologist, but also the critical changes in those parameters from year to year during the patient’s life span. Eventually, the cardiologist uses the data associated with the patient’s medical history to judge whether or not a surgical intervention is required, at which time an angiogram is taken to make direct observations of the heart to provide the necessary instantaneous information to the surgeon. When evaluating the historical data, the cardiologist is cognizant of the fact that the heart’s parameters are directly influenced by the evolving technologies used in taking the measurements, and the interpretation of the measurements is influenced by the changing medical knowledge and experience of the cardiologist.

Similarly, the reservoir geophysicist uses the 3-D seismic method combined with the direct observations made at well locations to monitor the reservoir conditions that are crucial for optimum development of the field. The objective in optimum field development is to lengthen the life span of the field, prevent water invasion, and recover as much hydrocarbon as possible. By recording 3-D seismic data over the field at various time intervals, which may be from months to years, we introduce the fourth dimension to the the analysis of the data — calendar time, thus the term 4-D seismic method to describe the time-lapse 3-D seismic exploration.

What type of reservoir parameters can we monitor using the 4-D seismic method? Known applications of the 4-D seismic method are summarized below [40] [41],[42] [43] [44].

  1. Monitoring the spatial extent of the steam front following in-situ combustion or steam injection used for thermal recovery,
  2. Monitoring the spatial extent of the injected water front used for secondary recovery,
  3. Imaging bypassed oil,
  4. Determining flow properties of sealing or leaking faults, and
  5. Detecting changes in oil-water contact.

A demonstrative example of the 4-D seismic method is shown in Figure 11.5-1 [44]. The two seismic sections along the same inline traverse have been extracted from the image volumes derived from 3-D poststack time migrations of two time-lapse 3-D seismic data. One survey was conducted in 1989 before the production commenced in the field and the other survey was conducted in 1998 sometime after the production was started. The oil-water contact (OWC) is distinctly visible in the 1989 section, while it is not apparent in the 1998-section. Additionally, the top-sand reflector is stronger in the 1998 section. This may be attributed to the increased impedance contrast between the overlying shales and the reservoir sands which have higher water saturation as a result of continuing production.

Figure 11.5-2 shows another application of the 4-D seismic method to steam flooding of a reservoir for thermal recovery. In this field study, six time-lapse 3-D seismic surveys were conducted [42]. The first survey at time T1 was conducted prior to steam injection, while the subsequent surveys were conducted after the steam injection was started. The red, near-circular feature on the time slices correspond to the spatial extent of the injected steam; note how the steam front expands further in the northwesterly direction with increasing calendar time.

Processing of 4-D seismic data

Just as in the case of the medical example given above, processing, inversion, and interpretation of 4-D seismic data are influenced by the evolving technologies in 3-D seismic exploration. The different vintages of 3-D seismic data that are used in a 4-D seismic project are often recorded with different vessels, source and cable geometries, and source and receiver types and arrays. In fact, some 4-D seismic projects may involve, say, two time-lapse 3-D surveys — one conducted using the conventional streamer cable and the other conducted using the ocean-bottom 4-C technique (next section). The 3-D surveys most likely would be conducted using different recording directions and bin sizes. Figure 11.5-3 shows the base maps for two time-lapse 3-D surveys with different recording directions and bin sizes [45]. The data associated with the 1979 survey and the 1991 survey were recorded with a 34-degree difference in grid orientation. Also, the bin size for the 1979-survey was 80 × 27.5 m, whereas the bin size for the 1991 survey was 12.5 × 12.5 m. Additionally, these 3-D seismic data sets most likely would be processed differently — not only the processing sequences would be different but also the processing parameters. Figure 11.5-4 [45] shows a section and a time slice from each of the two time-lapse 3-D surveys which are referred to in Figure 11.5-3. The 1991 survey data have produced a more accurate image of the salt flank.

Hence, the time-lapse 3-D data sets used in a 4-D seismic project need to be cross-equalized prior to the interpretation of the results. Cross-equalization involves the following steps [45].

  1. Align the grids of the time-lapse data to a common grid orientation. In many cases of 4-D seismic projects, grid alignment and subsequent steps in cross-equalization are applied to poststack data. As such, grid alignment may be achieved by remigrating the poststack data to a specified common grid orientation. If you have access to prestack data, one way to align the grids of the time-lapse data is by crossline migration (3-D prestack time migration), the output of which would be common-azimuth gathers.
  2. Apply spectral balancing to the time-lapse 3-D data to account for the differences in the spectral bandwidth and shape. Figure 11.5-5 shows the amplitude spectra computed from the 1979 and 1991 survey data shown in Figure 11.5-4. Note the significant differences in the shape and bandwidth of the spectra before cross-equalization. These differences have been minimized by cross-equalization.
  3. Derive amplitude gain curves from the time-lapse 3-D data based on trace envelopes, and apply the gain curves for amplitude balancing.
  4. Estimate static shifts between the time-lapse data traces and apply them to eliminate vertical time diffferences.
  5. Examine and determine differences in event positioning in the migrated data volumes associated with the time-lapse 3-D surveys. Eliminate the differences in event positioning by residual migration.
Figure 11.5-1  Time-lapse seismic data: (a) preproduction (1989), and (b) postproduction (1998). ([44]; figure courtesy MacLeod et al., Chevron and Schlumerger Geco-Prakla; data courtesy Chevron.)

Seismic reservoir monitoring

Figure 11.5-6 shows a difference section from the time-lapse data as in Figure 11.5-4 following cross-equalization. This difference section exhibits a strong amplitude anomaly at the reservoir level situated at the salt flank. Such an amplitude difference may be attributed to changes in the reservoir conditions as a result of production [45]. Because of a wide range of factors associated with acquisition and analysis of the 4-D data, in addition to the difference data volume, the individual data volumes themselves are also visualized and interpreted.

The example of cross-equalization shown in Figure 11.5-7 relate to a steam injection project. Note the differences in the time slices from the image volumes associated with the 1996 survey and 1997 survey before and after cross-equalization. The bubbles correspond to the location of the injection wells.

The 4-D seismic anomalies are characterized as differences between time-lapse 3-D data that are present after cross-equalization as exemplified by Figure 11.5-5. Calibration of these anomalies often is ambiguous, in that, they may be attributable to changes in one or more of the reservoir conditions, such as change in fluid saturation caused by water displacing oil, pore pressure change caused by injection, or a temperature change caused by steam injection [43].

Although significant progress has been made in the 4-D seismic method, its value in determining dynamic reservoir properties is just beginning to be demonstrated. The information regarding the dynamic reservoir properties much sought after by the production engineer includes changes in oil saturation, water saturation, and pore pressure. Future developments in seismically driven reservoir characterization and monitoring should contribute significantly to optimum management of oil and gas fields.

11.7 4-C seismic method

Throughout this textbook, we discussed processing, inversion, and interpretation of compressional or P-wave seismic data. Specifically, we had in mind a compressional seismic source and the reflected signal recorded by each receiver being a compressional seismic wave-field. Recall from analysis of amplitude variation with offset that for non-normal incidence at a layer boundary, an incident compressional plane wave is partitioned into not just reflected and transmitted compressional-wave components, but also reflected and transmitted shear-wave, or S-wave components. Hence, a fraction of the incident P-wave is converted into a reflected S-wave. The amplitudes of the individual components are described by the Zoeppritz equations (12). By way of prestack amplitude inversion of nonzero-offset P-wave data described in analysis of amplitude variation with offset, we are then able to make an estimate of the S-wave reflectivity as an AVO attribute (equations 37 and 49).

In a conventional marine seismic survey, we cannot record P-to-S converted-wave energy even if we deploy sensors that can register the shear-wave energy. This is because the upcoming converted-wave energy is not transmitted through the water column to reach the recording cable since fluids cannot support shear strain. Thus, to capture the converted-wave energy, we need to record it at the water bottom using an ocean-bottom cable (OBC). And to record it, we need to use geophones that register velocity of the particle motion that is perpendicular to the direction of the wave propagation. Since the upcoming wave is primarily in the vertical direction, we need to use a geophone that records the particle motion in the horizontal direction. In fact, for a good reason that will be obvious later in the section, we need to deploy not one, but two horizontal geophones that are oriented perpendicular to one another. To complement the recording of the pressure wave by a hydrophone, again for a reason that will be obvious later in the section, we might also wish to record the vertical component of the particle motion using a vertical geophone. Hence, an OBC recording is done using three geophones and one hydrophone for each receiver unit along the cable, making it a four-component (4-C) seismic survey. The final product from the analysis of a 4-C survey data is a pair of P-wave and S-wave image sections (in the case of a 2-D survey) or volumes (in the case of a 3-D survey). Strictly the P-wave data are associated with P-to-P reflections and S-wave data are associated with P-to-S converted waves. Heretofore, we shall refer to these two wave types as PP and PS, respectively, so as to explicitly indicate that they both are generated by a P-wave source. (An S-wave source would have generated SS and SP data.)

Much of the P-to-S conversion takes place, not at the water bottom, but at reflectors that correspond to layer boundaries with significant contrast in elastic properties [46]. This fortuitous phenomenon is caused by the very low speeds of shear waves in seabed sediments [47] [48].

In this section, we shall briefly review acquisition and analysis techniques specific to 4-C seismic data. Before we proceed, however, we need to ask the question why we would want to go through the expense of conducting a 4-C OBC survey. Is there any exploration or development objective that we cannot achieve using conventional P-wave data, but we can with S-wave data? To answer this crucial question, first, we make reference to one of the AVO interpretation strategies discussed in analysis of amplitude variation with offset. Specifically, by AVO inversion of prestack amplitudes, we wish to estimate P-wave and S-wave reflectivities and do a crossplot analysis to infer fluid types and saturation in reservoir rocks. If, instead of indirectly extracting these AVO attributes from conventional P-wave data, we record both PP and PS data by a 4-C survey, we may be able to broaden our understanding of properties of reservoir fluids and rocks.

Known potential applications of converted-wave data are summarized below [48] [49] [50].

  1. Imaging beneath gas plumes,
  2. Imaging beneath salt domes,
  3. Imaging beneath basalts,
  4. Delineating reservoir boundaries with a higher S-wave impedance contrast than P-wave impedance contrast,
  5. Differentiating sand from shale,
  6. Detection of fluid phase change from oil-bearing to water-bearing sands,
  7. Detection of vertical fracture orientation,
  8. Mapping hydrocarbon saturation, and
  9. Mapping oil-water contact.

We now refer to a few examples illustrating the use of PS data. Figure 11.6-1a shows portions of a dipole sonic log measured at a well from a producing field. The S-wave velocity curve shows a marked contrast at the top- (A) and base- (B) reservoir unit. It, however, does not show a significant contrast at oil-water contact (OWC). The P-wave velocity curve shows a difference in the gradients within the postreservoir and reservoir units, but it does not show a marked contrast at the top-reservoir boundary as does the S-wave velocity curve (A). This is because of a lack of acoustic impedance contrast between the shales of the postreservoir unit and the oil sands of the reservoir unit. The oil-water contact, on the other hand, corresponds to a significant contrast on the P-wave velocity curve. While the position and geometry of the oil-water contact are important in terms of the production history of the field, this is not the only strategic information that is needed for the development. Specifically, for production, it is the top-reservoir boundary that needs to be delineated accurately so as to place the horizontal well trajectory close to the top and avoid missing a significant vertical oil column.

Compare the PP section derived from the conventional streamer 3-D survey and the PS section derived from the 4-C OBC survey, both of the same vintage, shown in Figures 11.6-1b and 11.6-1c. The PP section shows a strong event at 2 s that corresponds to the strong contrast on the P velocity curve associated with the oil-water contact (Figure 11.6-1a). Nevertheless, the top reservoir is nearly impossible to identify on this section. The PS section, on the other hand, shows two strong events with irregular geometry at about 3.6 s and 3.8 s. These events correspond to the strong contrasts on the S velocity curve associated with the top- and base-reservoir boundaries labeled A and B in Figure 11.6-1a, respectively.

Figure 11.5-6  Difference section after cross-equalization of the time-lapse data as in Figure 11.5-4. ([45]; figure courtesy Rickett, Stanford Exploration Project, and 4th Wave Imaging; data courtesy Chevron.)

It is important to note that the PS section is not a replacement for the PP section; instead, they are complementary. The PP section provides the information about the oil-water contact while the PS section provides the information about the top-reservoir boundary. Both are needed for optimum development of the field.

Another case for converted-wave data is shown in Figure 11.6-2. The PP section shows a zone of gas plume associated with the underlying leaky reservoir. The gas-saturated formations cause amplitude and traveltime distortions of the P-wave that passes through the anomalous zone. This can make the geometry of the underlying reservoir unit difficult to delineate. The gas plume actually represents a complex overburden that gives rise to strong lateral velocity variations. As such, the complexity of the overburden can be resolved by earth modeling in depth and the underlying reservoir zone can, in some cases, be imaged by prestack depth migration with an acceptable accuracy. The S-wave, on the other hand, is relatively unscathed by the presence of gas within the overburden; hence, the PS data provide a more accurate time image of the reservoir zone as compared to the time image derived from the PP data. Note also from both the sections in Figure 11.6-2 and the 3-D image volumes shown in Figure 11.6-3 that the reflectors within the overburden above the gas plume zone are stronger in the PP data compared to the PS data. This is because only a fraction of the incident P-wave energy is converted to S-wave energy at layer boundaries. Once again, the PP and PS data complement each other, and it certainly is not one or the other. While the PS section provides a better image of the reservoir structure, the PP section clearly indicates the presence of a gas plume above. This anomalous pressure zone has to be taken seriously during well planning.

Recording of 4-C seismic data

Marine 4-C data are recorded by using ocean-bottom cables with receiver units, each containing one hydrophone to detect the pressure wavefield and three geophones to detect particle motions in a Cartesian system. The receivers used in marine multicomponent surveys are usually of gimballed type; as such, the vertical geophone is guaranteed to measure the vertical component of the particle motion. The two horizontal components measure the particle motions in two orthogonal directions, and they are intended to be oriented in such positions that one of them is aligned in the direction of the receiver cable. Shown in Figure 11.6-4 is a diagram of an ideal three-component geophone layout that would be possible to achieve in land surveys. The orientations of the three components coincide with a right-handed Cartesian coordinate system. This means that the vertical z-component is positive downward, while the inline x-component is defined to have positive direction when the crossline y-component is clockwise with respect to the x-component.

The geophone orientation of the layout shown in Figure 11.6-4 is not achievable in an ocean-bottom survey. Although the vertical geophone is indeed oriented in the vertical direction and it measures the particle motion as positive downward, the two horizontal geophones are not guaranteed to be in the inline and crossline directions. Instead, these two geophones position themselves at various different, still orthogonal, directions. As a result, the horizontal geophones associated with a common-shot record measure particle motions in arbitrary directions, rather than the desired common inline and crossline directions (Figure 11.6-5). This arbitrary horizontal geophone orientation is primarily caused by the seabed conditions such as currents, unconsolidated sediments, and the roughness of the seabed surface.

The sensor systems used in OBC surveys are of two types — the node systems and cable systems. A node is an individual unit that houses a single hydrophone and three geophones in the Cartesian orientation. The nodes are pressed into the seabed sediments by a remotely operated vehicle. Most nodal systems are derivatives of the SUMIC system pioneered by Statoil. The more widely used cable systems deploy different designs for housing the hydrophones and geophones. The cable usually is dragged a certain distance and then draped down to the seabed along the desired traverse.

The recording geometry for a 4-C OBC survey is an adaptation of a typical land 3-D survey geometry (3-D survey design and acquisition). As illustrated in Figure 11.6-5, two or more cables are laid down on the seabed parallel to each other, and data are recorded using a conventional seismic vessel with source locations aligned in the direction perpendicular to the receiver lines.

Shown in Figure 11.6-6 is a composite common-shot gather from a 4-C OBC survey made up of four records. The individual components are the hydrophone record, and the inline, crossline and vertical geophone records. This survey was conducted using two receiver cables laid in parallel, 300-m apart on the seabed. While the total cable length is 6000 m, maximum offset for some shots was up to 9000 m. Each cable carried 240 receiver units at 25-m interval; hence, each of the four records shown in Figure 11.6-6 contains 240 traces. The same common-shot gather is displayed in Figure 11.6-7 with AGC applied in order to better examine the signal and noise character in each record. Note that the hydrophone and vertical geophone records exhibit events with similar moveout since they both contain P-waves. The records associated with the two horizontal components exhibit events with much larger moveout since these records contain S-waves which travel with a slower velocity than P-waves, almost twice as slow in many rock types. Therefore, when identifying the same event on both a P-wave record and an S-wave record, keep in mind that the event time in the latter can be twice the event time in the former. This means that if you plan a 5-s record length for conventional P-wave data, you would need to record the 4-C data using a 10-s duration. The slower S-wave velocities would also result in a much larger moveout on the S-wave record compared to the moveout on the P-wave record. This means that spatial aliasing would be more serious when applying multichannel processing applications, such as f − k filtering or migration, to S-wave data than P-wave data.

Figure 11.6-8 shows a common-receiver gather which was created by sorting the common-shot gather data as in Figure 11.6-6. The same receiver gather with AGC applied is shown in Figure 11.6-9. The horizontal geophone records from both the common-shot gather (Figure 11.6-7) and the common-receiver gather (Figure 11.6-9) exhibit events with relatively more complex moveout than those observed on the hydrophone or vertical geophone records. You should not expect a one-to-one correspondence between the events on the two sets of records. Because layer boundaries give rise to differences in P-wave and S-wave impedance contrasts, there may be some events that appear on both records with different strengths, or some events may be present in one record set and absent in the other.

Figure 11.6-10 shows a close-up portion of the composite shot gather shown in Figure 11.6-7 and the spectra of the individual components. Note the polarity reversal from one side of the cable to the other on the inline record as a consequence of the right-handed recording convention (Figure 11.6-5). The frequency content of the hydrophone data appears to be broader than the vertical geophone component; this is because of the imperfect coupling of the geophones with the seabed sediments. This difference in bandwidth because of the coupling issue is observed also in the common-receiver gather shown in Figure 11.6-11.

Acquisition of 4-C OBC data is different from conventional streamer recording in respect of the receivers. In fact, it is like land acquisition at the seabed. When the seabed has irregular geometry, it gives rise to both long- and short-wavelength statics. Therefore, in processing 4-C data, receiver statics need to be calculated and applied to both PP and PS data.

Gaiser’s coupling analysis of geophone data

Variations in geophone coupling contaminate signal amplitudes registered by the geophone components, and need to be compensated for in a surface-consistent manner. Because of coupling problems, what is recorded by each one of the three geophones is not exactly the same as the ground motion at the seabed. A frequency-domain model equation that relates the recorded signal components {X′(ω), Y′(ω), Z′(ω)} by the three geophones in the inline, crossline, and vertical directions (x, y, z), respectively, and the actual ground motions in the three orthogonal directions {X(ω), Y(ω), Z(ω)} is given by [52]


(66)

where ω is the angular frequency, I is unity, and the nonzero elements Cy, Cz, Vy, and Vz describe the coupling response of the geophones.

Figure 11.6-3  3-D image volumes associated with (a) PP data and (b) PS data as in Figure 11.6-2. (Figure courtesy [51], and Schlumberger Geco-Prakla; data courtesy BP-Amoco.)

Note from equation (66) that X′(ω) = X(ω); this means that we assume that the inline geophone is perfectly coupled. Since the inline geophone is guided by the cable itself, this is considered a valid assumption in practice. Whereas the vertical and crossline geophones are not coupled completely — hence the nonzero elements of the coupling matrix. The imperfect coupling leads to vertical and crossline geophone signals mutually contaminating each other in a manner that can be modeled by equation (66).

We wish to estimate the ground motion vector {X(ω), Y(ω), Z(ω)}; this requires inverting equation (66) as given by Gaiser [52]


(67)

where D = VzCy − VyCz and is the determinant of the coupling matrix in equation (66).

From the matrix equation (67), write explicitly the recovered ground motions


(68a)

and


(68b)

The coupling compensation operators are estimated in a surface-consistent manner [53] with the constraint that, following the rotation, the energy of the transverse component is minimum.

Gaiser [52] reported a coupling experiment to verify the validity of the coupling theory described above. Figure 11.6-12a,b,c show the inline, crossline, and vertical geophone records obtained from an OBC survey. In the same figure, the record associated with the crossline geophone is shown after compenating for coupling (Figure 11.6-12d). To study the validity of the compensation based on equation (68), a diver firmly planted the receiver unit into the seabed and the recording was repeated. The resulting crossline record is shown in Figure 11.6-12e. If the coupling theory holds, then the records in Figures 11.6-12d,e should look very similar. Differences may be attributed to poor coupling of the planted receiver unit.

Figure 11.6-13 shows the result of surface-consistent coupling analysis. To apply the coupling corrections, scale the amplitudes in a given geophone record by the product of the source scalar and the associated component scalar. Figure 11.6-14 shows the common-shot gather as in Figure 11.6-6 after the application of surface-consistent amplitude corrections (Figure 11.6-12). The same shot gather with AGC is shown in Figure 11.6-15. The common-receiver gather as in Figure 11.6-8 after the application of surface-consistent amplitude corrections (Figure 11.6-12) is shown in Figure 11.6-16. The same receiver gather with AGC is shown in Figure 11.6-17. To examine the degree of compensation for differences in geophone coupling, refer to the close-up displays shown in Figures 11.6-18 and 11.6-19.

Processing of PP data

The next step in processing of the 4-C seismic data is to calibrate the vertical geophone component Z(t) and sum it with the hydrophone component P(t) to obtain the total PP data. This dual-sensor summation is done to remove water-column reverberations [54]. Calibration may involve just a single scalar applied to the entire geophone data prior to summation with the hydrophone data. More sophisticated calibration techniques include the application of surface-consistent scalars computed for each receiver location or application of surface-consistent scalars computed for each receiver location and for each frequency component [55] [56] [57].

The merged PP data are now ready for conventional processing. First, apply a vertical time shift that is equal to the water depth divided by water velocity to bring the receivers from the seabed to the same datum as the shots. If the water depth is greater than 100 m, the vertical shift may not be valid; instead, the datuming may have to be done by wave-equation datuming (layer replacement).

The remaining prestack processing sequence for the PP data is no different from the land data processing sequence and includes geometric spreading correction, deconvolution, refraction and residual statics corrections to account for the variations in the seabed geometry, velocity analysis, NMO and DMO corrections. The poststack processing sequence typically includes deconvolution, time-variant filtering, and migration. Shown in Figure 11.6-20 is a CMP gather associated with the PP data as in Figure 11.6-14. The CMP stack is shown in Figure 11.6-21. The average fold of coverage is 150.

Rotation of horizontal geophone components

Return to the OBC geometry shown in Figure 11.6-5. As discussed earlier in the section, by using gimballed geophones, the vertical geophone orientation can be maintained in a true vertical direction. The two horizontal geophones, however, cannot be oriented along the desired inline and crossline directions. Instead, they align themselves in arbitrary orientations from one receiver station to the next. We need to realign the horizontal geophones associated with one common-shot gather to a common orientation. One common orientation that is the source-centered Cartesian coordinates is shown in Figure 11.6-22. This means that the horizontal geophones of all receivers that contribute traces to the shot station with a circle around it are rotated from inline-crossline (x, y) coordinates (the acquisition coordinate system) to radial-transverse (r, t) coordinates relative to the source location (the processing coordinate system). As a result, the radial geophone response will represent the horizontal particle motion in the source-receiver plane and the transverse geophone response will represent the horizontal particle motion perpendicular to the radial response. Following the rotation, a common-shot or a common-receiver gather associated with the radial component will comprise traces with radial responses in the source-receiver azimuthal directions.

Figure 11.6-23 shows a common-receiver gather associated with the inline and crossline geophone components before and after rotation. The equations for coordinate transformation of the particle motions from inline-crossline (x, y) coordinates to radial-transverse (r, t) coordinates are [48]


(69)

where the rotation angle θ is as labeled in Figure 11.6-23c, Y(t) and X(t) are the inline and crossline components as recorded in the field following the compensation for variations in coupling, and R(t) and T(t) are the rotational and transverse components after rotation.

It is generally assumed that most significant signal components — reflections, diffractions, and converted waves, are polarized in the source-receiver direction [52]. This means that, following the rotation of the horizontal geophone components, the radial component would contain most of the PS energy while the transverse component would contain negligible PS energy (Figure 11.6-23). If the transverse component does contain anomalously high level of signal energy, it may be attributable to anisotropy that causes shear-wave splitting in fractured rocks. We shall review this pheonomenon briefly in seismic anisotropy.

Figure 11.6-24 shows the common-shot gather as in Figure 11.6-14 after rotation of the horizontal components. The hydrophone record (first record from left) and the vertical geophone record (fourth record from left) are the same in both figures. The second and third record from the left in Figure 11.6-14 represent the inline and crossline geophone components of the particle motion; whereas, the second and third record from the left in Figure 11.6-24 represent the radial and transverse geophone components of the particle motion. Note from the AGC applied version of the same shot gather in Figure 11.6-25 that the radial component, unlike the inline component shown in Figure 11.6-15, does not exhibit the polarity reversal at zero offset. Also, the transverse component in Figure 11.6-25 contains relatively weak signal energy when compared to the radial component. These observations are also verified in the common-receiver gather shown in Figure 11.6-26. Compare the AGC applied version of the receiver gather before (Figure 11.6-16) and after (Figure 11.6-27) rotation and note the switch in polarity and low signal level in the transverse component. The spectral analysis of the shot and receiver gathers are shown in Figures 11.6-28 and 11.6-29; compare these results with Figures 11.6-18 and 11.6-19.

Common-conversion-point binning

We learned in analysis of amplitude variation with offset that an incident P-wave is partitioned at a layer boundary into reflected and transmitted P- and S-wave components. Consider the raypath geometry in Figure 11.6-30a for an incident P-wave generated by the source S1 and a flat reflector. The reflection angle for the PP-wave is equal to the angle of incidence; however, the reflection angle for the PS-wave is smaller than the angle of incidence. As a result, the PP reflection will follow a symmetric raypath and be recorded at receiver location R2, while the PS reflection will follow an asymmetric raypath and be recorded at receiver location R1.

Now consider the common-midpoint (CMP) raypath geometry for a source-receiver pair S1R1 shown in Figure 11.6-30b. There are two reflection arrivals at the receiver location R1 associated with the PP and PS raypaths. The reflection point B at which the incident P-wave is converted to the S-wave is displaced in the lateral direction by some distance d away from the reflection point A at which the incident P-wave is reflected and recorded at the same receiver location R1 as the converted S-wave. This means that, for an earth model with flat layers, the PP-wave reflection points coincide with the midpoint location (Figure 11.6-31a); whereas, the PS conversion points do not (Figure 11.6-31b). As a direct consequence of this observation, the notion of a CMP gather based on sorting PP data from acquisition coordinates — source and receiver, to processing coordinates — midpoint and offset, such that traces in the gather have the same midpoint coordinate, is not applicable to PS data. Instead PS data need to be sorted into common-conversion-point (CCP) gathers such that traces in this gather have the same conversion point coordinate.

An important aspect of CCP sorting is that the asymmetric raypath associated with the PS reflection gives rise to a periodic variation in fold of the CCP gathers. As for the conventional P-wave data with variations in fold caused by irregular recording geometry, amplitudes of the stacked PS data are adversely affected by the variation in the CCP fold [58] [48]. Just as one resorts to flexible bin size in the processing of 3-D seismic data to accommodate variations in fold, the same strategy may be applied for the PS data processing.

Binning the PS data into CCP gathers requires knowledge of the conversion-point coordinate xP. Referring to Figure 11.6-31b, note that the conversion-point coordinate follows a trajectory indicated by the broken curve that, in general, depends on the reflector depth [59].

To derive an expression for xP, refer to the geometry of the PS-raypath shown in Figure 11.6-32. By Snell’s law, we know that


(70)

where α and β are the P-wave and S-wave velocities, respectively, and φ0 is the P-wave angle of incidence and ψ1 is the reflection angle for the converted S-wave.

From the geometry of Figure 11.6-32, note that


(71a)

and


(71b)

where xP and xS are the lateral distances from the CCP to the source and receiver locations, respectively. Substitute equations (71a,71b) into equation (70), square and rearrange the terms to get


(72a)

Apply some algebraic manipulation to solve equation (72a) for xS


(72b)

where γ = α/β. Finally, substitute the relation xS = x − xP, where x is the source-receiver offset, into equation (72b) to get the desired expression


(72c)

From Figure 11.6-31b, note that the CCP location moves closer to CMP location as the depth of the reflector increases. At infinite depth, the CCP location reaches an asymptotic conversion point (ACP) [60]. In the limit z → ∞, equation (72c) gives the ACP coordinate xP with respect to the source location


(73a)

Since β < α, the conversion point is closer to the receiver location than the source location (Figure 11.6-31b). The displacement d = xP − x/2 of the asymptotic conversion point from the midpoint is, by way of equation (73a),


(73b)

While CCP binning may be performed using the ACP coordinate given by equation (73a), more accurate binning techniques account for the depth-dependence of the CCP coordinate xP based on a solution to equation (72c) [59] [61]. Because of the quartic form of equation (72c) in terms of xP, an iterative solution may be preferred in practice [61] [62] [63]. The iteration may be started by substituting the asymptotic form of xP given by equation (73a) into the right-hand side of equation (72c). The new value of xP may then be back substituted into equation (72c) to continue with the iteration.

Whatever the estimation procedure, note from equation (72c) that xP depends both on depth to the reflector and the velocity ratio γ = α/β. Unless a value for the velocity ratio is assumed, it follows that CCP binning requires velocity analysis of PS data to determine the velocity ratio γ. Additionally, an accurate CCP binning strictly requires the knowledge of reflector depths; thus, the advocation of an implicit requirement that 4-C seismic data analysis should be done in the depth domain. This requirement may be waivered if we only consider a horizontally layered earth model as in the next subsection.

Velocity analysis of PS data

From the raypath geometry of Figure 11.6-32, it follows that


(74)

where t is the two-way PS reflection traveltime from the source to the conversion point to the receiver, and z is the reflector depth.

Set x = xP = 0 in equation (74) to see that the two-way zero-offset PS traveltime is given by


(75)

Substitute equation (75) into equation (74) for the depth variable z to get


(76)

where γ = α/β.

Equation (76) describes the PS-wave moveout observed on a CCP gather. Although it is derived for a single flat layer in a constant-velocity medium, this equation also is applicable to a horizontally layered earth model. In that case, α and β refer to the P- and S-wave rms velocities.

Note from equation (76) that the asymmetric raypath associated with the PS reflection shown in Figure 11.6-32 gives rise to a nonhyperbolic moveout even in the case of a flat reflector in a constant-velocity medium. A way to avoid dealing with nonhyperbolic moveout would be to make the small-spread approximation and consider the best-fit hyperbola


(77)

to the traveltime trajectory associated with a PS reflection on a CCP gather [59]. In equation (77), t and t0 mean the same traveltimes as in equation (76), and vNMO is the moveout velocity for PS-wave derived from the best-fit hyperbola; as such, it is neither the P-wave velocity α nor the S-wave velocity β. In fact, the Taylor expansion of equation (74) yields the relation [60]. By assuming a hyperbolic moveout described by equation (77), the PS data can be corrected for moveout and stacked in the same manner as for the PP data.

Practical experience, however, points to the unavoidable fact that the PS data exhibit strong nonhyperbolic moveout behavior. Shown in Figure 11.6-33 is a CCP gather after hyperbolic moveout correction using equation (77) and nonhyperbolic moveout correction using equation (76). Note the overcorrection at far offsets within 0-2.5 s.

The need for nonhyperbolic moveout correction for the PS data makes it compelling to conduct a multiparameter velocity analysis. Unlike velocity analysis using the hyperbolic moveout equation (77), where we only need to scan for one parameter, vNMO, velocity analysis using the nonhyperbolic moveout equation (76) suggests scanning for three parameters — the PP-wave velocity α, the velocity ratio γ = α/β, and the CCP displacement xP. But, in practice, we do not have to scan for all three parameters. Instead, the iterative procedure described below may be followed.

  1. To begin with, note that the PP-wave velocity α can be estimated directly by velocity analysis of the PP data set itself.
  2. We may assume an initial value for the velocity ratio γ = α/β and estimate an initial value for xP using equation (73a).
  3. Knowing α and xP, use the nonhyperbolic moveout equation (76) to scan for γ as a function of t0. Figure 11.6-34 shows a γ-spectrum computed from the CCP gather in Figure 11.6-33a.
  4. Pick a function γ(t0) at each CCP analysis location along the line over the survey area and derive a γ(x, t0)-section as shown in Figure 11.6-35.
  5. Use the γ(x, t0)-section and the PP-wave velocity α to calculate an updated value for xP(t0) from equation (72c).
  6. Substitute the updated xP(t0) and the estimated γ and α into equation (76) to perform the nonhyperbolic moveout correction (Figure 11.6-33c).
Figure 11.6-33  (a) A common-conversion-point (CCP) gather, (b) after hyperbolic moveout correction using equation (77), and (c) after nonhyperbolic moveout correction using equation (76). (Figure courtesy [64].)

Another strategy for velocity analysis of the PS-wave data is the direct estimation of the PS-wave velocity β, rather than estimating the velocity ratio γ. Return to equation (76) and rewrite it explicitly in terms of α and β as


(78)

By using the nonhyperbolic equation (78), follow the alternative procedure for PS-wave velocity analysis outlined below.

  1. Again, estimate the PP-wave velocity α as before using the PP data set itself.
  2. Also, assume an initial value for γ = α/β and estimate an initial value for xP using equation (73a).
  3. Knowing α and xP, use the nonhyperbolic moveout equation (78) to scan for β as a function of t0. Figure 11.6-36 shows a CCP gather and the computed β-spectrum. Compare with the α-spectrum in Figure 11.6-20 and note the difference in the velocity ranges in the two spectra.
  4. Pick a function β(t0) at each CCP analysis location along the line over the survey area and derive a β(x, t0)-section.
  5. Use the β(x, t0)-section and the PP-wave velocity α to calculate an updated value for xP(t0) from equation (72c).
  6. Substitute the updated xP(t0) and the estimated β and α into equation (78) to perform the non-hyperbolic moveout correction.

Figure 11.6-37 shows the PS-wave stack based on the alternative procedure described above. Note the differences between this section and the PP-wave stack shown in Figure 11.6-21. Specifically, the PS-wave stack shows some interesting reflector geometry between 4-5 s; this behavior is absent within the equivalent time window (2-2.5 s) in the PP-wave stack.

One subtle issue is related to the time at which a is specified in equations (76) and (78). The two-way zero-offset time t0 in these equations is associated with the PS-wave; whereas, α is estimated at two-way zero-offset time associated with the PP-wave. To distinguish the two zero-offset times, first, rewrite equation (75) for the PS two-way time


(79a)

For the same depth z, the PP two-way time is given by


(79b)

Now eliminate z between the two equations to get the relation between the PP and PS zero-offset times


(80a)

or in terms of γ = α/β


(80b)

So, α in equation (78) is specified at time given by equation (80a), and a in equation (76) is specified at time given by equation (80b).

Dip-moveout correction of PS data

In principles of dip-moveout correction we learned that the DMO impulse response is in the form of an ellipse (Figure 5.1-12). We also learned that the aperture of the DMO operator depends on the reflector dip, the source-receiver separation, the reflection time of moveout-corrected common-offset data, and the velocity above the reflector. Since shear-wave velocities are generally lower than compressional-wave velocities, given that all other factors are the same, the DMO impulse response associated with the PS data will be different in shape as compared to the DMO impulse response associated with the conventional PP data shown in Figure 5.1-12. Specifically, the slower the velocity the more the action of the DMO operator (principles of dip-moveout correction). This means that the symmetric form of the DMO ellipse associated with the PP data is replaced with an asymmetric form as shown in Figure 11.6-38 (Alfaraj, 1993; Alfaraj and Larner, 1993; [65]. Additionally, the resulting curve is shifted laterally to the CCP location indicated by the dotted trajectory in Figure 11.6-38b (Alfaraj, 1993).

The PP DMO impulse response for variable velocity v(z) increasing with depth can be formed by squeezing the PP DMO impulse response for constant velocity [66]. Similarly, the asymmetric shape of the PS DMO impulse response shown in Figure 11.6-38b may be formed by squeezing the PP DMO impulse response on the source side and stretching it on the receiver side (Alfaraj, 1993).

Another important difference between the DMO corrections of the PP and PS data is in respect to the reflection-point dispersal (Section E.1). We learned, again in principles of dip-moveout correction, that the DMO correction removes the reflection-point dispersal along a dipping reflector. With the PP data, reflection-point dispersal is an issue only in the case of a dipping reflector. However, with the PS data, reflection-point dispersal takes place even for the case of a flat reflector. This is illustrated in Figure 11.6-39 where the raypaths of Figure 11.6-32 have been adapted to the CMP geometry. Note that you need to apply DMO correction to the PS data even in the absence of dip to remove the reflection-point dispersal (Alfaraj, 1993). Expressed differently, implicit to the PS DMO correction is CCP binning that involves mapping of amplitudes along the dotted trajectory in Figure 11.6-38b.

Figure 11.6-36  A CMP gather from the PS data as in Figure 11.6-14 and the semblance spectrum for PS velocity analysis. Compare with the velocity analysis of the PP data shown in Figure 11.6-20. (Processing by Orhan Yilmaz, 1999).

Migration of PS data

We learned in migration principles that migration of the stacked PP data can be performed by summation of amplitudes (equation 4-5) along a hyperbolic traveltime trajectory (equation 4-4) associated with a coincident source-receiver pair at the surface and a point diffractor at the subsurface. The amplitude from the resulting summation is placed at the apex time of the diffraction hyperbola. Prior to summation (equation 4-5), amplitude and phase factors inferred by the Kirchhoff integral solution to the scalar wave equation are applied to the stacked data (migration principles).

Similarly, migration of the stacked PS data can be conceptualized as a summation along a traveltime trajectory associated with a coincident source-receiver pair at the surface and a point diffractor at the subsurface. The difference between the zero-offset PP and the PS diffraction summation trajectories is in the velocities. Set x = 0 in equation (76) to get the traveltime equation for the zero-offset PS diffraction summation trajectory as


(81)

where t is the time associated with the unmigrated PS-stack, and t0 is the time associated with the migrated PS-stack given by equation (75).

Equation (81) is of the same form as equation (4-4) that describes the zero-offset PP diffraction summation trajectory. In fact, setting γ = 1 in equation (81) yields equation (4-4) for the case of constant velocity. Note from equation (81) that the migration velocity for the PS data is given by α/(γ + 1), which is neither the P-wave velocity α nor the S-wave velocity β. The PS-velocity α/(γ + 1) can be obtained from the PP-velocity α and the velocity ratio γ, both estimated from the velocity analysis described earlier in this section. Once the velocity field is estimated, the PS-stack can, in principle, be migrated using any of the algorithms described in migration.

11.8 Seismic anisotropy

A medium is anisotropic if its intrinsic elastic properties, measured at the same location, change with direction, and it is isotropic if its properties do not change with direction [67]. Most of seismic data analysis is based on the assumption that the subsurface behaves seismically isotropic. A case of anisotropy that is of interest in exploration seismology is the change in velocity with direction.

A medium is transversely isotropic if its elastic properties do not change in any direction perpendicular to an axis of symmetry. The usual meaning of seismic anisotropy is variation of seismic velocity, which itself depends on the elastic properties of the medium, with the direction in which it is measured [7]. Anisotropic variation of seismic velocity must not be confused with the source-receiver azimuthal variation of moveout velocity for a dipping reflector in an isotropic medium (equation 7-3).

There are two cases of seismic anisotropy that we shall review in this section, both of which are special cases of transverse isotropy. For convenience, consider a horizontally layered earth model. First, is the vertical transverse isotropy, otherwise simply referred to as transverse isotropy, for which velocities do not vary from one lateral direction to another, but vary from one direction to another on a vertical plane that coincides with a given lateral direction. Horizontal bedding and fracturing parallel to the bedding produce transverse isotropy. Second, is the horizontal transverse isotropy, otherwise known as azimuthal anisotropy, for which velocities vary from one lateral direction to another. Fracturing in a direction other than the bedding direction gives rise to azimuthal anisotropy.

Transverse isotropy stems from the fact that, as a result of a depositional process, velocities in a layer are different in the direction of bedding and the direction perpendicular to the bedding. Specifically, fine layering of isotropic beds that constitute the depositional unit give rise to an anisotropy with its axis of symmetry perpendicular to the bedding plane. Azimuthal anisotropy can stem from the fact that, as a result of a tectonic process, rock material associated with a layer may have different rigidity in different azimuthal directions. A specific example is secondary fracturing in rocks whereby the velocity in the fracture direction is higher than the velocity in the orthogonal direction, again, giving rise to an anisotropy with its axis of symmetry parallel to the bedding plane.

Refer to Figure 11.7-1 to review the physical aspects of wave propagation in an anisotropic medium. The directional change of the velocity is illustrated by the skewed ellipse in Figure 11.7-1a, with the fast velocity in the direction of the major axis and the slow velocity in the direction of the minor axis. Now consider Huygens’ secondary sources situated along the wavefront A, say at time t. This wavefront actually coincides with an exploding reflector dipping at an angle ϕ. Because of the anisotropic behavior of the propagation medium, the Huygens’ sources do not emanate semicircular wavefronts; instead, the wavefronts are skewed in the direction of the fast velocity. These skewed wavefronts will form the plane wavefront B at a later time t + Δt. While the energy was transmitted along the raypath SP at group velocity, the wavefront that represents a constant phase actually traveled from position A to B along TP normal to the wavefront at phase velocity. Because the group velocity is associated with the raypath, it is sometimes referred to as the ray velocity. Similarly, because the phase velocity is associated with the wavefront, it is sometimes referred to as the wavefront velocity. Note that the wavefront angle θ associated with the phase velocity is different from the ray angle ϕ associated with the group velocity. Only if the medium were isotropic, Huygens’ secondary sources would produce semicircular wavefronts and the phase angle θ would coincide with the ray angle ϕ.

Note from Figure 11.7-1 that the zero-offset raypath SP does not make a right angle with the reflector; thus, in the case of anisotropy, the zero-offset ray is not a normal-incident ray as would be in the case of isotropy. This behavior can be better explained by sketching the isotropic and anisotropic wavefronts as shown in Figure 11.7-1b. Specifically, in the case of isotropic medium, the wavefront emanating from a point P is circular and the zero-offset ray is normal-incident to the reflector. Whereas, in the case of an anisotropic medium, the wavefront is skewed and the zero-offset ray impinges on the reflector at a non-normal incidence angle.

An example of how velocity changes with direction in an anisotropic medium is shown in Figure 11.7-2 [68]. Based on velocity characteristics of the Green River shale [69], this polar plot of group velocities shows that, as for most rock types, the P-wave velocity is nonelliptical, whereas the SH-wave velocity behaves elliptically anisotropic. Also note that the horizontal P-wave velocity is greater than the vertical P-wave velocity.

We shall associate seismic anisotropy primarily with seismic velocities. Aside from anisotropic velocity analysis, this means that we will need to review those processes that are intimately influenced by velocity anisotropy, such as dip-moveout correction, migration and AVO analysis.

Figure 11.6-39  (a) Common-midpoint (CMP) raypath geometry for PP-data, and (b) common-conversion-point (CCP) raypath geometry for PS-data. Adapted from [65].

The generalized form of Hooke’s law, which is the foundation of linear elastic theory, states that each stress component can be expressed as a linear combination of all the strain components [70]. Hooke’s law is based on the assumption that elastic deformations in solids are infinitesimally small. The stiffness matrix, otherwise known as the elastic modulus matrix [69], that relates the stress components to the strain components is (Section L.1)


(82a)

The elements of the stiffness matrix are the elastic constants of an elastic solid. Since this matrix is symmetric, cij = cji, there are 21 independent constants for an elastic medium.

For an isotropic solid, the elastic behavior of which is independent of direction at a point within the solid, the number of independent elastic constants is only two, known as Lamé’s constants, λ and μ. Consequently, the stiffness matrix given by equation (82a) reduces to the special form


(82b)

To describe the P- and S-wave propagation in isotropic solids, only two elastic parameters are needed — Lamé’s constants, λ and μ. Other elastic parameters — Young’s modulus E, Poisson’s ratio σ and bulk modulus κ, and the P- and S-wave velocities can all be expressed in terms of λ and μ (Figure 11.0-1 and Appendix L.1).

For a transversely isotropic solid, the elastic behavior of which is the same in two orthogonal directions but different in the third direction, the number of independent constants is five [70] [69] [7]. For the more specific case of vertically transverse isotropy (VTI), which has a vertical symmetry axis, the five independent elastic constants are c11, c13, c33, c44, and c66 [69].

To explicitly describe the effect of anisotropy in wave propagation, Thomsen [69] has elegantly redefined the five elastic constants for the VTI media — the vertical P- and S-wave velocities, α0 and β0, in the vertical direction,


(83a)


(83b)

and three constants that describe the degree of anisotropy, ε, γ, and δ, in terms of the five constants c11, c13, c33, c44, and c66


(83c)


(83d)

and


(83e)

For most sedimentary rocks, the parameters ε, δ, and γ are of the same order of magnitude and usually much less than 0.2 [69]. In fact, the Thomsen parameters given by equations (83c,83d,83e) relate to the case of weak anisotropy described by small values (≪ 1) of ε, γ, and δ. While applications of anisotropy in seismic data analysis are primarily based on the assumption of weak anisotropy, these parameters still are useful to describe the general case of transverse isotropy.

The special case of δ = ε is known as elliptical anisotropy [71]. The ellipticity is associated with the shape of the wavefront expanding from a point source. Albeit its underlying theory is simpler than the general theory of anisotropy, elliptical anisotropy occurs in nature only rarely. Figure 11.7-3 shows a crossplot of the two Thomsen parameters, ε and δ, for various types of sedimentary and crystalline rocks based on field and laboratory studies [69]. Note that a few of the rock samples closely satisfy the ellipticity condition, ε = δ. Also, for most rock types, the Thomsen parameters, ε and δ, are positive constants less than 0.2; thus, the supporting evidence for weak anisotropy theory.

Anisotropic velocity analysis

We shall confine the discussion on anisotropy primarily to the practical case of P-wave propagation in weakly anisotropic rocks. Consider the wavefront geometry shown in Figure 11.7-1. The P-wave phase velocity is given by [69]


(84)

Note that the P-wave velocity depends on the anisotropy parameters δ and ε, but not on the parameter γ. In fact, the SV-wave velocity also depends only on δ and ε, while the SH-wave velocity depends only on γ.

In the special case of vertical incidence, θ = 0, equation (84) gives the vertical P-wave velocity, α0. In the special case of horizontal incidence, θ = 90 degrees, equation (84) gives


(85)

Solving for ε, we obtain


(86)

This equation provides a physical meaning for the Thomsen parameter ε. Specifically, the parameter ε indicates the degree of anisotropic behavior of the rock, measured as the fractional difference between vertical P-wave velocity α0 and the horizontal P-wave velocity αh. Since for most rocks ε > 0, note that the horizontal P-wave velocity is greater than the vertical P-wave velocity.

The normal-moveout velocity vNMO(ϕ = 0), where ϕ is the dip angle, for a flat reflector in an anisotropic medium is given by Thomsen [69]


(87)

In the special case of δ = 0, the moveout velocity is the same as the velocity of an isotropic medium represented by α0 (normal moveout).

The P-wave traveltime equation for a flat reflector in an anisotropic medium is given by Tsvankin and Thomsen [72]


(88)

where t is the two-way time from source to reflector to receiver, t0 is the two-way zero-offset time, x is the source-receiver offset, and A2 and A4 are coefficients which are written below to acknowledge their complexity


(89a)


(89b)

and


(89c)

Just for comparison, we write the traveltime equation (3-11) for a flat reflector in an isotropic medium using the current notation


(90)

For isotropic velocity analysis using equation (90), we only need to scan one parameter — the velocity vNMO (velocity analysis). For anisotropic velocity analysis using equation (88), note that we have to do a multiparameter scan involving α0, αh, β0, δ, and ε — an impractical and unacceptable proposal.

As a way to circumvent this insurmountable problem of a multiparameter scan in velocity analysis, Alkhalifah and Tsvankin [73] define a new effective anisotropy parameter


(91)

By way of equations (85), (87), (91), and (89a,89b,89c), and neglecting the effect of β0, equation (88) takes a form that is suitable for practical implementation [73]


(92)

where vNMO is the flat-reflector moveout velocity.

Equation (92) indicates that the traveltime for a reflector in an anisotropic medium follows a nonhyperbolic trajectory. By setting η = 0, equation (92) reduces to the isotropic case of equation (90). Note that the effect of anisotropy on reflection traveltimes is most significant at far offsets. The derivation of the anisotropic moveout equation (92) is based on a single flat layer. Nevertheless, in practice, it also is applicable to the case of a horizontally layered earth model with vertical transverse isotropy [74].

Figure 11.7-4a shows CMP raypaths associated with a flat reflector in a transversely isotropic medium with a velocity behavior as shown in Figure 11.7-2. Transverse isotropy causes fundamental departures from the CMP raypath geometry associated with an isotropic medium:

  1. The zero-offset raypath is non-normal incident.
  2. A single common reflection point does not exist; instead, anisotropy has given rise to reflection point dispersal even for a flat reflector.
  3. The anisotropic moveout given by equation (92) in general is nonhyperbolic.
  4. The NMO velocity measured from the slope of the t2x2 curve indicates that anisotropy makes the velocity offset dependent (Figure 11.7-4b).
Figure 11.7-1  (a) Application of Huygens’ principle to anisotropic plane-wave propagation from an exploding reflector after [68]. (b) The isotropic and anisotropic wavefront associated with an expanding P-wave after [7]. See text for details.

Figure 11.7-5 shows a CMP gather and its velocity spectrum computed using the isotropic traveltime equation (90). Note that there are some overcorrected events between the two-way zero-offset times of 1.5-2 s and some undercorrected events above the two-way zero-offset time of 1.5 s. The undercorrected events are clearly multiples which produce distinct semblance peaks on the velocity spectrum. The overcorrected events result from one of the three possibilities — that the overcorrection is caused by the dip effect on moveout velocities, anisotropy, or a combination of the two phenomena. In Figure 11.7-5, a close-up display of the zone of interest between 1.5-2 s shows the overcorrected events with a distinct moveout behavior that is not quite the same as what a typical overcorrected event with hyperbolic moveout looks like. Specifically, we see that the events are nearly flat with no moveout within the near-offset range, while the moveout is like the end of a hockey-stick within the far-offset range.

Try flattening the events by experimenting with various velocity picks as shown in Figure 11.7-6. Note that, whatever the pick you make on the velocity spectrum, the hockey-stick events do not quite flatten; there is always some curvature left along the moveout-corrected event.

Equation (92) indicates that, for anisotropic velocity analysis, we need to scan two parameters, vNMO and η. Note that the parameter η is only present in the fourth-order moveout term that is significant at far offsets. This suggests a two-stage parameter scan as described below.

  1. Perform hyperbolic velocity analysis using only the first two terms on the right-hand side of equation (92) and the near offsets where there is no hockey-stick effect. Figure 11.7-7 shows the same CMP gather as in Figure 11.7-5, but with far-offset traces muted, and the associated velocity spectrum. This first-stage analysis would give an estimate of the moveout velocity vNMO as in equation (92).
  2. Next, insert the estimated vNMO function into equation (92) and compute the η spectrum as shown in Figure 11.7-8. After picking an η function in time, apply the fourth-order moveout correction given by equation 92 to the CMP gather.

Figures 11.7-9 and 11.7-10 show portions of a CMP stacked section with and without accounting for anisotropy in velocity analysis and moveout correction. The portion shown in Figure 11.7-9 is abundantly rich in diffractions and the portion shown in Figure 11.7-10 contains steeply dipping fault-plane reflections conflicting with gently dipping reflections. Panel (a) in both figures shows the full-fold stack and panel (b) shows the near-offset stack based on isotropic velocity analysis and moveout correction using the first two terms of equation (92). Note that the full-fold stack has been adversely affected by the hockey-stick effect on moveout at far offsets, while the near-offset stack has better preserved the diffractions and dipping events. Panel (c) of both Figures 11.7-9 and 11.7-10 shows the full-fold stack derived from anisotropic velocity analysis and moveout correction using both the second-order and fourth-order terms in equation (92). Note that the full-fold stacking based on anisotropic moveout correction has better preserved the diffractions and dipping events compared to the full-fold stacking based on isotropic moveout correction. These observations are more evident on the migrated sections shown in Figures 11.7-11 and 11.7-12.

Figure 11.7-2  Group velocity functions for P- and SH-waves in Green River shale (After [75] [69].

Anisotropic dip-moveout correction

Since anisotropy directly influences propagation velocity, it certainly should have an impact on dip-moveout correction. Recall from principles of dip-moveout correction that 2-D dip-moveout (DMO) correction removes the dip effect on moveout velocities, and from processing of 3-D seismic data that 3-D dip-moveout correction removes both the dip and source-receiver azimuth effects on moveout velocities. Rewrite the normalized moveout velocity relation for the 2-D case [76] using the current notation


(93)

and note that, from the standpoint of velocities, isotropic DMO correction maps the moveout velocity vNMO(ϕ) for a dipping reflector to the moveout velocity vNMO(0) = α0 for a flat reflector.

While the DMO process corrects for the dip effect on moveout velocities as described by equation (93), it also maps common-offset data, which have been moveout-corrected using the flat-event velocities vNMO(0), to zero offset. Recall from principles of dip-moveout correction that this mapping is achieved by an integral transform given by equation (5-14a) where the factor A of equation (5-5) in the current notation is given by


(94)

Here, tn is the event time after NMO correction using the flat-event velocity vNMO, h is the half-offset and ϕ is the reflector dip.

As discussed in principles of dip-moveout correction, the integral transform factor A given by equation (94) requires knowledge of the reflector dip ϕ to perform the DMO correction. To circumvent this requirement, we use the relation (Section D.1)


(95)

which states that the reflector dip ϕ can be expressed in terms of midpoint wavenumber ky and frequency ω0 — the Fourier duals of midpoint y0 and zero-offset event time τ0 associated with the DMO-corrected data, respectively. Substitute equation (95) into equation (94) to get


(96)

The frequency-wavenumber domain dip-moveout correction that transforms the normal-moveout-corrected prestack data with a specific offset 2h from yn − tn domain to y0τ0 domain is achieved by the integral of equation (5-14a) where A is given by equation (96).

We now examine the DMO mapping, specifically the transform factor A of equation (94) for the case of an anisotropic medium. Consider a dipping reflector situated in an anisotropic medium with its symmetry axis not coincident with the normal to the reflector. For this general case, the normal-moveout velocity is given by [73]


(97)

where θ is the phase angle measured from the vertical, ϕ is the reflector dip (Figure 11.7-2), and α is the phase velocity evaluated in the direction coincident with the reflector dip, θ = ϕ.

By using equation (84) for the phase velocity, Tsvankin [73] computes the derivatives in equation (97) and derives an expression for the moveout velocity associated with the case of weak anisotropy. Below, we write his equation in normalized and compact form


(98)

where


(99)

Equation (98) states that anisotropic DMO correction maps the moveout velocity vNMO(ϕ) for a dipping reflector to the moveout velocity vNMO(0) for a flat reflector. By setting the anisotropy parameters ε = δ = 0, or making the assumption that the anisotropy is elliptical (ε − δ = 0), note that equation (98) reduces to the isotropic case given by equation (93).

Compare equations (93) and (98) and note that, for the case of ε − δ > 0, the isotropic DMO correction factor is smaller than the anisotropic DMO correction factor. This means that an anisotropic DMO correction would have a larger aperture than the isotropic DMO correction. However, depending on the magnitude of the anisotropy parameters, ε and δ, and their signs, the opposite could also be true. In fact, setting ε = δ in equation (98), the case of elliptical anisotropy, the effect of anisotropy on dip-moveout cancels out altogether.

We are now ready to redefine the transform factor A of equation (94) for anisotropic DMO correction [78]. To use the anisotropic moveout velocity of equation (98), first, rewrite equation (94) as


(100a)

Then substitute equation (93) to get the desired expression


(100b)

By substituting equations (98) and (95) into equation (100b), and making the weak anisotropy assumption, we obtain the transform factor given by Anderson and Tsvankin [78], which we write below in compact form,


(101)

where


(102)

and sin ϕ is given by equation (95) in terms of vNMO(0), ky, and ω0.

To implement the frequency-wavenumber domain anisotropic dip-moveout correction, again, use the integral of equation (5-14a) where A is given by equation (101).

By setting the anisotropy parameters ε = δ = 0, or making the assumption that the anisotropy is elliptical (ε − δ = 0), note that equation (101) reduces to the isotropic case given by equation (96). The implementation of anisotropic DMO correction requires a simple modification to the isotropic implementation, given by the factor [1 + (ε − δ)C]. Note that both for the moveout velocity given by equation (98) and the DMO transform factor given by equation (101), it is the difference ε − δ, and not the individual anisotropy parameters, that dictates the anisotropic effect. What remains to be determined is the difference ε − δ. Combine equations (87) and (91) to get the relation


(103)

The effective anisotropy parameter η and the flat-reflector moveout velocity vNMO(0) are estimated by anisotropic velocity analysis described earlier in this section (equation 92). The scaling velocity α0 is associated with a vertically incident P-wave in an isotropic, horizontally layered earth. Finally, note that isotropic DMO correction (equation 96) is velocity-independent, whereas anisotropic DMO correction (equation 101) is velocity-dependent by way of equation (102).

The anisotropic DMO impulse response can vary in shape and depart from the familiar elliptical shape associated with the isotropic DMO impulse response. Figure 11.7-13 shows some of the anisotropic impulse responses created by Anderson and Tsvankin [78] using a frequency-wavenumber domain DMO correction based on the above formulation. For ε − δ = 0, the case of elliptical anisotropy, the impulse response kinematically is identical to the isotropic impulse response. For ε − δ > 0, the aperture of the anisotropic impulse response broadens while maintaining the upward curvature, and for ε − δ < 0, the curvature is reversed downward.

An interesting theoretical conjecture can be drawn from the behavior of isotropic and anisotropic DMO impulse responses shown in Figures 11.7-13a,b. While anisotropy causes the DMO aperture to be broadened, vertically increasing velocity causes the DMO aperture to be narrowed (principles of dip-moveout correction). As a result of these two counteracting effects, we may conclude that anisotropic DMO correction that accounts for vertically varying velocity may be equivalent to constant-velocity isotropic DMO correction. Does this mean that, for most field data cases where velocities vary vertically, it suffices to perform constant-velocity isotropic DMO correction even in the presence of anisotropy? This would implicitly account for any anisotropic behavior that might be present in the data. In fact, Levin [79] conducted numerical studies of reflection times from a dipping plane in a transversely isotropic media using elastic parameters associated with four different sandstone, shale, and limestone samples, and concluded that, when the symmetry axis of transverse anisotropy is perpendicular to the reflector, the isotropic moveout velocity given by equation (93) is largely valid. Larner [80] extended this work to the case of a medium with vertically varying velocity and reached a similar conclusion that, for most cases of weak anisotropy, isotropic DMO correction is largely valid.

We refer to the field data example shown in Figure 11.7-14. Following isotropic DMO correction, the CMP gather as in Figure 11.7-5 exhibits the hockey-stick behavior of the anisotropic moveout more distinctly. Whatever the velocity pick, the isotropic moveout correction using equation (90) fails to flatten the events within the time gate 1.5-2 s (Figure 11.7-15). By first applying isotropic moveout correction using the near-offset traces in the gather, we derive an estimate of the normal-moveout velocity vNMO(0) as shown in Figure 11.7-16. This then is followed by the analysis of the anisotropy parameter η and applying the fourth-order moveout correction using equation (92) (Figure 11.7-17). Note that, despite the DMO correction that did not account for anisotropy, the subsequent anisotropic moveout correction has been successful in flattening the events almost completely.

Scan the anisotropic velocity analyses shown in Figures 11.7-18 through 11.7-23 to examine the extent of anisotropy manifested by these DMO-corrected gathers and how well the anisotropy has been accounted for by the post-DMO velocity analysis. The failure in flattening the events completely in some gathers may stem from various sources of error associated with the assumptions made about anisotropy, and it may even be caused by not correcting for anisotropy during DMO correction that preceded the velocity analysis.

Figures 11.7-24 and 11.7-25 show portions of a stacked section with isotropic DMO correction which was followed by velocity analysis with and without anisotropy accounted for. Compare with the corresponding panels in Figures 11.7-9 and 11.7-10 and note that the steep fault-plane reflections and diffraction flanks have been preserved by DMO correction. In fact, much of the task of preserving steep dips with conflicting dips already has been achieved by isotropic dip-moveout processing as seen in panel (a) in both Figures 11.7-24 and 11.7-25 that shows the full-fold stack based on isotropic velocity analysis and moveout correction using the first two terms of equation (92). For comparison, panel (b) in both Figures 11.7-24 and 11.7-25 shows the near-offset stack based on isotropic processing. Panel (c) of both Figures 11.7-24 and 11.7-25 shows the full-fold DMO stack derived from post-DMO anisotropic velocity analysis and moveout correction using both the second-order and fourth-order terms in equation (92). Compared with panel (a) in the same figures, the difference between the isotropic and anisotropic processing may be viewed at best as marginal. Again, this inconclusive result is attributable to the fact that DMO correction was done without accounting for anisotropy. The marginal differences are also evident on the migrated sections shown in Figures 11.7-26 and 11.7-27. By accounting for anisotropy in DMO correction, however, the subsequent imaging of the steeply dipping fault-plane reflections can be improved [81].

Throughout the various stages in data analysis, we encounter different moveout types:

  1. Normal moveout associated with flat events,
  2. Dip moveout associated with dipping events,
  3. Azimuthal moveout associated with 3-D recording geometry of varying source-receiver directions, and
  4. Anisotropic moveout caused by directional changes in velocity.

In practice, the moveout associated with an event observed in a CMP gather actually is a combination of all the contributions of the moveout effects listed above. It makes sense in practice to attempt to break up the total moveout into individual components and correct for them one at a time. A pragmatic workflow for moveout correction is NMO correction (moveout type (a)), followed by DMO correction (combined moveout types (b) and (c)), and anisotropic moveout correction (moveout type (d)). The 2-D field data example shown in Figures 11.7-26 and 11.7-27 is based on this moveout correction sequence with the exclusion of moveout type (c) related to 3-D data.

Anisotropic migration

A convenient way to understand the effect of anisotropy on migration is to refer to the zero-offset dispersion relation of the one-way scalar wave equation that describes wave propagation associated with exploding reflectors (migration principles), using the present notation,


(104)

where kz and ky are the wavenumbers associated with the depth z and midpoint y axes, ω is the temporal frequency and α is the medium velocity used to migrate a zero-offset compressional wavefield. The dispersion relation given by equation (104) is the basis from which all zero-offset finite-difference and frequency-wavenumber migration algorithms are developed (Sections D.1, D.4, and D.7). Hence, it makes good sense to examine what form this dispersion relation takes in the case of an anisotropic medium.

Accompanying the dispersion relation is the equation for wave extrapolation used in finite-difference and frequency-wavenumber algorithms (migration principles). Given a zero-offset upcoming wavefield P(y, z = 0, t) at the surface z = 0, for which a stacked section is a good approximation, the objective is to extrapolate it downward at depth steps of Δz using the extrapolation equation


(105)

This is then followed by invoking the imaging principle which states that the migrated section P(x, z, t = 0) at each depth z corresponds to the summation over the frequency components of the extrapolated wavefield, or equivalently collecting the extrapolated wavefield at time t = 0.

To account for the effect of anisotropy in migration, the medium velocity α0 in the dispersion relation given by equation (104) needs to be replaced with the P-wave phase velocity given by equation (84). This means that velocity used in anisotropic migration is dependent on the phase angle of the propagating wavefront (Figure 11.7-1). We shall first rewrite the dispersion relation of equation (104) explicitly in terms of the velocity


(106)

and replace α0 with the phase velocity of equation (84) specified for the reflector dip ϕ


(107)

Next, we shall use the relations (Section D.1)


(108a)

and


(108b)

for substitution into equation (106b) to get


(109)

By setting the anisotropy parameters ε = δ = 0, note that the anisotropic dispersion relation given by equation (109) reduces to the isotropic case given by equation (106).

A straightforward practical implementation of the anisotropic dispersion relation of equation (109) is within the framework of a frequency-wavenumber algorithm such as Stolt or phase-shift migration. By using an anisotropic dispersion relation defined for time migration, Anderson [82] have modified Fowler DMO and prestack time migration (prestack time migration) to account for anisotropy. Nevertheless, various practical forms of the dispersion relation can also be used to include the effect of anisotropy in frequency-space finite-difference explicit or implicit schemes [83] [84]. Although it is only a theoretically interesting conjecture, elliptical anisotropy may be accounted for by a vertical stretching in depth that follows an isotropic time migration [85].

The effect of anisotropy on an expanding wavefront is seen in Figure 11.7-28. Depending on the values of the Thomsen parameters, the wavefront can take various shapes, sometimes quite warped. These wavefronts can be considered as the kinematic aspect of an anisotropic migration impulse response. Note that, at certain dips, more migration takes place if anisotropy is accounted for as compared to isotropic migration. On the other hand, at some other dips, less migration may take place, depending on the anisotropy parameters.

A widely recognized case of anisotropy is from offshore West Africa. Figure 11.7-29 shows a popular example of anisotropic processing [81]. The data were processed by applying both isotropic and anisotropic DMO correction and migration. Note the better imaging of the steep fault planes by accounting for anisotropy in processing.

Note from the anisotropic dispersion relation given by equation (109) that it is important to supply an anisotropic migration algorithm with accurate estimates of the anisotropy parameters ε and δ. It is most likely safer to do isotropic migration than anisotropic migration with incorrect values for ε and δ. Alkhalifah and Larner [86] conducted some numerical studies for migration error in transversely isotropic media and concluded that, fortuitously, isotropic migration in the presence of anisotropy yields fairly accurate results for dips up to 50 degrees. For steeper dips, the isotropic assumption no longer is valid. Nevertheless, Alkhalifah and Larner [86] also caution the practitioner of the undesirable consequences of anisotropic migration using poorly estimated ε and δ parameters.

Effect of anisotropy on AVO

The P-to-P reflection amplitude as a function of angle of incidence given by equation (16) can be modified to accommodate for transverse isotropy as follows [88] [89]


(110)

where Δε and Δδ are changes in anisotropy across the flat interface that separates the upper and lower anisotropic media. By setting Δε = Δδ = 0, equation (110) reduces to the isotropic form given by equation (16). Since the origin of equation (16) is the Aki-Richards equation (15), equation (110) also is based on the assumption that, in addition to the small changes in elastic parameters, the changes in anisotropy parameters ε and δ are also small across the interface.

Figure 11.7-32  Shear-wave splitting from a North Sea 3-D survey: (a) radial and (b) transverse components as a function of source receiver-azimuth in degrees. See text for details. (Figure courtesy [90], and Baker-Hughes Western Geophysical.)

The terms of equation (110) can be split into two parts, R(θ) = Ri(θ) + Ra(θ) — the isotropic component Ri(θ) identical to the terms of equation (16), and the anisotropic component Ra(θ) given by [89] [87]


(111)

Figure 11.7-30 shows P-to-P reflection amplitudes as a function of angle of incidence at two different interfaces associated with a horizontally layered earth model with transverse isotropy [87]. The curves labeled as Ri correspond to the isotropic component given by equation (16) for three different values of β/α ratios and the curves labeled as Ra correspond to the anisotropic component of the reflection amplitudes given by equation (111) for a specific combination of Δε and Δδ values. Since Ra(θ) does not depend on β/α, the Ra(θ) curve is common to each of the three cases of the Ri(θ) curves. To get the total response given by equation (110), add the curves associated with the two components.

There are two important effects of anisotropy on AVO:

  1. The polarity of P-to-P reflection amplitudes can reverse for some combination of the anisotropy parameters.
  2. The anisotropic effect on reflection amplitudes is most significant at large angles of incidence.

Shear-wave splitting in anisotropic media

In 4-C seismic method, we reviewed the 4-C seismic method based on the conversion of P-waves to S-waves at layer boundaries. The converted S-wave propagating in an anisotropic medium exhibits a unique behavior. Consider a rock layer that contains vertical fractures oriented in, say the north-south direction, as sketched in Figure 11.7-31 — a case of a horizontally transverse isotropy, or equivalently, a case of azimuthal anisotropy. The S-wave that enters this transversely isotropic layer from below emerges from above as split and polarized in two orthogonal directions — the shear-wave component that is polarized parallel to the fracture orientation and the shear-wave component that is polarized perpendicular to the fracture orientation. Shear-wave splitting in azimuthally anisotropic media is sometimes referred to as shear-wave birefringence.

Figure 11.E-1  The moveout-corrected CMP gather referred to in Exercise 11-6.

The shear-wave component polarized parallel to the fracture orientation travels with a speed faster than the shear-wave component in the perpendicular direction. In practice, we have the radial and transverse components; but not the slow and fast shear-wave components. The fast and slow shear-wave components, however, can be extracted from the radial and transverse geophone components (4-C seismic method) by a special rotation developed by Alford [91]. The Alford rotational analysis of the shear-wave data provides an estimate of fracture orientation. This is strategically important in developing reservoirs within fractured rocks [92] [93]; specifically, in planning trajectories for production wells.

The split shear waves interfere with one another in a way that depends on the source-receiver azimuth. Figure 11.7-32 shows radial and transverse components sorted with respect to source-receiver azimuth. Each trace in these gathers is associated with a radial or transverse component in a specific azimuthal direction. Note that the radial component is fairly consistent from one azimuthal direction to the next. The transverse component, on the other hand, shows polarity change every 90 degrees. From these gathers, the Alford rotation produces the fast and slow shear-wave components while minimizing the energy in the transverse component.

Exercises

Exercise 11-1. Derive equation (2a) for the Fresnel zone using the geometry of Figure 11.1-3.

Exercise 11-2. Refer to Figure 11.4-1a. Consider a surface multiple from the first reflecting interface. Trace the traveltime on the VSP diagram in Figure 11.4-1b. Multiples do not reach the downgoing wave path; thus, they can be eliminated by corridor stacking.

Exercise 11-3. Refer to Figure 11.4-1b. Should the slopes of the downgoing and upcoming waves associated with a layer be the same in magnitude?

Exercise 11-4. What procedure does CMP stacking correspond to in the f − k domain?

Exercise 11-5. Sketch the traveltime response for a point scatterer on a zero-offset VSP record.

Exercise 11-6. Consider a CMP gather with a single reflection event. Suppose you have applied hyperbolic moveout correction using equation (90) and discovered that the event is not flat for all offsets. Instead, the moveout-corrected event may have one of the three shapes shown in Figure 11.E-1. Match each of the curves A, B, and C with the following three possibilities:

  1. You have applied hyperbolic moveout correction using an erroneously low velocity in equation (90).
  2. You have applied a second-order moveout correction (equation 3-4b) to an event that has a fourth-order moveout behavior described by equation (3-5a).
  3. You have ignored anisotropy in moveout correction described by equation (92).

Appendix L: Mathematical foundation of elastic wave propagation

L.1 Stress-strain relation

Seismic waves induce elastic deformation along the propagation path in the subsurface. The equation of wave propagation in elastic solids are derived by using Hooke’s law and Newton’s second law of motion. We shall begin with the stress-strain relation for elastic solids.

Solid bodies, such as rocks are capable of propagating forces that act upon them. Imagine an infinitesimally small volume around a point within a solid body with dimensions (dx, dy, dz) as depicted in Figure L-1. Stress is defined as force per unit area. The stress acting upon one of the surfaces, say dy − dz, can in general be at some arbitrary direction. It can, however, be decomposed into three components: Pxx which is normal to the surface, and Pxy and Pxz which are tangential to the surface. The first subscript refers to the direction of the normal to the surface, and the second subscript refers to the direction of the stress component. The components of the stress that are normal to the surfaces associated with the infinitesimal volume shown in Figure L-1 are called the normal stress components and the components of the stress that are tangential to the surfaces are called the shear stress components. A normal stress component is tensional if it is positive, and compressional if it is negative. To retain the solid volume depicted in Figure L-1 in equilibrium requires nine stress components that make up the stress tensor


(L-1)

The diagonal elements are the normal stress components and the off-diagonal elements are the shear stress components.

As the dimensions of the volume surrounding a point inside a solid body are made infinitesimally smaller, the sum of the moments of all surface forces about any axis must vanish. This requirement causes Pxy = Pyx, Pxz = Pzx, and Pyz = Pzy, making the stress tensor (L-1) symmetrical.

Fluids cannot support shear stress. In a fluid medium, only one independent stress component remains, since Pxx = Pyy = Pzz = −p, the hydrostatic pressure (the negative sign signifies the compressive nature of the hydrostatic pressure), and Pxy = Pxz = Pyz = 0.

Stress induces deformation of solid bodies. Consider two points, P and Q, within a solid body as indicated in Figure L-2. Subject to a stress field, the solid is deformed in some manner and the particles at points P and Q are displaced to new locations P′ and Q′. The displacement components δu, δv, and δw can be expressed by the relation


(L-2)

To understand the physical meaning of the elements of the tensor containing the partial derivatives in equation (L-2), instead of arbitrary displacement as shown in Figure L-2, consider deformations of specific types. The simplest deformation is an extension in one direction as a result of a tensional stress as illustrated in Figure L-3a. (For simplicity, only one surface of the volume in Figure L-1 is considered.) The fractional change in length in the x-direction is δu/δx, which, in the limit as the volume becomes infinitesimally small, is defined as the normal strain component.

Strain is a dimensionless quantity. There are three normal strain components, exx, eyy, and ezz:


(L-3a)


(L-3b)


(L-3c)

The stress field away from the typical seismic source is so small that it does not cause any permanent deformation on rock particles along the propagation path. Hence, the strain induced by seismic waves is very small, usually around 10−6. A positive strain refers to an extension and a negative strain refers to a contraction.

Figure L-3  Deformations caused by stress acting on one surface of the volume in Figure L-1: (a) linear deformation that results in extension of the side AB in the x direction by an amount BB′; (b) shearing only; (c) rotation only; (d) combined angular deformation (α) and rotation (β). See text for details.

Other types of deformation are caused by shearing (Figure L-3b), rotation (Figure L-3c), and a combination of the two (Figure L-3d). The angular deformations ξ and ζ, in the limit as the volume becomes infinitesimally small, can be expressed in terms of displacement components δu and δw, by using the relations


(L-4a)

and


(L-4b)

which, by adding both sides, yield


(L-5a)

Similarly, we obtain


(L-5b)

and


(L-5c)

The quantities described by equations (L-3) and (L-5) make up the strain tensor


(L-6)

Note that, by way of equations (L-5a,L-5b,L-5c), the strain tensor (L-6) is symmetric.

Also from equations (L-4a,L-4b), we obtain by subtraction


(L-7a)

Similarly,


(L-7b)

and


(L-7c)

Note that by definition

The angular deformations given by equations (L-5a,L-5b,L-5c) are called shear strains since they result in a shearing of the volume around a point within a solid body (Figure L-3b). The quantities given by equations (L-7a,L-7b,L-7c) are associated with rotation without deformation (Figure L-3c).

The displacement tensor in equation (L-2) can be expressed as


(L-8)

The elements of the matrices in equation (L-8) can then be expressed in terms of the normal strain components exx, eyy, and ezz, given by equations (L-3a,L-3b,L-3c), the shear strain components, exz, exy, and eyz, given by equations (L-5a,L-5b,L-5c), and the rigid rotational components θxz, θxy, and θyz, given by equations (L-7a,L-7b,L-7c). All of these components are then substituted into equation (L-2) to get


(L-9)

Next, we need to establish the relation between the stress tensor (L-1) and the strain tensor (L-6). For an elastic solid, Hooke’s law states that the strain at any point is directly proportional to the stresses applied at that point. First, consider a volume (Figure L-1) under tensional principle stresses, Pxx, Pyy, and Pzz. The stress-strain relations are written as


(L-10a)


(L-10b)

and