User:Ageary/Yilmaz 9

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

1 9.0 Introduction
2 9.1 Models with horizontal layers
3 9.2 Model with low-relief structure
4 9.3 Model with complex overburden structure
- 4.1 Image gathers
- 4.2 Constant half-space velocity analysis
5 9.4 Model building
6 9.5 Model updating
7 Exercises
8 Appendix J: Data modeling by inversion
9 See also
10 References
11 External links

9.0 Introduction

Seismic representation of an earth model in depth usually is described by two sets of parameters — layer velocities and reflector geometries. Depth migration is the ultimate tool for delineation of reflector geometries (earth imaging in depth). If layer velocities are determined accurately, reflector geometries can be recovered by iterative depth migration (2-D poststack depth migration). Difficulties in estimating layer velocities with a required level of accuracy make the earth model estimation a challenging task for the geophysicist.

In this chapter, first, we shall discuss practical inversion methods to estimate layer velocities (Sections models with horizontal layers, model with low-relief structure, and model with complex overburden structure). Next, by using appropriate combinations of inversion methods, we shall develop strategies to build an earth model (model building). Finally, we shall use residual moveout analysis and reflection tomography to update an initial model (model updating) and thus obtain a final model that can be used to image the subsurface in depth as discussed in earth imaging in depth. By being equipped with methods and strategies for interval velocity estimation, model building, and updating (earth modeling in depth), and depth migration (earth imaging in depth), we shall then study in structural inversion complete case studies for earth modeling and imaging in depth.

Nearly all of the practical methods of layer velocity estimation are based on ray theory, and more specifically, on inversion of seismic traveltimes. Velocity estimation methods include Dix conversion of rms velocities, inversion of stacking velocities, coherency inversion, and analysis of image gathers from prestack depth migration.

As described in introduction to earth imaging in depth, velocity variations within the earth may be characterized in two ways — structure-dependent and structure-independent. A structure-dependent earth model comprises geological formations with interfaces that coincide with distinct velocity contrasts. We encounter structure-dependent earth models in areas with extensional and compressional tectonics, and especially in areas with salt and overthrust tectonics. A structure-dependent earth model usually requires a layer-by-layer estimation of layer velocities and delineation of reflector geometries that coincide with the layer boundaries themselves.

A structure-independent earth model comprises geologic formations with interfaces that do not necessarily coincide with distinct velocity contrast. Structure-independent earth models often are associated with low-relief structures or stratigraphic plays involving a depositional sequence with facies changes. A structure-independent earth model may be estimated initially by Dix conversion of rms velocities without requiring a layer-by-layer analysis.

The simplest method for estimating layer velocities is Dix conversion of rms velocities ^[1]. The method requires the rms velocities associated with the layer boundaries that are included in the earth model to be constructed. The rms velocities ideally are estimated by prestack time migration (prestack time migration). Alternatively, a smoothly varying form of stacking velocities estimated from dip-moveout corrected data may be a reasonable substitute for rms velocities (dip-moveout correction in practice). Less desirably, stacking velocities themselves with a fair degree of smoothing applied may be used in lieu of rms velocities. The Dix conversion formula (Section J.4) is valid for horizontally layered earth models with constant layer velocities and small offsets. For an earth model with dipping layer boundaries and layer velocities with vertical and lateral variations, more accurate methods are required such as stacking velocity inversion, coherency inversion and image-gather analysis.

Stacking velocity inversion (Thorson et al., 1985) requires time horizons picked from unmigrated CMP-stacked data and stacking velocities at analysis locations. Assume that a velocity-depth model already has been estimated for the first n − 1 layers, and that we want to estimate the layer velocity for the nth layer below a CMP location. For a trial constant velocity assigned to the nth layer, the method involves normal-incidence time-to-depth conversion of the time horizon associated with the base of the nth layer, then modeling of the nonzero-offset traveltimes associated with the CMP reflection event that corresponds to the base of the nth layer and determining the moveout velocity by fitting a hyperbola to the modeled traveltime trajectory. This procedure is repeated for a range of constant trial velocities, and the velocity that yields the minimum discrepancy between the actual stacking velocity and the modeled moveout velocity is assigned to the nth layer below the CMP location where the stacking velocity inversion is being performed.

Coherency inversion ^[2] also requires time horizons picked from unmigrated CMP-stacked data. However, in lieu of stacking velocities as for stacking velocity inversion, coherency inversion requires analyzing CMP gathers themselves. Again, assume that a velocity-depth model already has been estimated for the first n − 1 layers, and that we want to estimate the layer velocity for the nth layer below a CMP location. For a trial constant velocity assigned to the nth layer, coherency inversion involves normal-incidence time-to-depth conversion of the time horizon associated with the base of the nth layer, then modeling of the nonzero-offset traveltimes associated with the CMP reflection event that corresponds to the base of the nth layer, and computing the semblance (velocity analysis) within a CMP data window that follows the modeled traveltime trajectory. This procedure is repeated for a range of constant trial velocities and the velocity that yields the highest semblance value is assigned to the nth layer below the CMP location where the coherency inversion is being performed.

Time horizons used in normal-incidence time-to-depth conversion as part of the stacking velocity inversion and coherency inversion procedures are picked from unmigrated CMP-stacked data. Alternatively, time horizons interpreted from the time-migrated volume of data can be unmigrated to obtain the time horizons equivalent to the time horizons picked from the unmigrated data. We circumvent the picking of prestack reflection traveltimes in coherency inversion by measuring the discrepancy between the modeled and actual traveltimes by way of semblance. Similarly, we avoid the picking of prestack reflection traveltimes in stacking velocity inversion by measuring the discrepancy between the modeled and actual stacking velocities.

Stacking velocity inversion and coherency inversion both take into account vertical velocity gradients which may be available from sonic logs. The methods also honor ray bending at layer boundaries. While Dix conversion assumes a hyperbolic moveout for the reflection event that corresponds to the base of the layer under consideration, both interval velocity estimate from coherency inversion and stacking velocity inversion are based on nonhyperbolic CMP traveltime modeling.

Both stacking velocity inversion and coherency inversion can be considered accurate for velocity-depth models with smoothly varying reflector geometries and lateral velocity variations greater than the effective cable length associated with the layer boundary under consideration. As for conventional stacking velocity estimation (velocity analysis), the accuracy in interval velocity estimation from Dix conversion, stacking velocity inversion, and coherency inversion are all influenced by the reflector depth, magnitude of the velocity, and the cable length. Specifically, the deeper the reflector, the larger the layer velocity above, and the shorter the cable length, the less accurate is the interval velocity estimate.

To estimate, update and verify velocity-depth models for targets beneath complex overburden structures, such as those associated with overthrust and salt tectonics, ultimately, we have to do image-gather analysis ^[3] ^[4] ^[5], ^[6]. An image gather is the output from prestack depth migration and is a true (common depth-point) CDP gather at a surface location. Stacking of image gathers yields an earth image in depth. If the velocity-depth model is correct, then events on an image gather are flat. In this respect, an image gather can be considered like a moveout-corrected CMP gather, except the vertical axis on an image gather is in depth.

An event on an image gather with a moveout indicates an erroneously too low or too high velocity. By examining a panel of image gathers from the same location but with different constant trial velocities for the layer under consideration, one can pick the velocity that yields a flat event and assign it as the velocity of the layer above. Image gathers also can be used to make residual corrections to velocity estimates at analysis locations. This normally is done by first converting the gather to the time domain, performing residual moveout velocity analysis, and converting back to the depth domain. The resulting residual correction should favorably improve the power of the stack obtained from image gathers and yield an updated velocity-depth model. The model updating based on the image-gather analysis usually is repeated until residual moveouts on image gathers are reduced to a minimum ^[7] ^[8].

In the following sections, we shall review the velocity estimation methods described above by using three synthetic model data sets with varying complexity. These are models with horizontal layers, low-relief structure and complex structure. Emphasis will be on practical aspects of the velocity estimation techniques.

Inversion methods for data modeling

Historically, the term seismic inversion often has been used within the context of acoustic impedance estimation from a broad-band time-migrated CMP-stacked data. This narrow meaning commonly is referred to as trace inversion. In practice, however, seismic inversion has a broader scope of applications which can be grouped in two categories — data modeling and earth modeling.

What we do in seismic data processing described in fundamentals of signal processing through 3-D seismic exploration is based largely on data modeling. An observed seismic wavefield can be described in two parts — traveltimes and amplitudes. Seismic amplitudes are more prone to the detrimental effects of noise as compared to traveltimes. Hence, in seismic inversion, we almost always treat traveltimes and amplitudes separately. When modeling the observed data, we either model the traveltimes or amplitudes. When modeling the earth, again, we use the traveltimes, as in structural inversion (structural inversion), or amplitudes, as in stratigraphic inversion (reservoir geophysics).

It is most appropriate to provide a summary list of applications of seismic inversion for data modeling that are spread throughout fundamentals of signal processing to 3-D seismic exploration.

Deconvolution is based on modeling a one-dimensional (1-D) seismogram by optimum Wiener filtering for a minimum-phase estimate of the source wavelet, to predict multiples, and obtain an estimate of white reflectivity series (the convolutional model).
We model traveltime deviations on moveout-corrected CMP gathers to estimate surface-consistent shot and receiver residual statics (residual statics corrections).
We model refracted arrival times to estimate, again, surface-consistent shot and receiver intercept time anomalies, and thus obtain shot and receiver refraction statics (refraction statics corrections).
One type of formulation of the discrete Radon transform is by generalized linear inversion. The discrete Radon transform is used to model a CMP gather so as to attenuate multiples and random noise, while compensating for missing data and finite cable length in recording (the radon transform).
We model the seismic signal represented by reflection events assumed to be linear from trace to trace and attenuate random noise uncorrelated from trace to trace by using spatial prediction filters (linear uncorrelated noise attenuation).
Based on the same data modeling concept, we design spatial prediction filters to perform trace interpolation (processing of 3-D seismic data).
Data modeling also can be used in the design of a three-dimensional (3-D) dip-moveout correction operator which accounts for irregular spatial sampling and undersampling of recorded data (processing of 3-D seismic data).

Most data modeling applications are based on the theory of generalized linear inversion. Although each of the applications listed above is treated in the appendixes of the respective chapters, for completeness, a mathematical summary based on the generalized linear inversion theory is provided in Sections J.1, J.2 and J.3.

Inversion procedures for earth modeling

Practical methods for estimating layer velocities and delineating reflector geometries can be appropriately combined to form inversion procedures to construct earth models in depth from seismic data. Listed in Table 9-1 are four such combinations.

These combinations are ordered from top to bottom with an increasing level of accuracy. Also, for a given combination, the methods for layer velocity estimation and reflector geometry delineation are compatible. This means that, for instance, if you think you can afford the accuracy of stacking velocity inversion to estimate the layer velocities, it should suffice to perform image-ray depth conversion to define the reflector geometries. Nevertheless, in practice, you may wish to choose other combinations of the methods from the left- and right-hand columns. For instance, the combination of Dix conversion with image-ray depth conversion is another inversion procedure that is used widely. Also, you may be compelled to apply an inversion procedure that involves multiple combinations. For instance, in areas where salt tectonics has caused formation of diapiric structures, the earth model may be estimated in three parts — the overburden above the salt diapir, the salt diapir itself, and the substratum. You may then use coherency inversion combined with poststack depth migration to estimate the overburden model, and image-gather analysis combined with prestack depth migration to define the base-salt geometry and estimate the substratum model.

The following are the general guidelines for implementing the procedures listed in Table 9-1. For convenience, we shall refer to these procedures with keywords — vertical stretch, map migration, poststack depth migration, and prestack depth migration. The primary consideration in the choice for an inversion procedure is the degree of lateral velocity variations and the complexity of reflector geometries. A mild-to-moderate lateral velocity variation is associated with a zero-offset diffraction response that is represented by a skewed, but almost hyperbolic traveltime trajectory (Figure 8.0-12). A strong lateral velocity variation is associated with a zero-offset diffraction response that is represented by a distorted, nonhyperbolic traveltime trajectory (Figure 8.0-11). A severe lateral velocity variation is associated with a zero-offset diffraction response that is represented by a complex, multivalued traveltime trajectory (Figure 8.0-13).

Vertical Stretch: A combination of Dix conversion of stacking velocities to estimate layer velocities and vertical-ray time-to-depth conversion of time horizons picked from a time-migrated volume of data to delineate reflector geometries: This is a procedure appropriate for cases with negligible ray bending at layer boundaries, gentle dips, and lateral velocity variations judged to be within the bounds of time migration.
Map Migration: A combination of stacking velocity inversion to estimate layer velocities and image-ray time-to-depth conversion of time horizons picked from a time-migrated volume of data to delineate reflector geometries: This is a procedure appropriate for cases with moderate ray bending at layer boundaries, moderate vertical velocity gradients, and moderate lateral velocity variations.
Poststack Depth Migration: A combination of coherency inversion to estimate layer velocities and poststack depth migration to delineate reflector geometries: This is a procedure appropriate for cases with significant ray bending at layer boundaries and significant vertical velocity gradients, and strong lateral velocity variations with sharp changes in reflector curvatures.
Prestack Depth Migration: A combination of image-gather analysis to estimate and update layer velocities, and stacking of image gathers to delineate reflector geometries: This is a procedure appropriate for cases with significant ray bending at layer boundaries, and severe lateral velocity variations associated with salt and overthrust tectonics.

**Table 9-1.** A set of inversion procedures for earth modeling in depth to estimate layer velocities and delineate reflector geometries.
Layer Velocities	Reflector Geometries
Dix conversion of rms velocities	vertical-ray time-to-depth conversion (vertical stretch)
stacking velocity inversion	image-ray time-to-depth conversion (map migration)
coherency inversion	poststack depth migration
image-gather analysis	prestack depth migration

The inversion methods listed in Table 9-1 are used to estimate an initial earth model in depth. Seismic inversion also is used to update the estimated model (model updating). A common application of inversion to estimate the errors in the initial model parameters — layer velocities and reflector depths, is reflection traveltime tomography. Tomographic inversion involves perturbing the model parameters by a small amount so as to match the modeled reflection traveltimes with the observed traveltimes. Refraction traveltime tomography (Section C.9) and reflection traveltime tomography (Section J.6) both are based on the assumption that the perturbation required to update the model parameters is very small compared to the spatial variations in the model parameters themselves. In practice, tomography is best used strictly to touch-up a carefully estimated earth model based on some plausable geologic constraints; it should never be used by itself to estimate the model.

Velocity-depth ambiguity

A fundamental problem with inversion is velocity-depth ambiguity. This means that an error in layer velocity can be indistinguishable from an error in reflector geometry. To help resolve the velocity-depth ambiguity, we must use nonzero-offset data to estimate the layer velocities and reflector geometries (Section J.7). We can never claim to completely resolve the velocity-depth ambiguity and obtain the one and the only solution — the true velocity-depth model, from inversion of seismic data. But, we can hope to reduce the many possible solutions to a few geologically plausable velocity-depth models. Of course, if there is ample well control, we may be able to discard all but one of the few solutions and declare that as the solution.

Consider the four geologically plausable velocity-depth models of a salt diapir shown in Figure 9.0-1. Compare these models and note the differences in reflector geometries, which are subtle in some places and distinctively obvious in others. There also are differences in layer velocities from one model to the other.

Compute the zero-offset reflection traveltime responses by way of normal-incidence rays and display them in the form of zero-offset sections as shown in Figure 9.0-2. While there are differences between the velocity-depth models in Figure 9.0-1, differences in the zero-offset traveltime sections are less than marginal. Any differences that we observe in the traveltime sections may be attributed largely to the digitizing errors that have occurred in building the velocity-depth models. We may, within the bounds of ray theory associated with reflectors with reasonable curvature and shape, retain one of the four zero-offset traveltime sections and discard the other three sections shown in Figure 9.0-2. We may also conclude that the zero-offset traveltime section we have kept is consistent with all four of the velocity-depth models shown in Figure 9.0-1.

In reality, we record wavefields that comprise not just traveltimes but also amplitudes. Pretend that the zero-offset traveltime section which we retained in Figure 9.0-2 is equivalent to a zero-offset wavefield section as in Figure 8.2-1, insofar as traveltimes are concerned.

Now, consider the inverse problem. Suppose that the analyst is given the zero-offset section in Figure 8.2-1 and is asked to estimate an earth model in depth from it. As an answer, the analyst may provide any one of the earth models shown in Figure 9.0-1. They all are equally valid answers since they all are consistent with the input zero-offset section in Figure 8.2-1. In fact, there are infinitely many possible earth models, not just the four models in Figure 9.0-1, which are consistent with the zero-offset section. While there are infinite number of solutions, there should only be one valid answer — the geologic model shown in Figure 8.2-1. We will never know which one, if any, of the answers corresponds to the true velocity-depth model.

In practice, the analyst deals with a stacked section and not a zero-offset section. Shown in Figure 9.0-3 are three velocity-depth models estimated from the seismic data associated with the stacked section in the same figure by means of model estimation procedures described in this chapter. The zero-offset traveltimes associated with the layer boundaries included in each of the three models are computed by normal-incidence ray tracing. When modeled traveltimes are repeatedly overlayed on the same stacked section, we note that they all coincide with the observed traveltimes associated with the reflection events that correspond to the layer boundaries included in the velocity-depth models. Which one of the three models best represents the true subsurface structure — that is the question. And this question has a profound implication in practice; in the present case, an accurate delineation of the geometry for the target reflector represented by the deepest layer boundary has an impact on the exploration and development of a potential reservoir.

Mathematically, the nature of any inversion problem with a multiple number of solutions that are all consistent with observed data is characterized as nonuniqueness. This mathematical nonuniqueness is known to the exploration seismologist as velocity-depth ambiguity. Refer back to Figure 9.0-1 and note that in reality only Model 1 is the correct model — the same as the true geologic model shown in Figure 8.2-1. However, any combination of layer velocities and reflector geometries represented by the other three earth models is equally consistent with the input data. This means that errors in layer velocities and errors in reflector geometries are indistinguishable — a restatement of the velocity-depth ambiguity. A robust quantification of velocity-depth ambiguity is provided in Section J.7.

Figure 9.0-1 Four different earth models in depth associated with a diapiric structure. Note the differences in layer velocities and reflector geometries.
Figure 9.0-2 The zero-offset traveltime sections associated with the earth models in depth shown in Figure 9.0-1.

Since we can never obtain the true representation of an earth model from inversion of seismic data, the plausable strategy is to estimate an initial model and then update it to get a final model that may be considered an acceptable approximation to the true model. Model building strategies are discussed in model building, and model updating by reflection tomography is presented in model updating.

The important question in practice is, following the model update, how much velocity-depth ambiguity remains unresolved in a final velocity-depth model. As illustrated in Figures 9.0-1 and 9.0-2, the ambiguity with zero-offset data is infinite. By using nonzero-offset data, we can hope to resolve the velocity-depth ambiguity up to a certain theoretical limit. An important rule to keep in mind is that, for data with good quality, velocity-depth ambiguity for a reflector can be resolved with an acceptable degree of accuracy if the data used in inversion have been recorded with offsets greater than the reflector depth ^[9] ^[10].

While an earth model can only be estimated with an accuracy that is within the threshold of velocity-depth ambiguity, it does have to be consistent with the input data used in inversion to estimate the model. Consistency is a necessary condition for an earth model to be certified as an acceptable estimate of the true model. A quick way to check for consistency is by ray-theoretical forward modeling of zero-offset traveltimes associated with the reflector boundaries that are in the earth model itself, and then comparing them with the actual traveltimes picked from the stacked data. Any discrepancy between the modeled and actual traveltimes is an indication of errors in the earth model parameters — layer velocities and/or reflector geometries. By using nonzero-offset data, we can hope to reduce the many possible solutions to a few. Furthermore, by introducing constraints, we may be able to converge to a single solution provided the set of constraints are reliable. One set of constraints is the depth information at well locations. This well-top information can be used to calibrate results of inversion of surface seismic data and obtain a single earth model in depth that not only is consistent with the surface data set itself but also with the borehole data.

Aside from consistency, model verification has to include a test of flatness of events on image gathers derived from prestack depth migration. A correct model, again, within the limitations of velocity-depth ambiguity and accuracy of inversion methods, would yield an accurate image from prestack depth migration irrespective of the source-receiver offset. Thus, with the correct model, the resulting image gathers would have flat events. An erroneous earth model, on the other hand, would cause residual moveout on image gathers. The use of image gathers in model estimation, verification, and update is discussed in model with complex overburden structure and illustrated by some of the case studies presented in structural inversion.

Model representation and visualization

In earth modeling, a surface corresponding to a layer boundary is usually represented by a set of triangles, the size and shape of which vary depending on the complexity of the reflector geometry. Shown in Figure 9.0-4 are four surfaces associated with layer boundaries included in an earth model. Also included in Figure 9.0-4 are the representations of these surfaces by triangulation.

A velocity-depth model usually is represented either in the form of a gridded or tessellated volume. Gridding means dividing the whole volume into a set of 3-D cells of equal size with appropriate dimensions in the inline, crossline, and depth directions. Tessellation means dividing the volume associated with each layer into a set of tetrahedra, the size and shape of which depend on the geometry of layer boundaries (Figure 9.0-5). In a tessellated velocity-depth model, a velocity and a gradient, if available from sonic logs, are assigned to each corner of the tetrahedra.

A velocity-depth model is represented either in gridded or tessellated form depending on the application that needs it as input. For instance, ray tracing in coherency inversion and prestack depth migration may be performed using gridded or tessellated models. On the other hand, wave extrapolation in 3-D poststack depth migration based on finite-difference schemes is performed conveniently using gridded models.

Images of layer boundaries included in an earth model in depth can be converted to a physical model using various image construction techniques. Figure 9.0-6 shows the four surfaces that represent the layer boundaries from the earth model as in Figure 9.0-4 and the slabs that represent the layers of the earth model created by laser lithography.

In areas with complex structures, we often have to deal with mulitvalued depth surfaces. For instance, salt overhangs associated with diapirism and imbricate structures associated with overthrusting cause a surface to fold onto itself. Shown in Figure 9.0-7 is an earth model that comprises a complex diapiric structure. The top-salt boundary is represented by the multivalued yellow surface and the base-salt boundary is represented by the pink surface. The surfaces that represent the layer boundaries within the overburden must be attached to the top-salt boundary surface without any gaps so as to form individual volumes associated with each layer. Figure 9.0-8a shows the solid-model representation of the earth model in Figure 9.0-7. When expanded, the solid model shows the individual volumes associated with the layers included in the earth model (Figure 9.0-8b).

Figure 9.0-3 Three velocity-depth models (left column) all of which are consistent with the observed traveltimes on the stacked section (right column). Superimposed on the same stacked section are the modeled zero-offset traveltimes associated with the velocity-depth models in the left column.
Figure 9.0-4 Surfaces that represent the reflector geometries (left column) and the triangulated form of these surfaces (right column).
Figure 9.0-5 (a) The model volume that includes the surfaces in Figure 9.0-4, (b) the tessellated form of the model.
Figure 9.0-6 Physical modeling of an earth model — surfaces that represent layer boundaries (left) and slabs that represent the layers themselves (right). The physical model dimensions have an aspect ratio of one.

Isolate the salt mass from the model as shown in Figure 9.0-9 to examine the complexity of the salt diapir. The multivalued surface that represents the top-salt boundary is shown in Figure 9.0-10a. The triangulated mesh for this surface is shown in Figure 9.0-10b with an enlarged view in Figure 9.0-11.

9.1 Models with horizontal layers

In this section, we shall examine the accuracy of Dix conversion and coherency inversion to estimate layer velocities using two earth models with horizontal layers, but with lateral velocity variations. The near-surface layer is of constant velocity in the first model and its velocity varies laterally in the second model. Otherwise, both models are identical. Shown in Figure 9.1-1a are the velocity profiles for the six layers in the model. We shall refer to the layers by the horizon names corresponding to the base of each layer, H1 through H6. Listed in Table 9-2 are the layer velocities and depths to the base of each layer. When the layer velocity is not constant, the range is given in Table 9-1. The lateral velocity gradients in layers H3 and H4 are about 125 m/s and 200 m/s over one cable length, respectively. Figure 9.1-1b shows the velocity-depth model with the color bar on the right-hand margin.

A total of 384 shot records was modeled using the two-way acoustic wave equation. The simulated recording geometry consists of an off-end cable with 96 receivers and offset range 25-2400 m. Shot and receiver intervals are both 25 m, and the CMP interval is 12.5 m and CMP fold is 48. Figure 9.1-2a shows the CMP-stacked section with the picked time horizons that correspond to the layer boundaries H1 through H6 in Figure 9.1-1b. Compare with the velocity-depth model (Figure 9.1-1b) and note that flat horizons in depth correspond to curved horizons in time because of the lateral velocity variations (Figure 9.1-1a). The stacking velocity section is shown in Figure 9.1-2b with the color bar on the right-hand margin. The stacking velocity section was derived from the horizon-consistent stacking velocity profiles shown in Figure 9.1-3.

**Table 9-2.** Parameters of the model with horizontal layers and constant-velocity near-surface layer.
Layer	Velocity (m/s)	Depth (m)
H1	1500	100
H2	2000	1000
H3	2400 – 2700	1500
H4	3000 – 3500	1800
H5	4500	2250
H6	3000	2700

Dix conversion

The horizon-consistent stacking velocity profiles (Figure 9.1-3) at each of the layer boundaries are used to perform Dix conversion to derive the interval velocity profiles for each of the layers. Dix conversion is based on the formula

v_{n}={\sqrt {\frac {V_{n}^{2}\tau _{n}-V_{n-1}^{2}\tau _{n-1}}{\tau _{n}-\tau _{n-1}}}},

(1)

where v_n is the interval velocity within the layer bounded by the (n − 1)st layer boundary above and the nth layer boundary below, τ_n and τ_n−1 are the corresponding two-way zero-offset times, and V_n and V_n−1 are the corresponding rms velocities. Derivation of equation (1) is provided in Section J.4.

Equation (1) is based on the assumptions that the layer boundaries are flat and the offset range used in estimating the rms velocities V_n and V_n−1 corresponds to a small spread. Additionally, keep in mind that the rms velocities used in equation (1) are based on a straight-ray assumption; thus, ray bending at layer boundaries are not accounted for in Dix conversion.

The procedure for estimating the layer velocities and reflector depths using Dix conversion of stacking velocities includes the following steps:

For each of the layers in the model, pick the time of horizon on the unmigrated CMP-stacked data that corresponds to the base-layer boundary (Figure 9.1-2a). These times are used in lieu of the two-way zero-offset times in equation (1).
Extract the rms velocities at horizon times (Figure 9.1-3).
Use equation (1) to compute the interval velocities for each of the layers from the known quantities — rms velocities and times at top- and base-layer boundaries.
Use interval velocities and times at layer boundaries to compute depths at layer boundaries. If the input times are from an unmigrated stacked section as in Figure 9.1-2a, use normal-incidence rays for depth conversion. If the input times are from a migrated stacked section, use image rays for depth conversion.

Figure 9.0-7 Two different views of a complex earth model associated with a salt diapirism. (Modeling by Cyril Gregory; courtesy Paradigm Geophysical.)
Figure 9.0-8 (a) The solid-model representation of the earth model in Figure 9.0-7, (b) the exploded form of (a). (Modeling by Cyril Gregory; courtesy Paradigm Geophysical.)
Figure 9.0-9 Two different views of the salt mass (the pink solid) and the substratum isolated from the solid model shown in Figure 9.0-8. (Modeling by Cyril Gregory; courtesy Paradigm Geophysical.)
Figure 9.0-10 (a) The surface that represents the top-salt boundary, and (b) its representation by triangulation. (Modeling by Cyril Gregory; courtesy Paradigm Geophysical.)
Figure 9.0-11 A detailed portion of the triangulated surface shown in Figure 9.0-10b. The surface underneath the triangulated mesh is the base-salt boundary. (Modeling by Cyril Gregory.)

Interval velocity profiles derived from Dix conversion are shown in Figure 9.1-4a. The earth model can be constructed by combining the estimated interval velocity profiles and depth horizons (Figure 9.1-4b). Comparison with the true model shown in Figure 9.1-1b clearly demonstrates that the interval velocity estimation based on Dix conversion is not completely accurate. The interval velocity profiles derived from Dix conversion (Figure 9.1-4a) exhibit the sinusoidal oscillations caused by the swings in the stacking velocity profiles themselves (Figure 9.1-3).

The fundamental problem is that the stacking velocity estimation is based on fitting a hyperbola to CMP traveltimes associated with a laterally homogeneous earth model. If there are lateral velocity variations in layers above the layer under consideration, and if these variations are within a cable length, then stacking velocities would oscillate in a physically implausable manner ^[11] ^[12] ^[13]. As a consequence, the resulting interval velocity estimation based on Dix conversion is adversely affected. In the present case, Dix conversion has produced fairly accurate estimates for the interval velocities of the top three layers — H1, H2, and H3 as shown in Figure 9.1-4a. But the interval velocity estimates for layers H4 and H5 have been adversely affected by the laterally varying velocities within the layer above, H3.

The pragmatic approach would be to smooth out the oscillations in the stacking velocities before Dix conversion and smooth out the oscillations in the velocity profiles after Dix conversion (Figure 9.1-5a). Then, the resulting earth model is expected to be free of the adverse effects of stacking velocity anomalies (Figure 9.1-5b). A closer look at the central portions of the estimated models using Dix conversion is shown in Figure 9.1-6. Note that the model derived from the smoothed interval velocities is closer to the true model. We shall make an attempt in model updating to update this result by using tomography.

Coherency inversion

The procedure for estimating layer velocities from coherency inversion requires CMP gathers at analysis locations and horizon times picked from unmigrated stacked data. Alternatively, time horizons picked from time-migrated data can be forward-modeled to derive the zero-offset traveltimes needed for coherency inversion. As for Dix conversion, the velocity estimate from coherency inversion is local, independent of data away from the analysis location.

A procedure for velocity-depth model estimation that includes coherency inversion is conducted layer-by-layer starting from the surface. As for Dix conversion, consider the synthetic data set associated with the model with horizontal layers shown in Figure 9.1-1. We shall adopt the interval velocity profile for the first layer H1 estimated from Dix conversion and start the application of coherency inversion with layer H2. Assume that the velocity-depth model for the first n − 1 layers already has been estimated. For the nth layer, follow the steps below for coherency inversion:

For a trial constant velocity assigned to the nth layer, perform normal-incidence traveltime inversion to convert the time horizon corresponding to the base-layer boundary to a trial depth horizon.
Given the geometry of the CMP gather at the analysis location, assign a trial velocity to the half space that includes the layer yet to be determined and compute the CMP traveltimes using the known overburden velocity-depth model. The modeled CMP traveltime trajectory that corresponds to the base of the layer under consideration is in general nonhyperbolic, for the ray tracing used to compute the CMP traveltimes accounts for ray bending at layer boundaries and incorporates vertical velocity gradients within layers above. Shown in Figures 9.1-7a through 9.1-11a are the CMP ray-paths at an analysis location through the unknown and half space that includes the layer under consideration and the known overburden part of the model. The corresponding CMP traveltime trajectories are shown in color in Figures 9.1-7b through 9.1-11b overlayed on the CMP gather at the analysis location.
Extract data window along the modeled traveltime trajectory. Shown in Figures 9.1-7c through 9.1-11c are the data windows for seven trial velocities around the optimum velocity that best flattens the event within the data window.
Compute semblance using the data within the selected window as a measure of discrepancy between the modeled and the actual traveltimes. Figure 9.1-12 shows the semblance spectra computed from the CMP gather in Figures 9.1-7b through 9.1-11b from the reflection events associated with layers H2 through H6.
Repeat all the steps above for a range of constant velocities.
Pick the constant trial velocity as the layer velocity for which the semblance is the maximum.

Figure 9.1-1 An earth model that comprises six flat layers: (a) the interval velocity profiles for the six horizons H1-H6; (b) true velocity-depth model created from the profiles in (a).
Figure 9.1-2 (a) The CMP-stacked section derived from the modeled common-shot gathers using the earth model shown in Figure 9.1-1b; (b) the stacking velocity section.
Figure 9.1-3 Horizon-consistent stacking velocity semblance spectra computed from the CMP gathers of the synthetic data as in Figure 9.1-2a along the time horizons H1-H6.
Figure 9.1-4 (a) The interval velocity profiles derived from Dix conversion of the horizon-consistent stacking velocity profiles picked from the semblance spectra shown in Figure 9.1-3; (b) estimated velocity-depth model. Compare with the true velocity-depth model shown in Figure 9.1-1b.
Figure 9.1-5 (a) The interval velocity profiles as in Figure 9.1-4a displayed by the thick curves and smoothed interval velocity profiles displayed by the thin curves; (b) estimated velocity-depth model using the smoothed interval velocity profiles in (a). Compare with the true velocity-depth model shown in Figure 9.1-1b and the model derived from the unsmoothed interval velocities shown in Figure 9.1-4b.
Figure 9.1-6 Central portions of (a) the true velocity-depth model shown in Figure 9.1-1b, (b) the velocity-depth model shown in Figure 9.1-4b estimated using the unsmoothed interval velocity profiles shown in Figure 9.1-4a, and (c) the velocity-depth model shown in Figure 9.1-5b estimated using the smoothed interval velocity profiles shown in Figure 9.1-5a.

The results of coherency inversion are used to pick velocity nodes at analysis locations. Specifically, the layer velocity at a given location is selected based on the semblance curve and the data window along the modeled traveltime trajectory making sure that the estimated velocities are geologically plausable. Data windows along modeled traveltime trajectories can be examined for flatness criterion to pick an optimum velocity node. Specifically, the events in Figures 9.1-7c through 9.1-11c can be considered as events on a moveout-corrected CMP gather. A flat event on a moveout-corrected CMP gather suggests that the stacking velocity associated with that event is optimum. Likewise, a flat event in the data window panels in Figures 9.1-7c through 9.1-11c suggests that the layer velocity from coherency inversion is optimum. An erroneously too low or too high velocity causes residual moveout which can be observed on the event within the data window.

The semblance spectra in Figure 9.1-12 were computed using three different maximum offsets — 2400 m, 1400 m, and 400 m. Note that, for a given layer, the longer the cable length, the sharper the semblance peak — the higher the velocity resolution. Also, the deeper the event and the higher the velocity, the broader the semblance peak — the poorer the velocity resolution.

Figures 9.1-7 through 9.1-11 show the results of coherency inversion at one CMP location. In practice, for 2-D data, coherency inversion often is applied continuously along the line. As for horizon-consistent stacking velocity analysis, for each layer, a horizon-consistent semblance spectrum is computed using coherency inversion. For 3-D data, as for conventional velocity analysis, coherency inversion normally is applied at uniformly spaced grid points over the survey area.

Figure 9.1-7 Coherency inversion for layer H2 as in Figure 9.1-2a: (a) CMP raypaths, (b) the CMP gather at the analysis location, (c) the data window that includes the color-marked event in the CMP gather in (b) with moveouts associated with the different trial velocities assigned to layer H2. See text for details.
Figure 9.1-8 Coherency inversion for layer H3 as in Figure 9.1-2a: (a) CMP raypaths, (b) the CMP gather at the analysis location, (c) the data window that includes the color-marked event in the CMP gather in (b) with moveouts associated with the different trial velocities assigned to layer H3. See text for details.
Figure 9.1-9 Coherency inversion for layer H4 as in Figure 9.1-2a: (a) CMP raypaths, (b) the CMP gather at the analysis location, (c) the data window that includes the color-marked event in the CMP gather in (b) with moveouts associated with the different trial velocities assigned to layer H4. See text for details.
Figure 9.1-10 Coherency inversion for layer H5 as in Figure 9.1-2a: (a) CMP raypaths, (b) the CMP gather at the analysis location, (c) the data window that includes the color-marked event in the CMP gather in (b) with moveouts associated with the different trial velocities assigned to layer H5. See text for details.
Figure 9.1-11 Coherency inversion for layer H6 as in Figure 9.1-2a: (a) CMP raypaths, (b) the CMP gather at the analysis location, (c) the data window that includes the color-marked event in the CMP gather in (b) with moveouts associated with the different trial velocities assigned to layer H6. See text for details.
Figure 9.1-12 Semblance spectra derived from coherency inversion applied to the central CMP gather associated with the stacked data in Figure 9.1-1a using maximum offsets of 2400 m, 1400 m, and 400 m. See text for details.

Figure 9.1-13 shows the semblance spectra derived from coherency inversion for layers H2 through H6 of the model shown in Figure 9.1-1. The interval velocity profiles shown in Figure 9.1-14a have been picked by tracking the semblance peaks. Compare with the profiles associated with the true model shown in Figure 9.1-1a, and note that by coherency inversion the interval velocities for layers H2 and H3 have been estimated accurately. But the semblance spectra in Figure 9.1-13 for the underlying layers H4, H5, and H6 exhibit oscillations akin to those associated with the interval velocities derived from Dix conversion of stacking velocities (Figure 9.1-4). Specifically, when there are lateral velocity variations within a cable length in one layer, in this case layer H3, the interval velocity estimation for the layers below is adversely affected as shown in Figure 9.1-13. This phenomenon also was observed by Sorin ^[14].

Shown in Figure 9.1-14b is the velocity-depth model derived from coherency inversion. The model was constructed layer-by-layer starting from the top. Coherency inversion was used to estimate the layer velocities (Figure 9.1-14a), and normal-incidence traveltime inversion of the time horizons picked from the stacked section (Figure 9.1-1) was used to obtain the depth horizons. Note the distortions caused by the oscillations of the interval velocity profiles in the estimated velocity-depth model. As in Dix conversion, these oscillations must be smoothed out for one layer before estimating the interval velocity for the next layer. The resulting semblance spectra exhibit a reduced degree of oscillations (Figure 9.1-15). The velocity-depth model based on the revised interval velocity profiles (Figure 9.1-16a), which have been picked from the new set of semblance spectra, is shown in Figure 9.1-16b. A closer look at the central portion of the estimated model using coherency inversion is shown in Figure 9.1-17. Provided that the oscillations shorter than a cable length are eliminated from the interval velocity profiles, coherency inversion seems to produce an acceptable estimate of the layer velocities.

Near-surface layer with laterally varying velocities

We now examine the accuracy of Dix conversion and coherency inversion for the model with horizontal layers (Table 9-2) but with a near-surface layer H1 with laterally varying velocities between 800 m/s and 1500 m/s. The interval velocity profiles are shown in Figure 9.1-18a and the velocity-depth model is shown in Figure 9.1-18b.

As for the model with a constant-velocity near-surface layer (Table 9-2), a total of 384 shot records was modeled using the two-way acoustic wave equation with identical line geometry. Figure 9.1-19a shows the CMP-stacked section with the picked time horizons that correspond to the layer boundaries H1 through H6 in Figure 9.1-18b. The stacking velocity section is shown in Figure 9.1-19b with the color bar on the right-hand margin.

The stacking velocity section was derived from the horizon-consistent stacking velocity profiles shown in Figure 9.1-20. As a direct consequence of the lateral velocity variations in the near-surface layer H1, the stacking velocities oscillate violently for the layers below. If Dix conversion is done using the unedited stacking velocity profiles, the resulting interval velocity profiles exhibit geologically implausable variations (Figure 9.1-21a). It is imperative in practice to remove the oscillations from the stacking velocity profiles before Dix conversion. The resulting interval velocity profiles shown in Figure 9.1-21b are closer to the true profiles shown in Figure 9.1-18a.

By using the interval velocity profiles, convert the time horizons shown in Figure 9.1-19a to depth horizons. Then, combine the interval velocity profiles with the depth horizons to compose the velocity-depth model shown in Figure 9.1-21c. A comparison with the true model shown in Figure 9.1-19b, once again, clearly demonstrates that the interval velocity estimation based on Dix conversion is not completely accurate. We shall make an attempt in model updating to update the model in Figure 9.1-21c by using tomography.

Figure 9.1-23 shows the semblance spectra derived from coherency inversion for layers H2 through H6 of the model shown in Figure 9.1-18. The analysis was conducted layer by layer starting from the top. Based on the lessons learned from the model experiments shown in Figures 9.1-13 and 9.1-15, the oscillations in the semblance spectra were rejected while tracking the interval velocity profile from the semblance spectrum for the layer under consideration before moving down to the next layer. Despite the rejection of the oscillations, the semblance spectra still exhibit the influence of the near-surface layer with lateral velocity variations (Figure 9.1-23).

Shown in Figure 9.1-24b is the velocity-depth model using the interval velocity profiles (Figure 9.1-24a) derived from coherency inversion. A closer look at the central portion of the estimated model using coherency inversion is shown in Figure 9.1-25. As long as the oscillations shorter than a cable length are eliminated from the interval velocity profiles, coherency inversion seems to produce a better estimate of the velocity-depth model compared to Dix conversion (Figure 9.1-22). The difference between the two estimates is in the reflector geometries. Dix conversion has introduced spurious structures into the model (Figure 9.1-22), while coherency inversion has introduced a bulk shift in the reflector depths (Figure 9.1-25). In earth modeling, an error in the form of a distorted reflector geometry is worse than an error in the form of a bulk shift in the reflector depth. While the error in the form of a bulk shift can be corrected for by calibrating the estimated model to well tops, the error in the form of a distorted reflector geometry may require a serious revision of the estimated model.

Figure 9.1-13 Horizon-consistent coherency inversion semblance spectra computed from the CMP gathers associated with the stacked data in Figure 9.1-2a. No smoothing has been applied to the velocity profiles picked from these semblance spectra.
Figure 9.1-14 (a) The interval velocity profiles derived from coherency inversion semblance spectra shown in Figure 9.1-13; (b) estimated velocity-depth model. Compare with the true velocity-depth model shown in Figure 9.1-1b.
Figure 9.1-15 Horizon-consistent coherency inversion semblance spectra computed from the CMP gathers associated with the stacked data in Figure 9.1-2a. When computing the semblance spectrum for a given layer, smoothing has been applied to the velocity profiles picked from the semblance spectra associated with the layers above.
Figure 9.1-16 (a) The interval velocity profiles derived from coherency inversion semblance spectra shown in Figure 9.1-15; (b) estimated velocity-depth model. Compare with the true velocity-depth model shown in Figure 9.1-1b.
Figure 9.1-17 Central portions of (a) the true velocity-depth model shown in Figure 9.1-1b, and (b) the velocity-depth model shown in Figure 9.1-14b estimated using the smoothed interval velocity profiles shown in Figure 9.1-16a.

9.2 Model with low-relief structure

In this section, we shall review velocity estimation via stacking velocity inversion and coherency inversion in the presence of low-relief structures. Figure 9.2-1 shows a velocity-depth model with such characteristics. The model simulates a transgressive depositional sequence within the first 1-km depth, a deltaic sequence between 1.5-2 km, and a deeper depositional sequence between 2-2.5 km. Our goal is to detect the subtle lateral velocity variations within the individual sequences.

A total of 154 shot records was modeled using the two-way acoustic wave equation. The simulated recording geometry consists of a split-spread cable with 97 receivers and an offset range of 0-2350 m. Shot and receiver intervals are both 50 m. Figure 9.2-2 shows the CMP-stacked section with and without the interpreted time horizons. Assuming that the CMP-stacked section is largely equivalent to a zero-offset section, these time horizons correspond to two-way zero-offset time picks which are used in stacking velocity inversion and coherency inversion.

The color-coded true velocity-depth model is shown in Figure 9.2-3. Using coherency inversion and stacking velocity inversion to estimate layer velocities, the velocity-depth models also shown in Figure 9.2-3 are obtained. In both cases, normal-incidence traveltime inversion is used to delineate the reflector geometries. Note that both techniques are able to delineate the shallow transgressive sequence. The individual units within the deltaic sequence are not included in the velocity-depth model estimation. Instead, only the top (Horizon 5) and base (Horizon 6) of the sequence are included in the model. Nevertheless, both techniques are able to detect the existence of velocity variations from one unit to the next. The deltaic sequence, however, is delineated by coherency inversion more accurately. Finally, stacking velocity inversion fails to estimate the internal velocity distribution of the deepest sequence between 2-2.5 km, correctly. Coherency inversion, on the other hand, has at least been able to detect the relative magnitude of the velocity variations within this sequence with reasonable accuracy. In the following paragraphs, we shall discuss how the two velocity-depth model estimates in Figure 9.2-3 were made.

Figure 9.1-18 An earth model that comprises six flat layers: (a) the interval velocity profiles for the six horizons H1-H6; (b) true velocity-depth model created from the profiles in (a). This model is the same as the model in Figure 9.1-1 except that H1 has laterally varying velocities.
Figure 9.1-19 (a) The CMP-stacked section associated with the earth model shown in Figure 9.1-18; (b) stacking velocity section.
Figure 9.1-20 Horizon-consistent stacking velocity semblance spectra computed from the CMP gathers of the synthetic data as in Figure 9.1-19a.
Figure 9.1-21 (a) The interval velocity profiles derived from Dix conversion of the horizon-consistent stacking velocity profiles picked from the semblance spectra shown in Figure 9.1-20, (b) the interval velocity profiles of (a) after lateral smoothing, (c) estimated velocity-depth model. Compare with the true velocity-depth model shown in Figure 9.1-18b.
Figure 9.1-22 Central portions of (a) the true velocity-depth model shown in Figure 9.1-18b, (b) the velocity-depth model shown in Figure 9.1-21c estimated using the smoothed interval velocity profiles shown in Figure 9.1-21b.

Stacking velocity inversion

The procedure for estimating layer velocities from stacking velocity inversion requires the horizon times picked from unmigrated stacked data (Figure 9.2-2) and stacking velocities at time horizons that correspond to layer boundaries in the model (Figure 9.2-4). Horizon-consistent stacking velocities are estimated by computing semblance for a range of constant velocities continuously along the time horizon picked from the stacked data (Figure 9.2-2). The velocity spectrum for each time horizon (Figure 9.2-4) then is picked to derive a stacking velocity curve along the midpoint axis. As for coherency inversion, the velocity estimate from stacking velocity inversion is local, independent of data away from the analysis location.

A procedure for velocity-depth model estimation that includes stacking velocity estimation is conducted layer-by-layer starting from the surface. Assume that the velocity-depth model for the first n − 1 layers already has been estimated. For the nth layer, follow the steps below for stacking velocity inversion:

For a trial constant velocity assigned to the nth layer, perform normal-incidence traveltime inversion to convert the time horizon corresponding to the base-layer boundary to a trial depth horizon.
Given the geometry of the CMP gather at the analysis location (not the CMP gather itself, but only the source-receiver geometry associated with it), compute the CMP traveltimes. The modeled CMP traveltime trajectory that corresponds to the base of the layer under consideration is in general nonhyperbolic, because the ray tracing used to compute the CMP traveltimes accounts for ray bending at layer boundaries and incorporates vertical velocity gradients within layers above.
Compute the best-fit hyperbolic traveltime trajectory, and thus determine the modeled stacking velocity for the trial interval velocity.
Measure the discrepancy between the modeled and the actual stacking velocities by way of semblance.
Repeat all the steps above for a range of constant velocities.
Pick the constant trial velocity as the layer velocity for which the difference between the modeled and the actual stacking velocities is minimum or the semblance is maximum.

Figure 9.1-23 Horizon-consistent coherency inversion semblance spectra computed from the CMP gathers associated with the stacked data in Figure 9.1-19a. When computing the semblance spectrum for a given layer, smoothing has been applied to the velocity profiles picked from the semblance spectra associated with the layers above.
Figure 9.1-24 (a) The interval velocity profiles derived from coherency inversion semblance spectra shown in Figure 9.1-23; (b) estimated velocity-depth model. Compare with the true velocity-depth model shown in Figure 9.1-18b.
Figure 9.1-25 Central portions of (a) the true velocity-depth model shown in Figure 9.1-18b, and (b) the velocity-depth model shown in Figure 9.1-24b estimated using the smoothed interval velocity profiles shown in Figure 9.1-24a.

Following the steps described above, perform stacking velocity inversion to obtain the layer-velocity semblance spectra shown in Figure 9.2-5. Note that velocity estimates down to Horizon 3 are fairly accurate. The estimates for Horizons 4 and 5 show some departure from the true velocity but are largely acceptable. Note the variations in the estimate for Horizon 6 — the method is attempting to distinguish the units within the deltaic sequence with different velocities. The velocity estimate for Horizon 7a has significant departures from the true velocity of 3500 m/s. Finally, the method has failed to estimate the internal velocity variations within the deep sequence, correctly.

Figure 9.2-1 A velocity-depth model with low-relief structures.
Figure 9.2-2 The CMP-stacked section associated with the velocity-depth model in Figure 9.2-1 with the interpreted time horizons.
Figure 9.2-3 (a) True velocity-depth model associated with the stacked section in Figure 9.2-2; (b) result of stacking velocity inversion; (c) result of coherency inversion.
Figure 9.2-4 Horizon-consistent stacking velocity spectra for the model in Figure 9.2-1.
Figure 9.2-5 Stacking velocity inversion for the model in Figure 9.2-1. Input to inversion are horizon-consistent stacking velocities (Figure 9.2-4) and the horizon times picked from stacked data (Figure 9.2-2). Posted velocities are the true layer velocities as in Figure 9.2-1.
Figure 9.2-6 Coherency inversion for the model in Figure 9.2-1. Picked velocity profiles are overlayed on the semblance spectra. Input to inversion are CMP gathers and the horizon times picked from stacked data (Figure 9.2-2). Posted velocities are the true layer velocities as in Figure 9.2-1.
Figure 9.2-7 Coherency inversion for the model in Figure 9.2-1. (a) CMP gather at analysis location CMP 75 as in Figure 9.2-8; (b) semblance curve; (c) data window for three trial velocities. CMP raypaths that correspond to these trial velocities are shown in Figure 9.2-8.
Figure 9.2-8 Left column: Normal-incidence rays associated with Horizon 7a as in Figure 9.2-2 using three trial velocities; right column: CMP rays at the analysis location for the same trial velocities. The CMP gather and the data windows for Horizon 7a are shown in Figure 9.2-7.
Figure 9.2-9 Semblance curves from coherency inversion at CMP location 75 as in Figure 9.2-8 applied to the horizons as in Figure 9.2-2. The CMP gather is shown in Figure 9.2-7a and the data windows are shown in Figure 9.2-10.
Figure 9.2-10 Data windows from coherency inversion at CMP location 75 as in Figure 9.2-8 applied to the horizons as in Figure 9.2-2. The CMP gather is shown in Figure 9.2-7a and the semblance curves are shown in Figure 9.2-9. True layer velocities are posted at the bottom of each panel.
Figure 9.2-11 Data windows and semblance curves from coherency inversion for Horizon 7a at CMP location 75 as in Figure 9.2-8 using velocity increments of (a) 25 m/s, (b) 50 m/s, (c) 100 m/s, and (d) 200 m/s. The CMP gather is shown in Figure 9.2-7a. True layer velocity is 3500 m/s.

Coherency inversion

The results of coherency inversion for the model with low-relief structures (Figure 9.2-1) are shown in Figure 9.2-6. Note that, for shallow layers, the inversion yields accurate velocity estimates. It also has detected the velocity variations within the deltaic sequence over the full-fold CMP range, correctly — note the increase in velocity for Horizon 6 from around 3000 m/s on the left to 3200 m/s in the middle and back to 3000 m/s on the right, and compare with Figure 9.2-1. Nevertheless, note the very short-wavelength variations in the velocity curves derived from picking the maxima of the semblance plots. As was demonstrated by the coherency inversion tests applied to the model with horizontal layers (models with horizontal layers), this observation has an important practical implication with regard to evaluation and use of the results of velocity estimation. Specifically, lateral velocity variations of very short-wavelengths that are much less than a cable length should not be incorporated into a velocity-depth model. Instead, some lateral smoothing of velocity estimates is almost always needed.

Shown in Figure 9.2-7a is a CMP gather used in coherency inversion to obtain a velocity estimate for the layer above Horizon 7b. Following the procedure described in the previous section, use the normal-incidence rays to convert the time horizon picked from the stacked data (Horizon 7a in Figure 9.2-2) at the analysis location to a trial depth horizon. The normal-incidence depth conversions are shown in Figure 9.2-8 for the three trial velocities as in Figure 9.2-7. The normal-incidence rays in Figure 9.2-8 are overlayed on the true velocity-depth model. Note that, in case of 3000 m/s, the trial depth horizon is shallower than the actual layer boundary for Horizon 7a; it coincides with the actual boundary in case of 3500 m/s (the true velocity of the layer above Horizon 7a); and, it is deeper in the case of 4000 m/s velocity.

For each trial constant velocity assigned to the layer above Horizon 7a, perform modeling of CMP traveltimes at the analysis location. The raypaths associated with the CMP traveltimes for the three trial velocities are also displayed in Figure 9.2-8. Note that the velocity contrast at layer boundaries and the geometry of these boundaries within the overburden has caused reflection point smearing at Horizon 7a.

The CMP data window that includes the reflection event for Horizon 7a is displayed in Figure 9.2-7c along the modeled traveltime trajectories. The flatness criterion suggests that the optimum layer velocity at the analysis location is 3500 m/s. The maximum of the semblance curve derived from coherency inversion coincides with the optimum choice for the layer velocity (Figure 9.2-7b).

Velocity resolution

Figure 9.2-9 shows the semblance curves derived from coherency inversion at CMP location 75. (Horizon 3 is missing at this location.) The CMP gather itself is shown in Figure 9.2-7a. Note that the sharpness of the peak in the semblance curve, hence the velocity resolution, depends upon the depth of the layer boundary and the magnitude of the layer velocity. Also, recall from Figure 9.1-12 that the velocity resolution also depends on the effective cable length.

The sampling interval for the velocity axis in the semblance curves should be chosen by taking into consideration the velocity resolution that can be achieved. Figure 9.2-10 shows the CMP data windows with the reflections that correspond to the selected horizons. The horizontal axis spans the offset range used in coherency inversion. For shallow horizons (1, 2, and 4), the offset range is smaller because of muting. The CMP gather is shown in Figure 9.2-7a and the semblance curves are shown in Figure 9.2-9. For a given horizon, the data window with the flattest event is distinguishable when the velocity sampling is appropriate. Specifically, for shallow horizons with low velocity (e.g., Horizons 1 and 2), velocity increment needs to be small enough to pick a layer velocity, accurately. However, for deeper events with high velocity (e.g., Horizons 7a and 7e), the velocity increment does not have to be as small, since event curvature in the data windows becomes indistinguishable as seen in Figure 9.2-10. It is only with large velocity increments that we observe a marked difference in event curvature as demonstrated in Figure 9.2-11. A good rule of thumb in practice is that the sampling interval in velocity used in stacking velocity inversion and coherency inversion needs to be specified as small as 25 m/s for velocities as low as 1500 m/s and can be specified as large as 200 m/s for velocities as high as 5000 m/s.

9.3 Model with complex overburden structure

Figure 9.3-1 illustrates a complex overburden structure associated with overthrust tectonics. Such structures were formed as a result of the tectonic movements during the Lower Miocene and Upper Cretaceous, and are common in North America, South America, and the Middle East. The target horizon is the flat reflector at 2.5 km below the imbricate structures.

Figure 9.3-1 A velocity-depth model with complex overburden structure caused by overthrust tectonics.

The velocity-depth model comprises a shallow sequence with a relatively simple structure. Underneath this shallow sequence is a shale-marl sequence with a strong vertical velocity gradient (0.5 m/s/m). Then, we have the imbricated fault structures of the carbonate sequence, and finally the target level at 2.5 km characterized as the detachment zone that separates the incompetent rock layers above from the competent rock layers below.

A total of 154 shot records was modeled using the two-way acoustic wave equation. The simulated recording geometry consists of a split-spread cable with 97 receivers and an offset range of 0-2350 m. Shot and receiver intervals are both 50 m. Figure 9.3-2a shows the CMP-stacked section. Note the traveltime distortions along the deepest reflection caused by the severe ray bending within the overburden. Time migration (Figure 9.3-2b) — whether post- or prestack, will not resolve the deleterious effect of the complex overburden structure.

Image gathers

Suppose that we already have estimated the velocity-depth model in Figure 9.3-1 down to 1-km depth. We shall construct the remaining part of the model layer-by-layer — the shale-marl sequence, the imbricate structures and the substratum, by prestack depth migration. Our goal here is to examine image gathers from prestack depth migration for layer velocity estimation and the stack of the image gathers for reflector geometry delineation. The results of the layer-by-layer velocity-depth model estimation based on image-gather analysis are shown in Figures 9.3-3, 9.3-4, and 9.3-5.

Figure 9.3-3a shows the velocity-depth model with the known shallow sequence down to 1 km, and the unknown part that is represented by the half-space below. The velocity assigned to the half-space is that of the shale-marl sequence (Figure 9.3-1). Note that this layer has a vertical velocity gradient. Figure 9.3-3b shows the depth image from prestack depth migration using the velocity-depth model in Figure 9.3-3a. Superimposed on this section are the layer boundaries in the velocity-depth model.

Figure 9.3-6a shows selected image gathers from prestack depth migration. The event T associated with the top of the carbonate sequence exhibits a flat character on the image gathers. This indicates that the velocity field for the layer above is correct. Consequently, the depth image in Figure 9.3-3b, which was obtained by stacking the image gathers as in Figure 9.3-6a, yields the correct reflector geometry for the top of the carbonate sequence. Interpret the top of the carbonate sequence from this section and insert it as a layer boundary into the velocity-depth model (Figure 9.3-4a).

Next, assign the velocity for the carbonate sequence (5700 m/s) to the half-space below the top-carbonate boundary as shown in Figure 9.3-4a. Using this model, perform prestack depth migration and obtain the depth image shown in Figure 9.3-4b. The image gathers indicate flat events for the top T and base B of the carbonate sequence (Figure 9.3-6b), thus confirming that the velocity of the carbonate sequence is correct. Additionally, the depth image in Figure 9.3-4b, which was obtained by stacking the image gathers as in Figure 9.3-6b, yields the correct reflector geometry for the base of the carbonate sequence. Interpret the base of the carbonate sequence from this section and insert it as a layer boundary into the velocity-depth model (Figure 9.3-5a).

Finally, assign the substratum velocity (5000 m/s) to the half-space below the base-carbonate boundary as shown in Figure 9.3-5a. Using this model, perform prestack depth migration and obtain the depth image shown in Figure 9.3-5b. The image gathers indicate flat events for the top T and base B of the carbonate sequence, and the flat reflector F within the substratum (Figure 9.3-6c). The velocity-depth model in Figure 9.3-5a is the same as the true velocity-depth model in Figure 9.3-1. Additionally, the section in Figure 9.3-5b, which was obtained by stacking the image gathers as in Figure 9.3-6c, represents the correct image in depth. Interpret the flat reflector at 2.5 km from this section and insert it as a layer boundary into the velocity-depth model (Figure 9.3-5a).

Now examine further the image gathers in Figure 9.3-6. The flatness of an event on image gathers is an indication of the accuracy of the velocity field associated with the layer above the layer boundary that is represented by that event. Actually, to declare the layer velocity as correct and accurate, not only the event associated with the base-layer, but also all events above that event should be flat. Note that in Figure 9.3-6a all events down to and including the top-carbonate are fairly flat.

The "nonflatness" of an event on image gathers is detectable only if there is sufficient cable length. For instance, the shallow events on image gathers in Figure 9.3-6a extend to a narrow range of offsets as a result of muting. As a result, it is difficult, if not impossible, to infer how flat these events are.

The events associated with the base-carbonate and the underlying flat reflector exhibit residual moveout on image gathers in Figure 9.3-6a. The event curvature for both layer boundaries is upward; this suggests that the velocity field assigned to the half-space which includes these two layer boundaries (Figure 9.3-3a) is erroneously low.

Figure 9.3-2 (a) CMP-stacked section associated with the velocity-depth model in Figure 9.3-1; (b) poststack time migration.
Figure 9.3-3 (a) Velocity-depth model associated with the overthrust data in two parts — the known part down to 1 km and the unknown part represented by the half-space below; (b) prestack depth migration using the model in (a). Half-space velocity is that of the shale-marl sequence — 3200 m/s referenced to the top of the layer with 0.5 m/s/m gradient.
Figure 9.3-4 (a) Velocity-depth model associated with the overthrust data in two parts — the known part down to the top of the imbricate structures and the unknown part represented by the half-space below; (b) prestack depth migration using the model in (a). Half-space velocity is that of the carbonate sequence — 5700 m/s.
Figure 9.3-5 (a) Velocity-depth model associated with the overthrust data in two parts — the known part down to the base of the imbricate structures and the unknown part represented by the half-space below; (b) prestack depth migration using the model in (a). Half-space velocity is that of the substratum — 5000 m/s.
Figure 9.3-6 Image gathers from prestack depth migration using the velocity-depth models shown in Figures (a) 9.3-3a, (b) 9.3-4a, and (c) 9.3-5a, respectively. T, B, and F are the events associated with the top-carbonate, base-carbonate, and the underlying flat reflector, respectively.
Figure 9.3-7 (a) Velocity-depth model associated with the overthrust data in two parts — the known part down to the base of the imbricate structures and the unknown part represented by the half-space below; (b) prestack depth migration using the model in (a). Half-space velocity is 4500 m/s.
Figure 9.3-8 (a) Velocity-depth model associated with the overthrust data in two parts — the known part down to the base of the imbricate structures and the unknown part represented by the half-space below; (b) prestack depth migration using the model in (a). Half-space velocity is 5500 m/s.
Figure 9.3-9 Image gathers from prestack depth migration using the velocity-depth models shown in Figures (a) 9.3-7a, and (b) 9.3-8a, respectively. T, B, and F are the events associated with the top-carbonate, base-carbonate, and the underlying flat reflector, respectively.

The detectability of residual moveout on image gathers is possible, again, only if there is sufficient cable length. The smaller the effective cable length for an event, the less detectable is the residual moveout, thus the poorer the velocity resolution. In the limit of zero offset, velocity resolution becomes nill. The residual moveout also is influenced by the magnitude of the layer velocity and the depth of the layer boundary.

Using the correct velocity for the carbonate sequence (Figure 9.3-4a), the event associated with the base-carbonate becomes flat on image gathers (Figure 9.3-6b). Nevertheless, the event associated with the flat reflector now exhibits a downward curvature, suggesting that the half-space velocity (5700 m/s) in Figure 9.3-4a is erroneously low for the substratum.

Finally, using the correct velocity for the substratum (Figure 9.3-5a), the event associated with the flat reflector F exhibits a flat character on image gathers (Figure 9.3-6c). Not only is this event flat in Figure 9.3-6c, but so are all the other events above. Suppose that an image gather contains 10 events, and that all are flat except the sixth from the top. This does not imply that the velocity-depth model is correct for all the layers except for the sixth layer. Instead, it implies that the velocity-depth model is correct down to and including the fifth layer, and the deeper part of the model is incorrect.

Flatness of all the events on image gathers is a means of verifying the accuracy of the velocity-depth model used in prestack depth migration. This is a necessary but not sufficient condition for verifying the accuracy of a model. As stated earlier in introduction to earth modeling in depth, by using nonzero-offset data, we can hope to resolve the velocity-depth ambiguity for a reflector if the data used in inversion have been recorded with offsets greater than the reflector depth. The additional limitations in model verification based on image-gather analysis are the reflector depths and the magnitude of the layer velocities. Specifically, the deeper the reflector or the higher the velocity of the layer above, the less the detectability of a residual moveout on image gathers (Figure 9.3-6).

Constant half-space velocity analysis

Image gathers from prestack depth migration are used to estimate layer velocities in two ways — constant half-space velocity analysis and residual moveout analysis. Suppose that the velocity-depth model has been established for the first n − 1 layers, and that we want to estimate the layer velocity for the nth layer. Using the known overburden velocity-depth model for the first n − 1 layers, assign a constant velocity to the half-space below that includes the nth layer. Perform prestack depth migration and output image gathers at some appropriate interval along the line. The image gathers should exhibit flat character for the events associated with the n − 1 layers, but show a residual moveout for the event associated with the nth layer. Repeat the analysis using the same overburden model and a range of constant velocities assigned to the half-space.

Consider the velocity-depth model in Figure 9.3-5a. The substratum below the base-carbonate layer boundary is the half-space that contains the flat reflector at 2.5 km. The half-space velocity is 5000 m/s. Selected image gathers from prestack depth migration using this model are shown in Figure 9.3-6c, and the depth image is shown in Figure 9.3-5b. Keep the same overburden model as in Figure 9.3-5a, but assign a constant velocity of 4500 m/s to the half-space as shown in Figure 9.3-7a. The resulting image gathers are shown in Figure 9.3-9a and the depth image in Figure 9.3-7b. The event on the image gathers associated with the flat reflector exhibits a small, but still detectable residual moveout that suggests erroneously low velocity for the layer above the flat reflector. Repeat the analysis for a half-space velocity of 5500 m/s (Figure 9.3-8a). Selected image gathers are shown in Figure 9.3-9b and the stack of image gathers in Figure 9.3-8b.

The velocity-depth models (Figure 9.3-5a, 9.3-7a, 9.3-8a) corresponding to the image gathers in Figures 9.3-6c, 9.3-9a, and 9.3-9b, respectively, have the same overburden, but with three different constant half-space velocities — 5000, 4500, and 5500 m/s, assigned to the substratum below the carbonate sequence. First, examine the residual moveout on the image gathers for the event associated with the flat reflector F. Although the moveout differences are subtle, the flatness is achieved with the 5000-m/s half-space velocity. By analyzing the image gathers from constant half-space velocity scans, we can make optimum velocity picks at analysis locations that best satisfy the flatness criterion for the layer under consideration.

9.4 Model building

Although doing it right the first time is most desirable, there is never a situation where this is possible when estimating an earth model in depth. The velocity-depth ambiguity that is inherent to inversion makes it very difficult getting the right answer — the true geological model, let alone the first time. Limitations in the resolving power of the methods to estimate layer velocities that arise from the band-limited nature of the recorded data and finite cable length used in recording further compound the problem. Finally, traveltime picking that is needed for most velocity estimation techniques and time-to-depth conversion as well as picking depth horizons from depth-migrated data to delineate reflector geometries are all adversely affected by noise present in the data.

All things considered, we can only expect to do our best in estimating what may be called an initial model, and update this model to get an acceptable final model. In this section, we shall discuss ways to estimate an initial model, and in the next section, we shall discuss application of residual moveout analysis and reflection traveltime tomography to update the initial model.

We shall discuss two strategies applicable to both 2-D and 3-D seismic data for initial model building:

A time-to-depth conversion strategy based on interpretation in the time domain, and
A layer-by-layer inversion strategy based on interpretation in the depth domain.

Practical methods to estimate layer velocities and delineate reflector geometries used in implementing the two strategies are listed in Table 9-1. A widely used combination for time-to-depth conversion is Dix conversion to estimate layer velocities and image-ray depth conversion to delineate reflector geometries. Whereas for layer-by-layer inversion, a widely used combination is coherency inversion to estimate layer velocities and poststack depth migration to delineate reflector geometries.

Time-to-depth conversion

The time-to-depth conversion strategy involves the following steps:

Interpret a set of time horizons from an image volume derived from time migration; these time horizons are usually associated with layer boundaries with velocity contrast or geological formations of interest.
Intersect rms velocity functions picked at specified analysis locations over the survey area with the time horizons from step (a) to derive horizon-consistent rms velocity maps. The rms velocity functions are preferably picked from gathers derived from prestack time migration.
Perform Dix conversion of the rms velocity maps from step (b) to derive interval velocity maps.
Perform vertical-ray or image-ray depth conversion of the time horizons from step (a) using the interval velocity maps from step (c).

The combination of the interval velocity maps from step (c) with the depth horizons from step (d) constitutes the initial model derived from time-to-depth conversion. This initial model may then be calibrated to well data if the depth horizon maps are to be used as depth structure maps for well positioning. Alternatively, the estimated initial model may be updated and used to perform 3-D post- or prestack depth migration to derive an image volume in depth.

In reference to step (c) of the procedure outlined above, interval velocities may be estimated by methods other than Dix conversion (Table 9-1). Nevertheless, as part of the common strategy for time-to-depth conversion, interval velocities are derived from Dix conversion of rms velocities. In reference to step (d) of the procedure outlined above, depth conversion of time horizons may be performed by one of the following three strategies:

Most commonly applied strategy is based on a combination of Dix conversion of rms velocities to interval velocities and image-ray depth conversion of time horizons interpreted from the time-migrated volume of data. This is the usual implementation of map migration. Stacking velocity inversion sometimes may be substituted for Dix conversion to estimate interval velocities.
Alternatively, depth conversion may be performed using vertical rays. This is acceptable only if lateral mispositioning because of lateral velocity variations is negligible (introduction to earth imaging in depth). Again, the interval velocities are estimated by Dix conversion.
Albeit rarely, a third option is to use normal-incidence rays for depth conversion. Time horizons interpreted from the time-migrated volume of data may first be forward-modeled to derive 3-D zero-offset traveltimes, which are then depth-converted using normal-incidence rays. Dix conversion still is the robust method for interval velocity estimation.

We now describe the details of the steps involved in the usual implementation of map migration based on image-ray depth conversion to derive depth structure maps.

Time structure maps

Time horizons are picked from image volumes obtained from either 3-D post- or 3-D prestack time migration. Aside from improved imaging of conflicting dips with different stacking velocities, the latter offers the advantage of providing a 3-D rms velocity field that is associated with events in their migrated positions (3-D prestack time migration). Figure 9.4-1 shows a 3-D view of the image volume derived from 3-D prestack time migration of a marine 3-D data set.

The interpreter identifies the time horizons that are associated with depositional sequence boundaries and geologically and lithologically significant layer boundaries within some of the depositional units. Then, reflection times are picked by combining seed detection and line-based interpretation strategies (interpretation of 3-D seismic data). A simplified form of an interpretation session without explicit fault identification is given below.

Seed points are placed on the event that is being picked at locations with good signal-to-noise ratio. These are then used to drive a seed detection algorithm to pick patches of the time surface around each seed point. Figure 9.4-2 shows picks from six time horizons from the shallowest (H1) to the deepest (H6). Depending on the signal-to-noise ratio and complexity of geometry of the time horizon, the extent of the surface patches varies from horizon to horizon. Some horizons are almost entirely picked by seed detection (H2), some are covered by limited amounts of seed-detected surface patches (H4), and some are not eligible for seed detection (H6). Note that seed detection has failed especially in intensively faulted areas.
To ensure structural control and adequate coverage of the time horizons, additional picking along inlines and crosslines is required. Figure 9.4-2 shows the horizon strands derived from line-based picking.
The surface patches derived from seed detection and horizon strands derived from line-based picking are then combined to form the complete set of control points for each horizon as shown in Figure 9.4-2. At this stage, a comprehensive editing and repicking are required to ensure consistency in picking. The edited control points are then input to a surface fitting algorithm (Section J.5) to create grid points that define the surface by a map function t_n(x, y) at every inline and crossline intersection, where the function value t_n represents the reflection time at the (x, y) location on the nth surface. Figure 9.4-3 shows the gridded surfaces derived from the control points shown in Figure 9.4-2 for the six horizons. The red corresponds to the highs and the blue corresponds to the lows on each horizon. Each horizon has been color-coded independently.

Note the fault signatures especially evident on horizons H1, H2, and H3. Actually, the sharp boundary between the green-blues and the yellow-reds observed on all the maps is clear evidence of a major fault through the middle of the survey area parallel to the longer dimension. Adjoining this major fault are the oblique faults that can be detected on H1, H2, and H3. Since we use rays to perform time-to-depth conversion, it is important to ensure that ray tracing is made stable by applying a carefully measured amount of smoothing to the gridded surfaces (Figure 9.4-4). This smoothing also is needed to edit outliers among the control points that have inevitably corrupted the grid points. Finally, the gridded surfaces are usually displayed in the form of contour maps as shown in Figure 9.4-5.

Interval velocity maps

The Dix equation (1), which relates rms velocities to interval velocities, is used to derive interval velocity maps. RMS velocities, in principle, are most appropriately estimated from prestack time-migrated data (3-D prestack time migration). Recall from normal moveout that the type of velocity that can be most reliably estimated from CMP data is the velocity used to apply normal-moveout correction. To stack the data we also substitute NMO velocities for stacking velocities. The use of NMO velocities as stacking velocities is based on the small-spread hyperbola assumption. To further substitute stacking velocities for rms velocities is only allowed if the CMP data are associated with horizontally layered earth. To justify the use of stacking velocities as rms velocities, we first need to correct for the dip effect on stacking velocities by way of dip-moveout (DMO) correction (principles of dip-moveout correction). In the case of a 3-D survey, we also need to correct for the source-receiver azimuthal effects on stacking velocities by way of 3-D DMO correction (processing of 3-D seismic data). This means that it is the stacking velocity field derived from 3-D DMO-corrected data that should be considered as a plausable substitute for the rms velocity field. But then the DMO velocities are associated with CMP gathers in their unmigrated positions. Strictly, we need the moveout velocities not only corrected for dip and azimuth effects but also estimated from gathers in their migrated positions. This is because the rms velocities used in Dix equation (1) are defined for a horizontally layered earth model (Section J.4). Thus the desired strategy is that velocities derived from 3-D prestack time migration should be substituted for rms velocities.

Although prestack time migration velocities are most desired to substitute for rms velocities, the interpreter may be compelled to use whatever velocity functions that may be available. These may have been derived from velocity analysis applied to DMO-corrected data or even, although hardly desirable, to CMP data without DMO correction. Under those circumstances, the velocity functions picked at analysis locations need to be edited for any dip effect by either eliminating the suspect functions altogether or by smoothing.

Figure 9.4-1 An image volume of data derived from 3-D prestack time migration. The color-coded top surface is the water bottom with the acquisition footprint exhibited by the striations along the inline direction.
Figure 9.4-2 Results of interpretation of six time horizons from the image volume shown in Figure 9.4-1.
Figure 9.4-3 Gridded surfaces derived from the horizon picks shown in Figure 9.4-2.
Figure 9.4-4 Smoothed versions of the gridded surfaces shown in Figure 9.4-3.
Figure 9.4-5 Time contour maps associated with the smoothed surfaces shown in Figure 9.4-4.
Figure 9.4-6 Horizon-consistent rms velocity maps associated with the time horizons shown in Figure 9.4-5.
Figure 9.4-7 Interval velocity maps derived from the horizon-consistent rms velocity maps shown in Figure 9.4-6.

Whatever the source of information, the interpreter starts with a set of velocity functions, each made up of a set of time-velocity pairs and associated with analysis locations over the survey area. The analysis grid typically varies from 500 × 500 m to 2 × 2 km; hence, there may be as many as 400 velocity functions per 100 km-squared of the survey area. The gridded time horizons are intersected with the velocity functions and, for each horizon, velocity nodes are extracted from the velocity functions coincident with the horizon times at the locations of the velocity functions themselves. These velocity nodes are then used as control points input to a gridding algorithm to create horizon-consistent rms velocity maps (Figure 9.4-6). There may be a need for further editing and smoothing of the rms velocity grids.

Finally, the horizon-consistent rms velocity values and the horizon times at each grid point are used in Dix equation (1) to compute the interval velocity values, which are then used to create the horizon-consistent interval velocity maps (Figure 9.4-7). Once again, there may be further need for editing and smoothing of the interval velocity maps to remove any geologically implausable velocity variations.

Depth structure maps

We now have a set of time horizon maps (Figure 9.4-4) and a corresponding set of interval velocity maps (Figure 9.4-7). Recall from introduction to earth imaging in depth that, in principle, a time-migrated image can be converted to a depth section by mapping the amplitudes along image rays. This notion also can be employed to convert time horizons into depth horizons. The process is done layer by layer starting with the shallowest horizon.

A comprehensive mathematical discussion on image-ray tracing is given by Hubral and Krey ^[15]. Refer to Figure 9.4-8 for a brief description of image-ray tracing. Consider an image ray that departs the nth surface at point S_n with coordinates (x_n, y_n, z_n) and emerges at the right angle at the earth’s surface at point S₀ with coordinates (x₀, y₀, z = 0). Our goal is to determine the coordinates of the output point (x_n, y_n, z_n) on the depth map from image-ray depth conversion of the input point (x₀, y₀, t_n) on the time map. To achieve this goal, we want to trace the image ray from the point of emergence (x₀, y₀, 0) back to its point of departure (x_n, y_n, z_n).

Suppose that the first n − 1 horizons have already been converted to depth, and that next we want to convert the nth horizon to depth. Since the earth model is known for the first n − 1 layers, then we know the coordinates of the intersection point S₁ of the image ray with the first layer, (x₁, y₁, z₁), where by definition of the image ray, x₁ = x₀ and y₁ = y₀. By using Snell’s law, we can determine the direction of the ray as it departs point S₁ reaching point S₂ on the next surface. As the image ray moves from one surface to the next, we add up the time it takes to travel. When it reaches the (n − 1)st layer, the elapsed two-way time t_n−1 is

t_{n-1}=2\sum _{k}^{n-1}{\frac {\Delta s_{k}}{v_{k}}},

(2)

where v_k is the interval velocity of the kth layer, and Δs_k is the distance between the intersection points of the image ray, S_k−1 and S_k, on the (k − 1)st and kth surfaces given by

\Delta s_{k}=\left[(x_{k}-x_{k-1})^{2}+(y_{k}-y_{k-1})^{2}+(z_{k}-z_{k-1})^{2}\right]^{1/2}.

(3)

Now, we examine the situation when the image ray departs the (n − 1)st layer at point S_n−1 on the way to the nth surface. Again, by Snell’s law we know the direction of the ray. We also know the elapsed time t_n − t_n−1 from the (n − 1)st surface to the nth surface since we know the total elapsed time t_n−1 from equation (2) and the total elapsed time t_n from the input time horizon read at point (x₀, y₀, z = 0). Finally, we know the interval velocity v_n of the nth layer from the interval velocity map. Therefore, we can calculate the elapsed distance Δs_n along the raypath as it departs the point S_k−1 on the (k − 1)st surface in the direction dictated by Snell’s law. The quantity Δs_n is given by

\Delta s_{n}={\frac {v_{n}}{2}}(t_{n}-t_{n-1}).

(4)

Finally, the coordinates of the point S_n that we need to know to perform the time-to-depth conversion are given by

x_{n}=x_{n-1}+\Delta s_{n}\cos \ \alpha ,

(5a)

y_{n}=y_{n-1}+\Delta s_{n}\cos \ \beta ,

(5b)

z_{n}=z_{n-1}+\Delta s_{n}\cos \ \gamma ,

(5c)

where α, β, and γ are the directional cosines of the ray at point S_n−1. The directional cosines are known by the application of Snell’s law at point S_n−1 with known coordinates (x_n−1, y_n−1, z_n−1).

To summarize, given the depth and interval velocity maps for the first n − 1 horizons, and the time and interval velocity maps for the nth horizon, we can trace an image ray associated with the time t_n(x₀, y₀) on the time map and derive the depth value z_n(x_n, y_n) on the depth map.

Figure 9.4-8 Principles of image-ray depth conversion. See text for details.

Figure 9.4-9 shows the depth maps derived from image-ray depth conversion of the time maps shown in Figure 9.4-4 using the interval velocity maps shown in Figure 9.4-7. As for the time maps, the usual way to post the depth maps is to contour them (Figure 9.4-10).

The depth maps are compatible with the time maps; nevertheless, there can be subtle differences because of velocity variations that would give rise to the departure of image rays from the vertical. To quantify the differences between vertical-ray and image-ray trajectories, and thus to quantify the differences between the time maps and depth maps, we can calculate the modulus Δd_n of the lateral displacement vector between the points S₀ and S_n as

\Delta d_{n}={\sqrt {(x_{n}-x_{0})^{2}+(y_{n}-y_{0})^{2}}},

(6a)

and create the displacement modulus maps shown in Figure 9.4-11. Note that the most significant displacement between the vertical rays and image rays is at the fault zones.

The displacement vector also has a directional azimuth ϕ_n which is given by

\phi _{n}=\tan ^{-1}{\frac {y_{n}-y_{0}}{x_{n}-x_{0}}}

(6b)

as measured from the inline x direction. The displacement azimuth maps are shown in Figure 9.4-12. Again, note that the most significant azimuthal variations are along the fault zones.

Calibration to well tops

The depth structure maps derived from time-to-depth conversion or layer-by-layer inversion invariably will not match the well tops. The sources of discrepancy between the estimated reflector depths and the well tops include limitations in the methods for interval velocity estimation, mispicking of time horizons input to depth conversion, and limitations in the actual depth conversion itself within the context of ray tracing through an earth model that includes complex layer boundaries. For the depth structure maps to be usable in subsequent reservoir modeling and simulation, it is imperative to calibrate them to well tops.

Figure 9.4-9 Depth horizons derived from time-to-depth conversion of the time horizons shown in Figure 9.4-4 using the interval velocity maps shown in Figure 9.4-7.
Figure 9.4-10 Contour maps of the depth horizons shown in Figure 9.4-9.
Figure 9.4-11 Modulus maps for the displacement vectors that describe the lateral departure between vertical rays and image rays. See text for details.
Figure 9.4-12 Azimuth maps for the displacement vectors that describe the directional difference between vertical rays and image rays. See text for details.
Figure 9.4-13 (a) Depth horizons as in Figure 9.4-9 derived from image-ray depth conversion of the time horizons shown in Figure 9.4-4; (b) solid model created from the depth horizons in (a).
Figure 9.4-14 Explosion of the solid model shown in Figure 9.4-13b.

Consider a seismically derived depth structure map z_s(x, y) based on time-to-depth conversion, say, for the layer boundary associated with the top-reservoir. Also consider N_w well tops z_w (x_i, y_i) for this horizon at locations (x_i, y_i), i = 1, 2, …, N_w. Since the velocity-depth model derived from time-to-depth conversion is supposed to be consistent with the input data — the time structure map τ(x, y) created from the interpretation of the time-migrated volume of data, we have

\tau (x,y)=2{\frac {z_{s}(x,y)}{V_{s}(x,y)}},

(7a)

where, for the purpose of calibration, V_s(x, y) can be either the average or rms velocity map associated with the horizon z_s(x, y). Also, for simplicity in calibration, we consider vertical rays rather than image rays as in equation (7a).

There exists a calibration velocity V_c(x, y) such that, at a well location (x_i, y_i), it satisfies the relation

\tau (x_{i},y_{i})=2{\frac {z_{w}(x_{i},y_{i})}{V_{c}(x_{i},y_{i})}}.

(7b)

Combine equations (7a) and (7b) to get a relation that is satisfied at the well locations

{\frac {V_{c}(x_{i},y_{i})}{V_{s}(x_{i},y_{i})}}={\frac {z_{w}(x_{i},y_{i})}{z_{s}(x_{i},y_{i})}}.

(8)

From the knowledge of the well tops z_w(x_i, y_i) and the seismically derived reflector depths z_s(x_i, y_i) at the well locations, equation (8) gives a calibration factor c(x_i, y_i) = V_c(x_i, y_i)/V_s(x_i, y_i) computed at each of the well locations. Next, apply kriging or some other interpolation technique to the sparsely defined calibration factors c(x_i, y_i), i = 1, 2, …, N_w to derive a calibration factor map c(x, y) specified at all grid locations (x, y). Kriging is a statistical method of determining the best estimate for an unknown quantity such as c(x, y) at some location (x, y) using a sparse set of values such as c(x_i, y_i) specified at locations (x_i, y_i) ^[16] ^[17].

The final step in calibration is to scale the depth structure map z_s(x, y) by the calibration factor map c(x, y)

z_{c}(x,y)=c(x,y)\ z_{s}(x,y),

(9)

where z_c(x, y) is the calibrated depth structure map. Note that, by way of equations (8) and (9), the calibrated depth z_c coincides with the well top z_w at well location (x_i, y_i).

Calibration to well tops is done only after the completion of model building, and just before well planning and reservoir modeling. When estimating an earth model by following a layer-by-layer inversion procedure (next subsection), depth horizon associated with the (n − 1)st layer should not be calibrated before estimating the model for the next layer n. This is because seismically derived layer velocities almost never match with well velocities. The discrepancy between the two is attributable to several factors, including the limited resolution in velocities estimated from seismic data (Sections models with horizontal layers, model with low-relief structure, and model with complex overburden structure) and seismic anisotropy (seismic anisotropy). Additionally, the high-frequency variations in the well velocities are absent from the seismically derived velocities.

The calibrated depth maps can be used to create a solid model of the earth as illustrated in Figure 9.4-13. Each layer is represented by a solid (Figure 9.4-14) with its interior populated by specific layer parameters. These may include compressional- and shear-wave velocities, densities, and rock physics parameters such as porosity, permeability, pore pressure, and fluid saturation. When populated by the petrophysical parameters, the solid associated with the reservoir layer represents a reservoir model. For the purpose of reservoir modeling, the solid for the reservoir layer usually is downscaled in the vertical direction by dividing it into thin slices with a thickness as small as 1 m — much less than the threshold for vertical seismic resolution (seismic resolution). Additionally, the solid for the reservoir layer is upscaled in the lateral direction by dividing each thin slice into finite elements with a varying size of up to 250 m on one side. The reservoir model is eventually fed into a reservoir simulation scheme to predict the geometry of the fluid flow from the given reservoir parameters.

Layer-by-layer inversion

The second strategy for initial model building involves a layer-by-layer application of an appropriate combination of inversion methods listed in Table 9-1. Layer-by-layer inversion involves the following steps:

Interpret a set of time horizons from unmigrated data to be used in lieu of zero-offset reflection traveltimes required by coherency inversion. Alternatively, interpret a set of time horizons from time-migrated data and perform forward modeling to obtain the required zero-offset traveltimes.
Assume that interval velocities and reflector geometries for the first n − 1 layers have been estimated, and that we want to estimate the same for the nth layer. By using the time horizons from step (a), first apply one of the inversion methods listed in the left-hand column of Table 9-1 to estimate the interval velocity field.
Assign the estimated interval velocity field for the nth layer to the half space that includes the unknown part of the model — the nth layer and the layers below, and perform depth migration. Often, it is sufficient to do poststack depth migration.
Interpret the depth horizon associated with the base of the nth layer from the depth image and incorporate it into the velocity-depth model.

The layer-by-layer inversion strategy alternates between layer velocity estimation and reflector geometry delineation for each layer starting from the earth’s surface and moving down one layer at a time. The layer-by-layer estimation strategy facilitates checking the results of inversion for one layer before moving down to the next. As such model updating can be interleaved with model estimation to circumvent accumulation of errors in layer velocities and reflector geometries as the analysis proceeds from the top down.

We shall demonstrate the layer-by-layer inversion strategy by applying a combination of coherency inversion to estimate layer velocities and poststack depth migration to delineate reflector geometries to the marine data set shown in Figure 9.4-15. Superimposed on the stacked section are the segments of interpreted reflection traveltime horizons. Starting from the top, H1 is the water bottom whereas H6 is the top-salt boundary. The abundance of diffractions within the salt layer is associated with high-velocity anhydrite-dolomite rafts. Input to coherency inversion is the zero-offset reflection traveltimes associated with normal-incidence rays and unmigrated data as shown in Figure 9.4-15. Diffraction flanks present in the unmigrated data must be avoided when interpreting the reflection traveltime segments from unmigrated data.

Details of coherency inversion for interval velocity estimation are provided in models with horizontal layers. Results of the layer-by-layer analysis have been compiled in Figures 9.4-16 through 9.4-21. Begin with the trivial task of modeling the water layer by normal-incidence time-to-depth conversion of the time horizon H1 in Figure 9.4-15 using a constant layer velocity of 1500 m/s. For each layer H2 through H6, and one layer at a time, the analysis includes the following steps:

Apply coherency inversion to compute the horizon-consistent semblance spectrum (panel (a) in Figures 9.4-16 through 9.4-21), and
pick the interval velocity profile from the semblance spectrum by tracking the semblance peaks.
Assign the interval velocity profile to the half space that includes the unknown layer itself and use the known overburden layer velocities and reflector geometries to create the intermediate velocity-depth model (panel (b) in Figures 9.4-16 through 9.4-21).
Perform poststack depth migration using the intermediate velocity-depth model (panel (c) in Figures 9.4-16 through 9.4-21).
Interpret the depth horizon that corresponds to the base boundary of the layer under consideration.
Incorporate the interpreted depth horizon into the intermediate velocity-depth model (panel (d) in Figures 9.4-16 through 9.4-21).
Proceed to the next layer below and repeat the above steps.

Figure 9.4-15 An unmigrated CMP-stacked section with seven time horizons that correspond to layer boundaries with significant velocity contrast.

The coherency semblance spectra shown in panel (a) of Figures 9.4-16 through 9.4-21 have been compiled into a single panel as shown in Figure 9.4-22. First, note the gradual decrease in resolution as we move down to the deeper layers. Much like conventional stacking velocity analysis, the resolving power of coherency inversion is governed by the cable length, the reflector depth, the layer velocity, and bandwidth of the data. Second, observe the short-wavelength variations in the velocity profiles defined by the semblance peaks. As was demonstrated by the synthetic data experiments in models with horizontal layers, these are caused by lateral velocity variations in the overlying layers that are much less than a cable length (Figures 9.1-13 and 9.1-20). In the present case shown in Figure 9.4-22, note that the interval velocity profiles for layers H6 and H7 exhibit much more pronounced fluctuations compared to the layers above. This is caused by the strong lateral variations in the interval velocity profile for layer H5. In practice, the interval velocity profile picked from the semblance spectrum must exclude such rapid fluctuations. Otherwise, the reflector geometry associated with the base of the layer under consideration will be corrupted by geologically implausable variations. Also note from the semblance spectra that we have to interpolate the interval velocity profiles through the zones with missing reflection events.

The accuracy in velocity estimation by coherency inversion can be monitored closely by examining the modeled CMP traveltimes and the associated raypaths. Shown in Figure 9.4-23 are the CMP gathers at a selected analysis location for each of the layers H2 through H7. Superimposed on these gathers are the modeled CMP traveltime trajectories that correspond to the optimum layer velocities picked from the semblance spectra shown in Figure 9.4-24. The raypaths associated with the modeled CMP traveltimes shown in Figure 9.4-25 illustrate the ray bending at layer boundaries. As such, the modeled traveltimes in Figure 9.4-23 are in general nonhyperbolic. From the surface lateral extent of the raypath family for each layer in Figure 9.4-25 and the offset range of the modeled traveltimes in Figure 9.4-23, we note the resolution that can be attained from the semblance spectra. Specifically, note that the sharpness of the spectra in Figure 9.4-24 decreases with increasing layer velocity, decreasing cable length, and increasing reflector depth.

Figure 9.4-16 Velocity estimation and reflector geometry delineation for layer H2 as denoted in Figure 9.4-15. See text for details.
Figure 9.4-17 Velocity estimation and reflector geometry delineation for layer H3 as denoted in Figure 9.4-15. See text for details.
Figure 9.4-18 Velocity estimation and reflector geometry delineation for layer H4 as denoted in Figure 9.4-15. See text for details.
Figure 9.4-19 Velocity estimation and reflector geometry delineation for layer H5 as denoted in Figure 9.4-15. See text for details.
Figure 9.4-20 Velocity estimation and reflector geometry delineation for layer H6 as denoted in Figure 9.4-15. See text for details.
Figure 9.4-21 Velocity estimation and reflector geometry delineation for layer H7 as denoted in Figure 9.4-15. See text for details.

The results of the alternating steps of layer velocity estimation and reflector geometry delineation involved in layer-by-layer model building have been compiled in Figure 9.4-26. The left and the right columns correspond to the compilations of panels (b) and (c) of Figures 9.4-16 through 9.4-21, respectively. The analysis sequence illustrated by Figure 9.4-26 begins with panel (a) and ends with panel (l). As outlined above, perform coherency inversion to generate the interval velocity profile and, hence, the intermediate velocity-depth model for layer H2 (panel (a)). Next, perform poststack depth migration and delineate the reflector geometry associated with the base of layer H2 (panel (b)). Then move down to the next layer and repeat the same analysis (panels (c) and (d)), until you reach the last layer in the sequence (panels (k) and (l)).

Figure 9.4-22 Semblance spectra from coherency inversion with picked interval velocity profiles for layers H2 through H7 as denoted in Figure 9.4-15.

The final form of the velocity-depth model estimated by applying the layer-by-layer inversion strategy is shown in Figure 9.4-27a. Check the consistency of this velocity-depth model with the depth image derived from poststack depth migration (Figure 9.4-27b) and the stacked section (Figure 9.4-27c). Specifically, overlay the depth horizons from the velocity-depth model onto the depth image and note that the reflector geometries implied by the latter are in agreement with the depth horizons. Then, perform zero-offset normal-incident modeling of the traveltimes associated with the layer boundaries included in the velocity-depth model and overlay them onto the stacked section. Again, note that the actual reflection traveltimes on the stacked section are in good agreement with the modeled traveltimes. This exercise confirms the consistency of the estimated velocity-depth model with the input stacked data. Yet, the estimated model is but one of the many possible solutions. It can only be considered an initial model; thus, it needs to be verified by checking its consistency with prestack data and, finally, it needs to be updated (model updating).

We now examine the performance of the layer-by-layer inversion strategy applied to structurally complex data. Figure 9.4-28a shows a time-migrated stacked section with a complex overburden structure associated with salt tectonics. The top-salt boundary is represented by horizon H6. The section has been interpreted to identify the layer boundaries with significant velocity contrast. For coherency inversion to estimate interval velocities layer by layer, the required zero-offset traveltimes were modeled from the time horizons interpreted from the time-migrated section. The modeled zero-offset traveltimes are shown in Figure 9.4-28 superimposed on the unmigrated stacked section; note that they are consistent with the observed reflection times. The alternative to using the modeled zero-offset times for coherency inversion is clearly the reflection times picked directly from the unmigrated stacked data. While this alternative may be appropriate for cases of simple structures (Figure 9.4-15), in the present case of a complex structure it is easier and safer to interpret the time-migrated data.

Figure 9.4-23 CMP gathers associated with the stacked section in Figure 9.4-15 displayed with the modeled traveltimes overlayed. The corresponding CMP raypaths are shown in Figure 9.4-25.
Figure 9.4-24 Semblance spectra derived from coherency inversion to estimate velocities for layers H2 through H7 as denoted in Figure 9.4-15.
Figure 9.4-25 Part 1: The CMP raypaths associated with the modeled traveltimes shown in Figures 9.4-23.
Figure 9.4-25 Part 2: The CMP raypaths associated with the modeled traveltimes shown in Figures 9.4-23.
Figure 9.4-26 Model building layer by layer starting from the surface. Left column shows the velocity-depth models and the right column shows the image sections from poststack depth migration using the stack shown in Figure 9.4-15.

Results of the alternating steps of layer velocity estimation using coherency inversion and reflector geometry delineation using poststack depth migration involved in layer-by-layer model building have been compiled in Figure 9.4-29. The left column shows the intermediate models for layers H2 to H5 (H1 is the water layer), and the right column shows the intermediate depth images from poststack depth migration. Start with the intermediate velocity-depth model in panel (a) and perform poststack depth migration to get the depth image in panel (b). Interpret the depth horizon associated with the base of layer H2 and insert this horizon into the velocity-depth model in panel (c). Then, move down to the next layer and perform coherency inversion to derive the interval velocity profile for it. Assign this velocity profile to the half space below the base of layer H2. This constitutes the updated intermediate velocity-depth model (panel (c)) which is used to obtain the depth image in panel (d). Continue in this alternating manner for estimating the layer velocities and delineating the reflector geometries until you reach horizon H5 (panel (h)).

We do not intend to proceed further and complete the model building for all the layers. Instead, we shall examine the accuracy of velocity estimation for the next layer H6 just above the salt layer. Figure 9.4-30 shows the coherency semblance spectra for layers within the overburden from H2 to H6. We observe that the quality of the semblance peaks is very good for layer H2, while it begins to degrade almost immediately for the layers below. Specifically, in the case of layer H3, the semblance peaks lose their sharpness at the center portion of the line where layer H3 becomes very thin (Figure 9.4-28a) causing instability in ray tracing. Next, the lateral velocity variation in layer H4 has caused rapid ondulations in the semblance spectrum for layer H5. Finally, note that the quality of the semblance peaks has deteriorated significantly for layer H6.

Figure 9.4-27 (a) Final velocity-depth model derived from the layer-by-layer application of coherency inversion and poststack depth migration, (b) depth migration using the velocity-depth model in (a), (c) modeled zero-offset traveltimes using the velocity-depth model in (a) overlayed on the unmigrated stacked section as in Figure 9.4-15.
Figure 9.4-28 (a) A time-migrated stacked section with a set of interpreted horizons, (b) the unmigrated stacked section with the zero-offset traveltimes modeled from the time horizons in (a).
Figure 9.4-29 Model building layer by layer starting from the surface. Left column shows the velocity-depth models and the right column shows the image sections from poststack depth migration using the stack shown in Figure 9.4-28b. Layers are labeled as the horizons that correspond to the base labeled as in Figure 9.4-28a.
Figure 9.4-30 Semblance spectra derived from coherency inversion to estimate velocities for layers H2 through H6 as denoted in Figure 9.4-28a.

To closely examine this rapid degradation in the quality of semblance, we shall conduct velocity estimation for layer H6 at four selected locations. For each location, Figures 9.4-31 to 9.4-34 show the modeled CMP raypaths that correspond to the modeled CMP traveltime trajectory computed by assigning the velocity associated with the semblance peak to the half space below horizon H5. Where the overburden is relatively less complex, we note that the semblance peak is distinctive and thus the estimate velocity is not ambiguous (Figure 9.4-31). At CMP locations where the overburden is sufficiently complex to cause significant ray bending at layer boundaries, however, the semblance curves have poor quality (Figures 9.4-32 to 9.4-34). Note that, at these locations, reflection point dispersal at the vicinity of the normal-incidence point is large and raypaths are complex. As a result, the complex layer boundaries adversely affect the quality of interval velocity estimation using methods that rely on ray tracing. Also recall from the model experiments in models with horizontal layers that lateral velocity variations in the layers above cause rapid fluctuations in the semblance profile for the layer below. We conclude that all known practical velocity estimation techniques based on ray theory alone suffer from a degradation of lateral resolution in areas with complex overburden structures.

Structure-independent inversion

In areas with low-relief structures and moderate lateral velocity variations, a structure-independent inversion strategy can be used to circumvent interpretation of time horizons when deriving an initial estimate of the earth model. Compared to a layer-by-layer inversion strategy, it can prove to be robust and less labor-intensive. We shall outline a procedure for structure-independent earth model estimation using a field data example.

Shown in Figure 9.4-35a is a time-migrated CMP-stacked section that exhibits a highly developed deltaic depositional sequence. Note from the stacking velocity field shown in Figure 9.4-35b that the lateral velocity variations are mild to moderate. Because the dips are gentle and the structures have low reliefs, we may substitute the stacking velocity field in Figure 9.4-35b for the rms velocity field that we need for Dix conversion.

Consider a set of fictitious, flat time horizons as displayed in Figure 9.4-36a, and extract the rms velocity profiles along these horizons from the rms velocity section (Figure 9.4-36b).
Perform Dix conversion to generate interval velocity profiles from the rms velocity profiles (Figure 9.4-37a). Note that for deeper horizons, lateral velocity variations in the layers above have caused oscillations in the interval velocity profiles.
Apply lateral smoothing to remove these oscillations and use the edited interval velocity profiles to convert the flat time horizons (Figure 9.4-36a) to depth horizons.
Combine the interval velocity profiles (Figure 9.4-37a) with the depth horizons to build an initial interval velocity field as shown in Figure 9.4-37b. Note that lateral velocity variations have caused the flat time horizons (Figure 9.4-36a) to transform to nonflat depth horizons.
Perform poststack depth migration (Figure 9.4-38a) using the initial interval velocity field (Figure 9.4-37b) and overlay the depth horizons derived in step (c) onto the depth section. Note that the depth horizons do not conform to the geometry of the reflectors inferred by depth migration; thus, the term structure-independent model estimation.
Discard the structure-independent depth horizons and replace them with the depth horizons interpreted from the depth-migrated section (Figure 9.4-38b).
Overlay the depth horizons from step (f) onto the interval velocity section from step (d) (Figure 9.4-39a).
Extract the interval velocity profiles along the depth horizons from the interval velocity section (Figure 9.4-39b).
Eliminate the oscillations from these profiles and combine them with the depth horizons from step (f) to build a structurally consistent earth model in depth (Figure 9.4-40a).
Perform prestack depth migration and obtain the image section shown in Figure 9.4-40b from the image gathers as in Figure 9.4-41.

Events on image gathers, except for the multiples, are mostly flat. This means that the estimated earth model in depth (Figure 9.4-40a) is fairly accurate. In practice, to attain consistency of the estimated model with the input data, depth migration may have to be iterated a few times (2-D poststack depth migration). This then is followed by model updating with reflection tomography, which is discussed in the next section.

9.5 Model updating

Limitations in the techniques for velocity estimation and velocity-depth ambiguity inherent to seismic inversion are compelling reasons for the need to update an estimated earth model in depth, however it has been constructed. Unfortunately, model updating tools themselves also have limitations in terms of their ability to resolve lateral velocity variations and refine reflector geometries. Again, the cable length and reflector depth dictate the extent that model updating techniques can resolve the velocity-depth ambiguity. In this section, we shall review residual moveout corrections applied to image gathers as a local method and reflection tomography as a global method for model updating.

Figure 9.4-31 (a) Raypaths associated with the optimum layer velocity estimated from coherency inversion applied to the gather at CMP location 216, (b) the coherency semblance spectrum with its peak corresponding to the velocity assigned to the layer under consideration for computing the raypaths in (a), (c) enlarged view of the raypaths as in (a).
Figure 9.4-32 (a) Raypaths associated with the optimum layer velocity estimated from coherency inversion applied to the gather at CMP location 479, (b) the coherency semblance spectrum with its peak corresponding to the velocity assigned to the layer under consideration for computing the raypaths in (a), (c) enlarged view of the raypaths as in (a).
Figure 9.4-33 (a) Raypaths associated with the optimum layer velocity estimated from coherency inversion applied to the gather at CMP location 596, (b) the coherency semblance spectrum with its peak corresponding to the velocity assigned to the layer under consideration for computing the raypaths in (a), (c) enlarged view of the raypaths as in (a).
Figure 9.4-34 (a) Raypaths associated with the optimum layer velocity estimated from coherency inversion applied to the gather at CMP location 803, (b) the coherency semblance spectrum with its peak corresponding to the velocity assigned to the layer under consideration for computing the raypaths in (a), (c) enlarged view of the raypaths as in (a).
Figure 9.4-35 (a) A CMP-stacked section, and (b) the stacking velocity field.
Figure 9.4-36 (a) The stacking velocity field as in Figure 9.4-35b with fictitious flat time horizons superimposed, and (b) the rms velocity profiles extracted from the section in (a) along the fictitious flat time horizons.
Figure 9.4-37 (a) The interval velocity profiles derived from the rms velocity profiles in Figure 9.4-36b by way of Dix conversion, and (b) a structure-independent velocity-depth model created by combining the interval velocity profiles as in (a) with the fictitious depth horizons that correspond to the flat time horizons as in Figure 9.4-36a.
Figure 9.4-38 (a) Depth migration of the stacked section in Figure 9.4-35 using the velocity-depth model in Figure 9.4-37b with its associated structure-independent depth horizons superimposed, and (b) the same section as in (a) with a set of interpreted structurally consistent depth horizons.
Figure 9.4-39 (a) The velocity-depth model as in Figure 9.4-37b with a set of interpreted structurally consistent depth horizons as in Figure 9.4-38b, and (b) the interval velocity profiles extracted from the section in (a) along the structurally consistent depth horizons superimposed onto the same section.
Figure 9.4-40 (a) A structurally consistent velocity-depth model created by combining the interval velocity profiles as in Figure 9.4-39b with the depth horizons as in Figure 9.4-38b, and (b) prestack depth migration of the data associated with the stacked section in Figure 9.4-35.
Figure 9.4-41 Selected image gathers associated with the image section derived from prestack depth migration shown in Figure 9.4-40.

Residual moveout analysis

We begin with building an initial velocity-depth model from the data as in Figure 9.4-15 using the time-to-depth conversion strategy described in model building. Figure 9.5-1a shows the time-migrated stacked section with a set of interpreted horizons that correspond to layer boundaries with significant velocity contrast. Perform Dix conversion of the horizon-consistent rms velocity profiles (Figure 9.5-1b) to derive the interval velocity profiles as shown in Figure 9.5-1c. Clearly, you often are compelled to apply some smoothing to the rms velocity profiles before Dix conversion and even apply additional smoothing to the interval velocity profiles afterwards. Then, perform image-ray depth conversion of the time horizons interpreted from the time-migrated section (Figure 9.5-1a) to generate the depth horizons. Finally, combine the interval velocity profiles with the depth horizons to create the velocity-depth model shown in Figure 9.5-2a.

Complete the analysis by checking for consistency of this estimated initial velocity-depth model with the depth image derived from poststack depth migration and the stacked section. Note that the reflector geometries inferred by the depth image shown in Figure 9.5-2b are in agreement with the depth horizons. Additionally, observe that the actual reflection traveltimes on the stacked section are in good agreement with the modeled zero-offset traveltimes as shown in Figure 9.5-2c.

We now check for consistency of the initial velocity-depth model in Figure 9.5-2a with prestack data and update it by correcting for the residual moveout observed on image gathers derived from prestack depth migration. Image gathers are like moveout-corrected CMP gathers with vertical axis in depth. Unlike in CMP gathers, however, events in image gathers are in their migrated positions. Shown in Figure 9.5-3a is an image section derived from prestack depth migration of the data associated with the stacked section in Figure 9.5-2c using the initial velocity-depth model in Figure 9.5-2a. Superimpose the depth horizons from the velocity-depth model onto the image section and make some minor adjustments where necessary by re-interpreting the depth horizons (Figure 9.5-3b). Selected image gathers associated with the depth image from prestack depth migration are shown in Figure 9.5-4.

How do we make use of the image gathers to update the initial velocity-depth model? First, consider applying to a CMP gather conventional stacking velocity analysis. Compute the velocity spectrum (velocity analysis) and pick a velocity function. Following the normal-moveout correction of the CMP gather using this velocity function, events should look flat if the velocity function had been picked correctly. If the picking was done incorrectly, then you would observe events with residual moveout. In principle, this residual moveout can be computed and used to update the initially picked velocity function.

Analogous to the conventional stacking velocity analysis, if the initial velocity-depth model has been estimated with sufficient accuracy, then the image gathers derived from prestack depth migration using this model should exhibit flat events associated with the layer boundaries included in the model. Any errors in layer velocities and/or reflector geometries, on the other hand, should give rise to residual moveout along those events on the image gathers. Again, in principle, this residual moveout can be determined and used for model updating as illustrated in Figure 9.5-5. The image gather in Figure 9.5-5a exhibits events with residual moveout represented by the purple trajectories. Assume that the residual moveout is parabolic and compute the semblance spectrum as shown in Figure 9.5-5b. The horizontal axis of the semblance plane represents the depth error and the vertical axis represents the depth of the event. The semblance spectrum has two quadrants that correspond to positive and negative residual moveouts, or equivalently, to positive and negative depth errors. A flat event would yield a semblance peak that coincides with the vertical axis with zero depth error, whereas an event with residual moveout would yield a semblance peak situated either in the left or right quadrant depending on the sign of the depth error.

Much like picking a velocity function from a conventional stacking velocity semblance spectrum, a vertical function that represents the depth-dependent residual moveout can be picked from the semblance spectrum in Figure 9.5-5b. This function can then be used to correct for the residual moveout as shown in Figure 9.5-5c. The actual steps to make the residual moveout correction are as follows:

Extract the interval velocity function from the velocity-depth model at the image-gather location where the residual moveout analysis is to be done.
Convert the image gather from depth to time using the interval velocity function.
Assume that the residual moveout of events on the image gather in time is parabolic and compute the semblance spectrum for a range of negative and positive moveouts.
Apply the residual moveout correction to the image gather.
Compute a new rms velocity function from the results of the residual moveout analysis.
Compute a new interval velocity function (equation 1) from the updated rms velocity function at the image-gather location.
Convert the image gathers back to depth using the new interval velocity functions.
Finally, update the velocity-depth model using the new interval velocity functions.

Figure 9.5-1 (a) A poststack time-migrated stacked section with a set of interpreted horizons, (b) the horizon-consistent rms velocity profiles, (c) the interval profiles derived from Dix conversion of the rms velocity profiles.
Figure 9.5-2 (a) The velocity-depth model created by using a combination of Dix conversion to estimate layer velocities (Figure 9.5-1c) and image-ray depth conversion of the time horizons interpreted from the time-migrated stacked section (Figure 9.5-1a), (b) poststack depth migration with depth horizons as in (a), (c) modeled zero-offset traveltimes overlayed on the unmigrated stacked section.
Figure 9.5-3 (a) Prestack depth migration using the velocity-depth model in Figure 9.5-2a, (b) the same image section with depth horizons as in Figure 9.5-2a.
Figure 9.5-4 Selected image gathers associated with the prestack depth-migrated section in Figure 9.5-3b and residual moveout semblance spectra computed from these gathers.
Figure 9.5-5 (a) An image gather associated with the prestack depth-migrated section in Figure 9.5-3b, (b) residual moveout spectrum computed from (a), (c) the image gather in (a) after residual moveout correction, (d) the residual moveout section derived from the semblance spectra shown in Figure 9.5-6.
Figure 9.5-6 Residual moveout semblance spectra computed from the image gathers as in Figure 9.5-4 along the depth horizons shown in Figure 9.5-3b. The velocity-depth model used in prestack depth migration is shown in Figure 9.5-2a.

In principle, residual moveout analysis can be carried out for image gathers at some spatial interval. The residual moveout spectra computed from the image gathers displayed in Figure 9.5-4 indicate varying degrees of errors in the initial velocity-depth model shown in Figure 9.5-2a. An alternative way to measure the residual moveout is by computing it along the depth horizons themselves. Shown in Figure 9.5-6 are the horizon-consistent residual moveout semblance spectra. The vertical axis of the semblance plane represents the depth error which corresponds to either positive or negative residual moveout.

Pick the residual moveout profiles from the semblance spectra as shown in Figure 9.5-6 and combine them with the depth horizons to create the residual moveout section shown in Figure 9.5-5d. This section gives an impression of how much residual moveout, thus the range of errors in the initial velocity-depth model, is present on image gathers as a function of depth and distance along the line.

Figure 9.5-7a shows the interval velocity profiles before (as in Figure 9.5-1c) and after the residual moveout update, and Figure 9.5-7b shows the updated depth horizons displayed on top of the initial velocity-depth model as in Figure 9.5-2a. When combined, the updated interval velocity profiles and depth horizons yield the updated velocity-depth model shown in Figure 9.5-8a.

Following the update, the new model needs to be checked for consistency with the input seismic data. Figure 9.5-8b shows the image section derived from prestack depth migration using the updated model in Figure 9.5-8a. Note that the depth horizons associated with the updated model, when superimposed on the image section, coincide with the reflectors that correspond to the layer boundaries. The modeled zero-offset reflection traveltimes using the updated model also are in good agreement with the observed traveltimes of the events on the unmigrated stacked section that are associated with the layer boundaries included in the velocity-depth model of Figure 9.5-8a.

To further verify the accuracy of the updated model, compare the residual moveout semblance spectra computed from the image gathers at selected locations along the line after the model update (Figure 9.5-9) with those before the update (Figure 9.5-4). Note that most of the semblance peaks are now positioned along the vertical axis of the semblance spectra that corresponds to zero residual moveout. The horizon-consistent residual moveout spectra after the model update are shown in Figure 9.5-10. Superimposed on these spectra are the residual moveout profiles that were picked from the spectra before the update (Figure 9.5-6). Note that the residual moveout errors have been reduced significantly after the model update.

The residual moveout analysis of image gathers and the update of velocity-depth model should be performed iteratively until the velocity-depth model and the depth image are consistent ^[8]. Consistency may be achieved with only a few iterations for cases where residual moveouts are small. Sometimes, consistency is never achieved even after several iterations. This often occurs when the initial residual moveouts are caused largely by significant errors in the initial velocity-depth model. An erroneous initial estimate most likely is a result of rapid lateral velocity variations less than a spread length (models with horizontal layers). In general, updating velocity-depth models based on residual moveout analysis of image gathers yields acceptable results for moderately complex structures associated with compressional and extensional tectonics. However, it may not be suitable for complex overburden structures associated with overthrust or salt tectonics.

Reflection traveltime tomography

Reflection traveltime tomography is based on perturbing the initial model parameters by a small amount and then matching the change in traveltimes to the traveltime measurements made from residual moveout analysis of image gathers ^[18] ^[19]. A mathematical treatment of the subject is given in Section J.6. Here, we remind ourselves of the underlying assumptions, outline the theory, and examine the performance of the method with model experiments. We shall complement the discussion on tomography with a field data example.

We must do the best we can in building an accurate earth model in depth (model building) so that only small changes remain to be made to the model by tomography. Specifically, a tomographic update can be expected to work provided the changes to be made to the initial earth model parameters in terms of slownness and depths at layer boundaries are small compared to the model parameters themselves. Additionally, we shall assume that the initial model is made up of horizontal layers with laterally invariant model parameters.

In the usual implementation of reflection traveltime tomography, the model parameters are perturbed while preserving the offset values of the seismic data (Section J.6). The tomographic update Δp to the model parameters that comprise the changes in the slowness and depths to layer boundaries is given by the generalized linear inversion (GLI) solution (equation J-88 in Section J.6)

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

(10)

where Δt denotes the column vector that represents the residual moveout times measured from the image gathers, L is a sparse matrix — its elements are in terms of the slowness and depth parameters associated with the initial model, and T denotes matrix transposition.

Figure 9.5-7 Model updating by residual moveout corrections: (a) The interval velocity profiles as in Figure 9.5-1c before (thin curves) and after (thick curves) residual moveout corrections, (b) the depth horizons in Figure 9.5-2a after an update with the velocity-depth model before update as in the same figure shown in the background.
Figure 9.5-8 (a) The velocity-depth model in Figure 9.5-2a after model updating by residual moveout corrections, (b) prestack depth migration with depth horizons as in (a), (c) modeled zero-offset traveltimes overlayed on the unmigrated stacked section.
Figure 9.5-9 Selected image gathers associated with the prestack depth-migrated section in Figure 9.5-8b and residual moveout semblance spectra computed from these gathers.
Figure 9.5-10 Residual moveout semblance spectra computed from the image gathers as in Figure 9.5-9 along the depth horizons shown in Figure 9.5-8b. The velocity-depth model used in prestack depth migration is shown in Figure 9.5-8a.

Consider the earth model with horizontal layers shown in Figure 9.1-1. We shall make an attempt to update the initial estimate shown in Figure 9.1-5 using reflection tomography.

Generate a set of image gathers (Figure 9.5-4) from prestack depth migration (Figure 9.5-3) using an initial velocity-depth model (Figure 9.5-2a).
Convert the image gathers from depth to time using the interval velocity functions extracted from the initial velocity-depth model at the image-gather locations.
Compute the horizon-consistent residual moveout for all offsets along events on image gathers that correspond to the layer boundaries included in the model. Figure 9.5-6 shows the residual moveout spectra along the six horizons of the model in Figure 9.5-2a. The vertical axis represents the residual moveout measured at a reference offset, usually the maximum offset. The horizontal axis represents the CMP locations along the line. Since residual moveout can be either negative or positive, the vertical axis is in both the positive and negative directions.
Pick the residual moveout profiles for all the horizons by tracking the semblance peaks. Any departure from the horizontal axis indicates a nonzero value for residual moveout.
Build the traveltime error vector Δt using the residual moveout times. As an example, you may have 10 layers, 1000 CMPs with a fold of 30. This means that the length of the traveltime error vector Δt is 300 000.
Define the initial model by a set of slowness and depth parameters and construct the coefficient matrix L in equation (10) (Section J.6). As an example, you may have 10 layers and each layer may be defined by 50 pairs of slowness and depth values in the lateral direction. This means you would have 1000 parameters in your model space. It also means that for the example given in step (b), you have 300 000 equations to solve for 1000 parameters. The solution for this overdetermined system is given by equation (10).
Estimate the change in parameters vector Δp, by way of the GLI solution given by equation (10).
Update the parameter vector p + Δp. Figure 9.5-11a shows the interval velocity profiles before and after the tomographic update and Figure 9.5-11b shows the updated depth horizons superimposed on the initial velocity-depth model of Figure 9.5-2a. In a tomographic update, you may wish to perturb a subset of layer velocities and/or reflector geometries. This depends on your confidence in the initial model parameters for the layers and the quality of the residual moveout profiles to be used in inversion.

By combining the updated interval velocity profiles (Figure 9.5-11a) with the new depth horizons (Figure 9.5-11b), we obtain the velocity-depth model after the first iteration of tomographic update (Figure 9.5-12a). Following the update, we check for consistency of the new model with the input seismic data. Overlay the depth horizons from the updated model in Figure 9.5-12a onto the image section derived from prestack depth migration shown in Figure 9.5-12b and note that they coincide with the reflectors associated with the layer boundaries included in the model. Also, the modeled zero-offset reflection traveltimes using the updated model coincide with the observed traveltimes of the events on the unmigrated stacked section (Figure 9.5-12c) that are associated with the layer boundaries included in the velocity-depth model of Figure 9.5-12a.

Next, compute the residual moveout semblance spectra from the image gathers at selected locations along the line after the model update (Figure 9.5-13) and compare them with those before the update (Figure 9.5-4). While most of the semblance peaks are now aligned with the vertical axis of the spectra, the tomographic update may be repeated to further remove any remaining residual moveout errors. The horizon-consistent residual moveout semblance spectra after the first iteration of the tomographic model update are shown in Figure 9.5-14. The residual moveout profiles that were picked from the spectra before the update (Figure 9.5-6) have been overlayed on the spectra after the update to observe the extent of the removal of residual moveout errors by the tomographic update.

Proceed to a second iteration of tomographic update by picking a new set of residual moveout profiles from the spectra shown in Figure 9.5-14. Then, follow the updating procedure described above to obtain a new set of interval velocity profiles (Figure 9.5-15a) and depth horizons (Figure 9.5-15b). Note that the changes in interval velocities and depth horizons that result from the second iteration are smaller compared to those from the first iteration (Figure 9.5-11).

Combine the new set of interval velocity profiles and depth horizons to create the next update of the velocity-depth model shown in Figure 9.5-16a. Verify the updated model by performing prestack depth migration to generate an image section (Figure 9.5-16b) and zero-offset modeling of traveltimes associated with the layer boundaries included in the model (Figure 9.5-16c). Next, compute the residual moveout semblance spectra to examine any remaining moveout errors (Figure 9.5-17). Finally, compute the horizon-consistent residual moveout semblance spectra and overlay the residual moveout profiles from the first iteration (Figure 9.5-18). Note that the changes from the first to the second iteration are marginal. At this point, you may wish to end the iterations for tomographic update.

Figure 9.5-11 First iteration of model updating by reflection tomography: (a) The interval velocity profiles as in Figure 9.5-1c before (thin curves) and after (thick curves) tomographic updating, (b) the depth horizons in Figure 9.5-2a after the update, with the velocity-depth model before the update as in Figure 9.5-2a shown in the background.
Figure 9.5-12 (a) The velocity-depth model in Figure 9.5-2a after model updating by reflection tomography, (b) prestack depth migration with depth horizons as in (a), (c) modeled zero-offset traveltimes overlayed on the unmigrated stacked section.
Figure 9.5-13 Selected image gathers associated with the prestack depth-migrated section in Figure 9.5-12b and residual moveout semblance spectra computed from these gathers.
Figure 9.5-14 Residual moveout semblance spectra computed from the image gathers as in Figure 9.5-13 along the depth horizons shown in Figure 9.5-12b. The velocity-depth model used in prestack depth migration is shown in Figure 9.5-12a.

The extent to which an initial velocity-depth model is perturbed by a tomographic update depends on the accuracy of that initial model, which in turn, depends on how it has been estimated. Consider the two initial velocity-depth models shown in Figures 9.4-27a and 9.5-2a, which, for convenience, will be referred to as Models 1 and 2. Model 1 in Figure 9.4-27a was built by applying the layer-by-layer inversion strategy with a combination of coherency inversion and poststack depth migration to estimate the layer velocities and delineate the reflector geometries. Model 2 in Figure 9.5-2a was built by applying the time-to-depth conversion strategy with a combination of Dix conversion to estimate the layer velocities and image-ray depth conversion to delineate the reflector geometries.

As with the model in Figure 9.5-2a, perform prestack depth migration using Model 1 to obtain the image section shown in Figure 9.5-19a. Then, compute the residual moveout semblance spectra for selected image gathers (Figure 9.5-20). By comparing these spectra with those in Figure 9.5-4, it may not be obvious as to whether Model 1 or Model 2 requires more perturbation. However, a comparison of the horizon-consistent residual moveout semblance spectra shown in Figures 9.5-6 and 9.5-21 reveals that Model 1 gives rise to less residual moveout on image gathers than does Model 2. As such, the tomographic update of Model 1 yields changes in interval velocities and reflector geometries (Figure 9.5-22) that are smaller than those for Model 2 (Figure 9.5-7).

Irrespective of the strategy followed to derive the initial model, the results of model updating need to be verified for consistency with the input seismic data (Figure 9.5-23) and examined for any remaining residual moveouts (Figure 9.5-24) to decide whether or not to continue with the iterations of tomographic update.

Limitations in resolving velocity-depth ambiguity by tomography

Reflection traveltime tomography can be surprisingly successful in updating an initial model with significant errors. Consider, for instance, the velocity-depth model in Figure 9.5-25 with constant layer velocities. This model is indeed practically consistent with the input seismic data as verified by the depth image in Figure 9.5-25b and the modeled zero-offset traveltimes overlayed on the unmigrated stacked section in Figure 9.5-25c.

Perform prestack depth migration using the initial model to derive the image section in Figure 9.5-26a and compute the residual moveout semblance spectra of the selected image gathers in Figure 9.5-27. Note the significant moveout errors manifested by the semblance peaks which are off the zero moveout centerline on the spectra. Now, compute the horizon-consistent residual moveout semblance spectra along the depth horizons shown in Figure 9.5-26b and note the significant departures from the zero moveout line (Figure 9.5-28). Pick the residual moveout profiles from these spectra and perform a tomographic update to obtain a new set of interval velocity profiles and depth horizons as shown in Figure 9.5-29.

Combine the updated interval velocities and reflector geometries to create the updated velocity-depth model shown in Figure 9.5-30a. Note that this model is fairly close to the updated models derived from the applications of time-to-depth conversion (Figure 9.5-16a) and layer-by-layer inversion (Figure 9.5-23a) strategies. Model verification tests (Figures 9.5-30b,c) and residual moveout semblance spectra (Figures 9.5-31 and 9.5-32) demonstrate the feasibility and accuracy of the updated model in Figure 9.5-30a.

As demonstrated by the extreme case of an erroneous initial model (Figure 9.5-25), a tomographic update can steer the model toward an acceptable final model. How far, though, can we go with tomography to obtain a finely accurate model? In the following model experiments, we shall examine this fundamentally important question regarding tomography.

Figure 9.5-33 shows horizon-consistent residual moveout semblance spectra derived from image gathers that were created by prestack depth migration of the synthetic data as in Figure 9.1-2a using the initial velocity-depth model shown in Figure 9.1-5b. Shown in Figure 9.5-34 are the results of tomographic updating — the interval velocity profiles before (thin curves) and after (thick curves) the tomographic update, and the updated velocity-depth model itself. Compare with the initial velocity-depth model shown in Figure 9.1-5b and the true velocity-depth model shown in Figure 9.1-1b. Compare with the true model (Figure 9.1-1) and the initial model (Figure 9.1-5) and note that tomographic update has produced a result fairly close to the true model. A close-up of the center portions of the initial model, the final model, and the true model after the tomographic update is shown in Figure 9.5-35.

Figure 9.5-15 Second iteration of model updating by reflection tomography: (a) The interval velocity profiles as in Figure 9.5-11a before (thin curves) and after (thick curves) tomographic updating, (b) the depth horizons in Figure 9.5-12a after update, with the velocity-depth model before the update as in Figure 9.5-12a shown in the background.
Figure 9.5-16 (a) The velocity-depth model in Figure 9.5-12a after model updating by reflection tomography, (b) prestack depth migration with depth horizons as in (a), (c) modeled zero-offset traveltimes overlayed on the unmigrated stacked section.
Figure 9.5-17 Selected image gathers associated with the prestack depth-migrated section in Figure 9.5-16b and residual moveout semblance spectra computed from these gathers.
Figure 9.5-18 Residual moveout semblance spectra computed from the image gathers as in Figure 9.5-17 along the depth horizons shown in Figure 9.5-16b. The velocity-depth model used in prestack depth migration is shown in Figure 9.5-16a.
Figure 9.5-19 (a) Prestack depth migration using the velocity-depth model in Figure 9.4-27a, (b) the same image section with depth horizons as in Figure 9.4-27a.
Figure 9.5-20 Selected image gathers associated with the prestack depth-migrated section in Figure 9.5-19b and residual moveout semblance spectra computed from these gathers.
Figure 9.5-21 Residual moveout semblance spectra computed from the image gathers as in Figure 9.5-20 along the depth horizons shown in Figure 9.5-19b. The velocity-depth model used in prestack depth migration is shown in Figure 9.4-27a.
Figure 9.5-22 Model updating by reflection tomography: (a) The interval velocity profiles as in Figure 9.4-22 before (thin curves) and after (thick curves) tomographic update, (b) the depth horizons in Figure 9.4-27a after, with the velocity-depth model before the update as in Figure 9.5-27a shown in the background.
Figure 9.5-23 (a) The velocity-depth model in Figure 9.4-27a after model updating by reflection tomography, (b) prestack depth migration with depth horizons as in (a), (c) modeled zero-offset traveltimes overlayed on the unmigrated stacked section.
Figure 9.5-24 Residual moveout semblance spectra computed from the image gathers associated with the prestack depth-migrated section in Figure 9.5-23b along the depth horizons shown in the same figure. The velocity-depth model used in prestack depth migration is shown in Figure 9.5-23a.
Figure 9.5-25 (a) A velocity-depth model created by iterative poststack depth migration using constant layer velocities, (b) poststack depth migration with depth horizons as in (a), (c) modeled zero-offset traveltimes overlayed on the unmigrated stacked section.
Figure 9.5-26 (a) Prestack depth migration using the velocity-depth model in Figure 9.5-25a, (b) the same image section with depth horizons as in Figure 9.5-25a.
Figure 9.5-27 Selected image gathers associated with the prestack depth-migrated section in Figure 9.5-26b and residual moveout semblance spectra computed from these gathers.
Figure 9.5-28 Residual moveout semblance spectra computed from the image gathers as in Figure 9.5-27 along the depth horizons shown in Figure 9.5-26b. The velocity-depth model used in prestack depth migration is shown in Figure 9.5-26a.

Iterate the updating process until the residual moveout for all the layers is minimized. Specifically, perform prestack depth migration using the tomographically updated velocity-depth model in Figure 9.5-34 and generate a new set of image gathers. Compute a new set of residual moveout semblance spectra from the image gathers as shown in Figure 9.5-36. Compare with the residual moveout spectra in Figure 9.5-33 before the tomographic update and note that the residual moveout after the update has been reduced significantly.

We now test the tomographic update in the presence of large errors in the initial velocity-depth model. Consider the velocity-depth model shown in Figure 9.1-1 based on the parameters listed in Table 9-2. We shall deliberately introduce errors to this model by setting the velocities for layer H3 to a constant value of 2550 m/s and for layer H4 to a constant value of 3250 m/s. The erronenous velocity-depth model is shown in Figure 9.5-37.

By adopting the distorted model in Figure 9.5-37 as the initial velocity-depth model, perform prestack depth migration, generate image gathers and compute residual moveout semblance spectra as shown in Figure 9.5-38. Note the large moveouts caused by the velocity errors associated with layers H3 and H4. The tomographic update shown in Figure 9.5-39a has modified the velocities for layers H3 and H4, and for the underlying layers. Specifically, the constant-velocities assigned to layers H3 and H4 have been replaced with laterally varying velocity profiles that are consistent with the actual velocity profiles shown in Figure 9.1-1. Nevertheless, compare the profiles in Figures 9.1-1a and 9.5-39a and note that the velocity estimate for the underlying layer H5 has been adversely influenced by the tomographic update applied to the layers above.

A close-up of the central portions of true velocity-depth model (Figure 9.1-1b), the initial velocity-depth model (Figure 9.5-37b), and the updated velocity-depth model (Figure 9.5-39b) is shown in Figure 9.5-40. Although the errors in the initial velocity-depth model for the shallow layers have been reduced, note that, in general, starting with an initial velocity-depth model with large errors, tomographic updating does not always produce satisfactory results.

The fact that the tomographic update for one layer is influenced by the errors in the parameters for the layers above is best illustrated by the model in Figure 9.1-18. This is the model with horizontal layers (Table 9-2) but with a near-surface layer H1 with laterally varying velocities between 800 m/s and 1500 m/s. Results of earth model estimation using Dix conversion are shown in Figure 9.1-21c. Consider this model as the initial velocity-depth model and perform prestack depth migration of the data associated with the stacked section in Figure 9.1-19a, generate image gathers, and compute the residual moveout semblance spectra shown in Figure 9.5-41. Pick the spectra to obtain the residual moveout profiles by tracking the semblance peaks, then use these profiles in a tomographic update to get the new model shown in Figure 9.5-42. Use the resulting model as input to a second iteration for model updating. Compare the resulting residual moveout spectra shown in Figure 9.5-43 with those from the first iteration shown in Figure 9.5-41 and note that the second iteration has reduced the residual moveout, significantly. The updated model from the second iteration is shown in Figure 9.5-44. A close-up of the central portions of the intermediate velocity-depth model (Figure 9.5-42b), the velocity-depth model updated for the second time (Figure 9.5-44b), and the true velocity-depth model (Figure 9.1-18b) is shown in Figure 9.5-45. Compare these estimates with the close-up view (Figure 9.1-22b) of the initial velocity-depth model shown in Figure 9.1-21c. The lesson to learn from this experiment is not to start with an initial model (Figure 9.1-21c) with large errors and expect a tomographic update to correct for these errors.

Finally, we shall update the model shown in Figure 9.1-24 which was estimated by coherency inversion. Consider this model as the initial velocity-depth model and perform prestack depth migration, generate image gathers, and compute the residual moveout spectra shown in Figure 9.5-46. Pick the spectra, again, by tracking the semblance peaks to obtain the residual moveout profiles. Then use these profiles in a tomographic update to get the new model shown in Figure 9.5-47. A close-up of the central portions of the initial velocity-depth model (Figure 9.1-24b), the updated velocity-depth model (Figure 9.5-47b), and the true velocity-depth model (Figure 9.1-18b) is shown in Figure 9.5-48.

Based on the results of the experiments using Dix conversion and coherency inversion applied to the model with a near-surface layer that has lateral velocity variations (Figure 9.1-18), again, the lesson to learn is that, whatever the strategy used in the model building, if you start with an initial model with large errors, do not expect tomographic updating to correct for these errors.

Figure 9.5-29 Model updating by reflection tomography: (a) The interval velocity profiles before (thin curves) and after (thick curves) tomographic update, (b) the depth horizons in Figure 9.5-25a after update, with the velocity-depth model before the update as in Figure 9.5-25a shown in the background.
Figure 9.5-30 (a) The velocity-depth model in Figure 9.5-25a after model updating by reflection tomography, (b) prestack depth migration with depth horizons as in (a), (c) modeled zero-offset traveltimes overlayed on the unmigrated stacked section.
Figure 9.5-31 Selected image gathers associated with the prestack depth-migrated section in Figure 9.5-30b and residual moveout semblance spectra computed from these gathers.
Figure 9.5-32 Residual moveout semblance spectra computed from the image gathers as in Figure 9.5-31 along the depth horizons shown in Figure 9.5-30b. The velocity-depth model used in prestack depth migration is shown in Figure 9.5-30a.
Figure 9.5-33 Horizon-consistent residual moveout semblance spectra derived from image gathers that were created by prestack depth migration of the synthetic data as in Figure 9.1-2a using the initial velocity-depth model shown in Figure 9.1-5b.
Figure 9.5-34 (a) The interval velocity profiles before (thin curves) and after (thick curves) tomographic update; (b) the updated velocity-depth model. Compare with the initial velocity-depth model shown in Figure 9.1-5b and the true velocity-depth model shown in Figure 9.1-1b.
Figure 9.5-35 Central portions of (a) the estimated initial velocity-depth model shown in Figure 9.1-5b using the smoothed interval velocity profiles shown in Figure 9.1-5a, (b) the updated velocity-depth model shown in Figure 9.5-34b, and (c) the true velocity-depth model shown in Figure 9.1-1b.
Figure 9.5-36 Horizon-consistent residual moveout semblance spectra derived from image gathers that were created by prestack depth migration of the synthetic data as in Figure 9.1-2a using the updated velocity-depth model shown in Figure 9.5-34b.
Figure 9.5-37 (a) The velocity profiles as in Figure 9.1-1a except for horizons H3 and H4 which deliberately have been assigned wrong velocities; (b) the velocity-depth model created from the interval velocity profiles in (a).
Figure 9.5-38 Horizon-consistent residual moveout semblance spectra derived from image gathers that were created by prestack depth migration of the synthetic data as in Figure 9.1-2a using the initial velocity-depth model shown in Figure 9.5-37b.
Figure 9.5-39 (a) The interval velocity profiles before (thin curves) and after (thick curves) tomographic update; (b) the updated velocity-depth model. Compare with the initial velocity-depth model shown in Figure 9.5-37b and the true velocity-depth model shown in Figure 9.1-1b.
Figure 9.5-40 Central portions of (a) the true velocity-depth model shown in Figure 9.1-1b, (b) the initial velocity-depth model shown in Figure 9.5-37b, and (c) the updated velocity-depth model shown in Figure 9.5-39b.
Figure 9.5-41 Horizon-consistent residual moveout semblance spectra derived from image gathers that were created by prestack depth migration of the synthetic data as in Figure 9.1-19a using the initial velocity-depth model shown in Figure 9.1-21c.
Figure 9.5-42 (a) The interval velocity profiles before (thin curves) and after (thick curves) the first tomographic update; (b) the updated velocity-depth model. Compare with the initial velocity-depth model shown in Figure 9.1-21c and the true velocity-depth model shown in Figure 9.1-18b.
Figure 9.5-43 Horizon-consistent residual moveout semblance spectra derived from image gathers that were created by prestack depth migration of the synthetic data as in Figure 9.1-19a using the updated velocity-depth model shown in Figure 9.5-42b.
Figure 9.5-44 (a) The interval velocity profiles before (thin curves) and after (thick curves) the second tomographic update; (b) the velocity-depth model updated for the second time. Compare with the initial velocity-depth model shown in Figure 9.1-21c, the intermediate velocity-depth model shown in Figure 9.5-42b, and the true velocity-depth model shown in Figure 9.1-18b.
Figure 9.5-45 Central portions of (a) the intermediate velocity-depth model shown in Figure 9.5-42b, (b) the velocity-depth model updated for the second time shown in Figure 9.5-44b, and (c) the true velocity-depth model shown in Figure 9.1-18b.

We now examine reflection traveltime tomography for its ability to update an initial model in the presence of a complex overburden. Figure 9.5-49a shows the rms velocity profiles along the time horizons in Figure 9.4-28a. Dix conversion yields the interval velocity profiles shown in Figure 9.5-49b. Admittedly, some editing and smoothing before and after Dix conversion were done to tame the interval velocities to behave in a structurally sensible manner.

By using the interval velocity profiles shown in Figure 9.5-49b, the time horizons interpreted from the time-migrated section shown in Figure 9.4-28a were converted to depth by way of image rays. Then, these depth horizons were combined with the interval velocity profiles (Figure 9.5-49b) to create an initial velocity-depth model. This was followed by depth migration of the stacked section in Figure 9.4-28b using the initial velocity-depth model, which in turn, was revised based on the interpretation of the depth-migrated section. A second iteration of poststack depth migration (Figure 9.5-50a) produced a velocity-depth model (Figure 9.5-49c) that is consistent with the reflection times observed on the unmigrated stacked section (Figure 9.5-50b).

Actually, the now consistent velocity-depth model in Figure 9.5-49c is what we shall update by way of reflection tomography. To begin with, perform prestack depth migration and generate image gathers as in Figure 9.5-51. Compute the residual moveout semblance spectra and observe that there are semblance peaks that are off the center line, indicating the presence of residual moveout. Stacking of the image gathers yields the image section from prestack depth migration shown in Figure 9.5-52. The significant improvement in imaging the diapiric structure and the base-salt reflector (horizon H8 as in Figure 9.5-50a) becomes obvious when this section is compared with the poststack depth-migrated section in Figure 9.5-50a. Nevertheless, the fact that the semblance spectra (Figure 9.5-51) exhibit residual moveout along the events associated with the layer boundaries included in the velocity-depth model suggests that the image from prestack depth migration may be improved by updating the initial velocity-depth model.

Now here is a dilemma that we have to live with when updating velocity-depth models. As a result of the inevitable muting of the image gathers (Figure 9.5-51), there is often insufficient offset range to accurately measure residual moveouts for shallow horizons. Because of the very poor residual moveout spectra, shallow horizons can be unreliable in model updating. Errors in the model parameters for the shallow layers will of course adversely affect the estimation and updating of the model parameters for the deeper layers. At greater depths, although there is sufficient cable length, faster velocities prevent recognition of velocity or depth errors with adequate resolution.

There is of course some mid-range of depths and velocities for which residual moveout semblance spectra can be used for model updating. In the present example, residual moveout semblance spectra associated with horizons H5 to H8 (as denoted in Figure 9.5-50a) were computed and picked as input to the reflection traveltime tomography solver (Section J.6). When picking the residual moveout profiles from the semblance spectra in Figure 9.5-52, keep in mind not to track all the rapid fluctuations that are much less than a cable length. These are the manifestations of the lateral velocity variations in the layers above as was demonstrated by the experiments conducted using synthetic data in models with horizontal layers.

The workflow for tomographic updating — prestack depth migration to generate image gathers, computing residual moveout along events associated with the layer boundaries, picking residual moveout profiles and using them as input to the tomography to derive the updates for layer velocities (Figure 9.5-53a) and reflector geometries (Figure 9.5-53b), were repeated twice to produce the final model and image shown Figure 9.5-54. In each iteration, the first four layers were kept unchanged because of poor estimates for residual moveout along the horizons that correspond to the base of these layers (H1 through H4). Compare the images from prestack depth migrations shown in Figures 9.5-52a and 9.5-54a obtained from the initial and final velocity-depth models of Figures 9.5-52b and 9.5-54b, respectively, and note that the differences are marginal.

Iterative tomographic updating, as demonstrated by the model experiments in this section, does not necessarily result in convergence of an initial model to a final model with zero residual moveout. Instead, the solution may just wobble and never converge. This scenario is especially likely for cases of complex overburden structures as in the present example. In such cases, additional information; such as well control, check-shot velocities, and geologic constraints, is required to derive an acceptable final model.

Turning-ray tomography

Rapid discharge of the sediments by the Mississippi River has given rise to an accumulation of gas-charged mudflows within complex, meandering channels down the slopes of its delta. Velocities within these mudflows are significantly lower than the surrounding deltaic sediments and can be as low as 300 m/s. The underlying predeltaic Holocene sequence also is characterized by sediments with lateral velocity variations ^[20] ^[21].

Because of the absence of a near-surface refractor with a strong velocity contrast, first arrivals on shot gathers recorded over the Mississippi Delta often do not represent the typical refracted arrivals associated with a head wave. Therefore, neither refraction traveltime tomography (Section C.9) is applicable to estimate a near-surface model nor refraction statics (refraction statics corrections) may be a rigorous solution to resolving the near-surface complexity in the Mississippi Delta.

Figure 9.5-46 Horizon-consistent residual moveout semblance spectra derived from image gathers that were created by prestack depth migration of the synthetic data as in Figure 9.1-19a using the initial velocity-depth model shown in Figure 9.1-24b.
Figure 9.5-47 (a) The interval velocity profiles before (thin curves) and after (thick curves) the first tomographic update; (b) the updated velocity-depth model. Compare with the initial velocity-depth model shown in Figure 9.1-24b and the true velocity-depth model shown in Figure 9.1-18b.
Figure 9.5-48 Central portions of (a) the initial velocity-depth model shown in Figure 9.1-24b using the interval velocity profiles shown in Figure 9.1-24a, (b) the updated velocity-depth model shown in Figure 9.5-47b using the interval velocity profiles shown in Figure 9.5-47a, and (c) the true velocity-depth model shown in Figure 9.1-18b.
Figure 9.5-49 (a) RMS velocity profiles associated with the data shown in Figure 9.4-28a, (b) the interval velocity profiles computed from (a) by way of Dix conversion, and (c) an initial velocity-depth model derived from the combination of the interval velocity profiles in (b) with the depth horizons that were created by image-ray depth conversion of the time horizons shown in Figure 9.4-28a.
Figure 9.5-50 (a) Poststack depth migration of the stacked section in (b) using the initial velocity-depth model in Figure 9.4-49c, (b) zero-offset traveltimes modeled from the depth horizons in (a) and superimposed on the stacked section.
Figure 9.5-51 Part 1: Selected image gathers associated with the prestack depth-migrated section in Figure 9.5-52 and residual moveout semblance spectra computed from these gathers.
Figure 9.5-51 Part 2: Selected image gathers associated with the prestack depth-migrated section in Figure 9.5-52 and residual moveout semblance spectra computed from these gathers.
Figure 9.5-52 Prestack depth migration of the data associated with the stacked section shown in Figure 9.5-50b using the initial velocity-depth model shown in Figure 9.5-49c, and (b) the residual moveout semblance spectra computed from the image gathers associated with the prestack depth-migrated section in (a) along Horizons H5-H8 labeled as in Figure 9.5-50a.
Figure 9.5-53 (a) The interval velocity profiles before (thin curves) and after (thick curves) the second tomographic update, and (b) the depth horizons after the update superimposed on the intermediate velocity-depth model derived from the first iteration of tomographic update.
Figure 9.5-54 (a) The final velocity-depth model after the second tomographic update, and (b) prestack depth migration of the data associated with the stacked section shown in Figure 9.5-50b using the final velocity-depth model in (a).
Figure 9.5-55 (a) A map of statics derived from turning-ray tomography applied to a 3-D seismic data set, (b) an inline stack along the traverse indicated by the horizontal line in (a) with refraction statics corrections, and (c) the same inline stack with turning-ray tomographic statics corrections. ^[22]; courtesy Geosignal, Western Geophysical.

Instead, the first arrivals are often associated with diving waves through the deltaic and predeltaic Holocene sediments ^[23]. Because of the unusual velocity gradients within the near-surface layers, the downgoing incident wave rapidly turns around before being reflected and is recorded by the receiver as the first arrival.

Just as reflection traveltime tomography (Section J.6) can be used to update an initial estimate of a subsurface velocity-depth model, turning-ray tomography may be used to update an initial estimate of a near-surface velocity-depth model ^[23] ^[24].

Begin with the picking of the first arrivals that represent the diving waves through the near surface.
Define an initial velocity-depth model by a set of near-surface layers with constant velocities and thicknesses, and model the first-arrival times by ray tracing.
Compute the difference between the modeled and observed first-arrival times.
Estimate the change in parameters vector Δp of equation (10), by way of the GLI solution given by equation (J-88) by perturbing the velocities of the near-surface layers only.
Update the parameter vector p + Δp to obtain a new near-surface velocity-depth model.
Iterate steps (a) through (e) as necessary to minimize the discrepancy between the modeled and actual first-arrival times.

The final velocity-depth model resulting from the iterative application of turning-ray tomography is then used to compute the one-way traveltimes through the near-surface model along vertical raypaths. These are then used to apply the necessary source and receiver statics corrections to the prestack data. Figure 9.5-55a shows the statics solution derived from turning-ray tomography applied to a 3-D offshore seismic data set from the Mississippi Delta ^[22]. Note the mudflows characterized by the strings of negative statics shifts. Figure 9.5-55b shows a stacked section along an inline traverse with refraction statics corrections (refraction statics corrections), and Figure 9.5-55c shows the stacked section along the same traverse with statics corrections as in (a) based on turning-ray tomography. Note the significant improvement of event continuity in the central part of the section.

Exercises

Exercise 9-1. Consider the hypothetical case of a pure strike-slip fault with the same horizontally situated formation on both sides. Can such an earth model be identified by seismic imaging?

Exercise 9-2. There can be a multiple number of velocity-depth models such as shown in Figure 9.0-1 consistent with a single zero-offset traveltime section (Figure 9.0-2). Can these velocity-depth models also be consistent with a single nonzero-offset traveltime section?

Appendix J: Data modeling by inversion

J.1 The generalized linear inversion

We shall start by making the following formal statement:

Given an observed data set d, estimate a set of parameters p which are used to construct a model d′ of the observed data set d, such that the difference between the observed data set d and the modeled data set d′ is minimum based on a specific norm.

For most geophysical applications, the norm for minimization usually is chosen to be L₂, in which case the error energy in data modeling is made minimum in the least-squares sense. Another norm that is appropriate in some applications is L₁, in which case the magnitude of the error in data modeling is made minimum.

To form a mathematical framework for the objective stated above, we need a model equation that relates the modeled data with the model parameters to be estimated as

\mathbf {d} ^{\prime }=\mathbf {Lp} ,

(J-1)

where d′ is the modeled data vector, p is the model parameter vector, and L is the matrix that relates the modeled data vector to the model parameter vector.

The error vector e is defined as the difference between the modeled data vector and the observed data vector

\mathbf {e} =\mathbf {d} -\mathbf {d^{\prime }} .

(J-2)

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

\mathbf {e} =\mathbf {d} -\mathbf {Lp} .

(J-3)

Following Lines and Treitel ^[25], the least-squares solution for equation (J-3) can be determined. First, the cumulative squared error S is expressed as

S=\mathbf {e^{T}e} ,

(J-4a)

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

S=(\mathbf {d} -\mathbf {Lp} )^{\mathbf {T} }(\mathbf {d} -\mathbf {Lp} ).

(J-4b)

Expand the right-hand side to get

S=\mathbf {d^{T}d} -\mathbf {p^{T}L^{T}d} -\mathbf {d^{T}Lp} +\mathbf {p^{T}L^{T}Lp} .

(J-4c)

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

-\mathbf {d^{T}L} +\mathbf {p^{T}L^{T}L} =\mathbf {0} .

(J-5a)

Apply matrix transpose and rearrange the terms

(\mathbf {L^{T}L} )\mathbf {p} =\mathbf {L^{T}d} ,

(J-5b)

which yields the desired least-squares solution

\mathbf {p} =(\mathbf {L^{T}L} )^{\mathbf {-} 1}\mathbf {L^{T}d} ,

(J-6a)

where L^TL is the covariance matrix and (L^TL)⁻¹L^T is the least-squares (also called generalized linear) inverse of L.

Equation (J-6a) represents the generalized linear inverse (GLI) solution to the parameter vector p. This solution is widely used in many stages of seismic data analysis. Examples include deconvolution (optimum wiener filters), residual statics corrections (residual statics corrections), refraction statics corrections (refraction statics corrections), and the discrete Radon transform (the radon transform). In most applications, the solution needs to be constrained since the covariance matrix L^TL needs to be inverted. The constrained solution is given by

\mathbf {p} =(\mathbf {L^{T}L} +\beta \mathbf {I} )^{\mathbf {-} 1}\mathbf {L^{T}d} ,

(J-6b)

where β is called the damping factor and I is the identity matrix.

In some applications, the generalized linear inverse problem is formulated in the frequency domain. Then, the unconstrained solution is given by

\mathbf {p} =(\mathbf {L^{T\ast }L} )^{\mathbf {-} 1}\mathbf {L^{T\ast }d} ,

(J-7a)

and the constrained solution is given by

\mathbf {p} =(\mathbf {L^{T\ast }L} +\beta \mathbf {I} )^{\mathbf {-} 1}\mathbf {L^{T\ast }d} ,

(J-7b)

where the asterisk denotes complex conjugate, and p, d, and L are complex.

The computationally suitable technique to solve for the parameter vector p is dictated by the properties of the coefficient matrix L; that in turn is influenced by how the parameter vector p is defined. The coefficient matrix L may be sparse with many of its elements being zero as for the residual and refraction statics. It may be near singular with the elements of two rows being very similar in size as for the discrete Radon transform. The covariance matrix L^TL may have the Toeplitz structure with its elements being diagonally symmetric as for deconvolution and the discrete Radon transform. It may also be diagonally dominant as for deconvolution. In geophysical applications, techniques to solve for the parameter vector p in equations (J-6a,J-6b) or (J-7a,J-7b) include Levinson recursion, conjugate gradient, Gauss-Seidel and singular-value decomposition.

J.2 The GLI formalism of deconvolution

Deconvolution is fundamentally a data modeling technique. Specifically, we model a 1-D seismogram for a minimum-phase estimate of the source wavelet, to predict multiples, and ultimately obtain an estimate of white reflectivity series (optimum wiener filters). The mathematical foundation of deconvolution is provided in Appendix B. Here, we shall apply the theory of the generalized linear inversion outlined in the previous section to describe deconvolution within the framework of 1-D seismic inversion. Inversion techniques at higher dimensions will be summarized in the next section.

We shall consider designing a least-squares inverse filter f(t) that converts a wavelet w(t) to a desired form d(t) such that the difference e(t) between the actual output y(t) and the desired output d(t) is minimum in the least-squares sense. The zero-delay unit spike is a special case of the desired output d(t). Other forms of d(t) can also be considered, such as a zero-phase band-limited wavelet. The model equation for deconvolution is given by

y(t)=w(t)\ast f(t)

(J-8)

Consider the discrete form of equation (J-8), with w(t) represented by the m–length time series (w₀, w₁, w₂, …, w_m−1), and f(t) represented by the n–length time series (f₀, f₁, f₂, …, f_n−1).

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

{\begin{pmatrix}y_{0}\\y_{1}\\y_{2}\\\vdots \\y_{m+n-1}\end{pmatrix}}={\begin{pmatrix}w_{0}&0&\ldots &0\\w_{1}&w_{0}&&\\w_{2}&w_{1}&w_{0}&\\\vdots &w_{2}&w_{1}&\ddots \\w_{m-1}&\vdots &w_{2}&\ddots \\0&w_{m-1}&\vdots &\\\vdots &&w_{m-1}&\\0&&&w_{m-1}\end{pmatrix}}{\begin{pmatrix}f_{0}\\f_{1}\\f_{2}\\\vdots \\f_{n-1}\end{pmatrix}}.

(J-9)

This is the equation of complete transient convolution. Define the output vector on the left-hand side by y, the coefficient matrix on the right-hand side by L, and the filter vector by f. Equation (J-9) takes the compact form

\mathbf {y} =\mathbf {Lf} .

(J-10)

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

(\mathbf {L^{T}L} )\mathbf {f} =\mathbf {L^{T}d} .

(J-11)

Consider the special case of a three-point wavelet (w₀, w₁, w₂). Set up the L matrix of equation (J-9) for this special case

\mathbf {L} ={\begin{pmatrix}w_{0}&0&0\\w_{1}&w_{0}&0\\w_{2}&w_{1}&w_{0}\\0&w_{2}&w_{1}\\0&0&w_{2}\end{pmatrix}},

(J-12a)

with its transpose L^T given by

\mathbf {L^{T}} ={\begin{pmatrix}w_{0}&w_{1}&w_{2}&0&0\\0&w_{0}&w_{1}&w_{2}&0\\0&0&w_{0}&w_{1}&w_{2}\\\end{pmatrix}}.

(J-12b)

Now multiply the two matrices to get

\mathbf {L^{T}L} ={\begin{pmatrix}w_{0}^{2}+w_{1}^{2}+w_{2}^{2}&w_{1}w_{0}+w_{2}w_{1}&w_{2}w_{0}\\w_{0}w_{1}+w_{1}w_{2}&w_{0}^{2}+w_{1}^{2}+w_{2}^{2}&w_{1}w_{0}+w_{2}w_{1}\\w_{0}w_{2}&w_{0}w_{1}+w_{1}w_{2}&w_{0}^{2}+w_{1}^{2}+w_{2}^{2}\end{pmatrix}}.

(J-12c)

Compute the first three lags of the autocorrelation (r₀, r₁, r₂) of the wavelet (w₀, w₁, w₂):

{\begin{array}{l}r_{0}=w_{0}^{2}+w_{1}^{2}+w_{2}^{2}\\r_{1}=w_{0}w_{1}+w_{1}w_{2}\\r_{2}=w_{0}w_{2},\end{array}}

(J-13)

and note that the elements of the covariance matrix L^TL given by equation (J-12c) are the first three autocorrelation lags of the wavelet (w₀, w₁, w₂) given by equations (J-13). Thus, for the general case, we note that

\mathbf {L^{T}L} ={\begin{pmatrix}r_{0}&r_{1}&r_{2}&\ldots &r_{n-1}\\r_{1}&r_{0}&r_{1}&\ldots &r_{n-2}\\r_{2}&r_{1}&r_{0}&\ldots &r_{n-3}\\\vdots &\vdots &\vdots &\ddots &\vdots \\r_{n-1}&r_{n-2}&r_{n-3}&\ldots &r_{0}\end{pmatrix}},

(J-14)

where (r₀, r₁, r₂, …, r_n−1) are the first n autocorrelation lags of the input wavelet series (w₀, w₁, w₂, w_m−1).

For the special case of a desired output vector d

\mathbf {d} ={\begin{pmatrix}d_{0}\\d_{1}\\d_{2}\\d_{3}\\d_{4}\\\end{pmatrix}},

(J-15a)

obtain the matrix product L^Td

\mathbf {L^{T}d} ={\begin{pmatrix}w_{0}d_{0}+w_{1}d_{1}+w_{2}d_{2}\\w_{0}d_{1}+w_{1}d_{2}+w_{2}d_{3}\\w_{0}d_{2}+w_{1}d_{3}+w_{2}d_{4}\\\end{pmatrix}}.

(J-15b)

Now compute the first three lags of the crosscorrelation (g₀, g₁, g₂) of the desired output (d₀, d₁, d₂, d₃, d₄) with the input wavelet (w₀, w₁, w₂)

g_{0}=w_{0}d_{0}+w_{1}d_{1}+w_{2}d_{2},

(J-16a)

g_{1}=w_{0}d_{1}+w_{1}d_{2}+w_{2}d_{3},

(J-16b)

g_{2}=w_{0}d_{2}+w_{1}d_{3}+w_{2}d_{4},

(J-16c)

and note that the elements of the matrix L^Td given by equation (J-15b) are the first three lags of the crosscorrelation of the desired output (d₀, d₁, d₂, d₃, d₄) with the input wavelet (w₀, w₁, w₂) given by equations (J-16a, J-16b, J-16c). Thus, for the general case, we note that

\mathbf {L^{T}d} ={\begin{pmatrix}g_{0}\\g_{1}\\g_{2}\\\vdots \\g_{n-1}\end{pmatrix}},

(J-17)

where (g₀, g₁, g₂, …, g_n−1) are the first n lags of the crosscorrelation of the desired output d with the input wavelet w.

By substituting equations (J-14) and (J-17) into equation (J-11), we get

{\begin{pmatrix}r_{0}&r_{1}&r_{2}&\ldots &r_{n-1}\\r_{1}&r_{0}&r_{1}&\ldots &r_{n-2}\\r_{2}&r_{1}&r_{0}&\ldots &r_{n-3}\\\vdots &\vdots &\vdots &\ddots &\vdots \\r_{n-1}&r_{n-2}&r_{n-3}&\ldots &r_{0}\end{pmatrix}}{\begin{pmatrix}f_{0}\\f_{1}\\f_{2}\\\vdots \\f_{n-1}\end{pmatrix}}={\begin{pmatrix}g_{0}\\g_{1}\\g_{2}\\\vdots \\g_{n-1}\end{pmatrix}}.

(J-18)

This is the discrete form of the classic Wiener-Hopf integral equation to estimate the least-squares shaping filter f that converts an input wavelet w into a desired form d.

In general, the autocorrelation matrix L^TL in equation (J-11) is diagonally dominant since r₀ is the largest value of the autocorrelation lags (r₀, r₁, r₂, …, r_n−1). Nevertheless, practical implementations of equation (J-11) often require adding a small fraction ε of the zero-lag of the autocorrelation to the diagonal elements of the matrix L^TL:

(\mathbf {L^{T}L} +\epsilon \mathbf {I} )\mathbf {f} =\mathbf {L^{T}d} ,

(J-19)

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

Note that the autocorrelation matrix L^TL given by equation (J-14) is of Toeplitz form. A typical length for the deconvolution filter f in equation (J-18) is between 40 and 80 samples. This makes the size of the autocorrelation matrix L^TL to be between 40 × 40 and 80 × 80. If the autocorrelogram is computed from an input seismogram represented by a single trace, then the length of the input data vector is typically 1000 samples. The ratio of the length of the input seismogram used in computing the autocorrelation lags to the filter length should be no less than 8.

Equation (J-11) may be solved for the filter vector f using Levinson recursion (Section B.6). It may also be solved by the conjugate gradient method ^[26] ^[27] ^[28]. The conjugate gradient method by Hestenes ^[29] described and exemplified by Wang and Treitel ^[26] is outlined below.

Refer to equation (J-11) and define

\mathbf {A} =\mathbf {L^{T}L}

(J-20a)

and

\mathbf {B} =\mathbf {L^{T}d} ,

(J-20b)

so that

\mathbf {Af} =\mathbf {B} .

(J-21)

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

\mathbf {R} _{i}=\mathbf {B} -\mathbf {AX} _{i}

(J-22)

becomes a null vector itself. Here, X_i is the parameter vector estimate after the ith iteration. If the autocorrelation matrix A has dimensions n × n, the conjugate gradient method yields the solution f after m < n iterations.

The recursive scheme requires the following initial values

c_{-1}=1,

(J-23a)

\mathbf {P} _{-1}=0,

(J-23b)

and

\mathbf {R} _{0}=\mathbf {B} -\mathbf {A\ X} _{0}.

(J-23c)

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

c_{i}=\mathbf {R} _{i}^{T}\mathbf {R} _{i},\quad i=0,1,2,\cdots ,m.

(J-24a)

To begin the recursion, set i = 0 and compute the residual error c₀ using equation (J-24a). Then, insert the values for c₋₁ from equation (J-23a) and c₀ from equation (J-24a) into ^[26]

b_{i-1}=c_{i}/c_{i-1}

(J-24b)

to get the value for b₋₁. Next insert the values for P₋₁ from equation (J-23b), R₀ from equation (J-23c), and b₋₁ from equation (J-24b) into

\mathbf {P} _{i}=\mathbf {R} _{i}+b_{i-1}\mathbf {P} _{i-1}

(J-24c)

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

\mathbf {Q} _{i}=\mathbf {AP} _{i},\quad m\leq n,

(J-24d)

d_{i}=\mathbf {P} _{i}^{T}\mathbf {Q} _{i},

(J-24e)

and

a_{i}=c_{i}/d_{i}

(J-24f)

to get an estimate X₁ of the filter vector f given by

\mathbf {X} _{i+1}=\mathbf {X} _{i}+a_{i}\mathbf {P} _{i}

(J-24g)

and the associated new residual vector R₁ given by

\mathbf {R} _{i+1}=\mathbf {R} _{i}-a_{i}\mathbf {Q} _{i}.

(J-24h)

Return to the beginning of the recursion and repeat the computations given by equations (J-24a) through (J-24h). Stop the iteration after a specified value of m < n.

I have tested this algorithm and compared the results with those obtained from Levinson recursion. For a filter length n, conjugate gradient gives the exact solution derived from Levinson recursion after m = n iterations. Nevertheless, an approximate solution that is acceptable without much compromise can be achieved by m < n. This means that for a large filter length n, conjugate gradient can be more efficient than Levinson recursion.

J.3 Applications of the GLI technique

In this section, we shall review selected applications of generalized linear inversion (GLI) for data modeling as part of a seismic data processing sequence. Specifically, we shall compile a summary of the theoretical discussions in previous appendixes on deconvolution (Sections B.5 and B.8), residual statics corrections (Section C.4), refraction statics corrections (Section C.8), and the discrete Radon transform (Section F.3) within the context of the GLI solution given by equations (J-6a,J-6b) or (J-7a, J-7b).

For each application, we shall make reference to the model equation, input data vector d, the coefficient matrix L, the parameter vector p and their specific sizes, and the technique to solve for the parameter vector.

The model equation for surface-consistent deconvolution expressed in the frequency domain is given by (Section B.8):

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

(J-25)

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

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

{\begin{pmatrix}\ \\\ \\\ \\\vdots &\\{\tilde {X}}_{ij}^{\prime }&\\\vdots \\\ \\\ \\\ \\\end{pmatrix}}={\begin{pmatrix}\ \\\ \\\ \\\cdots 1\cdots \ \ 1\cdots \ \ 1\cdots \ \ 1\cdots \\\ \\\ \\\ \\\end{pmatrix}}{\begin{pmatrix}\vdots \\{\tilde {S}}_{j}\\\vdots \\\vdots \\{\tilde {H}}_{l}\\\vdots \\\vdots \\{\tilde {E}}_{k}\\\vdots \\\vdots \\{\tilde {G}}_{i}\\\vdots \end{pmatrix}}.

(J-26)

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

We want to estimate for each frequency ω the model parameters p such that the difference between the actual spectral component $\mathbf {\tilde {X}}$ and the modeled spectral component $\mathbf {{\tilde {X}}^{\prime }}$ is minimum in the least-squares sense. Follow the steps from equation (J-2) to (J-7a) for the complex case and derive the GLI solution for surface-consistent deconvolution, $\mathbf {p} =(\mathbf {L^{T\ast }L} )^{\mathbf {-1} }\mathbf {L^{T\ast }{\tilde {X}}} .$ The parameter vector p that contains the spectral components ${\tilde {S}}_{j},{\tilde {G}}_{i},{\tilde {H}}_{l}\ {\text{and}}\ {\tilde {E}}_{k},$ which are associated with the source and receiver locations, offset dependency, and earth’s impulse response, is solved for each frequency component ω. Results from all frequency components are then combined to obtain the terms in equation (J-25). The surface-consistent spiking deconvolution operator to be applied to each trace in the data set is then the minimum-phase inverse of s_j(t) * g_i(t) * h_l(t).

The size of the matrix L^TL is (n_s + n_h + n_e + n_r) × (n_s + n_h + n_e + n_r). Consider 500 frequency components for which equation (J-28) is to be solved for. For the example specified above, you would have to solve a matrix equation $\mathbf {p} =(\mathbf {L^{T\ast }L} )^{\mathbf {-1} }\mathbf {L^{T\ast }{\tilde {X}}}$ of the size 4060 × 4060 for each of the 500 frequency components. A practical scheme for solving this equation is based on the Gauss-Seidel iterative scheme described in Section B.8.

To estimate surface-consistent shot and receiver residual statics, we model traveltime deviations on moveout-corrected CMP gathers. The model equation for residual statics estimation is given by (Section C.4)

t_{ij}^{\prime }=s_{j}+r_{i}+G_{k}+M_{k}x_{ij}^{2},

(J-27)

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

Consider a data set with n_s shot locations and n_c channels, so that the maximum number of picks is n_s × n_c. For n_s × n_c picks of t_ij, and n_s shot locations, n_r receiver locations, n_G midpoint locations, we have the following set of equations:

{\begin{pmatrix}\vdots \\\vdots \\\vdots \\t_{ij}^{\prime }\\\vdots \\\vdots \\\vdots \\\end{pmatrix}}={\begin{pmatrix}\ \\\ \\\ \\\ \\\cdots 1\quad \cdots 1\quad \cdots 1\quad \cdots x_{ij}^{2}\cdots \ \\\ \\\ \\\ \\\ \\\ \\\end{pmatrix}}{\begin{pmatrix}\vdots \\s_{j}\\\vdots \\r_{i}\\\vdots \\G_{k}\\\vdots \\M_{k}\\\vdots \end{pmatrix}},

(J-28)

which in compact matrix notation is given by t′ = Lp as in equation (J-1). Here t′ is the column vector of (n_s × n_c)-length in equation (J-28), L is the sparse matrix in the same equation with dimensions (n_s × n_c) × (n_s + n_r + n_G + n_G), and p is the column vector of (n_s + n_r + n_G + n_G)-length on the right-hand side of the same equation. Except for the three elements in each row, the L matrix contains zeros. Consider an example of n_s = 1000, n_c = 240, n_r = 1000, and n_G = 2000. You would then have 240 000 equations to estimate 6000 parameters.

Follow the steps from equation (J-2) to (J-6) and derive the GLI solution for surface-consistent residual statics, p = (L^TL)⁻¹L^Tt. Here t denotes the column vector of (n_s × n_c)-length that represents the traveltime deviations picked from moveout-corrected CMP gathers. The size of the matrix L^TL is (n_s + n_r + n_G + n_G) × (n_s + n_r + n_G + n_G). For the example specified above, you would have to solve a matrix equation p = (L^TL)⁻¹L^Tt of the size 6000 × 6000. A practical scheme for solving this equation is based on the Gauss-Seidel method (Section C.4).

To estimate surface-consistent shot and receiver intercept time anomalies, and thus obtain shot and receiver refraction statics, we model refracted arrival times derived from picking first breaks (refraction statics corrections). The model equation for the refracted arrivals is given by (Section C.8)

t_{ij}^{\prime }=T_{j}+T_{i}+s_{b}x_{ij},

(J-29)

where T_j and T_i are the intercept time anomalies at shot and receiver locations, respectively, and s_b is the bedrock slowness. Here, we shall consider the special case of surface shots. Hence, for n shot/receiver stations the parameter vector is p : (T₁, T₂, …, T_n; s_b).

Consider a data set with n shot locations and n_c channels, so that the total number of traces is n × n_c. Assume that you will pick the first breaks from half the number of the traces. For n × n_c/2 picks of t_ij, and n + 1 parameters p : (T₁, T₂, …, T_n; s_b), we have the following set of equations:

{\begin{pmatrix}\vdots \\\vdots \\\vdots \\t_{ij}^{\prime }\\\vdots \\\vdots \\\vdots \\\end{pmatrix}}={\begin{pmatrix}\ \\\ \\\ \\\cdots \quad 1\quad \cdots \quad 1\quad \cdots \quad x_{ij}\\\ \\\ \\\ \\\end{pmatrix}}{\begin{pmatrix}\vdots \\T_{j}\\\vdots \\T_{i}\\\vdots \\s_{b}\end{pmatrix}},

(J-30)

which in matrix notation is t′ = Lp as in equation (J-1). Here t′ is the column vector of (n × n_c/2)-length in equation (J-30), L is the sparse matrix in the same equation with dimensions (n × n_c/2) × (n + 1), and p is the column vector of (n + 1)-length on the right-hand side of the same equation. Except for the three elements in each row, the L matrix contains zeros. Consider an example of n = 1000 and n_c = 240. You would then have 120 000 equations to estimate 1001 parameters.

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

In some applications, instead of directly modeling the observed data, such as refracted arrivals, it may be preferable to perturb the parameter vector p by a small amount Δp and compute the change in the model vector Δt′. Perturbation of the parameter vector associated with the earth is the basis for tomographic inversion. Since it is the change in the parameter vector and not the data vector that we estimate, tomographic application of the GLI technique may be appropriately considered an example of earth modeling rather than data modeling. In Section J.6, we review reflection traveltime tomography for updating earth model parameters. Here, we will refer to the model equation for refraction tomography (Section C.9)

\Delta t_{ij}^{\prime }=Z_{j}\Delta s_{wj}+Z_{i}\Delta s_{wi}+X_{ij}\Delta s_{b},

(J-31)

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

Examine the structure of equation (J-31) and note that, instead of modeling refraction traveltimes by way of equation (J-29), we model the change in traveltimes by way of equation (J-31) and thus estimate the near-surface parameters. Hence, for n shot/receiver stations, the parameter vector is Δp : (Δs_w1, Δs_w2, …, Δs_wn, Δs_b). We refer to the scheme based on equation (J-31) as the iterative GLI solution.

Consider a data set with n shot locations and n_c channels, so that the total number of traces is n × n_c. Again, assume that you will pick the first breaks from half the number of the traces. For n × n_c/2 picks of t_ij and n + 1 parameters Δp : (Δs_w1, Δs_w2, …, Δs_wn, Δs_b), we have the following set of equations:

{\begin{pmatrix}\vdots \\\vdots \\\vdots \\\Delta t_{ij}^{\prime }\\\vdots \\\vdots \\\vdots \\\end{pmatrix}}={\begin{pmatrix}\ \\\ \\\ \\\cdots \quad Z_{j}\quad \cdots \quad Z_{i}\quad \cdots \quad X_{ij}\\\ \\\ \\\ \\\end{pmatrix}}{\begin{pmatrix}\vdots \\\Delta s_{wj}\\\vdots \\\Delta s_{wi}\\\vdots \\\Delta s_{b}\end{pmatrix}}.

(J-32)

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

{\boldsymbol {\Delta }}\mathbf {t^{\prime }} =\mathbf {L} {\boldsymbol {\Delta }}\mathbf {p} ,

(J-33)

where Δt′ is the column vector of (n + n_c/2)-length on the left-hand side of equation (J-32), L is the sparse matrix in the same equation with dimensions (n + n_c/2) × (n + 1), and p is the column vector of (n + 1)-length on the right-hand side of equation (J-32). Except for the three elements in each row, the L matrix contains zeros. For the example of n = 1000 and n_c = 240, you would have 120 000 equations to estimate 1001 parameters.

The GLI solution to equation (J-34) satisfies the requirement that the energy of the error vector e = Δt − Δt′ is minimum and is given by

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

(J-34)

where Δt denotes the column vector of n + n_c/2-length that represents the difference between the observed (picked) refracted arrival times and the initial estimate of the modeled times in equation (J-32). The size of the matrix L^TL is (n + 1) × (n + 1). For the example specified above, you would have to solve a matrix equation p = (L^TL)⁻¹L^Tt of the size 1001 × 1001. A practical scheme for solving this equation is based on the Gauss-Seidel method (Section C.4).

We now review data modeling by the discrete Radon transform. Let d(h, t) represent a CMP gather in the offset space (the plane of offset versus two-way traveltime) and u(v, τ) be the transformed data in the velocity space (the plane of stacking velocity versus two-way zero-offset time). The model equation for the hyperbolic Radon transform in the Fourier transform domain is given by (Section F.3)

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

(J-35)

where ω′ is the Fourier dual of t′, with t′ = t², t being the two-way traveltime. Also, v is the hyperbolic moveout velocity and h is the (half) offset.

Equation (J-33) can be written in the matrix form d′ = Lu as in equation (J-1) for each ω′ component of d′(h, ω′) and u(v, ω′), where L now is a complex matrix of dimensions m × p:

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

(J-36)

and d′ and u are complex vectors of lengths m and p, respectively, where m and p are the number of offsets and the number of velocities used in computing the discrete Radon transform. Note that the elements of the L matrix only depend on the geometry of the input CMP gather and the range of velocities used in constructing the velocity-stack gather u(v, τ).

Follow the steps described by equations (J-2) through (J-7a) to obtain the GLI solution for each Fourier component of u(v, ω′) given by u = (L^T*L)⁻¹L^T*d, where the asterisk denotes complex conjugate. The size of the matrix L^TL is m × p. For m = 120, p = 120, and 500 Fourier components, you would have to solve a matrix equation p = (L^TL)⁻¹L^Tt of the size 120 × 120 for each of the 500 Fourier components. The elements of the matrix L given by equation (J-36) may have similar magnitudes in two different rows because of the small separation of offsets associated with the input CMP gather. As such, the matrix can be near singular and requires a solution based on singular-value decomposition (Section F.3).

Finally, slant-stack transform is a special case of the generalized discrete Radon transform based on linear moveout equation. For details of the generalized linear inversion formulation of the slant-stack transform, see Section F.3.

J.4 Dix conversion

We shall develop the traveltime equation for a CMP recording geometry and an earth model that comprises horizontal isovelocity layers (Figure J-1). Then, we shall derive the Dix equation for interval velocities computed from rms velocities associated with the horizontal layers. From the geometry of Figure J-1, note that the half-offset h is given by

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

or

h=\sum _{k=1}^{n}\tan \ \theta _{k}\Delta z_{k}.

(J-37a)

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

p={\frac {\sin \theta _{k}}{v_{k}}},

(J-37b)

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

\tan \ \theta _{k}={\frac {pv_{k}}{\sqrt {1-p^{2}v_{k}^{2}}}}.

(J-37c)

Also note that

\Delta z_{k}={\frac {v_{k}\Delta \tau _{k}}{2}},

(J-37d)

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

x=\sum _{k=1}^{n}{\frac {pv_{k}^{2}\Delta \tau _{k}}{\sqrt {1-p^{2}v_{k}^{2}}}}.

(J-38)

Figure J-1 Geometry of a horizontally layered earth model with constant velocities to develop the theory for the normal moveout and Dix equations (Section J.4).

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

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

or

t_{n}(x)=2\sum _{k=1}^{n}{\frac {\Delta z_{k}}{v_{k}\cos \theta _{k}}}.

(J-39a)

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

\cos \ \theta _{k}={\sqrt {1-p^{2}v_{k}^{2}}}.

(J-39b)

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

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

(J-40)

Equations (J-38) and (J-40) are the raypath equations for an earth model with horizontal isovelocity layers. Based on source-receiver reciprocity, the traveltime $t_{n}^{2}(x)$ can be represented by a curve that is symmetric with respect to offset x (Taner and Koehler, 1969)

t_{n}^{2}(x)=C_{0}+C_{1}x^{2}+C_{2}x^{4}+\ldots ,

(J-41)

which also is referred to in Section C.1. We need to expand the square-root term in equations (J-38) and (J-40) into the Taylor series to derive the expressions for t²(x) and the even powers of x displayed in equation (J-41) ^[15]. The series expansions are

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

(J-42a)

and

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

(J-42b)

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

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

(J-43a)

and

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

(J-43b)

In equation (J-43a), the terms (a₁, a₂, a₃, …) are the fractional coefficients (1, 1/2, 3/8, …) displayed in equation (J-42a). Similarly, in equation (J-43b), the terms (b₀, b₁, b₂, …) are the fractional coefficients (1, 1/2, 3/8, …) displayed in equation (J-42b). Interchange the order of the summations in equations (J-43a) and (J-43b) to obtain

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

(J-44a)

and

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

(J-44b)

Define two terms, A_j and B_j

A_{j}=\sum _{k=1}^{n}a_{j}v_{k}^{2j}\Delta \tau _{k},

(J-45a)

B_{j}=\sum _{k=1}^{n}b_{j}v_{k}^{2j}\Delta \tau _{k},

(J-45b)

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

x=\sum _{j=1}^{\infty }p^{2j-1}A_{j}

(J-46a)

and

t_{n}(x)=\sum _{j=0}^{\infty }p^{2j}B_{j}.

(J-46b)

The even powers of x that we need to substitute into equation (J-41) can also be expressed as a double summation of the form given by equation (J-46a)

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

(J-47a)

where m = 1, 2, 3, …, and A_{j, 2m} is a recursive combination of A_j of equation (J-45a) ^[15].

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

t_{n}^{2}(x)=\sum _{j=0}^{\infty }p^{2j}B_{j,2},

(J-47b)

where B_{j, 2} is a recursive combination of B_j of equation (J-45b) ^[15].

By using equations (J-47a) and (J-47b) in equation (J-41), and setting the terms on both sides of the equation with like powers of p equal to another, the coefficients, C₀, C₁, C₂, …, can be determined in terms of the model parameters — interval velocities v_k and layer thicknesses defined in terms of two-way vertical times Δτ_k. While details can be found in Hubral and Krey ^[15], here, we shall consider a truncated form of equation (J-41) to derive the Dix equation.

Specifically, we shall consider the following form of equation (J-41)

t_{n}^{2}(x)\approx C_{0}+C_{1}x^{2}.

(J-48)

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

(x)\approx p\sum _{k=1}^{n}v_{k}^{2}\Delta \tau _{k}

(J-49a)

and

t_{n}(x)\approx \sum _{k=1}^{n}\Delta \tau _{k}+{\frac {p^{2}}{2}}\sum _{k=1}^{n}v_{k}^{2}\Delta \tau _{k}.

(J-49b)

Refer to equation (J-37b) and note that by retaining only the terms up to the second power of p, we are making a small-spread approximation as in equation (J-48).

Define

\tau _{n}=\sum _{k=1}^{n}\Delta \tau _{k}

(J-50a)

and

V_{n}^{2}={\frac {1}{\tau _{n}}}\sum _{k=1}^{n}v_{k}^{2}\Delta \tau _{k},

(J-50b)

where τ_n is the total two-way zero-offset elapsed time along the vertical raypath from the surface down to the nth interface, and V_n is the rms velocity at the nth interface measured along the vertical raypath (Figure J-1).

Use the definitions given by equations (J-50a) and (J-50b) to rewrite equations (J-49a) and (J-49b) to yield

x\approx pV_{n}^{2}\tau _{n}

(J-51a)

and

t_{n}(x)\approx \tau _{n}+{\frac {p^{2}V_{n}^{2}\tau _{n}}{2}}.

(J-51b)

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

x^{2}\approx p^{2}V_{n}^{4}\tau _{n}^{2}

(J-52a)

and

t_{n}^{2}(x)\approx \tau _{n}^{2}+p^{2}V_{n}^{2}\tau _{n}^{2}.

(J-52b)

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

\tau _{n}^{2}+p^{2}V_{n}^{2}\tau _{n}^{2}=C_{0}+C_{1}p^{2}V_{n}^{4}\tau _{n}^{2},

(J-53)

and set the terms with like powers of p on both sides of the equation equal to one another to determine the coefficients C₀ and C₁

C_{0}=\tau _{n}^{2}

(J-54a)

and

C_{1}={\frac {1}{V_{n}^{2}}}.

(J-54b)

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

t_{n}^{2}(x)\approx \tau _{n}^{2}+{\frac {x^{2}}{V_{n}^{2}}}.

(J-55)

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

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

(J-56a)

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

V_{n-1}^{2}={\frac {1}{\tau _{n-1}}}\sum _{k=1}^{n-1}v_{k}^{2}\Delta \tau _{k}.

(J-56b)

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

V_{n}^{2}={\frac {V_{n-1}^{2}\tau _{n-1}+v_{n}^{2}\Delta \tau _{n}}{\tau _{n}}}.

(J-56c)

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

v_{n}={\sqrt {\frac {V_{n}^{2}\tau _{n}-V_{n-1}^{2}\tau _{n-1}}{\tau _{n}-\tau _{n-1}}}},

(J-57)

which is equation (1) of the main text. Equation (J-57) is the basis for Dix conversion of rms velocities to interval velocities. Given the rms velocities measured from the surface down to the (n − 1) and nth layer boundaries along the vertical raypaths and the associated two-way zero-offset times t_n−1 and t_n, the layer velocity within the interval between the two layer boundaries can be computed using equation (J-57). Keep in mind the underlying assumption that the earth model comprises horizontal isovelocity layers. Although reflector dip can be incorporated into the Dix equation ^[15], in practice, it is most desirable to work with rms velocities estimated from gathers generated by prestack time migration. Additionally, keep in mind that the rms velocities used in equation (J-57) are based on straight-ray assumption; thus, ray bending at layer boundaries is not accounted for in Dix conversion. Finally, since equation (J-48) is an approximation to equation (J-41), the offset range used in estimating the rms velocities V_n and V_n−1 corresponds to a small spread.

J.5 Map processing

A map is defined as a 2-D surface g(x, y). Depending on the quantity being mapped, g(x, y) may have many types of units; for example, gravitational attraction (mGal), magnetic intensity (gamma), elevation, or even times picked along marker horizons from seismic data. Here, we set the convention that the positive x-axis points eastward and the positive y-axis points northward. A discrete map function is represented by a grid of mesh points over the x − y plane. These mesh points are spaced commonly at equal intervals in the x and y directions. For many types of mapping, g(x, y) is a smooth function and such a map can be analyzed in the Fourier transform domain. However, there are situations (as for isochron and structure maps) in which the map function has discontinuities that represent faulting.

Maps usually are created from irregularly spaced observation values. Thus, the map function at a particular grid point must be computed by some fitting procedure. A local plane surface fitting procedure is described later in this section.

Figure J-2 shows the contour map of the test surface used in this study. Void grid points have been filled with the arithmetic mean of the map function. Before map creation, some correction may be applied to observed data, followed by various types of editing. The main topic of this section is the next processing stage, which involves various 2-D processing techniques, each of which is aimed at achieving a particular interpretive goal. To clarify the objectives in digital processing, consider the properties of the quantity represented by a map function. In particular, consider the gravity problem. Ultra-long wavelength anomalies generally are associated with variations in crustal thickness. Moderately long wavelengths usually are caused by variations in the basement topography. Medium to short wavelength anomalies are related closely to local tectonic disturbances such as folding. Finally, very short wavelength anomalies have a variety of sources; some are because of small, near-surface features, while others are just noise of various types. Note that gradation in wavelength variations is related somewhat to the depth of the source that causes the anomaly. The longer the wavelength, the deeper the anomaly; the shorter the wavelength, the shallower the anomaly.

Similar physical interpretations are possible with other map functions, such as the seismic (time) horizon map in Figure J-2. Some of the very short wavelengths here result from near-surface effects that are manifested as residual statics. Also note in Figure J-2 that some moderately long wavelengths correspond to structural undulations that exist in the area. The most important observation made from this map is that most of the features with different wavelengths are not spatially isolated, but are superimposed. This characteristic is common for all types of map functions, whether gravity, magnetic, elevation, or time. To separate the effects of different features from each other, we must analyze them in terms of wavelength. This wavelength analysis also is a way to discriminate the depth of the source for some types of data such as potential-field measurements.

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

The basic motivation behind digital processing of maps is to separate anomalies. Simple 2-D smoothing and wavelength filtering are techniques for separating anomalies. Vertical derivatives and analytic continuation also are useful for enhancing certain anomalies so that they appear more pronounced on the map. Some techniques for separating anomalies by their wavelengths are described next.

Before applying digital processing techniques, it is best to study the map in terms of wavelength composition. The 2-D amplitude spectrum of a map is an excellent tool for recognizing not only the wavelength content, but also the orientation of various components. The most useful display is the color contour plot of the amplitude spectrum (Figure J-3) from which various bands of wavelengths are distinguished clearly. Pink represents long, beige represents moderate, and yellow represents short wavelength anomalies. For a 2-D real function, such as a map, the amplitude spectrum is antisymmetric. Thus, only a pair of quadrants (the first and second) of the amplitude spectrum needs to be displayed.

A simple 2-D smoothing operation is the easiest way to obtain a map that represents the regional anomaly. Figure J-4 shows the map in Figure J-2 after 2-D smoothing. Most of the very short wavelength features present in the center of the original map have been eliminated. Depending on the interpreter’s goal, this output might be a satisfactory regional map. (However, a smoother regional map often is desired.)

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

Smoothing basically is done by computing the average value of the grid points that fall onto a ring of some desired radius. The center of the ring coincides with the output point. There is no reason why more than one ring should not be considered. For n concentric rings with m_i points over the ith ring, the average value g_i of the quantity g_ij that is being mapped is

{\bar {g}}_{i}={\frac {1}{m_{i}}}\sum _{j=1}^{m_{i}}g_{ij}.

(J-58a)

Thus, the cumulative average g over n rings is

{\bar {g}}={\frac {1}{n}}\sum _{i=1}^{n}{\bar {g}}_{i}.

(J-58b)

Weighting factors, which depend on the distance from the center of the rings, often are used in smoothing algorithms. In these cases, equation (J-58b) becomes

{\bar {g}}=\sum _{i=1}^{n}w_{i}{\bar {g}}_{i},

(J-58c)

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

g=g_{0}-{\bar {g}},

(J-59)

where g₀ is the grid value at the center of the rings. Figure J-5 shows the regional anomaly obtained by using 15 rings. In general, the more rings, the more smoothing of the data. Note that the output of the ring method is relative. Depending on the nature of the objectives of map processing, Figure J-5 may or may not be a good estimate of the regional anomaly. In most cases, determining what is regional and what is residual is a matter of interpretive judgment.

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

The hypothetical map in Figure J-6a has ten different anomalies of various shapes and orientations. Figure J-6b is the 2-D amplitude spectrum of this map with (k_x, k_y) as the Fourier duals of (x, y). Anomaly 3, which has infinite wavelength in the north-south direction (y-axis) and a wavelength of 2AA′ in the east-west direction (x-axis), lies on the k_x-axis in the transform domain. Anomaly 7 also lies on the k_x-axis, but is farther from the origin since it has a shorter wavelength component (2BB′) in the x direction. Anomaly 8 tilted counterclockwise by an angle θ degrees from the y direction has an apparent wavelength component 2CC′ in the y direction and 2DD′ in the x direction.

Note the following relations:

\tan \theta ={\frac {\lambda _{x}}{\lambda _{y}}},

(J-60a)

where λ_x = DD′ and λ_y = CC′ are the wavelength components. By definition

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

and

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

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

\tan \theta ={\frac {k_{y}}{k_{x}}}.

(J-60b)

Figure J-5 A regional anomaly map based on equation (J-58c). The input is the contour map in Figure J-2.

Equations (J-60a) and (J-60b) imply that map trend MM′ is orthogonal to transform trend TT′.

Unlike the simple, linear anomalies 3, 7, and 8, rounded anomalies 1, 2, and 10 map onto the entire transform plane. Moreover, because of their equal size, anomalies 1 and 2 cannot be distinguished on the plane of 2-D amplitude spectrum, even though they are isolated in the space domain. However, the spatial extent of these two anomalies is much wider than that of anomaly 10. Thus, they are confined to lower wavenumbers. Finally, elongated features like anomalies 4, 5, 6, and 9 map onto the transform plane at certain dominant directions that are orthogonal to their respective spatial trends. The dots on the k_x − k_y plane in Figure J-6b depict the maximum energy associated with the anomalies.

Gridding involves fitting a locally plane surface to a set of control points around each grid output point. Consider a plane-surface fit in the least-squares sense:

{\tilde {g}}(x,y)=a_{0}+a_{1}x+a_{2}y.

(J-61)

The least-squares error is

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

(J-62)

where g is the observed value at the grid point (x, y), and M is the number of observations at and around that grid point. For local plane fitting, M usually is set to 9 points. We want to find a set of (a₀, a₁, a₂) for which L is minimum;

{\frac {\partial L}{\partial a_{0}}}={\frac {\partial L}{\partial a_{1}}}={\frac {\partial L}{\partial a_{2}}}=0.

(J-63)

Figure J-6 Anomalies of various shapes in (a) the space and (b) the wavenumber domains.

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

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

(J-64)

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

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

When put into matrix form, we obtain

{\begin{pmatrix}M&\sum x&\sum y\\\sum x&\sum x^{2}&\sum xy\\\sum y&\sum xy&\sum y^{2}\\\end{pmatrix}}{\begin{pmatrix}a_{0}\\a_{1}\\a_{2}\end{pmatrix}}={\begin{pmatrix}\sum g\\\sum xg\\\sum yg\end{pmatrix}}.

(J-65)

Equation (J-65) is solved for the set of coefficients (a₀, a₁, a₂).

Unlike the hypothetical anomalies in Figure J-6a, real data consist of anomalies of various shapes and orientations that are superimposed on each other. However, in the transform domain, the real anomalies can be separated in terms of their wavelength contents and orientations; something that cannot be achieved in the space domain. Thus, the transform domain provides a way to apply various filtering operations to a map. A band of wavelengths can be passed-regardless of orientation, as shown by the radial filter transfer function in Figure J-7a. Anomalies that are oriented in the range (θ₁, θ₂) from the northerly direction can be passed regardless of their size, using a directional filter whose transfer function is given by Figure J-7b. Finally, a band of wavelengths with a directional orientation can be passed. This is achieved by a filter with a transfer function as shown in Figure J-7c. In practice, as in other filter design techniques, the transfer functions shown must be tapered in the neighborhood of cutoff wavelengths for optimal performance.

Once the transfer function is designed to suit the purpose, the filter itself can be applied to the map in the transform or the space domain. In the transform domain, the transfer function is multiplied with the 2-D Fourier transform of the map. Subsequent inverse transformation yields the filtered map. To apply the filtering in the space domain, first inverse Fourier transform the filter’s transfer function to obtain its 2-D impulse response. Two-dimensional convolution of this impulse response with the map yields the filtered map. Note that filters whose transfer functions are depicted in Figure J-7 do not cause any phase shift; i.e., anomalies are not displaced in the space domain after filtering.

Before applying the 2-D filters, the possible bands of the wavelengths present on the amplitude spectrum of the map to be filtered must be determined. Four regions were distinguished from the color plot (Figure J-3) of the amplitude spectrum of the test surface. We have the regional component down to 21-km wavelengths. A subregional part of the spectrum between the 21- and 12-km wavelengths is defined. A moderate range of wavelengths is between 12 and 8 km, and the residual components are beyond the 8-km cutoff wavelength. This residual region can be subdivided; anything less than 4 km corresponds to very short wavelength anomalies. These anomalies may be caused by various types of noise.

The bands that were just determined now are used as cutoff wavelengths of the radial filter transfer function. A given map can be scanned with a suite of low-pass filters and several filtered maps can be produced, each with a potentially unique interpretational value. A result from such a filtering operation is shown in Figure J-8. Note the resemblance between Figure J-8 and the regional map (Figure J-5) that was obtained from the ring method described in this section. Although not demonstrated here, note that as the bandwidth of the filter is increased, more and more of the short wavelength anomalies are included, making the output less and less regional in character.

A filter scan is not limited to low-pass filtering only. High-pass, band-pass, and band-reject filters can be applied to maps to achieve a particular interpretational goal. For example, we may want to obtain a residual map that is free of the very short wavelength anomalies that usually are considered noise. This is accomplished by properly selecting the band-pass filter. A high-pass filter is the proper choice to retain all of the short wavelength region of the spectrum.

Directional filters scan the map of interest at various angles to emphasize a particular trend that may exist in the data. Figure J-9 shows the result of one of the several directional filters applied to the test surface in Figure J-2. A low-pass radial filter was cascaded with each directional filter so that any regional trend could be delineated in the area. The results of directional filtering indicate the existence of a prominent regional trend approximately N30°W (Figure J-9). In some cases, a certain band of wavelengths may have one dominant trend that is different from that of another band of wavelengths. This situation may imply a change in the tectonic setting over the geologic history in the area.

J.6 Reflection traveltime tomography

We want to update the parameters of an earth model that has already been estimated by a suitable combination of inversion methods discussed in this chapter. Reflection traveltime tomography is an inversion method for estimating the earth model parameters from the reflection traveltimes associated with the observed seismic data. The reflection traveltime from a source at the surface to a reflection point at the subsurface and back to a receiver at the surface is represented by an integral of the traveltime segments along the raypath that depend on the earth model parameters themselves. This makes the direct inversion of the traveltimes to estimate the earth model parameters a nonlinear problem. Nevertheless, based on a formal linearization of the eikonal equation, which is a ray-theoretical approximation to the scalar wave equation (Section H.2), it can be shown that small changes in reflection traveltimes are linearly related to small changes in earth model parameters ^[30]. We are thus motivated to use reflection traveltime tomography not to estimate an initial earth model, but to update an already estimated model.

Figure J-7 Transfer functions for (a) band-pass, (b) directional, and (c) directional band-pass wavenumber filters.
Figure J-8 The contour map in Figure J-2 after low-pass filtering.
Figure J-9 The contour map in Figure J-2 after directional low-pass filtering.

To develop a theory for the model update, we shall set the following framework for our strategy:

Define the earth model parameters in terms of slowness functions and depths at the boundaries of the layers included in the initial model.
Assume that the initial model is made up of horizontal layers with laterally invariant model parameters.
Use the reflection times from CMP data associated with the layer boundaries included in the model and update the earth model such that the discrepancy between the modeled reflection times and the actual reflection times is minimum in the least-squares sense.
Estimate the changes in the model parameters, rather than the parameters themselves.
Perturb the initial model in two parts — slowness perturbation and depth perturbation.
Perturb the model parameters while preserving the offset values of the seismic data.

Consider the raypath geometry in a horizontally layered earth model shown in Figure J-10. We want to derive the reflection traveltime equation for the raypath from the reflection point R_i at the interface z_n in the subsurface to the receiver location G_j at the surface z = 0. This one-way raypath is associated with a CMP gather at midpoint location y and half-offset h. Our earth model consists of n layers above the reflection point R_i.

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

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

t_{AC}^{\prime }={\frac {AC}{v_{k}}}

(J-66a)

or

t_{AC}^{\prime }=(z_{k}-z_{k-1})s_{k}\sec \theta _{k},

(J-66b)

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

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

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

(J-67)

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

BC=AB\tan \theta _{k}

(J-68a)

or

BC=(z_{k}-z_{k-1})\tan \theta _{k}.

(J-68b)

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

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

(J-69)

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

s=s_{k}\sin \theta _{k},

(J-70)

which is the horizontal component of the slowness.

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

\left[t_{n}^{\prime }\right]_{modeled}=\left[t_{n}^{\prime }\right]_{initial}+\left[{\frac {\partial t_{n}^{\prime }}{\partial p}}\right]_{modeled}\Delta p.

(J-71)

The error in modeling the traveltimes is given by

e_{n}=\left[t_{n}\right]_{observed}-\left[t_{n}^{\prime }\right]_{modeled}.

(J-72)

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

e_{n}=\left[t_{n}\right]_{observed}-\left[t_{n}^{\prime }\right]_{initial}-\left[{\frac {\partial t_{n}^{\prime }}{\partial p}}\right]_{modeled}\Delta p.

(J-73)

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

e_{n}=\Delta t_{n}-{\frac {\partial t_{n}^{\prime }}{\partial p}}\Delta p.

(J-74a)

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

e_{n}=\Delta t_{n}-\Delta t_{n}^{\prime },

(J-74b)

where

\Delta t_{n}^{\prime }={\frac {\partial t_{n}^{\prime }}{\partial p}}\Delta p.

(J-75a)

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

\Delta t_{n}^{\prime }=\sum _{m=1}^{n}\Delta t_{n}^{\prime }(s_{m})+\sum _{m=1}^{n}\Delta t_{n}^{\prime }(z_{m}),

(J-75b)

where $\Delta t_{n}^{\prime }(s_{m})\ {\text{and}}\ \Delta t_{n}^{\prime }(z_{m})$ are the contributions of the slowness perturbation and depth perturbation, respectively.

To compute the traveltime perturbations caused by the slowness and depth perturbations, we shall follow the elegant formulation by Sherwood ^[18]. The slowness perturbation Δs_m in layer m, where 1 ≤ m ≤ n, yields a traveltime perturbation $\Delta t_{n}^{\prime }(s_{m})$ given by

\Delta t_{n}^{\prime }(s_{m})={\frac {\partial t_{n}^{\prime }}{\partial s_{m}}}\Delta s_{m}+\sum _{k=1}^{n}{\frac {\partial t_{n}^{\prime }}{\partial \theta _{k}}}\Delta \theta _{k}.

(J-76a)

Since the offset of the input data is preserved in the perturbation, note from Figure J-10 and equation (J-70) that, while perturbing the slowness, the propagation angle also is perturbed. That is why we include the second partial differential term on the right-hand side of equation (J-76a). In fact, slowness perturbation of a single layer m causes perturbation of propagation angles within all of the n layers along the traveltime path. And that is why we have the summation in the second term on the right-hand side of equation (J-76a).

Compute the partial derivatives in equation (J-76a) using equation (J-67) to yield

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

(J-76b)

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

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

(J-76c)

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

\Delta h_{n}(s_{m})={\frac {\partial h_{n}}{\partial s_{m}}}\Delta s_{m}+\sum _{k=1}^{n}{\frac {\partial h_{n}}{\partial \theta _{k}}}\Delta \theta _{k}.

(J-77a)

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

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

(J-77b)

Note that this summation is the same as the summation in the second term on the right-hand side of equation (J-76c). A direct consequence of the rule for preserving the offset value during perturbation leads to setting Δh_n(s_m) = 0 in equation (J-77b). This means that the second term on the right-hand side of equation (J-76c) vanishes, giving us the final expression for the slowness perturbation:

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

(J-78)

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

\Delta t_{n}^{\prime }(z_{m})={\frac {\partial t_{n}^{\prime }}{\partial z_{m}}}\Delta z_{m}+\sum _{k=1}^{n}{\frac {\partial t_{n}^{\prime }}{\partial \theta _{k}}}\Delta \theta _{k}.

(J-79a)

Since the offset of the input data is preserved in the perturbation, note from Figure J-1 that, while perturbing the depth, the propagation angle also is perturbed. Again, that is why we include the second partial differential term on the right-hand side of equation (J-79a). Since depth perturbation of a single layer m causes perturbation of propagation angles within all of the n layers along the traveltime path, we have the summation in the second term on the right-hand side of equation (J-79a).

Compute the partial derivatives in equation (J-79a) using equation (J-67) to yield

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

(J-79b)

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

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

(J-79c)

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

\Delta h_{n}(z_{m})={\frac {\partial h_{n}}{\partial z_{m}}}\Delta z_{m}+\sum _{k=1}^{n}{\frac {\partial h_{n}}{\partial \theta _{k}}}\Delta \theta _{k}.

(J-80a)

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

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

(J-80b)

The two parts of the first term of the right-hand side of this equation (tan θ_m) and (tan θ_m+1) are justified as for equation (J-79b).

A direct consequence of the rule for preserving the offset value during perturbation leads to setting Δh_n(z_m) = 0 in equation (J-80b). This then leads to

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

(J-81a)

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

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

(J-81b)

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

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

(J-81c)

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

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

(J-82)

We now combine, using equation (J-75b), the slowness perturbation given by equation (J-78) and the depth perturbation given by equation (J-82) and obtain

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

(J-83)

Examine the structure of equation (J-83) and note that, instead of modeling reflection traveltimes, we model the change in traveltimes, and thus estimate the change in the model parameters. The parameter vector is Δp : (Δs_m, …, Δz_m, …), where 1 ≤ m ≤ n.

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

{\begin{pmatrix}\vdots \\\vdots \\\vdots \\\Delta t_{n}^{\prime }\\\vdots \\\vdots \\\vdots \\\end{pmatrix}}={\begin{pmatrix}\\\ldots \quad Z_{m}\quad \ldots \quad S_{m}\quad \ldots \\\\\end{pmatrix}}{\begin{pmatrix}\vdots \\\Delta s_{m}\\\vdots \\\Delta z_{m}\\\vdots \\\end{pmatrix}},

(J-84)

where

Z_{m}=(z_{m}-z_{m-1})\ \sec \theta _{m}

(J-85a)

and

S_{m}=s_{m}\ \cos \theta _{m}-s_{m+1}\cos \theta _{m+1}.

(J-85b)

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

{\boldsymbol {\Delta }}\mathbf {t^{\prime }} =\mathbf {L} {\boldsymbol {\Delta }}\mathbf {p} ,

(J-86)

where Δt′ is the column vector on the left-hand side of equation (J-84), L is the sparse matrix with nonzero elements Z_m and S_m given by equations (J-85a,J-85b), and Δp is the column vector on the right-hand side of equation (J-84). Except for the two elements in each row, the L matrix contains zeros.

The generalized linear inversion (GLI) solution to equation (J-86) satisfies the requirement that the energy of the error vector

\mathbf {e} ={\boldsymbol {\Delta }}\mathbf {t} -{\boldsymbol {\Delta }}\mathbf {t^{\prime }}

(J-87)

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

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

(J-88)

where Δt denotes the column vector that represents the difference between the observed reflection traveltimes and the initial estimate of the modeled times, and T denotes matrix transposition.

We now outline the procedure for a tomographic model update based on the theory described in this section:

Perform prestack depth migration using the initial model and generate a set of image gathers.
Compute the residual moveout for all offsets along events on image gathers that correspond to the layer boundaries included in the model, and thus build the traveltime error vector Δt. For instance, you may have 10 layers, 1000 CMPs with a fold of 30. This means that the length of the traveltime error vector Δt is 300 000.
Define the initial model by a set of slowness and depth parameters, and construct the coefficient matrix L in equation (J-88) by computing the nonzero matrix elements Z_m and S_m given by equations (J-85a) and (J-85b). For instance, you may have 10 layers and each layer may be defined by 50 pairs of slowness and depth values in the lateral direction. This means you would have 1000 parameters in your model space. It also means that for the example given in step (b), you have 300 000 equations to solve for 1000 parameters. The solution for this overdetermined system is given by equation (J-88).
Estimate the change in parameters vector Δp, by way of the GLI solution given by equation (J-88).
Update the parameter vector p + Δp.
Iterate steps (a) through (e) as necessary to minimize the discrepancy between the modeled and actual traveltimes.

J.7 Threshold for velocity-depth ambiguity

Velocity-depth ambiguity states that an error in depth is indistinguishable from an error in velocity. Specifically, we shall investigate under what circumstances two earth models defined by layer velocities and reflector geometries in depth are indistinguishable. We shall adopt the derivation by Lines ^[10]. Albeit the underlying theory to be presented is somewhat simplistic compared to the comprehensive analysis by Bickel ^[9], it does provide a practical insight to important aspects of velocity-depth ambiguity.

Consider a single, horizontal reflector at depth z in a medium with constant velocity v. The reflection traveltime associated with a CMP recording geometry of a fixed offset 2h is given by the hyperbolic moveout equation

t(v,z)={\frac {2{\sqrt {z^{2}+h^{2}}}}{v}}.

(J-89)

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

\Delta t={\frac {\partial t}{\partial v}}\Delta v+{\frac {\partial t}{\partial z}}\Delta z.

(J-90)

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

{\frac {\partial t}{\partial v}}=-2{\sqrt {z^{2}+h^{2}}}{\frac {1}{v^{2}}},

(J-91a)

and simplify by way of equation (J-89)

{\frac {\partial t}{\partial v}}=-{\frac {t}{v}}.

(J-91b)

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

{\frac {\partial t}{\partial z}}={\frac {2z}{v{\sqrt {z^{2}+h^{2}}}}},

(J-92a)

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

{\frac {\partial t}{\partial z}}={\frac {4z}{v^{2}t}}.

(J-92b)

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

{\frac {\Delta t}{t}}=-{\frac {\Delta v}{v}}+{\frac {4z\Delta z}{v^{2}t^{2}}}.

(J-93)

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

{\frac {\Delta t}{t}}=-{\frac {\Delta v}{v}}+{\frac {\Delta z}{z}}.

(J-94)

Equation (J-94) states a rule of fundamental importance: for a zero-offset case, when the perturbation in velocity Δv/v is the same as the perturbation in depth Δz/z, no change occurs in traveltime; thus, velocity-depth ambiguity is infinite.

Now consider two earth models — an unperturbed model defined by the variables (z, v) and a perturbed model defined by the variables (z + Δz, v + Δv), such that the two models are indistinguishable at zero offset. This means that the zero-offset times t₀ associated with the two models are identical. Rewrite equation (J-89) in terms of the constant zero-offset time t₀ = 2z/v

t(v)=t_{0}{\sqrt {1+{\frac {4h^{2}}{t_{0}^{2}v^{2}}}.}}

(J-95)

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

\Delta t={\frac {\partial t}{\partial v}}\Delta v.

(J-96)

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

\Delta t=-{\frac {4h^{2}\Delta v}{v^{3}t_{0}}}{\frac {1}{\sqrt {1+{\dfrac {4h^{2}}{t_{0}^{2}v^{2}}}}}}.

(J-96)

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

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

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

{\frac {\Delta v}{v}}={\frac {z^{2}}{h^{2}}}{\frac {\Delta t}{t_{0}}},

(J-97)

where z = vt₀/2. Equation (J-97) states another rule for velocity-depth ambiguity: for a reflector at depth z, velocity ambiguity defined by Δv/v can be made smaller by increasing the offset h and decreasing the error in the time pick defined by Δt/t₀ at offset h. Note that for the zero-offset case, h = 0, the velocity ambiguity is infinite. This conclusion also was reached earlier by way of equation (J-94).^[31] ^[32]

References

↑ Dix, 1955, Dix, C. H., 1955, Seismic velocities from surface measurements: Geophysics, 20, 68–86.
↑ Landa et al., 1988, Landa, E., Kosloff, D., Keydar, S., Koren, Z., and Reshef, M., 1988, A method for determination of velocity and depth from seismic data: Geophys. Prosp., 36, 223–243.
↑ Faye and Jeannaut, 1986, Faye, J-P. and Jeannaut, J-P., 1986, Prestack migration velocities from focusing depth analysis: 56th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 438–440.
↑ Reshef and Kessler, 1989, Reshef, M. and Kessler, D., 1989, Practical implementation of three-dimensional poststack depth migration: Geophysics, 54, 309–318.
↑ Reshef, 1997, Reshef, M., 1997, The use of 3-D prestack depth imaging to estimate layer velocities and reflector positions: Geophysics, 62, 206–210.
↑ Reshef, 2001, Reshef, M., 2001, Some aspects of interval velocity analysis using 3-D depth migrated gathers: Geophysics (Scheduled for publication in Jan.-Feb. 2001 issue).
↑ Al-Yahya, 1989, Al-Yahya, K., 1989, Velocity analysis by iterative profile migration: Geophysics, 54, 718–729.
↑ ^8.0 ^8.1 Deregowski, 1990, Deregowski, 1990, Common-offset migrations and velocity analysis: First Break, 8, 225–234.
↑ ^9.0 ^9.1 Bickel, 1990, Bickel, S. H., 1990, Velocity-depth ambiguity of reflection traveltimes: Geophysics, 55, 266–276.
↑ ^10.0 ^10.1 Lines, 1993, Lines, L., 1993, Ambiguity in analysis of velocity and depth: Geophysics, 58. 596–597.
↑ Lynn and Claerbout, 1982, Lynn, W. S. and Claerbout, J. F., 1982, Velocity estimation in laterally varying media: Geophysics, 47, 884–897.
↑ Loinger, 1983, Loinger, E., 1983, A linear model for velocity anomalies: Geophys. Prosp., 31, 98–118.
↑ Rocca and Toldi, 1983, Rocca, F. and Toldi, J., 1983, Lateral velocity anomalies: 53rd Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 572–574.
↑ Sorin et al. (1996), Sorin, V., 1996, Velocity estimation in homogenous media: 56th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1515–1517.
↑ ^15.0 ^15.1 ^15.2 ^15.3 ^15.4 ^15.5 Hubral and Krey (1980), Hubral, P. and Krey, T., 1980, Interval velocities from seismic reflection time measurements: Soc. Expl. Geophys.
↑ Sheriff, 1991, Sheriff, R. E., 1991, Encyclopedic dictionary of exploration geophysics: Soc. Expl. Geophys.
↑ David, 1987, David, M., 1987, Geostatistics: in Encyclopedia of Science and Technology, 6, 141–144, Academic Press.
↑ ^18.0 ^18.1 Sherwood et al., 1986, Sherwood, J. W. C., Chen, K. C., and Wood, M., 1986, Depths and interval velocities from seismic reflection data for low-relief structures: Proc. Offshore Tech. Conf., 103–110.
↑ Kosloff et al., 1996, Kosloff, D., Sherwood, J. W. C., Koren, Z., Machet, E., and Falkovitz, Y., 1996, Velocity and interface depth determination by tomography of depth migrated gathers: Geophysics, 61, 1511–1523.
↑ Coleman et al., 1980, Coleman, J. M., Prior, D. B., and Garrison, L. E., 1980, Subaqueous sediment instabilities in the offshore Mississippi River Delta: USGS Open-File Report 80-01, 60.
↑ May et al., 1988, May, J. A., Meeder, C. A., Tinkle, A. R., and Wener, K. R., 1988, Seismic no-data zone, Offshore Mississippi Delta: Part II: using geologic information to predict acoustic properites: Proc. Offshore Tech. Conf., 75–84.
↑ ^22.0 ^22.1 Kim and Bell, 2000, Kim, H. S. and Bell, M. L., 2000, 3-D turning-ray tomography and its application to Mississippi Delta: 70th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 593–596.
↑ ^23.0 ^23.1 Zhou et al., 1992, Zhou, X., Sixta, D. P. and Angstman, B. G., 1992, Tomostatics: turning-ray tomography and statics corrections: The Leading Edge, 11, No. 12, 15–23.
↑ Bell et al., 1994, Bell, M. L., Lara, R., and Gray, W. C., 1994, Application of turning-ray tomography to the offshore Mississippi Delta: 64th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1509–1512.
↑ Lines and Treitel (1984), Lines, L. R. and Treitel, S., 1984, Tutorial: A review of least-squares inversion and its application to geophysical problems: Geophys. Prosp., 32, 159–186.
↑ ^26.0 ^26.1 ^26.2 Wang and Treitel, 1973, Wang, R. J. and Treitel, S., 1973, The determination of digital Wiener filters by means of gradient methods: Geophysics, 38, 310–326.
↑ Koehler and Taner, 1985, Koehler, F. and Taner, M. T., 1985, The use of the conjugate-gradient algorithm in the computation of predictive deconvolution operators: Geophysics, 50, 2752–2758.
↑ Claerbout, 1992, Claerbout, J. F., 1992, Earth soundings analysis — processing versus inversion: Blackwell Scientific Publications.
↑ Hestenes (1956), Hestenes, M. R., 1956, The conjugate gradient method for solving linear systems: Proc. Sympo. in Applied Math., VI, McGraw-Hill Book Co.
↑ Aldridge, 1994, Aldridge, D. F., 1994, Linearization of the eikonal equation: Geophysics, 59, 1631–1632.
↑ Landa, E., Thore, P., Sorin, V., and Koren, Z., 1991, Interpretation of velocity estimates from coherency inversion: Geophysics, 56, 1377–1383.
↑ Thorson, J. R., Gever, D. H., Swanger, H. J., Hadley, D. M., and Apsel, R. J., 1987, A model-based approach to interval velocity analysis: 57th Ann. internat. Mtg., Expanded Abstracts, 458–460.

External links

find literature about
Ageary

[ch09r9-1] Dix, 1955, Dix, C. H., 1955, Seismic velocities from surface measurements: Geophysics, 20, 68–86.

[ch09r16-2] Landa et al., 1988, Landa, E., Kosloff, D., Keydar, S., Koren, Z., and Reshef, M., 1988, A method for determination of velocity and depth from seismic data: Geophys. Prosp., 36, 223–243.

[ch09r10-3] Faye and Jeannaut, 1986, Faye, J-P. and Jeannaut, J-P., 1986, Prestack migration velocities from focusing depth analysis: 56th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 438–440.

[ch09r25-4] Reshef and Kessler, 1989, Reshef, M. and Kessler, D., 1989, Practical implementation of three-dimensional poststack depth migration: Geophysics, 54, 309–318.

[ch09r23-5] Reshef, 1997, Reshef, M., 1997, The use of 3-D prestack depth imaging to estimate layer velocities and reflector positions: Geophysics, 62, 206–210.

[ch09r24-6] Reshef, 2001, Reshef, M., 2001, Some aspects of interval velocity analysis using 3-D depth migrated gathers: Geophysics (Scheduled for publication in Jan.-Feb. 2001 issue).

[ch09r2-7] Al-Yahya, 1989, Al-Yahya, K., 1989, Velocity analysis by iterative profile migration: Geophysics, 54, 718–729.

[ch09r8-8] 8.0 ^8.1 Deregowski, 1990, Deregowski, 1990, Common-offset migrations and velocity analysis: First Break, 8, 225–234.

[ch09r4-9] 9.0 ^9.1 Bickel, 1990, Bickel, S. H., 1990, Velocity-depth ambiguity of reflection traveltimes: Geophysics, 55, 266–276.

[ch09r18-10] 10.0 ^10.1 Lines, 1993, Lines, L., 1993, Ambiguity in analysis of velocity and depth: Geophysics, 58. 596–597.

[ch09r21-11] Lynn and Claerbout, 1982, Lynn, W. S. and Claerbout, J. F., 1982, Velocity estimation in laterally varying media: Geophysics, 47, 884–897.

[ch09r20-12] Loinger, 1983, Loinger, E., 1983, A linear model for velocity anomalies: Geophys. Prosp., 31, 98–118.

[ch09r26-13] Rocca and Toldi, 1983, Rocca, F. and Toldi, J., 1983, Lateral velocity anomalies: 53rd Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 572–574.

[ch09r29-14] Sorin et al. (1996), Sorin, V., 1996, Velocity estimation in homogenous media: 56th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1515–1517.

[ch09r12-15] 15.0 ^15.1 ^15.2 ^15.3 ^15.4 ^15.5 Hubral and Krey (1980), Hubral, P. and Krey, T., 1980, Interval velocities from seismic reflection time measurements: Soc. Expl. Geophys.

[ch09r27-16] Sheriff, 1991, Sheriff, R. E., 1991, Encyclopedic dictionary of exploration geophysics: Soc. Expl. Geophys.

[ch09r7-17] David, 1987, David, M., 1987, Geostatistics: in Encyclopedia of Science and Technology, 6, 141–144, Academic Press.

[ch09r28-18] 18.0 ^18.1 Sherwood et al., 1986, Sherwood, J. W. C., Chen, K. C., and Wood, M., 1986, Depths and interval velocities from seismic reflection data for low-relief structures: Proc. Offshore Tech. Conf., 103–110.

[ch09r15-19] Kosloff et al., 1996, Kosloff, D., Sherwood, J. W. C., Koren, Z., Machet, E., and Falkovitz, Y., 1996, Velocity and interface depth determination by tomography of depth migrated gathers: Geophysics, 61, 1511–1523.

[ch09r6-20] Coleman et al., 1980, Coleman, J. M., Prior, D. B., and Garrison, L. E., 1980, Subaqueous sediment instabilities in the offshore Mississippi River Delta: USGS Open-File Report 80-01, 60.

[ch09r22-21] May et al., 1988, May, J. A., Meeder, C. A., Tinkle, A. R., and Wener, K. R., 1988, Seismic no-data zone, Offshore Mississippi Delta: Part II: using geologic information to predict acoustic properites: Proc. Offshore Tech. Conf., 75–84.

[ch09r13-22] 22.0 ^22.1 Kim and Bell, 2000, Kim, H. S. and Bell, M. L., 2000, 3-D turning-ray tomography and its application to Mississippi Delta: 70th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 593–596.

[ch09r32-23] 23.0 ^23.1 Zhou et al., 1992, Zhou, X., Sixta, D. P. and Angstman, B. G., 1992, Tomostatics: turning-ray tomography and statics corrections: The Leading Edge, 11, No. 12, 15–23.

[ch09r3-24] Bell et al., 1994, Bell, M. L., Lara, R., and Gray, W. C., 1994, Application of turning-ray tomography to the offshore Mississippi Delta: 64th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1509–1512.

[ch09r19-25] Lines and Treitel (1984), Lines, L. R. and Treitel, S., 1984, Tutorial: A review of least-squares inversion and its application to geophysical problems: Geophys. Prosp., 32, 159–186.

[ch09r31-26] 26.0 ^26.1 ^26.2 Wang and Treitel, 1973, Wang, R. J. and Treitel, S., 1973, The determination of digital Wiener filters by means of gradient methods: Geophysics, 38, 310–326.

[ch09r14-27] Koehler and Taner, 1985, Koehler, F. and Taner, M. T., 1985, The use of the conjugate-gradient algorithm in the computation of predictive deconvolution operators: Geophysics, 50, 2752–2758.

[ch09r5-28] Claerbout, 1992, Claerbout, J. F., 1992, Earth soundings analysis — processing versus inversion: Blackwell Scientific Publications.

[ch09r11-29] Hestenes (1956), Hestenes, M. R., 1956, The conjugate gradient method for solving linear systems: Proc. Sympo. in Applied Math., VI, McGraw-Hill Book Co.

[ch09r1-30] Aldridge, 1994, Aldridge, D. F., 1994, Linearization of the eikonal equation: Geophysics, 59, 1631–1632.

[ch09r17-31] Landa, E., Thore, P., Sorin, V., and Koren, Z., 1991, Interpretation of velocity estimates from coherency inversion: Geophysics, 56, 1377–1383.

[ch09r30-32] Thorson, J. R., Gever, D. H., Swanger, H. J., Hadley, D. M., and Apsel, R. J., 1987, A model-based approach to interval velocity analysis: 57th Ann. internat. Mtg., Expanded Abstracts, 458–460.

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15]

[16]

[17]

[18]

[19]

[20]

[21]

[22]

[23]

[24]

[25]

[26]

[27]

[28]

[29]

[30]

[31]

[32]

User:Ageary/Yilmaz 9

Contents

9.0 Introduction

Inversion methods for data modeling

Inversion procedures for earth modeling

Velocity-depth ambiguity

Model representation and visualization

9.1 Models with horizontal layers

Dix conversion

Coherency inversion

Near-surface layer with laterally varying velocities

9.2 Model with low-relief structure

Stacking velocity inversion

Coherency inversion

Velocity resolution

9.3 Model with complex overburden structure

Image gathers

Constant half-space velocity analysis

9.4 Model building

Time-to-depth conversion

Time structure maps

Interval velocity maps

Depth structure maps

Calibration to well tops

Layer-by-layer inversion

Structure-independent inversion

9.5 Model updating

Residual moveout analysis

Reflection traveltime tomography

Limitations in resolving velocity-depth ambiguity by tomography

Turning-ray tomography

Exercises

Appendix J: Data modeling by inversion

J.1 The generalized linear inversion

J.2 The GLI formalism of deconvolution

J.3 Applications of the GLI technique

J.4 Dix conversion

J.5 Map processing

J.6 Reflection traveltime tomography

J.7 Threshold for velocity-depth ambiguity

See also

References

External links

Navigation menu

Search